%Analyze lab ac-s data of dead sea waters
clear all
close all
%Data from sullivan's paper
tmp = xlsread('Sullivan_etal_2006_instrumentspecific.xls');


opts = optimset('fminsearch');      
opts = optimset(opts,'MaxIter',20000000); 
opts = optimset(opts,'MaxFunEvals',20000);
opts = optimset(opts,'TolX',1e-8);
opts = optimset(opts,'TolFun',1e-8);

D=dir('*_DI.dat');
numfiles=length(D);

%reading in the different calibration files
for i=1:numfiles
    filename=D(i).name;
    [t,dat,wla,wlc]=rd_wetview_acs(filename);
    c_cal=median(dat(end-240:end,[1:88]));
    a_cal=median(dat(end-240:end,[89:176]));
    c_cal_std=(prctile(dat(end-240:end,[1:88]),84)-prctile(dat(end-240:end,[1:88]),16))/2;
    a_cal_std=(prctile(dat(end-240:end,[89:176]),84)-prctile(dat(end-240:end,[89:176]),16))/2;   
    phi_T_c=interp1(tmp(:,1),tmp(:,2),wlc,'linear');
    phi_T_a=interp1(tmp(:,1),tmp(:,2),wla,'linear');
    dt(i)=datenum(filename(8:21),'yyyymmddHHMMSS')

    I = find(wlc>=700 & wlc<=743);  % spectral range for optimization (710 to 750nm)
    delT = 0;
    [xc] = fminsearch(@f_T, [delT,0], opts, a_cal(I), phi_T_c(I), wlc(I));

    I = find(wla>=700 & wla<=743);  % spectral range for optimization (710 to 750nm)
    delT = 0;
    [xa] = fminsearch(@f_T, [delT,0], opts, a_cal(I), phi_T_a(I), wla(I));
    plot(wlc,c_cal,'b:',wlc,c_cal-xc(1)*phi_T_c,'b-',wla,a_cal,'k:',wla,a_cal-xa(1)*phi_T_a,'k-');
    hold on
    A(i,:)=a_cal-xa(1)*phi_T_a;
    C(i,:)=c_cal-xc(1)*phi_T_c;
    A_std(i,:)=a_cal_std;
    C_std(i,:)=c_cal_std;
end
hold off
legend('c-cal','c-cal T-cor','a-cal','a-cal T-cor','Location','NorthWest')
 
data=[dt',A,C,A_std,C_std]

save NAAMES03_ACS_DIW_cals filename data wla wlc
hold off    
   
    
figure;
plot(A_std')