%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script integrates the counts under element-specific energy intervals 
% of the NGR spectra. Concentrations are calculated by comparing these
% counts to counts in a spectrum of a standard of known concentration.
%
% The algorithm is designed to quantify K, U and Th concentrations Hole by
% Hole. Therefore, we strongly suggest not to mix different Holes from one
% Site.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Potassium
clear; clc; close all;
load('$PATH/energylevels_K.mat') % !!! USER INPUT NEEDED !!! --> Specify where the algorithm finds the counts under the K peak in the NGR spectrum of the K standard
load('$PATH/results.mat') % !!! USER INPUT NEEDED !!! --> Specify where the algorithm finds the extracted spectra.

%%%% In this loop, we look for the channels that correspond to the 1270 -
%%%% 1570 keV energy range that we are interested in, i.e. the energy range 
%%%% of the 40K peak.  
cal_files=dir('background_calibration/*alibration*.Spe');
for j=1:length(cal_files)
    fid = fopen(strcat('background_calibration/',cal_files(j).name));
    tline = fgetl(fid);
                    
    while ~strcmp('$ENER_FIT:',tline)
        tline = fgetl(fid);
         % disp(tline)
    end
    tline = fgetl(fid);
    ENER_FIT(j,:)=cell2mat(textscan(tline,'%f'));
    
    Ener(j,:)=ENER_FIT(j,1)+[0:1023]*ENER_FIT(j,2);
    idx1=find(Ener(j,:)>=1270);
    idx2=find(Ener(j,:)<=1570);
    eval(strcat('K_peak_ind', num2str(j),'= intersect(idx1,idx2);'));
end

%%%% In this loop, the area under the 40K peak is integrated for each
%%%% NGR measurement. And then the ratio in peak area between
%%%% sample and K standard is calculated. The area under the 40K
%%%% peak for the K standard is given in the variables K_counts1 -
%%%% K_counts8.
for j=1:length(data_nobg_ec(1,:))
    detector=data_info(3,j);
    eval(strcat('data_K_counts(j)=sum(data_nobg_ec(K_peak_ind', num2str(detector),',j));')); % K_peak_ind1, K_peak_ind2, ... , K_peak_ind8 contain the indeces of the channels that fall within the 40K peak (1260 - 1560 keV)
    % In the next line, we divide the total counts under the peak by the
    % measurement time of the sample, to express it in cps. The same is
    % done for the K standard (measurement time = 1500 s). Next, the
    % area of both peaks (both in cps) is compared and expressed as the
    % ratio of the area of the sample and the area of the standard. 
    eval(strcat('data_K_counts(j)=(data_K_counts(j)/', num2str(t_measurement_sample(j)),')/(K_counts',num2str(detector),'/1500);')); 
end

%%%% The density of the sediment needs to be taken into account before an
%%%% etimate of the elemental concentration can be made. Therefore, we
%%%% convolve a Gaussian function that represents the detector response
%%%% function with GRA bulk densities.
[NUM,TXT,RAW]=xlsread('$PATH/GRA_356_U1463B.xlsx'); % !!! USER INPUT NEEDED !!! --> Specify here where the algorithm can find the .xlsx file that contains GRA bulk density. This data can be downloaded from LIMS. Make sure only to select the WRMSL data (not SHMSL). 
GRA=NUM(:,[4,6,8,9,11]);

for j=1:length(data_info)
    idx1=find(GRA(:,1)==data_info(1,j)); % Find the right core
    idx2=find(GRA(:,2)==data_info(2,j)); % Find the right section
    idx3=find(abs(GRA(:,3)-data_info(6,j))<20); % Find GRA measurements within +- 20 cm from the NGR measurement
    idx=intersect(intersect(idx1,idx2),idx3);
    
    if isempty(idx)==1
        GRA_i(j)=NaN;
        mCSFA(j)=NaN;
    else    
        GRA_ti=GRA(idx,5);
        GRA_weights=normpdf((GRA(idx,3)-data_info(6,j)),0,6.37);
        GRA_i(j)=sum(GRA_ti.*GRA_weights)/sum(GRA_weights); 

        idx3=find(abs(GRA(:,3)-data_info(6,j))<=20);
        idx=intersect(intersect(idx1,idx2),idx3);
        mCSFA(j)=GRA(idx(1),4)-(GRA(idx(1),3)-data_info(6,j))/100;
    end
end

rho_STND=2.07;
K_perc_STND=6.43; % Potassium standard density and concentration

% data_K_counts contains the ratio between the area of the sample's K-peak
% and the standard's K-peak. The following lines converts this ratio to the sample's
% concentration by taking into account the standard's density and
% concentration and the sample's density.
for j=1:length(data_info) 
data_K_perc(j)=data_K_counts(j)*rho_STND/GRA_i(j)*K_perc_STND;
end

data_K_perc=horzcat(mCSFA',data_K_perc');
data_K_perc=horzcat((data_info([1,2,6],:))',data_K_perc);
data_K_perc=horzcat(data_K_perc,GRA_i');

data_K_perc=sortrows(data_K_perc,4);
figure; plot(data_K_perc(:,4),data_K_perc(:,5));
xlabel('Depth m CSF-A'); ylabel('K (wt.%)')

%column headings in the output file are Core, Section, Offset (cm), Depth
%(m CSF-A), K (wt.%), GRA bulk density (g/cm3)
csvwrite('U9999A_K_concentration.csv',data_K_perc) % !!! USER INPUT NEEDED !!! --> Choose an appropriate file name.

%% Thorium
load('$PATH/energylevels_Th.mat') % !!! USER INPUT NEEDED !!! --> Specify here where the algorithm can find the area under the Thorium peak in the NGR spectrum of the Thorium standard

%%%% In this loop, we look for the channels that correspond to the 2410 -
%%%% 2710 keV energy range that we are interested in, i.e. the energy range 
%%%% of the Th peak.  
cal_files=dir('background_calibration/*alibration*.Spe');
for j=1:length(cal_files)
    fid = fopen(strcat('background_calibration/',cal_files(j).name));
    tline = fgetl(fid);
                    
    while ~strcmp('$ENER_FIT:',tline)
        tline = fgetl(fid);
         % disp(tline)
    end
    tline = fgetl(fid);
    ENER_FIT(j,:)=cell2mat(textscan(tline,'%f'));
    
    Ener(j,:)=ENER_FIT(j,1)+[0:1023]*ENER_FIT(j,2);
    idx1=find(Ener(j,:)>=2410);
    idx2=find(Ener(j,:)<=2710);
    eval(strcat('Th_peak_ind', num2str(j),'= intersect(idx1,idx2);'));
end

%%%% In this loop, the area under the 208 Tl peak is integrated for each
%%%% NGR measurement. And then the ratio in peak area between
%%%% sample and Th standard is calculated.
for j=1:length(data_nobg_ec(1,:))
    detector=data_info(3,j);
    eval(strcat('data_Th_counts(j)=sum(data_nobg_ec(Th_peak_ind', num2str(detector),',j));')); % Th_peak_ind1, Th_peak_ind2, ..., Th_peak_ind8 contain the indeces of the channels that fall within the 208Tl peak (2410 - 2710 keV)
    % In the next line, we divide the total counts under the peak by the
    % measurement time of the sample, to express it in cps. The same is
    % done for the Th standard (measurement time = 1500 s). Next, the
    % area of both peaks (both in cps) is compared and expressed as the
    % ratio of the area of the sample and the area of the standard. 
    eval(strcat('data_Th_counts(j)=(data_Th_counts(j)/', num2str(t_measurement_sample(j)),')/(Th_counts',num2str(detector),'/1500);'));  
end

rho_STND=2.156;
Th_perc_STND=20.7; % Thorium standard concentration and density

% data_Th_counts contains the ratio between the area of the sample's Th-peak
% and the standard's Th-peak. The following lines converts this ratio to the sample's
% concentration by taking into account the standard's density and
% concentration and the sample's density.
for j=1:length(data_info) 
data_Th_perc(j)=data_Th_counts(j)*rho_STND/GRA_i(j)*Th_perc_STND;
end

data_Th_perc=horzcat(mCSFA',data_Th_perc');
data_Th_perc=horzcat((data_info([1,2,6],:))',data_Th_perc);
data_Th_perc=horzcat(data_Th_perc,GRA_i');

data_Th_perc=sortrows(data_Th_perc,4);
figure; plot(data_Th_perc(:,4),data_Th_perc(:,5))
xlabel('Depth m CSF-A'); ylabel('Th (ppm)')

%column headings in the output file are Core, Section, Offset (cm), Depth
%(m CSF-A), Th (ppm), GRA bulk density (g/cm3)
csvwrite('U999A_Th_ppm.csv',data_Th_perc) % !!! USER INPUT NEEDED !!! --> Choose an appropriate file name.

%% Uranium

load('$PATH/energylevels_U.mat') % !!! USER INPUT NEEDED !!! --> Specify here where the algorithm can find the area under the Thorium peak in the NGR spectrum of the Thorium standard
load('results.mat')

%%%% In this loop, we look for the channels that correspond to the 1660 -
%%%% 1860 keV energy range that we are interested in, i.e. the energy range 
%%%% of the U peak.  
cal_files=dir('background_calibration/*alibration*.Spe');
for j=1:length(cal_files)
    fid = fopen(strcat('background_calibration/',cal_files(j).name));
    tline = fgetl(fid);
                    
    while ~strcmp('$ENER_FIT:',tline)
        tline = fgetl(fid);
    end
    tline = fgetl(fid);
    ENER_FIT(j,:)=cell2mat(textscan(tline,'%f'));
    
    Ener(j,:)=ENER_FIT(j,1)+[0:1023]*ENER_FIT(j,2);
    idx1=find(Ener(j,:)>=1660);
    idx2=find(Ener(j,:)<=1860);
    eval(strcat('U_peak_ind', num2str(j),'= intersect(idx1,idx2);'));
end

%%%% In this loop, the area under the 214 Bi peak is integrated for each
%%%% NGR measurement. And then the ratio in peak area between
%%%% sample and U standard is calculated
for j=1:length(data_nobg_ec(1,:))
    detector=data_info(3,j);
    eval(strcat('data_U_counts(j)=sum(data_nobg_ec(U_peak_ind', num2str(detector),',j));')); % U_peak_ind1, U_peak_ind2, ..., U_peak_ind8 contain the indeces of the channels that fall within the 208Tl peak (2410 - 2710 keV)
    % In the next line, we divide the total counts under the peak by the
    % measurement time of the sample, to express it in cps. The same is
    % done for the Th standard (measurement time = 1800 s). Next, the
    % area of both peaks (both in cps) is compared and expressed as the
    % ratio of the area of the sample and the area of the standard. 
    eval(strcat('data_U_counts(j)=(data_U_counts(j)/', num2str(t_measurement_sample(j)),')/(U_counts',num2str(detector),'/1800);'));  
end

rho_STND=2.23;
U_ppm_STND=0.1869; % U standard concentration and density

% data_U_counts contains the ratio between the area of the sample's U-peak
% and the standard's U-peak. The following lines converts this ratio to the sample's
% concentration by taking into account the standard's density and
% concentration and the sample's density.
for j=1:length(data_info) 
data_U_ppm(j)=data_U_counts(j)*rho_STND/GRA_i(j)*U_ppm_STND;
end

data_U_ppm=horzcat(mCSFA',data_U_ppm');
data_U_ppm=horzcat((data_info([1,2,6],:))',data_U_ppm);
data_U_ppm=horzcat(data_U_ppm,GRA_i');

data_U_ppm=sortrows(data_U_ppm,4);
figure; plot(data_U_ppm(:,4),data_U_ppm(:,5))
xlabel('Depth m CSF-A'); ylabel('U (ppm)')

%column headings in the output file are Core, Section, Offset (cm), Depth
%(m CSF-A), U (ppm), GRA bulk density (g/cm3)
csvwrite('U9999A_U_ppm.csv',data_U_ppm) % !!! USER INPUT NEEDED !!! --> Choose an appropriate file name.