close all
clear all
load('perims_msemean.mat')
%load('GATE_cloud_geometry_all_times.mat');
%load('GATE_Cloud_Geo_it23_yesHoles_yesChild');
%load('GATE_Cloud_Geo_it23_noHoles_yesChild');
%load('GATE_Cloud_Geo_it23_noHoles_noChild');
%error('yo')
%load('PowerLawGategeometryTall.mat')
%phmean: h data
%phz: height data
%perims: perimeter data

%%%% Renaming for new file %%%
%phmean = perim_mse_avg;
%phstd = perim_mse_std;
%phz = z(slice_level);

fontname = 'Helvetica';
set(0,'defaultaxesfontname',fontname);
fontsize = 12;
set(0,'defaultaxesfontsize',fontsize);

corrfac = 1.004;
phmean = phmean*corrfac; %kJ/kg
phstd = phstd*corrfac; %kJ/kg

hBinWidth = 2;
hbins = 331.5 : hBinWidth : 343.5;
stdthresh = 1; % kJ/kg
g = 9.8; % m/s2

nh = length(hbins) -1;
np = length(pbins) -1;
hbar = hbins(1:nh) + diff(hbins)/2;
pbinskm = pbins/1e3;%Perimeter bins in km
pkmbar = pbinskm(1:np) + diff(pbinskm)/2; % average perimeter in each bin
pkmbarlog = exp(log(pbinskm(1:np)) + diff(log(pbinskm))/2); % logarithmic average perimeter in each bin
perimskm = perims/1e3;%Perimeters in km

dP = diff(pbinskm);
dlogP = diff(log(pbinskm)); %Bin widths should be identical in log space
PbinMin = exp(log(pbinskm) - dlogP(1)/2); %Min perimeter bin
PbinMax = exp(log(pbinskm) + dlogP(1)/2); %Max perimeter bin
PbinWidth = PbinMax - PbinMin; %Perimeter bin width in linear space
Areakm = 204.8*204.8;%Model domain is 200 km each side. 
Heightkm = quantile(phz,0.99)/1e3; % Height in km below which 99% of the cloud perimeters are contained
Volumekm = Areakm*Heightkm;%Volume of the cloudy domain
densitykm = 0.63e9; %average density estimate for perturbations within 2K (5K is 0.7) kg/km3
mass = densitykm*Volumekm; %mass of domain for 2K

DensityFactor = (Volumekm); % To convert frequency distriubtions to #/km3
dPmat = repmat(diff(pbinskm),nh,1);


ModelScale = 100; %m
FractalDimension = 4/3; %From Lovejoy
Epsilon = 2.2e-4; %m3/s2 from model TKE dissipation for 12h to 24 h
Diffusivity = 0.211e-4*(263./273).^1.94.*(1e5./7e4); %m2/s molecular diffusivity estimate
Diffusivitykm = Diffusivity*1e-6;
Viscosity = 0.174e-4*(263./273).^1.94.*(1e5./7e4); %m2/s molecular diffusivity estimate
KolmogorovScale = (Viscosity^3/Epsilon)^(1/4); %m
%TRY SCALING EPSILON BY CLOUD AREA
KolmogorovScale = 1.1e-3; %m Estimated from GigaLES_stats
Jcharsatup = 0.02*1e6; %kg/km2/s Estimated from GigaLES_stats. good = find(saturatedhPert < 3); mean(saturatedMassFlux(good))
Jcharcld = 0.015*1e6; %kg/km2/s Estimated from GigaLES_stats. good = find(saturatedhPert < 3); mean(cldMassFlux(good))


TurbDiffkm = Diffusivitykm*(ModelScale/KolmogorovScale)^FractalDimension; %Krueger paper
EddyDiffModkm = 0.1207e-6; %From Model



%Adjust for different bin widths

low = find(phz < 1200);
lowmid = find(phz < 2400 & phz >= 1200);
highmid = find(phz < 5200 & phz >= 2400);
high = find(phz >= 5200);

%Randomly sample from oversampled data at low and mid altitudes to match
%high altitudes

lowfactor = 50/100;
lowmidfactor = mean(gradient(z(find(z < 2400 & z >= 1200))))/100;
highmidfactor = mean(gradient(z(find(z < 5200 & z >= 2400))))/100;

[perimskmlow, Ilow] = datasample(perimskm(low),round(length(low)*lowfactor));
[perimskmlowmid, Ilowmid] = datasample(perimskm(lowmid),round(length(lowmid)*lowmidfactor));
[perimskmhighmid, Ihighmid] = datasample(perimskm(highmid),round(length(highmid)*highmidfactor));
perimskmhigh = perimskm(high);

phmeanlow = phmean(Ilow);
phmeanlowmid = phmean(Ilowmid);
phmeanhighmid = phmean(Ihighmid);
phmeanhigh = phmean(high);

phzlow = phz(Ilow);
phzlowmid = phz(Ilowmid);
phzhighmid = phz(Ihighmid);
phzhigh = phz(high);

phstdlow = phstd(Ilow);
phstdlowmid = phstd(Ilowmid);
phstdhighmid = phstd(Ihighmid);
phstdhigh = phstd(high);

perimskm = [perimskmlow, perimskmlowmid, perimskmhighmid, perimskmhigh];
phmean = [phmeanlow, phmeanlowmid, phmeanhighmid phmeanhigh];
phzkm = [phzlow, phzlowmid, phzhighmid, phzhigh]/1000;
phstd = [phstdlow, phstdlowmid, phstdhighmid, phstdhigh];

phmeanz = phmean(find(phstd < stdthresh)); %vertical layers h
perimskmz = perimskm(find(phstd < stdthresh)); %vertical layers perimeters
phzz = phz(find(phstd < stdthresh)); % vertical layers heights


X = [phmeanz', perimskmz']; %h and perims matrix. Converted from K to kJ/kg

EDGES{1} = hbins;
EDGES{2} = pbinskm; %perimeter bins in km

fgrid = hist3(X, 'Edges', EDGES);%Frequency distribution in space of T and Perims.
fgrid = fgrid(1:nh,1:np); % For some reason fgrid is has a dimension
                                 % that is not one less than EDGES!
pfgrid = fgrid.*repmat(pkmbarlog,nh,1); %N*P frequency distribution at each T and P (km)

totnum_h = sum(fgrid,2)'; %Total number at each T
totnum_P = sum(fgrid,1); %Total number at each P
totP_hP  = pfgrid; %N*P frequency distribution at each T and P (km)
totP_h = sum(pfgrid,2)'; %Total perimeter at each T (km)
totP_P = sum(pfgrid,1)'; %Total perimeter at each P (km)
totP = sum(totP_h)'; %Total perimeter (km). Same as sum(totP_P)
totnum = sum(totnum_h); %Total number 

Number = fgrid; %Frequency distribution versus perimeter and temperature 
NumberAll_P = totnum_P; %Frequency distribution versus perimeter  
NumberAll_h = totnum_h; %Frequency distribution versus temperature  
TotalPerimeter_hP = totP_hP; %Total perimeter versus temperature and perimeter (km)
TotalPerimeter_h = totP_h; %Total perimeter versus temperature (km)
TotalPerimeter_P = totP_P; %Total perimeter versus perimeter (km)
TotalPerimeter = sum(TotalPerimeter_h'); % Total perimeter km
TotalNumber = sum(NumberAll_P'); %Total number
Probability = Number/totnum;

meanPerimeter = TotalPerimeter./TotalNumber;
%inbins = find(phmeanz >= min(hbins) & phmeanz <= max(hbins)); 
inbins = find(phmeanz >= 334.5 & phmeanz <= 340.5); 
meanh = sum(phmeanz(inbins))/length(phmeanz(inbins));
%h deviation
hdeviation = abs(hbar - meanh);
%Find variance in h
varh = sum((phmeanz(inbins) - meanh).^2)/length(phmeanz(inbins));
%varh = sum(perimskmz(inbins).*(phmeanz(inbins) - meanh).^2)/sum(perimskmz(inbins));
stdh = sqrt(varh);
%meanhdeviation = sum(TotalPerimeter_h.*hdeviation)./sum(TotalPerimeter_h);
meanhdeviation = sum(abs(phmeanz(inbins) - meanh))/length(phmeanz(inbins));

scaleh = find(phmeanz >= meanh - stdh & phmeanz <= meanh + stdh);


%Find maximum perimeter
maxPerimRat = 1/exp(1); %Ratio where max perimeter is defined
meanFirstFew = mean(TotalPerimeter_P(5:10));
for i = 1:nh
%Take mean of first few bins
    meanFirstFew_h(i) = mean(TotalPerimeter_hP(i,5:10));
    for j = 5:np
        perimeterRat_h = TotalPerimeter_hP(i,j)./meanFirstFew_h(i);
        perimeterRat = TotalPerimeter_P(j)./meanFirstFew;
        if perimeterRat_h < maxPerimRat
            maxPerimeter_h(i) = pkmbarlog(j);
            break
        end
    end
    
    for j = 5:np
        perimeterRat = TotalPerimeter_P(j)./meanFirstFew;
        if perimeterRat < maxPerimRat
            maxPerimeter = pkmbarlog(j);
            break
        end
    end
end

%initial number in power law
Number0_h = Number(:,5)'; %
Number0 = sum(Number0_h)*(1-exp(-hBinWidth/(stdh)));
%initial perimeter in exponential
TotalPerimeter0 = TotalPerimeter*(1-exp(-hBinWidth/(stdh)));
maxPerimeter0 = maxPerimeter*(1-exp(-hBinWidth/(stdh)));

%total mass Flux estimate
Jh = TurbDiffkm*densitykm*(hdeviation/2)./meanh.*TotalPerimeter0.*exp(-(hdeviation/2)./stdh); %kg/s
%mass flux versus delta h estimate
Jtot = TurbDiffkm*densitykm/meanh*stdh*TotalPerimeter; %kg/s
entrainmentspeed = Jtot/mass*2000;% m/s



%Average perimeter with respect to the temperature deviation
MeanPerimeter_h = sum(TotalPerimeter_h.*hdeviation)./sum(hdeviation);

%Average temperature deviation and variance for the perimeter
MeanhDeviation = sum(TotalPerimeter_h.*hdeviation)./sum(TotalPerimeter_h);
MeanhVariance = sum(TotalPerimeter_h.*hdeviation.^2)./sum(TotalPerimeter_h);

NumberDensity = fgrid./DensityFactor; %Frequency distribution versus perimeter and temperature 
NumberDensityAll = totnum_P./DensityFactor; %Frequency distribution versus perimeter  
PerimeterDensity = totP_h/DensityFactor; %Total perimeter density versus temperature (km/km^3)
TotalPerimeterDensity = sum(PerimeterDensity'); % Total perimeter density km/km3
TotalNumberDensity = sum(NumberDensityAll'); %Total number of /km3
MeanPerimeterDensity_h = MeanPerimeter_h/DensityFactor; %Mean Perimeter Density 

TimeScale_h = 1./(PerimeterDensity.*TurbDiffkm); %Estimate of dissipation timescale
TimeScale = 1./(TotalPerimeterDensity.*TurbDiffkm); %Estimate of dissipation timescale

TimeScale_h_Eddy = 1./(PerimeterDensity.*EddyDiffModkm); %Estimate of dissipation timescale
TimeScale_Eddy = 1./(TotalPerimeterDensity.*EddyDiffModkm); %Estimate of dissipation timescale

%Area estimate calculations
coeff = 1;
AreaDensity = coeff.*PerimeterDensity.^(2/3);
TotalAreaDensity = sum(AreaDensity');

%%%Region that satisfies power law
exids = [1:4,42:49];
exlocs = true(size(pkmbarlog)); exlocs(exids) = false; exlocs = exlocs & NumberDensityAll~=0;


%%% Calculate total as logarithm
NumberAll_P0 = NumberAll_P(5).*pkmbarlog(5)./pbinskm(5);
logNumberAll = log(totnum_P); %log Frequency distribution versus perimeter  
dlogNumberAll = gradient(log(totnum_P)); %dlog Frequency distribution versus perimeter  
dlogpkmbar = log(PbinMin(1) + PbinWidth(1)) - log(PbinMin(1)); %dlog Perimeter
logDist = dlogNumberAll./dlogpkmbar; % dlnn/dlnP
TotalPerimeter_log = sum(logDist(exlocs(2:end)).*NumberAll_P(exlocs(2:end)).*pkmbarlog(exlocs(2:end)).*dlogpkmbar);

%%%%%Plotting!!!
%load perims_spec_heights.mat

%gf = sum(sum(global_freqsdry8,3),2); %Sums over all times and temperatures for number distribution
%pbinskm = [pbinskm(:)]; %Perimeter bins

% plot perims at just one time
% include a -1 slope example
% fit line for 330 and 340
fig1 = figure('Color',[1,1,1],'Position',[50 50 700 550],'PaperPositionMode', 'auto','PaperType','<custom>','PaperSize',[9.5 11]);
%fig1 = figure('Color',[1,1,1]);
axes1 = axes('Parent',fig1,'YScale','log','YMinorTick','on','YMinorGrid','off','YGrid','on',...
    'XScale','log','XMinorTick','on','XMinorGrid','off','XGrid','on',...
    'FontWeight','bold',...
    'FontSize',14);
box(axes1,'on');
hold(axes1,'all');
%ylim(axes1,[1e-6 1.5]);
%xlim(axes1,[1e-1 2e4]);

% calculate fit line for all temperature values
NumberDensityAll = NumberDensityAll;
% define exclueded locations
[Pfit, Pstats] = polyfit(log10(pkmbarlog(exlocs)), log10(NumberDensityAll(exlocs)), 1);
seP = 2*sqrt(sum(inv(Pstats.R).^2,2))*Pstats.normr/sqrt(Pstats.df); %95% confidence bounds
yall = polyval(Pfit, log10(pkmbarlog(exlocs)));
plot(pkmbarlog(exlocs), 10.^yall, '-r','LineWidth',2)

%cmap = jet(length(Tbar));
cmap = parula(length(hbar));

plot(pkmbarlog, NumberDensityAll,'.','LineWidth',2,'Color',[0 0 0],'MarkerSize',20)

for ip =1:length(hbar)
    plot(pkmbarlog, NumberDensity(ip,:),'o','LineWidth',2,'Color',cmap(ip,:))
end

% plot a -1 line that goes through 1e7 for the first bin
y1p = -1*log10(pkmbarlog(exlocs)) - 3.2;
y1 = 10.^y1p;
plot(pkmbarlog(exlocs), y1, '--k','LineWidth',2,'DisplayName','Slope -2')
%ylab = strcat('\bf{Perimeter number density }', '$\bf{n}$',' ${\rm({km}^{-3}})$');
%xlab = strcat('\bf{Perimeter}  ', '$\bf{\,\,\Lambda}$',' $\rm{(km)}$');
ylab = strcat('\bf{Number density in logarithmically spaced bins n/V}',' {\rm({km}^{-3}})');
xlab = strcat('\bf{Perimeter at 100 m resolution}  ', '\bf{ \lambda}',' \rm{(km)}');

k = xlabel(xlab);
l = ylabel(ylab);

%set(k,'Interpreter','Latex','FontSize',16,'Color',cmap(1,:));
%set(l,'Interpreter','Latex','FontSize',16,'Color',cmap(1,:));
%set(k,'FontSize',16,'Color',cmap(1,:));
%set(l,'FontSize',16,'Color',cmap(1,:));
set(k,'FontSize',16,'Color',[0.1 0.1 0.8]);
set(l,'FontSize',16,'Color',[0.1 0.1 0.8]);

c1 = strcat('Simulation slope', '$ = $', num2str(round(100*Pfit(1))/100),...
    '$ \pm $', num2str(round(100*seP(1))/100));
t1 = text(0.6,4,c1);

set(t1,'FontWeight','bold','FontSize',18,'Color','red');

set(t1,'Interpreter','Latex')
%text(2,2,['Slope = ' sprintf('%5.2f',Pall(1)) ...
%    ],'FontWeight','bold','FontSize',14);

t2 = text(0.5,4e-6,['{Theoretical slope = -1}'  ...
    ],'FontWeight','bold','FontSize',18);

set(t2,'Interpreter','Latex')


set(axes1,'Clim',[min(hbins) max(hbins)]);
colormap(cmap)
cb = colorbar('peer',axes1,'FontWeight','bold','FontSize',14,'Ticks',hbins);
lcb = ylabel(cb,'$\rm{h^\star\;\left(kJ\,kg^{-1}\right)}$','FontWeight','bold','FontSize',14)
set([lcb],'FontWeight','bold','FontSize',18);
set([lcb],'Interpreter','Latex')

%%%% Create histogram
axes2 = axes('position',[0.548 0.65 0.23 0.23],'Parent',fig1,'YScale','lin','YMinorTick','off',...
    'YMinorGrid','off','YGrid','on','YTickLabel','','XTickLabel','',...
    'XScale','lin','XMinorTick','off','XMinorGrid','off','XGrid','on',...
    'FontWeight','bold',...
    'FontSize',12);


for i = 30:35
    rectangle('position',[PbinMin(i) 0  PbinWidth(i) NumberDensityAll(i)])
end
%axis([100 2000 0 2e-7])
 

% Create rectangle
annotation(fig1,'rectangle',...
    [0.538571428571429 0.538181818181818 0.0871428571428571 0.0890909090909091]);

% Create line
annotation(fig1,'line',[0.624285714285714 0.781428571428571],...
    [0.536363636363636 0.649090909090909]);

% Create line
annotation(fig1,'line',[0.54 0.547142857142857],...
    [0.628090909090909 0.649090909090909]);

% Create textbox
annotation(fig1,'textbox',...
    [0.618571428571429 0.814030870768575 0.23 0.0571428571428572],...
    'String',{'n\Delta\lambda is scale invariant'},...
    'LineStyle','none',...
    'FontSize',14);
% Create textbox
annotation(fig1,'textbox',...
    [0.618571428571429 0.734030870768575 0.23 0.0571428571428572],...
    'String',{'inset in a linear scale'},...
    'LineStyle','none',...
    'FontSize',14);
%    [0.618571428571429 0.814030870768575 0.135152664349352 0.0571428571428572],...

fname=['h_perims_all'];
%cd ../export_fig
%eval(['export_fig ' fname ]);
%print -dpdf -loose ./h_perims_all
export_fig(sprintf(fname), '-pdf', '-transparent', '-nocrop');
%export_fig(sprintf(fname), '-pdf');

%cd ../tim_perimeters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


fig2 = figure('Color',[1,1,1],'Position',[50 50 700 550],'PaperPositionMode', 'auto','PaperType','<custom>','PaperSize',[9.5 11]);
color_spectrum = 'parula';

colorgrid = eval([color_spectrum, '(nh)']);
%%%% put temperature deviation on top
axes7 = axes('Position', [0.13         0.13        0.611        0.75],'color','none');
%box(axes7,'off');
hold(axes7,'all');

for ip = 1 : nh
    
    plot(hdeviation(ip), TotalPerimeter_h(ip)/TotalPerimeter0, '.', 'color', colorgrid(ip,:),...
        'MarkerSize',40,'Parent',axes7) 
    
end

l1 = xlabel('$$\rm{Deviation\;from\;mean\;}\delta h^\star\;\left({kJ\,kg^{-1}}\right)$$');
l2 = ylabel('$$\Lambda_i/\Lambda_0$$');
set([l1],'FontWeight','bold','FontSize',22);
set([l1],'Interpreter','Latex')
set([l2],'FontWeight','bold','FontSize',22);
set([l2],'Interpreter','Latex')
% Create ylabel
%ylabel('Perimeter density (km^{-2})');

%box(axes6,'on');
% Set the remaining axes properties
set(axes7,'FontSize',14,'FontWeight','bold','LineStyleOrderIndex',3,...
    'XMinorTick','on','XScale','lin','YMinorTick','on','YScale','log','XLim',[0 6],'YLim',[0.0091 3.1]);
% Create legend
%leg1 = legend(axes6,'show');
%set(leg1, 'Position',[0.85 0.17 0.12 0.78]);

[Pfit, Pstats] = polyfit(hdeviation,log(TotalPerimeter_h/TotalPerimeter0),1);
seP = 2*sqrt(sum(inv(Pstats.R).^2,2))*Pstats.normr/sqrt(Pstats.df); %95% confidence bounds
yall = polyval(Pfit, hdeviation);
semilogy(axes7, hdeviation, exp(yall), '-r','LineWidth',2)
yall2 = polyval(Pfit, linspace(0,max(hdeviation)));
semilogy(axes7, linspace(0,max(hdeviation)), exp(yall2), '--r','LineWidth',2)

%Functional expectation
semilogy(axes7, 0:max(hdeviation), exp(-(0:max(hdeviation))/stdh), '--k','LineWidth',2)


sigmafit =  -round(100*(Pfit(1)))/100;
sigmaunc = -round(100*seP(1)./Pfit(1)*sigmafit)/100;
[logintfit,logintunc] = polyconf(Pfit,0,Pstats);
intfit = exp(logintfit);
intunc = exp(logintunc);
%c1 = strcat('$\rm{Slope}/\beta$', '\,=\,', num2str(round(100*(sigmafit*stdh))/100),...
%    '$ \pm $', num2str(round(100*(sigmaunc*stdh))/100));
c1 = strcat('$\rm{Simulated/Theoretical\;Slope}$', '\,=\,', num2str(round(100*(sigmafit*stdh))/100),...
    '$ \pm $', num2str(round(100*(sigmaunc*stdh))/100));
c2 = strcat('$$\Lambda_{i}$$', '\,=\,', '$$\Lambda_{0}\exp\left(-\beta{\delta{h}_i^\star}\right)$$');
%c3 = strcat('$\Lambda\left(0\right)/\Lambda_0$', '\,=\,', num2str(intfit),...
%    '$ \pm $', num2str(intunc));
%c3 = strcat('$\rm{Intercept}/\Lambda_{0}$', '\,=\,', num2str(round(100*(intfit))/100),  '$\pm $', num2str(round(100*(intunc))/100));
c3 = strcat('$\rm{Simulated/Theoretical\;Intercept}$', '\,=\,', num2str(round(100*(intfit))/100),  '$\pm $', num2str(round(100*(intunc))/100));
c4 = strcat('$\Lambda_{0}$','\,=\,','{\boldmath$\Lambda_{\rm{tot}}$}','$\left(1 - \exp\left(-{\beta}\Delta h\right)\right)$');
c5 = strcat('{${\beta}$}','\,=\,','$1/$','{\boldmath$\left<\delta h^\star\right>$}','$=$',num2str(round(100/stdh)/100),'$\rm{\,kg\,kJ^{-1}}$');
c6 = strcat('$\Delta h$','\,=\,',num2str(hBinWidth),'$\,\rm{\,kJ\,kg^{-1}}$');
c7 = strcat('Simulation');
c8 = strcat('Theoretical');

t1 = text(0.1,0.15e-1,c1);
t2 = text(1.5,2.8,c2);
t3 = text(0.1,0.11e-1,c3);
t4 = text(1.5,2.8/1.4,c4);
t5 = text(1.5,2.8/1.4^2,c5);
t6 = text(1.5,2.8/1.4^3,c6);
t7 = text(1.5, 0.4, c7);
t8 = text(0.6, 0.15, c8);

set([t1],'FontWeight','bold','FontSize',18);
set([t1],'Interpreter','Latex')
set([t2],'FontWeight','bold','FontSize',18);
set([t2],'Interpreter','Latex')
set([t3],'FontWeight','bold','FontSize',18);
set([t3],'Interpreter','Latex')
set([t4],'FontWeight','bold','FontSize',18);
set([t4],'Interpreter','Latex')
set([t5],'FontWeight','bold','FontSize',18);
set([t5],'Interpreter','Latex')
set([t6],'FontWeight','bold','FontSize',18);
set([t6],'Interpreter','Latex')
set([t7],'FontWeight','bold','FontSize',18);
set([t7],'Interpreter','Latex','Color','red')
set([t8],'FontWeight','bold','FontSize',18);
set([t8],'Interpreter','Latex')

cmap = parula(nh);
set(axes7,'Clim',[min(hbins) max(hbins)]);
colormap(cmap)
cb = colorbar('peer',axes7,'FontWeight','bold','FontSize',14,'Ticks',hbins);
lcb = ylabel(cb,'$\rm{h^\star\;\left(kJ\,kg^{-1}\right)}$','FontWeight','bold','FontSize',22);
set([lcb],'FontWeight','bold','FontSize',22);
set([lcb],'Interpreter','Latex');

% Create arrow
%annotation('arrow',[0.21 0.21],...
%    [0.775 0.866]);

% Create arrow
%annotation('arrow',[0.333333333333333 0.333333333333333],...
%    [0.279 0.16]);

fname=['Negative_exponents'];
export_fig(sprintf(fname), '-pdf', '-transparent', '-nocrop');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fig3 = figure('Color',[1,1,1],'Position',[50 50 700 550],'PaperPositionMode', 'auto','PaperType','<custom>','PaperSize',[9.5 11]);
color_spectrum = 'parula';

colorgrid = eval([color_spectrum, '(nh)']);
%%%% put temperature deviation on top
axes7 = axes('Position', [0.13         0.13        0.611        0.75],'color','none');
%box(axes7,'off');
hold(axes7,'all');

for ip = 1 : nh
    
    plot(hdeviation(ip), (maxPerimeter_h(ip)./maxPerimeter0), '.', 'color', colorgrid(ip,:),...
        'MarkerSize',40,'Parent',axes7) 
    
end


l1 = xlabel('$$\rm{Deviation\;from\;mean\;}\delta h^\star\;\left({kJ\,kg^{-1}}\right)$$');
l2 = ylabel('$$\lambda_{i,\rm{max}}/\lambda_{\rm{0},\rm{max}}$$');
set([l1],'FontWeight','bold','FontSize',22);
set([l1],'Interpreter','Latex')
set([l2],'FontWeight','bold','FontSize',22);
set([l2],'Interpreter','Latex')

% Create ylabel
%ylabel('Perimeter density (km^{-2})');

%box(axes6,'on');
% Set the remaining axes properties
set(axes7,'FontSize',14,'FontWeight','bold','LineStyleOrderIndex',3,...
    'XMinorTick','on','XScale','lin','YMinorTick','on','YScale','log','XLim',[0 6],'YLim',[0.0091 3.1]);
% Create legend
%leg1 = legend(axes6,'show');
%set(leg1, 'Position',[0.85 0.17 0.12 0.78]);

[Pfit, Pstats] = polyfit(hdeviation,(log(maxPerimeter_h/maxPerimeter0)),1);
seP = 2*sqrt(sum(inv(Pstats.R).^2,2))*Pstats.normr/sqrt(Pstats.df); %95% confidence bounds
yall = polyval(Pfit, hdeviation);
semilogy(axes7, hdeviation, exp(yall), '-r','LineWidth',2)
yall2 = polyval(Pfit, linspace(0,max(hdeviation)));
semilogy(axes7, linspace(0,max(hdeviation)), exp(yall2), '--r','LineWidth',2)

%Functional expectation
semilogy(axes7, 0:max(hdeviation), exp(-(0:max(hdeviation))/stdh), '--k','LineWidth',2)

sigmafit =  -round(100*(Pfit(1)))/100;
sigmaunc = -round(100*seP(1)./Pfit(1)*sigmafit)/100;
[logintfit,logintunc] = polyconf(Pfit,0,Pstats);
intfit = exp(logintfit);
intunc = exp(logintunc);
%c1 = strcat('$\rm{Slope}/\beta$', '\,=\,', num2str(round(100*(sigmafit*stdh))/100),...
%    '$ \pm $', num2str(round(100*(sigmaunc*stdh))/100));
c1 = strcat('$\rm{Simulated/Theoretical\;Slope}$', '\,=\,', num2str(round(100*(sigmafit*stdh))/100),...
    '$ \pm $', num2str(round(100*(sigmaunc*stdh))/100));
c2 = strcat('$$\Lambda_{i}$$', '\,=\,', '$$\Lambda_{0}\exp\left(-\beta{\delta{h}_i^\star}\right)$$');

c2 = strcat('$$\lambda_{i,\rm{max}}$$', '\,=\,', '$$\lambda_{0,\rm{max}}\exp\left(-\beta{\delta{h}_i^\star}\right)$$');
%c3 = strcat('$\Lambda_{max}\left(0\right)/\Lambda_{0,{max}}$', '\,=\,', num2str(intfit),...
%    '$ \pm $', num2str(intunc));
c3 = strcat('$\rm{Intercept}/\lambda_{0,\rm{max}}$', '\,=\,', num2str(round(100*(intfit))/100),  '$\pm $', num2str(round(100*(intunc))/100));
c3 = strcat('$\rm{Simulated/Theoretical\;Intercept}$', '\,=\,', num2str(round(100*(intfit))/100),  '$\pm $', num2str(round(100*(intunc))/100));
c4 = strcat('$\lambda_{0,\rm{max}}$','\,=\,','{\boldmath$\lambda_{\rm{tot},\rm{max}}$}','$\left(1 - \exp\left(-{\beta}\Delta h\right)\right)$');
c5 = strcat('{${\beta}$}','\,=\,','$1/$','{\boldmath$\left<\delta h^\star\right>$}','$=$',num2str(round(100/stdh)/100),'$\rm{\,kg\,kJ^{-1}}$');
c6 = strcat('$\Delta h$','\,=\,',num2str(hBinWidth),'$\,\rm{\,kJ\,kg^{-1}}$');
c7 = strcat('Simulation');
c8 = strcat('Theoretical');

t1 = text(0.1,0.15e-1,c1);
t2 = text(1.5,2.8,c2);
t3 = text(0.1,0.11e-1,c3);
t4 = text(1.5,2.8/1.4,c4);
t5 = text(1.5,2.8/1.4^2,c5);
t6 = text(1.5,2.8/1.4^3,c6);
t7 = text(2.7, 0.25, c7);
t8 = text(0.6, 0.15, c8);

set([t1],'FontWeight','bold','FontSize',18);
set([t1],'Interpreter','Latex')
set([t2],'FontWeight','bold','FontSize',18);
set([t2],'Interpreter','Latex')
set([t3],'FontWeight','bold','FontSize',18);
set([t3],'Interpreter','Latex')
set([t4],'FontWeight','bold','FontSize',18);
set([t4],'Interpreter','Latex')
set([t5],'FontWeight','bold','FontSize',18);
set([t5],'Interpreter','Latex')
set([t6],'FontWeight','bold','FontSize',18);
set([t6],'Interpreter','Latex')
set([t7],'FontWeight','bold','FontSize',18);
set([t7],'Interpreter','Latex','Color','red')
set([t8],'FontWeight','bold','FontSize',18);
set([t8],'Interpreter','Latex')

cmap = parula(nh);
set(axes7,'Clim',[min(hbins) max(hbins)]);
colormap(cmap)
cb = colorbar('peer',axes7,'FontWeight','bold','FontSize',14,'Ticks',hbins);
lcb = ylabel(cb,'$\rm{h^\star\;\left(kJ\,kg^{-1}\right)}$','FontWeight','bold','FontSize',22);
set([lcb],'FontWeight','bold','FontSize',22);
set([lcb],'Interpreter','Latex')

% Create arrow
%annotation('arrow',[0.21 0.21],...
%    [0.775 0.866]);

% Create arrow
%annotation('arrow',[0.333333333333333 0.333333333333333],...
%    [0.279 0.16]);

fname=['Negative_exponents_Lambda_max'];
export_fig(sprintf(fname), '-pdf', '-transparent', '-nocrop');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fig4 = figure('Color',[1,1,1],'Position',[50 50 700 550],'PaperPositionMode', 'auto','PaperType','<custom>','PaperSize',[9.5 11]);
color_spectrum = 'parula';

colorgrid = eval([color_spectrum, '(nh)']);
%%%% put temperature deviation on top
axes7 = axes('Position', [0.13         0.13        0.611        0.75],'color','none');
%box(axes7,'off');
hold(axes7,'all');

for ip = 1 : nh
    
    plot(hdeviation(ip), Number0_h(ip)./Number0, '.', 'color', colorgrid(ip,:),...
        'MarkerSize',40,'Parent',axes7) 
    
end
l1 = xlabel('$$\rm{Deviation\;from\;mean\;}\delta h^\star\;\left({kJ\,kg^{-1}}\right)$$');
l2 = ylabel('$$n_{i,\rm{min}}/n_{0,\rm{min}}$$');
set([l1],'FontWeight','bold','FontSize',22);
set([l1],'Interpreter','Latex')
set([l2],'FontWeight','bold','FontSize',22);
set([l2],'Interpreter','Latex')
% Create ylabel
%ylabel('Perimeter density (km^{-2})');

%box(axes6,'on');
% Set the remaining axes properties
set(axes7,'FontSize',14,'FontWeight','bold','LineStyleOrderIndex',3,...
    'XMinorTick','on','XScale','lin','YMinorTick','on','YScale','log','XLim',[0 6],'YLim',[0.0091 3.1]);
% Create legend
%leg1 = legend(axes6,'show');
%set(leg1, 'Position',[0.85 0.17 0.12 0.78]);

[Pfit, Pstats] = polyfit(hdeviation,log(Number0_h./Number0),1);
seP = 2*sqrt(sum(inv(Pstats.R).^2,2))*Pstats.normr/sqrt(Pstats.df); %95% confidence bounds
yall = polyval(Pfit, hdeviation);
semilogy(axes7, hdeviation, exp(yall), '-r','LineWidth',2)
yall2 = polyval(Pfit, linspace(0,max(hdeviation)));
semilogy(axes7, linspace(0,max(hdeviation)), exp(yall2), '--r','LineWidth',2)

%Functional expectation
semilogy(axes7, 0:max(hdeviation), exp(-(0:max(hdeviation))/stdh), '--k','LineWidth',2)

sigmafit =  -round(100*(Pfit(1)))/100;
sigmaunc = -round(100*seP(1)./Pfit(1)*sigmafit)/100;
[logintfit,logintunc] = polyconf(Pfit,0,Pstats);
intfit = exp(logintfit);
intunc = exp(logintunc);

%c1 = strcat('$\rm{Slope}/\beta$', '\,=\,', num2str(round(100*(sigmafit*stdh))/100),...
%    '$ \pm $', num2str(round(100*(sigmaunc*stdh))/100));
c1 = strcat('$\rm{Simulated/Theoretical\;Slope}$', '\,=\,', num2str(round(100*(sigmafit*stdh))/100),...
    '$ \pm $', num2str(round(100*(sigmaunc*stdh))/100));
c2 = strcat('$$\Lambda_{i}$$', '\,=\,', '$$\Lambda_{0}\exp\left(-\beta{\delta{h}_i^\star}\right)$$');

c2 = strcat('$$n_{i,min}$$', '\,=\,', '$$n_{0,min}\exp\left(-\beta{\delta{h}_i^\star}\right)$$');
c3 = strcat('$\rm{Simulated/Theoretical\;Intercept}$', '\,=\,', num2str(round(100*(intfit))/100),  '$\pm $', num2str(round(100*(intunc))/100));
c4 = strcat('$n_{0,min}$','\,=\,','{\boldmath$n_{\rm{tot},\rm{min}}$}','$\left(1 - \exp\left(-{\beta}\Delta h\right)\right)$');
c5 = strcat('{${\beta}$}','\,=\,','$1/$','{\boldmath$\left<\delta h^\star\right>$}','$=$',num2str(round(100/stdh)/100),'$\rm{\,kg\,kJ^{-1}}$');
c6 = strcat('$\Delta h$','\,=\,',num2str(hBinWidth),'$\,\rm{\,kJ\,kg^{-1}}$');
c7 = strcat('Simulation');
c8 = strcat('Theoretical');
t1 = text(0.1,0.15e-1,c1);
t2 = text(1.5,2.8,c2);
t3 = text(0.1,0.11e-1,c3);
t4 = text(1.5,2.8/1.4,c4);
t5 = text(1.5,2.8/1.4^2,c5);
t6 = text(1.5,2.8/1.4^3,c6);
t7 = text(1.5, 0.4, c7);
t8 = text(0.6, 0.15, c8);


set([t1],'FontWeight','bold','FontSize',18);
set([t1],'Interpreter','Latex')
set([t2],'FontWeight','bold','FontSize',18);
set([t2],'Interpreter','Latex')
set([t3],'FontWeight','bold','FontSize',18);
set([t3],'Interpreter','Latex')
set([t4],'FontWeight','bold','FontSize',18);
set([t4],'Interpreter','Latex')
set([t5],'FontWeight','bold','FontSize',18);
set([t5],'Interpreter','Latex')
set([t6],'FontWeight','bold','FontSize',18);
set([t6],'Interpreter','Latex')
set([t7],'FontWeight','bold','FontSize',18);
set([t7],'Interpreter','Latex','Color','red')
set([t8],'FontWeight','bold','FontSize',18);
set([t8],'Interpreter','Latex')

cmap = parula(nh);
set(axes7,'Clim',[min(hbins) max(hbins)]);
colormap(cmap)
cb = colorbar('peer',axes7,'FontWeight','bold','FontSize',14,'Ticks',hbins);
lcb = ylabel(cb,'$\rm{h^\star\;\left(kJ\,kg^{-1}\right)}$','FontWeight','bold','FontSize',22);
set([lcb],'FontWeight','bold','FontSize',22);
set([lcb],'Interpreter','Latex')

% Create arrow
%annotation('arrow',[0.21 0.21],...
%    [0.775 0.866]);

% Create arrow
%annotation('arrow',[0.333333333333333 0.333333333333333],...
%    [0.279 0.16]);

fname=['Negative_exponents_n'];
export_fig(sprintf(fname), '-pdf', '-transparent', '-nocrop');


