
clear all
close all
%setting options for fmisearch
opts = optimset('fminsearch');      
opts = optimset(opts,'MaxIter',4000); 
opts = optimset(opts,'MaxFunEvals',2000);   % usually 100*number of params
opts = optimset(opts,'TolFun',1e-9);

fit1=NaN*ones(154,1);
fit2=NaN*ones(154,1);
fit3=NaN*ones(154,1);

%make a large matrix with all CNR data
for i=[100:143,151:154]%1:99%80:88,90:91%50:58,60:91%2:58,60:91 %21:49, 60:62;%2:20, 63:91;%
    path='/Volumes/Data2/TaraEuropa/Partner_data/CNR-P/Absorption_spectra_2024/';
    filename=[path,'ST_ag',num2str(i),'.mat'];
    load(filename)
    SR=find(wl>=350 & wl<=500)
    wl_ag=wl(SR)';
    B=ag_ave(SR);
    C=ag_std(SR);

    %opts = optimset('LevenbergMarquardt','on');
    %guess for paramters (amplitude at 532 and slope)
    J=find(wl_ag-400 == min(wl_ag-400));
    x0=[B(J), 0.017];
    %minimization routine
    x1 = fminsearch(@exp_fit,x0,opts,B,wl_ag,C)
    x0=[B(J), 0.017,1];
    x2 = fminsearch(@stretched_exp_fit,x0,opts,B,wl_ag,C)
    x0=[B(J), 0.017,0];
    x3 = fminsearch(@exp_fit_offset,x0,opts,B,wl_ag,C)
    %plot data and fit
    plot(wl_ag, B,'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')
    %pause
    fit1(i,1:2)=x1;
    fit2(i,1:3)=x2;
    fit3(i,1:3)=x3;
    clear B C
end

plot([100:154],fit1(100:154,2),'LineWidth',2)
xlabel('station')
ylabel('spectral slope')
fontsize(gcf,16,"points")

% M = readtable("/Volumes/Data2/TaraEuropa/Partner_data/CNR-P/Station_lat_lon_2023.csv") 
% load coastlines
% plot(coastlon,coastlat,'-')
% hold on
% plot(M.Var3,M.Var2,'.')
% for i=[100:143,151:154]
%     text(M.Var3(i),M.Var2(i),num2str(M.Var1(i)),FontSize=16)
% end
% hold off
% axis([min(M.Var3)-0.5,max(M.Var3)+0.5,min(M.Var2)-0.5,max(M.Var2)+0.5 ])
   

function y = exp_fit(x0,spec,l,unc)
y=sum(((spec-x0(1).*exp(-x0(2).*(l-400)))./unc).^2);
end

function y = exp_fit_offset(x0,spec,l,unc)
y=sum(((spec-(x0(3)+x0(1).*exp(-x0(2).*(l-400))))./unc).^2);
end

function y = stretched_exp_fit(x0,spec,l,unc);
y=sum(((spec-x0(1).*exp(-x0(2).*(l-400).^x0(3)))./unc).^2);
end

