function [t, T, Tcal, wla, wlc, a_uncorr, c_uncorr, a_tscorr, c_tscorr]=process_acs_data(Sal)
%read and process the acs data, then plot time series and spectra of the ac data
%Select a directory to start
%Input:
%   Sal: salinity
%Outputs:
%   t: time vector, sec
%   Tcal: factory calibration temperature, oC
%   T: temperature of sample, oC
%   wla, wlc: wavelength vectors for absorption and attenuation, respectively
%   a_uncorr, c_uncorr: uncorrected absorption and attenuation data
%   a_tscorr, c_tscorr: absorption and attenuation data after T&S correction
%
%How to run the code:
%   1.Place 'Sullivan_etal_2006_instrumentspecific.xls' in the same folder
%   where you code stays
%   2.factory calibration file *.dev needs to be in the same folder where
%   you data files stay
%   3.copy 'process_acs_data(Sal)' into you Matlab command window, but
%   replace the Sal with the actual salinity of your sample, then press
%   Enter
%   4.it runs in batch. So just wait for a little bit, then you can find
%   the .mat files in the foler you just selected
%
%Written by Linhai Li, July 9, 2013, Updated by Linhai Li, July 12, 2013


%Read T&S correction coefficient
tmp = xlsread('Sullivan_etal_2006_instrumentspecific.xls');

path=uigetdir;
current_folder=cd(path);

filelist=dir('*.dat');
%[filename, path]=uigetfile({'*.dat'}, 'Select a ACS file', 'MultiSelect', 'off');

%Get the factory calibration file
fac_cal_file=dir(fullfile(path, '*.dev'));
fac_fid=fopen(fullfile(path, fac_cal_file.name));
for i=1:3
    junk=fgetl(fac_fid);
end
fgets(fac_fid, 6);
Tcal=fscanf(fac_fid, '%f', 1);
fclose(fac_fid);

for nf=1:numel(filelist)
    %Get data from actual data files of sample
    filename=filelist(nf).name;
    
    temp_ind=strfind(filename, 'p');
    for j=1:numel(temp_ind);
        if isempty(str2num(filename(temp_ind(j)-1))) || isempty(str2num(filename(temp_ind(j)+1)))
            continue;
        end
        
        if (0<=str2num(filename(temp_ind(j)-1))<=9) && (0<=str2num(filename(temp_ind(j)+1))<=9)
            T=str2num(filename(temp_ind(j)-2:temp_ind(j)-1))+str2num(filename(temp_ind(j)+1))/10;
        else
            continue;
        end
    end
    
    fid=fopen(fullfile(path,filename),'r');
    
    str_token='1			; acquisition binsize';
    str_token(str_token==' ')='';
    
    %skip non-data lines
    line=fgetl(fid);
    line(line==' ')='';
    while strcmp(line, str_token)==0
        line=fgetl(fid);
        line(line==' ')='';
    end
    
    %get the wavelength for C and A
    line=fgetl(fid);
    wl=regexp(line, '	', 'split');
    wl=wl(2:end)';
    wl=char(wl);
    wl=str2num(wl(:,2:end));
    
    %number of wavelength
    nwl=size(wl,1);
    
    wlc=wl(1:nwl/2)'; %first half is wavelength for C
    wla=wl(nwl/2+1:end)'; %second half is wavelength for A
    
    %T&S correction coefficient
    phi_Ta=interp1(tmp(:,1),tmp(:,2),wla,'linear');
    phi_Tc=interp1(tmp(:,1),tmp(:,2),wla,'linear');
    phi_Sc=interp1(tmp(:,1),tmp(:,4),wlc,'linear');
    phi_Sa=interp1(tmp(:,1),tmp(:,6),wla,'linear');
    
    i=1;
    while(feof(fid)==0)
        line=fgetl(fid);
        data=str2num(line);
        if size(data)<nwl+1
            break
        end
        t(i,1)=data(1)/1000; %first column of junk is time
        c_uncorr(i,:)=data(2:nwl/2+1)'; %second column to nwl/2+1 column is C data
        a_uncorr(i,:)=data(nwl/2+2:nwl+1)'; %nwl/2+2 column to nwl+1 column is A data
        i=i+1;
    end
    fclose(fid);
    
    for i=1:numel(t)
        a_tscorr(i,:)=a_uncorr(i,:)-(T-Tcal)*phi_Ta-Sal*phi_Sa;
        c_tscorr(i,:)=c_uncorr(i,:)-(T-Tcal)*phi_Tc-Sal*phi_Sc;
    end
    
    %save data to mat files
    data_struct=genvarname(filename(1:end-4));
    eval([data_struct '.time=t;']);
    eval([data_struct '.wla=wla;']);
    eval([data_struct '.wlc=wlc;']);
    eval([data_struct '.a_uncorr=a_uncorr;']);
    eval([data_struct '.c_uncorr=c_uncorr;']);
    eval([data_struct '.a_tscorr=a_tscorr;']);
    eval([data_struct '.c_tscorr=c_tscorr;']);
    save([filename(1:end-4) '.mat'], data_struct);
    
    % %plot time series
 
    figure('Name', ['Time Series_' filename], 'NumberTitle', 'off');
    subplot(2,1,1);
    plot(t, a_uncorr(:,:));
    ylabel('a_uncorr (1/m)', 'interpreter', 'none');
    
    subplot(2,1,2);
    plot(t, c_uncorr(:,:));
    ylabel('c_uncorr (1/m)', 'interpreter', 'none');
    xlabel('Time (sec)');
    
    figure('Name', ['Spectra_', filename], 'NumberTitle', 'off');
    subplot(2,1,1);
    hp=plot(wla, a_uncorr(1:10:end, :), '-', wla, a_tscorr(1:10:end, :), '--');
    ylabel('a_tscorr (1/m)', 'interpreter', 'none');
    legend(hp([1,end]), {'Uncorr', 'Tscorr'});
    
    subplot(2,1,2);
    hp=plot(wlc, c_uncorr(1:10:end, :), '-', wlc, c_tscorr(1:10:end, :), '--');
    xlabel('Wavelength (nm)');
    ylabel('c_tscorr (1/m)', 'interpreter', 'none');
    legend(hp([1,end]), {'Uncorr', 'Tscorr'});

    %clear data in case contaminate next file
    t=[];
    a_uncorr=[];
    c_uncorr=[];
    a_tscorr=[];
    c_tscorr=[];
end

cd(current_folder);


