
%% edit these!

clear
% directory with RAMS output files
ttlname='U10_CONTROL'; ttlname2=ttlname; ttlname2(ttlname=='_')='-';
ftype='L'; % "A" for analysis or "L" for lite files
iread=1; % read files? 1=read, 0=load matlab file

% set RAMSdir to the base directory where the output is located 
RAMSdir = ['[path_to_output_directory]/',ttlname,'/NOBAK/'];
SaveDir = [ttlname,'/',ftype,'files/'];

% grid setup
dx = 1.; % horiz grid spacing in km
nx = 3001; ny = 201; % # horizontal grid points

% grid point indices for reading and plotting
% x=1 and x=2 are the same; x=end-1 and x=end are the same, and same
% as x=1 and 2 -- so, 3:end-1 are unique x's; same for y
x1=3; x2=nx-1;
y1=3; y2=ny-1;
z1=2; z2=150;
tskip = 3; % skip interval to read; 3 = once per hour
tend = 9999; % last file to read

% reading / plotting settings
% frequency (days) to print out time display for read loops
frqdisp = .25;
% day interval for smoothing time series (e.g. 5-day running mean)
daymean = .25;
daystart = 20; % day start - in case first file is not from day 0


%% set up some arrays/values based on above settings
x = [0:nx-1]*dx; y = [0:ny-1]*dx; % km
zt = read_zlevs_hfile( RAMSdir, 't' )/1000.; % km
zm = read_zlevs_hfile( RAMSdir, 'm' )/1000.;
nz = length(zt);
if z2>nz; z2=nz; end
if ny==1; y1=1; y2=1; end % for 2D

% spatial arrays for reading and plotting (km)
xp = x(x1:x2); nxp=length(xp);
yp = y(y1:y2); nyp=length(yp);
zp = zt(z1:z2); nzp=length(zp);

% files
files = dir([RAMSdir,'/*-',ftype,'-*h5']);

% time info
nt = min(tend,length(files)); % # output times
%for t=1:tskip:nt; timestr(t,:)=files(t).name(5:21); end % yyyy-mm-dd-HHMMSS
for t=1:nt; timestr(t,:)=files(t).name(findstr(files(t).name,['-',ftype,'-'])+[3:19]); end % yyyy-mm-dd-HHMMSS
timevec = datevec(timestr,'yyyy-mm-dd-HHMMSS'); % Date Vector array
time1 = repmat(timevec(1,:),[size(timevec,1) 1]);
tp = etime(timevec,time1) / 3600/24. + daystart; % time plotting array (days)
ntp=length(tp);
% dt represents time in days between current file and the file prior; dt(1)=dt(2)
dt = zeros(size(tp)); dt(2:end)=diff(tp); dt(1)=dt(2);
clear timestr timevec time1

% [daymean, ex: 5]-day moving averages for the smallest output interval; this means that
% if output interval changes, the [5]-day moving average will actually be
% more than a [5]-day moving average where dt is coarser
avgdt = ceil( daymean/min(dt) );
avgdt = avgdt + mod(avgdt+1,2); % this makes sure avgdt is an odd number


% make output dir if it doesn't exist
if exist(SaveDir,'dir')~=7; mkdir(SaveDir); end

% save dims
save([SaveDir,'dims'],'xp','yp','zp','tp')


%% read section

if iread

    ipviewfig=1;
    if ipviewfig
        % create a figure handle
        ifig=figure; set_figsize(gcf,[15 10]); om=[.08 .1 .03 .06]; im=[.05 .08];
    end
    if exist([SaveDir,'vintcond_analysis'],'dir')~=7; mkdir([SaveDir,'vintcond_analysis']); end
    
    % initialize empty structure array regstatsfull
    regstatsall=struct([]);

    for t=1:tskip:nt
        
        cp=1004.; Rd=287.;
        if mod(tp(t),frqdisp)==0; disp(['day ',num2str(tp(t))]); end
        
        %% read in vintcond
        % vertically integrated condensate - need totcond, rho0, zzm
        totcond = (  h5read( [RAMSdir,files(t).name],'/RTP', [x1 y1 z1], [nxp nyp nzp] )...
            - h5read( [RAMSdir,files(t).name],'/RV',  [x1 y1 z1], [nxp nyp nzp] ) )*1000.; % g/kg
        if exist('rho0','var')==0
            % get base state density from level 2 to top-1
            fn0 = [RAMSdir,files(1).name]; % file name
            press0 = h5read( fn0,'/PI', [x1 y1 2], [nxp nyp nz-2] )/cp; % read from 2:nz-1 - first just exner function (non-dimensional)
            Tvirt0 = ( press0 .* h5read( fn0,'/THETA',[x1 y1 2],[nxp nyp nz-2] ) ) ... % press0 is exner0 here
                .* ( 1.+0.61*h5read( fn0,'/RV',[x1 y1 2],[nxp nyp nz-2] ) );
            press0 = 100000.*press0.^(cp/Rd); % pressure (Pa)
            rho0 = press0/Rd./Tvirt0; clear press0 Tvirt0 fn0
        end
        if exist('zzm','var')==0
            % get 3D height grid
            [~,~,zzm]=meshgrid(yp,xp,zm(1:nz-1)); % zm=1:nz-1 gives PW for zt=2:nz-1
        end
        % vertintcond = int(sfc->top){totcond*rho0*dz} [kg/m2=mm]; exclude last level in totcond
        vintcond = sum( totcond(:,:,1:end-1)/1000. .* rho0 .* diff(zzm,1,3)*1000, 3 );
        
        
        % also read / calculate density temperature at z=2, lowest level above ground
        % thetarho = theta*(1+.61*rv-rcond) = theta*(1+.61*rv-(rt-rv))
        %          = theta*(1 - rt + 1.61*rv)
        thrhopert = h5read([RAMSdir,files(t).name],'/THETA',[x1 y1 2],[nxp nyp 1]) ...
                        .* ( 1. - h5read([RAMSdir,files(t).name],'/RTP',[x1 y1 2],[nxp nyp 1]) ...
                                + 1.61* h5read([RAMSdir,files(t).name],'/RV',[x1 y1 2],[nxp nyp 1]) );
        % perturbation thetarho: subtract off median (less sensitive to
        % variations in cold pool strength in CTL than the mean is)
        thrhopert = thrhopert - median(thrhopert(:));
        
        % read in v-wind and save absolute value
        absv = abs( h5read([RAMSdir,files(t).name],'/VP',[x1 y1 2],[nxp nyp 1]) );
        
        % calculate cloud top height
        for i=1:nxp; for j=1:nyp
                idx=find( squeeze(totcond(i,j,:))>=0.01, 1, 'last' ); % 0.01 g/kg threshold for cloud top
                % 0.01 g/kg produces good match to OLR field based on one snapshot from U2_20m at time 1...
                if isempty(idx); CT_hgt(i,j)=0; % clear cloumn
                else CT_hgt(i,j) = zp( idx );
                end
        end; end

        
        %% smooth vertically integrated condensate field
        
        % expand vintcond field to repeat itself 2x in y
        field_rep2y = [vintcond vintcond];
        yp_rep2y = [yp yp+(yp(end)-yp(1)+1)];
        % also expand thrhopert, v, and CT_hgt
        thrhopert_rep2y = [thrhopert thrhopert];
        absv_rep2y = [absv absv];
        CT_hgt_rep2y = [CT_hgt CT_hgt];
        
        % create time string for plotting
        timestr=sprintf('%2.0f d %4.2f h',floor(tp(t)),(tp(t)-floor(tp(t)))*24);
        
        if ipviewfig
            % plot unsmoothed vintcond in panel 1
            figure(ifig); clf; multipan(ifig,3,1,1,'om',om,'im',im);
            pcolor(xp,yp_rep2y,field_rep2y'); shading flat; colorbar
            set(gca,'clim',[0 10]); ylabel('y (km)'); set(gca,'tickdir','out'); set(gca,'xticklabel',[]);
            title(['unsmoothed, ',ttlname2,', ',timestr])
        end
        
        % smooth 2D vintcond field using Gaussian filter with sigma 5km
        % create the 2D filter using multivariate normal function
        sigmax=5; % km
        sigmay=5;
        mux=0; % 0 for centered
        muy=0;
        
        % create array where half-width is 2x sigma
        [arrayy,arrayx] = meshgrid(-2*sigmay:dx:2*sigmay, -2*sigmax:dx:2*sigmax);
        arrlengthx=size(arrayx,1); arrlengthy=size(arrayy,2);
        
        % from functions online: need to reshape the 2D x and y arrays
        % also pass in mu, and a 1x2 array of sigma
        % (this means diagonal elements are sigma and off-diagonal elements are 0)
        gauss2d = mvnpdf([arrayy(:) arrayx(:)],[muy mux],[sigmay sigmax]);
        % reshape
        gauss2d = reshape( gauss2d, arrlengthx, arrlengthy );
        
        % checks on filter
        %figure; pcolor(gauss2d); axis equal
        %sum(gauss2d(:))
        
        smoothedfield = conv2( field_rep2y, gauss2d, 'same' );
        
        if ipviewfig
            % plot the smoothed field
            figure(ifig); multipan(ifig,3,1,2,'om',om,'im',im);
            pcolor(xp,yp_rep2y,smoothedfield'); shading flat; colorbar
            set(gca,'clim',[0 10]); ylabel('y (km)'); set(gca,'tickdir','out'); set(gca,'xticklabel',[]);
            title(['smoothed'])
        end
        
        % create binary image from smoothed field with values above threshold
        thresh=1.; % mm
        binary=zeros(size(smoothedfield));
        binary(smoothedfield>=thresh)=1;
        
        if ipviewfig
            % plot the binary field
            figure(ifig); multipan(ifig,3,1,3,'om',om,'im',im);
            pcolor(xp,yp_rep2y,binary'); shading flat; colorbar; hold all
            colormap(gca,parula(2)); xlabel('x (km)'); ylabel('y (km)'); set(gca,'tickdir','out');
            title(['binary >=',num2str(thresh),'mm'])
            linkaxes
        end
        
        
        %% connected regions stats
        
        % find connected regions with bwconncomp (uses default connectivity of 8 for 2D images)
        CC = bwconncomp(binary);
        % calculate stats, and pass in original vintcond field to get max and mean intensity
        regionstats = regionprops(CC,field_rep2y,'area','centroid','eccentricity','majoraxislength','minoraxislength','orientation','pixellist',...
            'maxintensity','meanintensity','weightedcentroid');
        % pass in thrhopert to get min and mean thtrho' intensity; same for absv and CT_hgt
        regionstats2 = regionprops(CC,thrhopert_rep2y,'minintensity','meanintensity');
        regionstats3 = regionprops(CC,absv_rep2y,'maxintensity','meanintensity');
        regionstats4 = regionprops(CC,CT_hgt_rep2y,'maxintensity','meanintensity');
        % add thrhopert values to regionstats from regionstats2; same for absv and regionstats3, and CThgt and regionstats4
        [regionstats.MeanThrhopert] = regionstats2.MeanIntensity;
        [regionstats.MinThrhopert] = regionstats2.MinIntensity; clear regionstats2
        [regionstats.MeanAbsv] = regionstats3.MeanIntensity;
        [regionstats.MaxAbsv] = regionstats3.MaxIntensity; clear regionstats3
        [regionstats.MeanCThgt] = regionstats4.MeanIntensity;
        [regionstats.MaxCThgt] = regionstats4.MaxIntensity; clear regionstats4
        % get rid of regions that do not have one axis at least 10 km in length
        % use arrayfun to define operations on the structure array
        %regionstats=regionstats( arrayfun( @(x) x.MajorAxisLength>=10||x.MinorAxisLength>=10, regionstats ) );
        % modification: get rid of regions that are smaller than 100 km2
        regionstats=regionstats( arrayfun( @(x) x.Area*dx*dx>=100., regionstats ) );
        % get rid of non-unique regions - first find any that have any border pixels and get rid of them
        idx=arrayfun( @(x) any(x.PixelList(:,1)==nyp*2)||any(x.PixelList(:)==1)||any(x.PixelList(:,2)==nxp), regionstats );
        regionstats=regionstats(~idx);
        % next save some stats for uniqueness and check for unique rows; use 'stable' to keep indices in the same order
        % use max intensity to define uniqueness
        maxintensity=arrayfun(@(x) x.MaxIntensity, regionstats);
        [~,idx,~]=unique(maxintensity,'stable');
        regionstats=regionstats(idx);
        % add day field to the structure
        [regionstats(:).day]=deal(tp(t));
        nregions=length(regionstats);
        
        if ipviewfig
            figure(ifig);
            % plot the centers of regions on the binary image
            arrayfun(@(x) plot(x.Centroid(2)*dx,x.Centroid(1)*dx,'or'), regionstats)
            
            if 1
            % plot the ellipses
            tt=linspace(0,2*pi,200); % 0-2pi
            for r=1:length(regionstats)
                theta=90+regionstats(r).Orientation; % array angle: angle from x-axis instead of y-axis
                a = regionstats(r).MajorAxisLength;
                b = regionstats(r).MinorAxisLength;
                a=a/2; b=b/2; % radius not diameter
                center = fliplr(regionstats(r).Centroid); % centroid is flipped around, change to (x,y)
                xtt = a*cos(tt); % parametric equations for x and y, from x-axis
                ytt = b*sin(tt);
                % now rotate the ellipse
                cosangle=cosd(theta); sinangle=sind(theta);
                xx = xtt*cosangle - ytt*sinangle + center(1);
                yy = xtt*sinangle + ytt*cosangle + center(2);
                plot(xx,yy,'r-','linewidth',3)
            end
            end
            print('-dpng',[SaveDir,sprintf('vintcond_analysis/3pan_smoothed_binary_day%06.2f.png',tp(t))])
        end
        
        % save the data
        if isempty(regstatsall)
            regstatsall=regionstats;
        else
            [regstatsall(end+[1:nregions])]=regionstats;
        end
        
    end
    % save regstatsall
    save([SaveDir,'vintcond_stats'],'regstatsall')
else
    load([SaveDir,'vintcond_stats'])
end


%% stats plots

%% set up for plots


if exist([SaveDir,'vintcond_feature_stats'],'dir')~=7; mkdir([SaveDir,'vintcond_feature_stats']); end

% save some arrays from the structure for plotting
majoraxlen=arrayfun(@(x) x.MajorAxisLength*dx, regstatsall); % major axis length (km)
minoraxlen=arrayfun(@(x) x.MinorAxisLength*dx, regstatsall); % minor axis length (km)
eccen=arrayfun(@(x) x.Eccentricity, regstatsall); % eccentricity (0: circle, 1: line segment) => how does this relate to axis ratio?
  % Xu Rutledge 2015 classified linear as major:minor aspect ratio > 5:2
orient=90.-abs(arrayfun(@(x) x.Orientation, regstatsall)); % angle from x-axis of feature: here since fields need to be transposed
  % for plotting, I think it is actually angle from the y-axis. So 0 = N/S oriented, and 90 = E/W oriented.
  % Transpose this to be angle from shear-vector: 90-orientation; 0=E/W; 90=N/S
farea=arrayfun(@(x) x.Area*dx*dx, regstatsall); % area (km^2) (= number of pixels in the region)
day=arrayfun(@(x) x.day, regstatsall); % day (20-25)
maxintensity=arrayfun(@(x) x.MaxIntensity, regstatsall); % max intensity (mm)
minthrhopert=arrayfun(@(x) x.MinThrhopert, regstatsall); % min theta-rho pert at z=2 (K)
meanthrhopert=arrayfun(@(x) x.MeanThrhopert,regstatsall); % mean theta-rho pert at z=2 (K)
maxabsv=arrayfun(@(x) x.MaxAbsv,regstatsall); % max |v| (m/s)
meanabsv=arrayfun(@(x) x.MeanAbsv,regstatsall); % mean |v| (m/s)
vararray=arrayfun(@(x) x.MaxCThgt,regstatsall); % max cloup top height (km)
%vararray=arrayfun(@(x) x.MeanCThgt,regstatsall); % mean cloud top height (km)

idays=day>=22; % indices for idays
imajax50km=majoraxlen>=50; % indices for major axis length >= 50 km
iarea500km2=farea>=500.; % indices for area >= 500 sq km
iarea1000km2=farea>=1000.;
iarea100km2=farea>=100.;

% define anonymous functions
% slope and intercept of linear fits
calc_a_b = @( xpts, ypts )  polyfit( xpts, ypts, 1 );
% calculate linear fit data given x, a, and b
calc_lfit = @( xpts, a_b_coefs )  xpts * a_b_coefs(1) + a_b_coefs(2);
% R^2 for variance explained by linear fit
calc_Rsqd = @( ypts,lfit )  1. - sum( ( ypts-lfit       ).^2 ) ...
                               ./sum( ( ypts-mean(ypts) ).^2 );

%% plots

% scatter plot of axis ratio vs orientation, colored by max vintcond 
figure; hold all
scatter(minoraxlen(idays)./majoraxlen(idays),abs(orient(idays)),farea(idays)/20,maxintensity(idays),'filled'); colorbar
plot([.4 .4],[0 90],'k:','linewidth',2)
set(gca,'clim',[0 80]); 
xlim([0 1]); ylim([0 90]);
xlabel('Minor:Major Axis Length')
ylabel('Angle from shear vector (\circ)')
title({'Convective feature stats';'(size~area, color=max integrated condensate, mm)'})
print('-dpng',[SaveDir,'vintcond_feature_stats/scatter_min-maj_orient_area_maxintensity'])


% scatter plot of axis ratio vs orientation, colored by time: all days
figure; hold all
% do for area >= 100 km^2
idx=iarea100km2;
scatter( minoraxlen(idx)./majoraxlen(idx), abs(orient(idx)), ...
    farea(idx)/20, day(idx), 'filled' ); colorbar
colormap(parula(10))
plot([.4 .4],[0 90],'k:','linewidth',2)
set(gca,'clim',[20 25]); xlim([0 1]); ylim([0 90]);
xlabel('Minor:Major Axis Length')
ylabel('Angle from shear vector (\circ)')
%title({'Convective feature stats';'(size~area, color=day; major ax. length >= 50 km)'})
%print('-dpng',[SaveDir,'vintcond_feature_stats/scatter_min-maj_orient_area_day_majaxlen50'])
title({'Convective feature stats';'(size~area, color=day; area >= 100 km^2)'})
print('-dpng',[SaveDir,'vintcond_feature_stats/scatter_min-maj_orient_area_day_area100'])


% heat map of axis ratio vs orientation, normalized by # objects
figure; hold all
axratio=minoraxlen./majoraxlen; absorient=abs(orient);
% do only ones with area >= 100 km^2
idx=iarea100km2;
axratio=axratio(idx); absorient=absorient(idx);
% calculate 2D PDF of # objects in each bin using histcounts2
% histcounts2( X, Y, Xedges, Yedges )
% define bins for x and y
bins_axratio = 0:.05:1; bins_orient = 0:5:90;
heatmap = histcounts2( axratio, absorient, bins_axratio, bins_orient );
% normalize by total # counts and convert to %
heatmap = heatmap / sum(heatmap(:)) * 100;
%heatmap(heatmap==0)=NaN;
% repeat last row and column for pcolor, which uses the bottom left point
% for shading information
heatmap(:,end+1)=heatmap(:,end); heatmap(end+1,:)=heatmap(end,:);
pcolor( bins_axratio, bins_orient, heatmap' ); shading flat; colorbar
set(gca,'clim',[0 5]); colormap(parula(25))
plot([.4 .4],[0 90],'k:','linewidth',2)
xlim([0 1]); ylim([0 90])
xlabel('Minor:Major Axis Length')
ylabel('Angle from shear vector (\circ)')
%title({'2D histogram of convective feature stats';'(% of total counts; major ax. length >= 50 km)'})
%print('-dpng',[SaveDir,'vintcond_feature_stats/heatmap_min-maj_orient_majaxlen50'])
title({'2D histogram of convective feature stats';'(% of total counts; area >= 100 km^2)'})
print('-dpng',[SaveDir,'vintcond_feature_stats/heatmap_min-maj_orient_area100'])


%% box plot - all data overtop one another

% set up arrays
r=0;
r=r+1; fullname{r}='U2_CONTROL';
r=r+1; fullname{r}='U2_NO-EVAP';
r=r+1; fullname{r}='U5_CONTROL';
r=r+1; fullname{r}='U5_NO-EVAP';
r=r+1; fullname{r}='U10_CONTROL';
r=r+1; fullname{r}='U10_NO-EVAP';
nr=r;
labelname=fullname;

boxdataall=[];
for r=1:nr
    
    % load data file
    matfile2read = [fullname{r},'/',ftype,'files/vintcond_stats.mat'];
    load(matfile2read);
    
    % pull out orientation
    boxdata=abs( arrayfun(@(x) x.Orientation, regstatsall) ); % angle from y-axis; 0 = N/S oriented, and 90 = E/W oriented
    boxdata=90-boxdata; % change to angle from E/W/shear vector: 0=E/W, 90=N/S
    % pull out aspect ratio, area, and day
    aspectratio=arrayfun(@(x) x.MinorAxisLength./x.MajorAxisLength, regstatsall); % minor:major axis ratio
    area=arrayfun(@(x) x.Area*dx*dx, regstatsall); % area of feature (km^2)
    day=arrayfun(@(x) x.day, regstatsall); % day (20-25)
    
    % only save stats for features with area >= 100 km2
    iarea100km2=area>=100.; % indices for area >= 100 sq km
    boxdata = boxdata(iarea100km2);
    aspectratio = aspectratio(iarea100km2);
    area = area(iarea100km2);
    day = day(iarea100km2);
    
    % diagnostic: save total # systems
    nsystems_total(r)=length(boxdata);
    % diagnostic: save total # linear systems
    boxdatatemp=boxdata(aspectratio<=.4); % temp array for below line of code
    nlinear_systems_total(r)=length(boxdatatemp); clear boxdatatemp
    
    % calculate box-data, aspect ratio, and area for only days 22-25
    boxdata=boxdata(day>=22);
    aspectratio=aspectratio(day>=22);
    area=area(day>=22);
    day=day(day>=22);
    
    % diagnostic: save # of systems from days 22-25
    nsystems_days22_25(r)=length(boxdata);
    % diagnostic: mean area of systems for days 22-25
    meanarea_days22_25(r)=mean(area);
    
    % save only linear systems with aspect ratio <= 2:5
    boxdata=boxdata(aspectratio<=.4);
    area=area(aspectratio<=.4);
    day=day(aspectratio<=.4);
    aspectratio=aspectratio(aspectratio<=.4);
    
    % diagnostic: save # of linear systems from days 22-25
    nlinear_systems_days22_25(r)=length(boxdata);
    % diagnostic: mean area of linear systems for days 22-25
    meanarea_linear_days22_25(r)=mean(area);
    
    if isempty(boxdataall)
        boxdataall=boxdata; % save box data in vector
        labeldata=zeros(size(boxdata))+r; % 1-6 for the different simulations
    else
        boxdataall(end+[1:length(boxdata)])=boxdata;
        labeldata(end+[1:length(boxdata)])=zeros(size(boxdata))+r;
    end
    
    % save 5 and 95 percentiles so I can change what the whiskers show
    pctiles_5_95(r,1:2) = prctile(boxdata,[5 95]);
    
end
clear boxdata

% box plot
figure; hold all
boxplot(boxdataall,labeldata,'colors','krkrkr','symbol','k')
set(gca,'xticklabel',labelname,'xtickLabelRotation',45)
ylim([0 90])
ylabel('Angle from shear vector')

% get all lines with tags to change data
UpWhisk=findall(gcf,'tag','Upper Whisker');
UpAdj=findall(gcf,'tag','Upper Adjacent Value');
LoWhisk=findall(gcf,'tag','Lower Whisker');
LoAdj=findall(gcf,'tag','Lower Adjacent Value');
% loop through runs and change y-data to be 5 and 95 percentiles
for r=1:nr
    % change upper whisker to 95 percentile
    irun = UpWhisk(r).XData(1); % NOTE: Assume plot ordering is the same for all sets of line types
    UpWhisk(r).YData(2) = pctiles_5_95(irun,2);
    UpAdj(r).YData(:) = pctiles_5_95(irun,2);
    LoWhisk(r).YData(1) = pctiles_5_95(irun,1);
    LoAdj(r).YData(:) = pctiles_5_95(irun,1);
end

% increase line widths
ln=findall(gcf,'type','line');
for i=1:length(ln)
        ln(i).LineWidth=3;
        if strcmp(ln(i).LineStyle,'--')
            ln(i).LineStyle='-';
        end
end
box off
title('Orientations of Linear Systems, days 22-25')
print('-dpng','boxplots_vintcond_stats_linsystems_day22-25_area100km2')


% pull the median, 25, and 75 lines
medline=findall(gcf,'tag','Median');
boxline=findall(gcf,'tag','Box');
for r=1:nr;
    median(r)=medline(r).YData(1);
    pct25(r)=boxline(r).YData(1);
    pct75(r)=boxline(r).YData(2);
end
% flip since i=1 in the plot is for u10_noevap etc.
median=fliplr(median); pct25=fliplr(pct25); pct75=fliplr(pct75);

%% PDFs of cloud top heights, comparing sets of runs

% set up arrays
r=0;
r=r+1; fullname{r}='U2_CONTROL'; labelname{r}='U2\_CONTROL';
r=r+1; fullname{r}='U2_NO-EVAP'; labelname{r}='U2\_NO-EVAP';
r=r+1; fullname{r}='U5_CONTROL'; labelname{r}='U5\_CONTROL';
r=r+1; fullname{r}='U5_NO-EVAP'; labelname{r}='U5\_NO-EVAP';
r=r+1; fullname{r}='U10_CONTROL'; labelname{r}='U10\_CONTROL';
r=r+1; fullname{r}='U10_NO-EVAP'; labelname{r}='U10\_NO-EVAP';
nr=r;

idays=1; % 0 = days 20-25, 1 = only 22-25

% variable to plot
var2plot='MaxCThgt'; bins=[3:.5:17]; 
  ttl='Max Cloud Top Heights'; xlab='Cloud Top Height (km)'; ylabunit='km^{-1}';
% var2plot='MeanCThgt'; bins=[3:.5:17];
%   ttl='Mean Cloud Top Heights'; xlab='Cloud Top Height (km)'; ylabunit='km^{-1}';
% var2plot='MaxAbsv'; bins=[0:1:20]; 
%   ttl='Max |v|'; xlab='|v| (m s^{-1})'; ylabunit='(m/s)^{-1}';
% var2plot='MeanAbsv'; bins=[0:.25:6]; 
%   ttl='Mean |v|'; xlab='|v| (m s^{-1})'; ylabunit='(m/s)^{-1}';
% var2plot='MaxIntensity'; bins=[0:5:100]; 
%   ttl='Max Integrated Condensate'; xlab='mm'; ylabunit='mm^{-1}';
% var2plot='MeanIntensity'; bins=[0:.25:10]; 
%   ttl='Mean Integrated Condensate'; xlab='mm'; ylabunit='mm^{-1}';
% var2plot='MajorAxisLength'; bins=[0:5:150]; 
%   ttl='Major Axis length'; xlab='km'; ylabunit='km^{-1}';
% var2plot='day'; bins=[20:.25:25]; 
%   ttl='day'; xlab='day'; ylabunit='day^{-1}';
% var2plot='Area'; bins=[0:200:5000]; 
%   ttl='Area'; xlab='km^2'; ylabunit='km^{-2}';
% var2plot='AspectRatio'; bins=[0:.05:1]; 
%   ttl='Aspect Ratio'; xlab='Aspect Ratio'; ylabunit='-';
% var2plot='MinThrhopert'; bins=[-8:.25:0]; 
%   ttl='Min \theta_{\rho}'''; xlab='\theta_{\rho}'''; ylabunit='K^{-1}';
% var2plot='MeanThrhopert'; bins=[-4:.1:0]; 
%   ttl='Mean \theta_{\rho}'''; xlab='\theta_{\rho}'''; ylabunit='K^{-1}';


% make plot directory if doesn't exist
plotdir=[var2plot,'_PDFs/'];
if ~exist(plotdir,'dir'); mkdir(plotdir); end

% info for idays
if idays==1
    ttlstr2 = ', days 22-25';
    prtstr2 = '_day22-25';
else
    ttlstr2 = '';
    prtstr2 = '_day20-25';
end

% figure handles
for r=1:nr
    
    idx=mod(r,2); if idx==0; idx=2; end % idx for run: either 1 or 2
    % clear figure if necessary
    if idx==1 % odd number: a new pair of simulations
        ifig=figure; clf;
        clear lgdstr lgd
    end
    
    % load file
    matfile2read = [fullname{r},'/',ftype,'files/vintcond_stats.mat'];
    load(matfile2read);
    % get variables
    if strcmp(var2plot,'AspectRatio')
        vararray=arrayfun(@(x) x.MinorAxisLength./x.MajorAxisLength, regstatsall);
    else
        vararray=arrayfun(@(x) x.(var2plot), regstatsall);
    end
    
    % read in area and only save features wtih area >= 100km2
    area = arrayfun( @(x) x.Area*dx*dx, regstatsall );
    iarea100km2=area>=100.; % indices for area >= 100 sq km
    vararray = vararray( iarea100km2 );
    area = area( iarea100km2 );
     
    % if idays=1, also get days and use data only for days>=22
    if idays==1
        day=arrayfun(@(x) x.day, regstatsall); % day (20-25)
        day = day( iarea100km2 );
        vararray = vararray( day>=22 );
        day = day( day>=22 );
        area = area( day>=22 );
    end
    
    % plot variable 1: max cloud top height
    figure(ifig);
    histogram(vararray,bins,'facealpha',.5,'normalization','pdf'); hold all;
    % legend
    lgdstr{idx} = labelname{r};
    % finish plot details
    if idx==2
        title([ttl,ttlstr2])
        lgd=legend(lgdstr,'location','northwest'); lgd.Box='off';
        xlabel(xlab)
        ylabel(['PDF (',ylabunit,')'])
        plotname=[var2plot,'_objects_area100km2_',labelname{r}(1:find(labelname{r}=='\',1)-1),prtstr2];
        print('-dpng',[plotdir,plotname])
        %saveas(gcf,[plotdir,plotname],'fig')
    end
    
end
