clear all
close all
%make a large matrix with all CNR data
for i=[1:99] 
    path='/Volumes/Data2/TaraEuropa/Partner_data/CNR-P/Absorption_spectra/';
    filename=[path,'ST_ag',num2str(i),'.mat'];
    load(filename)
    wl_ag=wl;
    B(i,:)=ag_ave;
    C(i,:)=ag_std;
end

%load SUVF data
opts = delimitedTextImportOptions("NumVariables", 19);
% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ",";
% Specify column names and types
opts.VariableNames = ["St", "time_start", "time_end", "lat_min", "lat_max", "lon_min", "lon_max", "Temp", "Temp_var", "Sal", "Sal_var", "Chl", "Chl_var", "POC", "POC_var", "gamma", "gamma_var", "FDOM", "FDOM_var"];
opts.VariableTypes = ["double", "datetime", "datetime", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"];
% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% Specify variable properties
opts = setvaropts(opts, "time_start", "InputFormat", "dd-MMM-yyyy HH:mm:ss");
opts = setvaropts(opts, "time_end", "InputFormat", "dd-MMM-yyyy HH:mm:ss");
% Import the data
A = readtable("/Volumes/Data2/TaraEuropa/Partner_data/Merge_data_discrete_inline/Summary_output_stations_within_15min_from_CTD.csv", opts);
%% Clear temporary variables
clear opts

j=1
for i=1:length(A.St)
    if A.FDOM(i) > 0
        D.FDOM(j)=A.FDOM(i);
        D.FDOM_var(j)=A.FDOM_var(i);
        D.ag(j,:)=B(A.St(i),:);
        D.ag_var(j,:)=C(A.St(i),:);
        j=j+1;
    end
end


for i=1:length(wl)
    AA=corrcoef(D.ag(:,i),D.FDOM);
    wl(i);
    plot(D.FDOM,D.ag(:,i),'*')
    R(i)=AA(1,2);
    p(i,:)= robustfit(D.FDOM,D.ag(:,i));
end

save FDOM_fit_coefficients p wl
AAA=[wl p];
writetable(AAA,'./parameter_to_convert_FDOM_to_ag_TaraEuropa.csv')

I=find(wl<550 & wl>400);
x=sum((p(I,2)*D.FDOM+p(I,1)*ones(1,length(D.FDOM))-D.ag(:,I)')./D.ag(:,I)')/149*100;
plot([1:93],x,'*')


%for wl=430nm
% figure;
% plot(Fdom_med,B(:,271),'*')
% hold on
% for i=1:length(Fdom_95)
% plot([Fdom_med(i),Fdom_med(i)],[B(i,271)-C(i,271),B(i,271)+C(i,271)],'-b')
% plot([Fdom_5(i),Fdom_95(i)],[B(i,271),B(i,271)],'-b')
% end
% hold off

% opts = optimset('fminsearch');      
% opts = optimset(opts,'MaxIter',4000); 
% opts = optimset(opts,'MaxFunEvals',2000);   % usually 100*number of params
% opts = optimset(opts,'TolFun',1e-9);
% 
% wl_ag=wl(151:301)
% x0=[p(301,1), 0.017];
% %minimization routine
% x1 = fminsearch(@exp_fit,x0,opts,p(151:301,1),wl_ag)
% x0=[p(301,1), 0.017,1];
% x2 = fminsearch(@stretched_exp_fit,x0,opts,p(151:301,1),wl_ag)
% x0=[p(301,1), 0.017,0];
% x3 = fminsearch(@exp_fit_offset,x0,opts,p(151:301,1),wl_ag)
% %plot data and fit
% plot(wl_ag, p(151:301,1),'k', wl_ag,x1(1)*exp(-x1(2)*(wl_ag-400)),'b',wl_ag,x2(1)*exp(-x2(2)*(wl_ag-400).^x2(3)),':b',wl_ag,x3(3)+x3(1)*exp(-x3(2)*(wl_ag-400)),'r')
% 
% 
%  plot(wl_ag,p(151:301,2),'b',LineWidth=2')
% % hold on;
% % plot(wl,1.35*exp(-0.026*(wl-400).^0.9),'k',LineWidth=2')
% % axis([400 700 0 1.5])
% % hold off
% 
% for i=1:length(Fdom_med)
%     figure(1)
%     plot(wl_ag,p(151:301,2)+Fdom_med(i).*x1(1).*exp(-x1(2).*(wl_ag-400)),'b',wl_ag,p(151:301,2)+Fdom_5(i).*x1(1).*exp(-x1(2).*(wl_ag-400)),':b',wl_ag,p(151:301,2)+Fdom_95(i).*x1(1).*exp(-x1(2).*(wl_ag-400)),':b','LineWidth',2);
%     hold on
%     plot(wl_ag,B(i,151:301),'k',wl_ag,B(i,151:301)+C(i,151:301),':k',wl_ag,B(i,151:301)-C(i,151:301),':k','LineWidth',2);
%     hold off
%     figure(2)
%     plot(wl_ag,(B(i,151:301)-(p(151:301,2)+Fdom_med(i).*x1(1).*exp(-x1(2)))')./B(i,151:301));
%     pause
% end
% 
% function y = exp_fit(x0,spec,l)
% y=sum(((spec-x0(1).*exp(-x0(2).*(l-400)))).^2);
% end
% 
% function y = exp_fit_offset(x0,spec,l)
% y=sum(((spec-(x0(3)+x0(1).*exp(-x0(2).*(l-400))))).^2);
% end
% 
% function y = stretched_exp_fit(x0,spec,l);
% y=sum(((spec-x0(1).*exp(-x0(2).*(l-400).^x0(3)))).^2);
% end

 
% pause
% 
% for i=1:length(I)
%     plot(wl_ag,B(I(i),:),'k',wl_ag,Fdom_med(I(i)).*0.21*exp(-0.0142*(wl_ag-500))+0.024,'r');
%     axis([400 700 0 B(I(i),300)+0.1])
% end


%0.23*exp(-0.013*(wl-500) +0.006 - Baltic (21:49)
%0.215*exp(-0.011*(wl-500) + 0.02 channel
