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','','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','','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','','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','','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');