%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script extracts the NGR energy spectra from the zip files,         
% subtracts the background spectrum, and carries out edge corrections.    
%
% 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; clc; close all;
cd('$PATH') % !!! USER INPUT NEEDED !!! --> Specify here the directory to which you downloaded the NGR .zip files
z_files=dir('*.zip'); % list of zipped files in the specified directory. File names should look like: ngr_sample_123-U1234A-1H-1-W_SHLF7001271_2015912172852.zip. In this example, this sample comes from Exp. 123, Site 1234, Hole A, Core 1, Type H, Section 1, Working half. 

for i=1:length(z_files)
    unzip(z_files(i).name,'unzipped'); % unzip all files and move them to the newly created folder 'unzipped'
end

movefile('unzipped/ngr_sample_STND*BACK*','background_calibration/') % Move background files to different folder
movefile('unzipped/*alibration*','background_calibration/') % Move calibration files to different folder
movefile('unzipped/ngr_edge*','background_calibration/') % Move edge correction file to different folder
spe_files=dir('unzipped/ngr_sample_*.spe'); % list of spe-files in the unzipped directory

data=zeros(1024,length(spe_files)); % This matrix will contain all NGR spectra
data_info=zeros(6,length(spe_files)); % This matrix will contain information on the spectrogram measurement 

%In this loop we fill in the NGR spectra in data, the information on
%the measured section in data_info, and the time of measurement.
for j=1:length(spe_files)
    fid = fopen(strcat('unzipped/',spe_files(j).name));
    tline = fgetl(fid);
    
    info=regexp(spe_files(j).name,'\d+','match');
    info=str2double(info);
    if length(info)==7 %In case we measured a Core Catcher, the variable "info" will have length = 7
        data_info(1,j)=info(3);
        data_info(2,j)=-999; % CC's get a section number -999
        data_info(3:4,j)=info([6:7]);
    else % In case we measured a 'normal' Section
    data_info(1:4,j)=info([3:4,7:8]); %Core, Section, Det, Position 1 or 2, length of section, offset (last 2 informations are filled later in script) 
    end
    
    while ~strcmp('$MEAS_TIM:',tline)
        tline = fgetl(fid);
    end
    tline = fgetl(fid);
    temp=cell2mat(textscan(tline,'%f'));
    t_measurement_sample(j)=temp(1);
     
    while ~strcmp('$DATA:',tline)
        tline = fgetl(fid);
    end
    tline = fgetl(fid);
    i=0;
    while i<1024
        i=i+1;
        tline = fgetl(fid);
        data(i,j)=cell2mat(textscan(tline,'%f')); % For the considered NGR measurement, we fill the matrix 'data(:,j)' with the NGR counts at the 1024 energy levels 
    end
    fclose(fid);
    j
end

%%% In this loop, we fill in the two last bits of information in data_info.
%%% These are the length of the section and the offset of the NGR reading
%%% within the section.
for j=1:length(spe_files)
    fid = fopen(strcat('unzipped/',spe_files(j).name));
    tline = fgetl(fid);
    for k=1:7
        tline = fgetl(fid);
    end
    tline = fgetl(fid);
    info=regexp(tline,'\d+','match');
    info=str2double(info);
    if isempty(info)==1
        info=150;
        warning(['Length of section was not given for Core ' num2str(data_info(1,j)) ' Section ' num2str(data_info(2,j)) '. We therefore assume a length of 150 cm'])
    end 
    data_info(5,j)=info(1); % Length of the section
    tline = fgetl(fid);
    info=regexp(tline,'\d+','match');
    info=str2double(info); 
    data_info(6,j)=info(1); % Offset within section
    fclose(fid);
end

%%% When a section < 150 cm, some detectors were not covered by core during
%%% the measurement. These are deleted from the dataset.
indeces=find(data_info(5,:)<data_info(6,:));
data(:,indeces)=[];
data_info(:,indeces)=[];
t_measurement_sample(indeces)=[];

%%% Extract background spectra
bg_files=dir('background_calibration/ngr_sample_STND*.spe'); % list of spe-files in the unzipped directory
bg=zeros(1024,length(bg_files));
for j=1:length(bg_files)
    fid = fopen(strcat('background_calibration/',bg_files(j).name));
    tline = fgetl(fid);
    
    if j==1
    while ~strcmp('$MEAS_TIM:',tline)
        tline = fgetl(fid);
    end
    tline = fgetl(fid);
    temp=cell2mat(textscan(tline,'%f'));
    t_measurement_backgr=temp(1);
    end
    
    while ~strcmp('$DATA:',tline)
        tline = fgetl(fid);
    end
    tline = fgetl(fid);
    i=0;
    while i<1024
        i=i+1;
        tline = fgetl(fid);
        bg(i,j)=cell2mat(textscan(tline,'%f'));
        bg(i,j)=bg(i,j)/t_measurement_backgr;
    end
    fclose(fid);
    j
end

%%%%%% Subtract background spectra from data spectra
if length(bg(1,:))==16
    for j=1:length(data)
        if data_info(4,j)==1 % Measurement position 1
            data_nobg(:,j)=data(:,j)-(bg(:,data_info(3,j))*t_measurement_sample(j));
        elseif data_info(4,j)==2 % Measurement position 2
            data_nobg(:,j)=data(:,j)-(bg(:,data_info(3,j)+8)*t_measurement_sample(j));
        end
    end
elseif length(bg(1,:))==8 % In case NGR measurements were only taken at one measurement position
    for j=1:length(data(1,:))
        data_nobg(:,j)=data(:,j)-bg(:,data_info(3,j));
    end
end

%%%%%% Make edge corrections
edge=[20	1
19	1
18	1.01
17	1.01
16	1.02
15	1.02
14	1.03
13	1.03
12	1.05
11	1.08
10	1.09
9	1.14
8	1.17
7	1.21
6	1.26
5	1.33
4	1.42
3	1.51
2	1.62
1	1.76
0	2]; % Table with edge correction coefficients
data_nobg_ec=data_nobg;
for j=1:length(data(1,:))
    if data_info(6,j)<20
        index=find(edge(:,1)==data_info(6,j));
        data_nobg_ec(:,j)=data_nobg(:,j).*edge(index,2);
    elseif data_info(5,j)-data_info(6,j)<20
        index=find(edge(:,1)==(data_info(5,j)-data_info(6,j)));
        data_nobg_ec(:,j)=data_nobg(:,j).*edge(index,2);
    end
end
    
%%%% plot example
k=100;
plot(data_nobg_ec(:,k));
title(strcat('Core ',num2str(data_info(1,k)),' Section ',num2str(data_info(2,k)),' Offset ',num2str(data_info(6,k))),'Interpreter', 'none');

%%%% Save the extracted and corrected NGR Spectra.
save('results.mat','spe_files','data','data_info','data_nobg_ec','t_measurement_backgr','t_measurement_sample');