clear
% close all

%% load data
dat0 = load([pwd,'\NEP_0_sesn2mV_5000msTC_40acgain_4umstep_v1']);
dat10 = load([pwd,'\NEP_10_sesn500uV_5000msTC_40acgain_4umstep_v1']);
dat20 = load([pwd,'\NEP_20_sesn200uV_5000msTC_40acgain_4umstep_v1']);
dat30 = load([pwd,'\NEP_30_sesn20uV_5000msTC_40acgain_4umstep_v1']);
dat40 = load([pwd,'\NEP_40_sesn10uV_5000msTC_40acgain_4umstep_v1']);
dat50 = load([pwd,'\NEP_50_sesn5uV_5000msTC_40acgain_4umstep_v1']);
dat60 = load([pwd,'\NEP_60_sesn2uV_5000msTC_40acgain_4umstep_v1']);
dat70 = load([pwd,'\NEP_70_sesn1uV_5000msTC_40acgain_4umstep_v1']);
dat80 = load([pwd,'\NEP_80_sesn1uV_10000msTC_40acgain_4umstep_v1']);
dat76 = load([pwd,'\NEP_76_sesn1uV_10000msTC_40acgain_4umstep_v1']);
dat90 = load([pwd,'\NEP_90_sesn1uV_10000msTC_40acgain_4umstep_v1']);
dat100 = load([pwd,'\NEP_blocked_sesn1uV_10000msTC_40acgain_4umstep_v1']);

datall = [dat0 dat10 dat20 dat30 dat40 dat50 dat60 dat70 dat76 dat80 dat90];
attendb = [0 10 20 30 40 50 60 70 76 80 100];

fielddb = attendb/2;
%% initial vals
f0 = 3.4e12;
lambda = 3e8/f0;
C = 0.5;
alpha = 0;
L0 = 0.48;
V0 = 0;
dL = 0;

%% fit SM fringes
for x = 1:size(datall,2)/2
    % load current data
    clear curdat
    curdat = datall(:,x*2-1:x*2);
    % scale to um and mV (AC coupled)
    curdat(:,1) = (curdat(:,1)-42.2)*1000;
    curdat(:,2) = (curdat(:,2)-mean(curdat(:,2)))*1000;
    % calculate initial params
    beta = (max(curdat(:,2))-min(curdat(:,2)))/2;
    data = [curdat(:,1)*1e-6 curdat(:,2)];
    curfft = abs(fft(curdat(:,2)));
    fftamp(x) = max(curfft(4:6))/2;
    % set C and alpha for second data onwards
    if x ~= 1
        C = Cstore(x);
        alpha = alphastore(x);
    end
    %% perform intial dL & lambda fits
    params = dL;
    vars = [beta V0 C alpha L0 lambda];
    options = optimset('TolX',1e-6,'MaxFunEvals',100,'MaxIter',100);
    dL = fminsearch('SM_fit_dL',params,options,data,vars);

    params = lambda;
    vars = [beta V0 C alpha L0 dL];
    options = optimset('TolX',1e-6,'MaxFunEvals',100,'MaxIter',100);
    lambda = fminsearch('SM_fit_lambda',params,options,data,vars);
    
    %% Perform full fit
    if x==1
        alpha = -1.92;
        params = [C alpha dL beta];
        vars = [V0 L0 lambda];
        options = optimset('TolX',1e-7,'MaxFunEvals',10000,'MaxIter',10000);
        [params,~] = fminsearch('SM_fit_4params',params,options,data,vars);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Store fitted params
        betastore(1) = params(4);
        Cstore(1) = params(1);
        alphastore(1) = params(2);
        dLstore(1) = params(3);
        V0store(1) = V0;
        lambdastore(1) = lambda;
        % calculate C for other attenuations and store alpha
        for i = 2:size(datall,2)/2
            Cstore(i) = Cstore(1)*10^-(fielddb(i)/10);
            alphastore(i) = alphastore(1);
        end
    else
        %% perform fitting for second dataset onwards using C and alpha calculated from first
        params = [dL beta];
        vars = [V0 L0 lambda alpha C];
        options = optimset('TolX',1e-7,'MaxFunEvals',10000,'MaxIter',10000);
        [params,~] = fminsearch('SM_fit_dL_beta',params,options,data,vars);
        
        betastore(x) = params(2);
        dLstore(x) = params(1);
        V0store(x) = V0;
        lambdastore(x) = lambda;
    end
end
%% correct individual values (see Fig3.m code)
dLstore(1) = 5.0114e-4;
betastore(1) = 0.5511;
Cstore(1) = 0.3927;
alphastore(1) = -1.8457;
dLstore(10) = 5.22e-4;
lambdastore(10) = 8.8373e-5;

%% plot fig 2
[plotstack2] = createStackPlot(3);
i=1;
for x = [1,6,10]
    curdat = datall(:,x*2-1:x*2);
    curdat(:,1) = (curdat(:,1)-42.2)*1000;
    curdat(:,2) = (curdat(:,2)-mean(curdat(:,2)))*1000;
    % calculate Vsm
    V_fit = zeros(1,length(curdat))';
    for loop = 1:length(curdat)
        z_inst = curdat(loop,1)*1e-6;
        L = L0 + dLstore(x) + z_inst;
        V_fit(loop) = selmixpower(betastore(x),Cstore(x),L,lambdastore(x),alphastore(x));
        V_fit(loop) = V_fit(loop) + V0store(x);
    end
    % scale 80dB data to uV
    if x == 10
        curdat(:,2) = curdat(:,2)*1000;
        V_fit = V_fit*1000;
    end
    % plot and set y labels
    plot(plotstack2(i),curdat(:,1),curdat(:,2),'bo','linewidth',1.2);
    hold(plotstack2(i),'on');
    plot(plotstack2(i),curdat(:,1),V_fit,'r','linewidth',1.2);
    
    if x~=10
        set(get(plotstack2(i),'ylabel'),'string','{\it{V}}_{SM} (mV)');
    else
        set(get(plotstack2(i),'ylabel'),'string','{\it{V}}_{SM} (\muV)');
    end
    set(plotstack2(i),'LineWidth',1);
    i=i+1;
end
% set figure styles
fsize = 14;
set(plotstack2(1:end-1),'xticklabel','');
set(get(plotstack2(end),'xlabel'),'string','Cavity Extension (\mum)','fontsize',fsize);
set(plotstack2(1),'ylim',[-1.2 1.2],'fontsize',fsize);
set(plotstack2(2),'ylim',[-13e-3 13e-3],'ytick',-10e-3:10e-3:10e-3,'fontsize',fsize);
set(plotstack2(3),'ylim',[-3e-1 3e-1],'fontsize',fsize);
axes(plotstack2(3));
for o = 1:3
    ylab = get(plotstack2(o),'ylabel');
    ylabpos(1:3,o) = get(ylab,'position');
end
for o = 1:3
    ylab = get(plotstack2(o),'ylabel');
    set(ylab,'position',[min(ylabpos(1,:)),0,-1]);
    set(get(plotstack2(o),'ylabel'),'fontsize',fsize);
end
t2 = text(3,14e-1,'0dB Attenuation','fontsize',fsize);
t3 = text(3,8e-1,'50dB Attenuation','fontsize',fsize);
t4 = text(3,2e-1,'80dB Attenuation','fontsize',fsize);
axes(plotstack2(3));
for x = 1:3
    pos = get(plotstack2(x),'position');
    set(plotstack2(x),'position',[pos(1)+0.03,pos(2)+0.05,pos(3:4)]);
end
clear pos
pos = get(t2,'position');
set(t2,'position',[pos(1)+0.03,pos(2)+0.05]);
pos = get(t3,'position');
set(t3,'position',[pos(1)+0.03,pos(2)+0.05]);
pos = get(t4,'position');
set(t4,'position',[pos(1)+0.03,pos(2)+0.05]);

%% plot fig 3
atten = sqrt(10.^(-attendb/10));
fftblocked = fftamp(end);
betastore = betastore(1:end-1);
atten = atten(1:end-1);

% % %% convet to log for fit
attenlog2 = log10(atten);
atten2 = sqrt(10.^(-attendb/10));
attenlogfit = log10(atten2);
betalog2 = log10(betastore/1000);
cf = fit(attenlog2',betalog2','poly1');
ydat = cf(attenlogfit');
fully = 10.^ydat;
betastore(1:4) = [0.5511,0.3505,0.1718,0.0437]; %% from other code

figure
loglog(atten.^2,betastore/1000,'ko','linewidth',1.2);
hold on
plot(atten2.^2,fully,'b','linewidth',1.2);
set(gca,'linewidth',1);
plot([1 1e-10],[fftblocked/1000 fftblocked/1000],'r-','linewidth',1.2); %% double check
set(gca,'ylim',[10^-8 10^-2],'xlim',[10^-10 10^0],'fontsize',14,'ytick',[10^-8, 10^-7, 10^-6, 10^-5, 10^-4, 10^-3, 10^-2],'xtick',[10^-10, 10^-8, 10^-6, 10^-4, 10^-2, 10^0]);
set(get(gca,'ylabel'),'string','{\it{V}}_{PkPk} (V)','fontsize',14);
set(get(gca,'xlabel'),'string','Power Attenuation {\it{\epsilon}}_a','fontsize',14);
% SaveForThesis(gcf,pwd,'FieldAtten');