
figure('color', 'w');

L=0.2;
f_0=3.517257185*10^8;% [MHz]
m=(132.90545/6.02e23)/1000; % [kg]
c=299792458; % [m/s]
k=1.3806488*10^(-23); % [m^2*kg/s^2]
hbar=1.05457173*10^(-34); % [m^2*kg/s]
dip=3.8014*10^(-29); % [C*m]
Av=6.02*10^23;
gamma_N=2*pi*5.234; % [MHz]
e_psi=8.854187817*10^(-12); % [F/m]
f_1=3.517257185*10^8;%


h(1)=subplot(2,1,1)


% Create data
 load('fit_data\smallLIAD150803VarGammaNEW.mat')
FIT3=FIT2;
sizeOf=size(FIT3,2);
t=linspace(0,10800,sizeOf);
gamma_D=FIT3(2,:);
gamma_N=5.234+1.16e-19*(FIT3(1,:).*FIT3(2,:)*2*c*hbar*e_psi)/(L*dip^2*2*pi*f_1*(pi*log(2))^0.5);
density=(FIT3(1,:).*FIT3(2,:)*2*c*hbar*e_psi)/(L*dip^2*2*pi*f_1*(pi*log(2))^0.5)*1e-6;
d=FIT3(1,:).*FIT3(2,:)./(gamma_N*(pi*log(2))^0.5);

for i=1070:length(d)-1000 %% correct peaks due to weird width fit
    if d(i)>1.1*d(i-1)||d(i)<0.9*d(i-1)
        d(i)=d(i-1);
    end
    
end

% Place axes at (0.1,0.1) with width and height of 0.8
%handaxes1 = axes('Position', [0.12 0.12 0.8 0.8]);

% Main plot
plot(t,smooth(d,1));
xlabel('Time (s)')
ylabel('OD_{experiment}') % left y-axis
%set(handaxes1, 'Box', 'off')
xlim([0 10800]);
 ylim([min(d) max(d)])

% Adjust XY label font
handxlabel1 = get(gca, 'XLabel');
set(handxlabel1, 'FontSize', 16, 'FontWeight', 'bold')
handylabel1 = get(gca, 'ylabel');
set(handylabel1, 'FontSize', 16, 'FontWeight', 'bold')

% Place second set of axes on same plot
handaxes2 = axes('Position', [0.69 0.66 0.2 0.2]);
plot(t-10800, smooth(d,1), 'r')
set(handaxes2, 'Box', 'off')
%set(handaxes2,'YTickLabel',[]);
% set(handaxes2,'XTickLabel',[]);
 xlabel('t (s)')
 xlim([-150 0]);

% Adjust XY label font
set(get(handaxes2, 'XLabel'), 'FontName', 'Times')
set(get(handaxes2, 'YLabel'), 'FontName', 'Times')

% Add another set of axes
handaxes3 = axes('Position', [0.22 0.66 0.2 0.2]);
plot(t, smooth(d,1), 'm')
set(handaxes3, 'Box', 'off')
%set(handaxes3,'YTickLabel',[]);
% set(handaxes3,'XTickLabel',[]);
 xlabel('t (s)')
xlim([0 30]);




h(2)=subplot(2,1,2)
% Generate data

T_total=10800;
Tsteps=sizeOf;


t=linspace(0,T_total,Tsteps);

TLIAD_switch=FIT3(6,:)-FIT3(6,1);

% TLIAD_switch=zeros(1,Tsteps);
% for i=1:Tsteps
%     if t(i)>=1 && t(i)<=10710
%         TLIAD_switch(i)=1;
%     else
%         TLIAD_switch(i)=0;
%     end
%     
% end


% N_tot=1.3e6;
% TLIAD_mag=0.1;
% eq_vapour=0.0001;
% eq_cluster=0.4;
% eq_reservoir=10000;
% R_vapour=0.15;
% R_cluster=0;%1e-8;
% R_reservoir=80e-8;
% R_desorb_layer=0.1;
% R_desorb_cluster=6;%10;
N_tot=2.3e6;
TLIAD_mag=0.1;
eq_vapour=0.0001;
eq_cluster=0.4;
eq_reservoir=10000;
R_vapour=0.15;
R_cluster=0;%1e-8;
R_reservoir=16e-7;
R_desorb_layer=0.06;
R_desorb_cluster=6;%10;


init=[N_tot/(1+eq_vapour+eq_cluster+eq_vapour*eq_reservoir),N_tot*eq_cluster/(1+eq_vapour+eq_cluster+eq_vapour*eq_reservoir),N_tot*eq_vapour/(1+eq_vapour+eq_cluster+eq_vapour*eq_reservoir),N_tot*eq_reservoir*eq_vapour/(1+eq_vapour+eq_cluster+eq_vapour*eq_reservoir)];
% [N_layer,N_cluster,N_vapour,N_reservoir]

M=zeros(length(init),length(init));
R=zeros(length(init),length(t));
LIAD=TLIAD_mag*TLIAD_switch;
R(:,1)=init';

for i=2:length(t)
    R_desorb_l=R_desorb_layer*LIAD(i-1);
    R_desorb_c=R_desorb_cluster*LIAD(i-1);
    M=[-R_vapour*eq_vapour-R_cluster*eq_cluster-R_desorb_l,R_cluster,R_vapour,0;
        R_cluster*eq_cluster,-R_cluster-R_desorb_c,0,0;
        R_vapour*eq_vapour+R_desorb_l,R_desorb_c,-R_vapour-R_reservoir*eq_reservoir,R_reservoir;
        0,0,R_reservoir*eq_reservoir,-R_reservoir];
    R(:,i)=expm(M*T_total/Tsteps)*R(:,i-1);   
end


% Main plot
plot(t,R(3,:));
xlabel('Time (s)')
ylabel('OD_{model}') % left y-axis
%set(handaxes1, 'Box', 'off')
xlim([0 10800]);
 ylim([R(3,end) max(R(3,:))])
 
% Adjust XY label font
handxlabel1 = get(gca, 'XLabel');
set(handxlabel1, 'FontSize', 16, 'FontWeight', 'bold')
handylabel1 = get(gca, 'ylabel');
set(handylabel1, 'FontSize', 16, 'FontWeight', 'bold')

% Place second set of axes on same plot
handaxes2 = axes('Position', [0.69 0.25 0.2 0.2]);
plot(t-10800, smooth(R(3,:),1), 'r')
set(handaxes2, 'Box', 'off')
%set(handaxes2,'YTickLabel',[]);
% set(handaxes2,'XTickLabel',[]);
 xlabel('t (s)')
 xlim([-150 0]);
%  ylim([R(3,end) R(3,sizeOf-1000)]);

% Adjust XY label font
set(get(handaxes2, 'XLabel'), 'FontName', 'Times')
set(get(handaxes2, 'YLabel'), 'FontName', 'Times')

% Add another set of axes
handaxes3 = axes('Position', [0.22 0.25 0.2 0.2]);
plot(t, smooth(R(3,:),1), 'm')
set(handaxes3, 'Box', 'off')
%set(handaxes3,'YTickLabel',[]);
% set(handaxes3,'XTickLabel',[]);
 xlabel('t (s)')
xlim([0 30]);



set(h(1),'xticklabel',[]);
pos=get(h,'position');
bottom=pos{2}(2);
top=pos{1}(2)+pos{1}(4);
plotspace=top-bottom;
pos{2}(4)=plotspace/2;
pos{1}(4)=plotspace/2;
pos{1}(2)=bottom+plotspace/2;

set(h(1),'position',pos{1});
set(h(2),'position',pos{2});

linkaxes(h,'x');

export_fig('file_name','-pdf','-q50')