%distinct colors
colors=[192 192 192;
    47 79 79;
    139 69 19;
    0 100 0;
    128 128 0;
    72 61 139;
    0 0 128;
    154 205 50;
    139 0 139;
    255 69 0;
    255 165 0;
    255 255 0;
    127 255 0;
    0 250 154;
    138 43 226;
    220 20 60;
    0 255 255;
    0 191 255;
    0 0 255;
    240 128 128;
    255 0 255;
    30 144 255;
    240 230 140;
    221 160 221;
    255 20 147]/256;
    
    
%%
%Sort simulations by error in c038-r038
field='c';
nvec=.999:-.001:.01;
proc='cc';
testset='../AMP_45_NP/output_3M_gamma/cc';
testsetb=testset;
vars{1}='cM8M8type1';
vars{2}='bin';

read_nonparametric
imom=3;
nonparametric_error_evolution

lastreal = find(~isnan(data(end,:,1)),1,'last');
datadiff=nan(1,ntest);
for ii=1:lastreal
    errinds = find(data(:,ii,2)>=0.5,1,'last');
    datadiff(ii)=data(errinds,ii,1)-data(errinds,ii,2);
end
[datasort,idatasort]=sort(datadiff);

%% AMP-F 6M
nu=1:2:15;
N=50*2.^[0:5];
M=1:5;

nkr=33;%# of bins
col=log(2)/3; %bin width
load masses_radii
diams=radii*2*1e-2; %convert cm to m

xl=3.14/6.*1000*diams.^3; %mass of one droplet in each bin

delimiter = ' ';

testset='../AMP_45_NP/output_6M_gamma/cc';
formatSpec2 = '%f%f%f%f%f%f%[^\n\r]';

vars{1}='cparamsM2467type0';
momz=6;
itest=0;
params=nan(6,1800,240);
for g=1:length(nu)
    for n=1:length(N)
        for m=1:length(M)            
            itest=itest+1;
            %Read data
            folder=[testset,'/nu',num2str(g),'_N',num2str(n),'_M',num2str(m),'/'];            
            filename = [folder,vars{1},'.txt'];            
            fileID = fopen(filename,'r');            
            dataArray = textscan(fileID, formatSpec2, 'Delimiter', delimiter, ...
                'MultipleDelimsAsOne', true,  'ReturnOnError', false);
            fclose(fileID);
            
            dataArray=[dataArray{1:6}];
            params(:,1:length(dataArray),itest) = dataArray';
            
            for itime=1:1800
                %Get distribution parameters
                nu1=params(3,itime,itest);
                nu2=params(5,itime,itest);
                dn1=params(4,itime,itest);
                dn2=params(6,itime,itest);
                
                %Create mass distributions
                for kr=1:nkr
                    gterm1a(kr)=-diams(kr)/dn1+(nu1+3)*log(diams(kr)/dn1);
                    gterm2a(kr)=-diams(kr)/dn2+(nu2+3)*log(diams(kr)/dn2);
                end
                
                gterm1b = max(gterm1a)*-1.;
                gterm2b = max(gterm2a)*-1.;
                
                for kr=1:nkr
                    md1(kr)=exp(gterm1a(kr)+gterm1b);
                    md2(kr)=exp(gterm2a(kr)+gterm2b);
                end
                
                %Calculate 3rd, 0th, and 1st moments of each mode
                m31=sum(md1./xl.*diams.^3)*col*1000.;
                m32=sum(md2./xl.*diams.^3)*col*1000.;
                
                m01=sum(md1./xl)*col*1000.;
                m02=sum(md2./xl)*col*1000.;
                
                m11=sum(md1./xl.*diams)*col*1000.;
                m12=sum(md2./xl.*diams)*col*1000.;
                
                %Calculate mean diameter of each mode
                nm1=m11/m01;
                nm2=m12/m02;
                
                %calculate standard deviation of diameter
                sq1=sqrt(sum((diams-nm1).^2.*md1./xl*col*1000)/m01);
                sq2=sqrt(sum((diams-nm2).^2.*md2./xl*col*1000)/m02);
                
                %relative dispersion
                rd1=sq1/nm1;
                rd2=sq2/nm2;
           
                %Figure out which is the smaller mode
                if nm1>nm2
                    params(1,itime,itest)=nm2;
                    params(2,itime,itest)=nm1;
                    params(3,itime,itest)=rd2;
                    params(4,itime,itest)=rd1;
                    params(5,itime,itest)=nu2;
                    params(6,itime,itest)=nu1;
                else 
                    params(1,itime,itest)=nm1;
                    params(2,itime,itest)=nm2;
                    params(3,itime,itest)=rd1;
                    params(4,itime,itest)=rd2;
                    params(5,itime,itest)=nu1;
                    params(6,itime,itest)=nu2;
                end
                
            end
            info.nu(itest)=g;
            info.N(itest)=n;
            info.M(itest)=m;
            
        end
    end
end
ntest=itest;

meand1=squeeze(params(1,:,:))*1e6;
meand2=squeeze(params(2,:,:))*1e6;
nu1=squeeze(params(5,:,:));
nu2=squeeze(params(6,:,:));
nu1(abs(nu1)>100)=NaN;
nu2(abs(nu2)>100)=NaN;
Nu1=nu1;Nu2=nu2;
Meand1=meand1;
Meand2=meand2;

%Do some data cleaning. Want to remove spurious time steps when the
%parameters are very different than parameters before or after. Set DSD
%properties at these time steps to NaN. Thresholds used here were based on
%inspection of the data.
dd=[zeros(1,240);diff(meand1)];
dn=[zeros(1,240);diff(nu1)];

dd2=[zeros(1,240);diff(meand2)];
dn2=[zeros(1,240);diff(nu2)];

nu1=Nu1;nu2=Nu2;
meand1=Meand1;
meand2=Meand2;

mask=abs(dd)>.1 | abs(dn)> .1;
nu1(mask==1)=NaN;
meand1(mask==1)=NaN;
meand1(nu1<1)=NaN;
nu1(nu1<1)=NaN;
nu1(meand1<8)=NaN;
meand1(meand1<8)=NaN;

mask=abs(dd2)>5 | abs(dn2)> 1;
nu2(mask==1)=NaN;
meand2(mask==1)=NaN;
meand2(nu2<0)=NaN;
nu2(nu2<0)=NaN;
nu2(meand2>251)=NaN;
meand2(meand2>251)=NaN;

ee=0.9;
nexs=25;

%Mode 1
f1=figure;
subplot(2,2,1)
hold on
for i=1:nexs
    %Use the sorted simulations now. Want the figure to include highlight
    %simulations that were bad to mediocre. Skipping every 4th sorted
    %simulation:
    il=idatasort(213-4*i);
    ist=indstart(il);
    ied=indend(il);
    ied=round(ee*(ied-ist))+ist;

    %confusing indexing
    il=dataIT(il);
    if ist>0 && ied>0
         nu1p=nu1(ist:ied,il);
         x=ist:ied;
         %Interpolate over the NaNs before plotting
         xp=x(~isnan(nu1p))';
         nu1pp=nu1p(~isnan(nu1p));
         nu1p=interp1(xp,nu1pp,x,'linear',NaN);
         
         meand1p=meand1(ist:ied,il);
         x=ist:ied;
         xp=x(~isnan(meand1p))';
         meand1pp=meand1p(~isnan(meand1p));
         meand1p=interp1(xp,meand1pp,x,'linear',NaN);
         
         ind0=find(~isnan(meand1p),1,'first');
         indf=find(~isnan(meand1p),1,'last');
   
         plot(meand1p,nu1p,'-','linewidth',3,'color',colors(i,:))
         plot(meand1p(ind0),nu1p(ind0),'o','linewidth',3,'color',colors(i,:))
         plot(meand1p(indf),nu1p(indf),'x','linewidth',3,'color',colors(i,:),'MarkerSize',8)
    end
end
xlabel('Mean Diam. (\mum)')
ylabel('Shape Parameter')
title('(a) Mode 1')
set(gca,'xlim',[5 35],'ylim',[-2 15])

%Mode 2
subplot(2,2,2)
hold on
for i=1:nexs
    il=idatasort(213-4*i);
    ist=indstart(il);
    ied=indend(il);
    ied=round(ee*(ied-ist))+ist;

    %confusing indexing
    il=dataIT(il);
    if ist>0 && ied>0
         nu1p=nu2(ist:ied,il);
         x=ist:ied;
         xp=x(~isnan(nu1p))';
         nu1pp=nu1p(~isnan(nu1p));
         nu1p=interp1(xp,nu1pp,x,'linear',NaN);
         
         meand1p=meand2(ist:ied,il);
         x=ist:ied;
         xp=x(~isnan(meand1p))';
         meand1pp=meand1p(~isnan(meand1p));
         meand1p=interp1(xp,meand1pp,x,'linear',NaN);
         
         ind0=find(~isnan(meand1p),1,'first');
         indf=find(~isnan(meand1p),1,'last');
   
         plot(meand1p,nu1p,'-','linewidth',3,'color',colors(i,:))
         plot(meand1p(ind0),nu1p(ind0),'o','linewidth',3,'color',colors(i,:))
         plot(meand1p(indf),nu1p(indf),'x','linewidth',3,'color',colors(i,:),'MarkerSize',8)
    end
end
xlabel('Mean Diam. (\mum)')
ylabel('Shape Parameter')
title('(b) Mode 2')
set(gca,'xlim',[0 250],'ylim',[-2 15])

%% AMP-CR c038-r038
nu=1:2:15;
N=50*2.^[0:5];
M=1:5;

nkr=33;%# of bins
col=log(2)/3; %bin width
load masses_radii
diams=radii*2*1e-2; %convert cm to m

xl=3.14/6.*1000*diams.^3; %mass of one droplet in each bin

delimiter = ' ';

testset='../AMP_45_NP/output_3M_gamma/cc';
formatSpec2 = '%f%f%f%f%f%f%[^\n\r]';

%Read gamma parameter values from a file
vars{1}='cparamsM8M8type1';
momz=6;
itest=0;
params=nan(6,1800,240);
for g=1:length(nu)
    for n=1:length(N)
        for m=1:length(M)
            
            itest=itest+1;
            
            folder=[testset,'/nu',num2str(g),'_N',num2str(n),'_M',num2str(m),'/'];            
            filename = [folder,vars{1},'.txt'];            
            fileID = fopen(filename,'r');            
            dataArray = textscan(fileID, formatSpec2, 'Delimiter', delimiter, ...
                'MultipleDelimsAsOne', true,  'ReturnOnError', false);
            fclose(fileID);
            
            dataArray=[dataArray{1:6}];
            params(:,1:length(dataArray),itest) = dataArray';
            
            for itime=1:1800
                nu1=params(2,itime,itest);
                nu2=params(5,itime,itest);
                dn1=params(3,itime,itest);
                dn2=params(6,itime,itest);
                
                %Some gamma math
                for kr=1:13
                    gterm1a(kr)=-diams(kr)/dn1+(nu1+3)*log(diams(kr)/dn1);
                end
                for kr=14:33
                    gterm2a(kr)=-diams(kr)/dn2+(nu2+3)*log(diams(kr)/dn2);
                end
                
                gterm1b = max(gterm1a)*-1.;
                gterm2b = max(gterm2a)*-1.;
                
                %Create mass distributions
                md1(1:nkr)=0;
                md2(1:nkr)=0;
                for kr=1:13
                    md1(kr)=exp(gterm1a(kr)+gterm1b);
                end
                for kr=14:nkr
                    md2(kr)=exp(gterm2a(kr)+gterm2b);
                end
                
                %3rd moment for modes 1 and 2. Corresponds to all 33 bins
                %to make it comparable to AMP-F in which modes could
                %overlap
                m31=sum(md1./xl.*diams.^3)*col*1000.;
                m32=sum(md2./xl.*diams.^3)*col*1000.;
                
                %0th moment for modes 1 and 2. 
                m01=sum(md1./xl)*col*1000.;
                m02=sum(md2./xl)*col*1000.;
                
                %1st moment for modes 1 and 2.
                m11=sum(md1./xl.*diams)*col*1000.;
                m12=sum(md2./xl.*diams)*col*1000.;
                
                %Mean diameter: M1/M0
                nm1=m11/m01;
                nm2=m12/m02;
                
                %Standard deviation
                sq1=sqrt(sum((diams-nm1).^2.*md1./xl*col*1000)/m01);
                sq2=sqrt(sum((diams-nm2).^2.*md2./xl*col*1000)/m02);
                
                %Relative dispersion
                rd1=sq1/nm1;
                rd2=sq2/nm2;
                
                %Mode 1 is definitely the cloud mode this time. No IF
                %statement needed
                params(1,itime,itest)=nm1;
                params(2,itime,itest)=nm2;
                params(3,itime,itest)=rd1;
                params(4,itime,itest)=rd2;
                params(5,itime,itest)=nu1;
                params(6,itime,itest)=nu2;
                
            end
            info.nu(itest)=g;
            info.N(itest)=n;
            info.M(itest)=m;
            
        end
    end
end
ntest=itest;

meand1=squeeze(params(1,:,:))*1e6;
meand2=squeeze(params(2,:,:))*1e6;
nu1=squeeze(params(5,:,:));
nu2=squeeze(params(6,:,:));
nu1(abs(nu1)>100)=NaN;
nu2(abs(nu2)>100)=NaN;

%Don't need to do any data cleaning this time. Just plot it.
figure(f1)
subplot(2,2,3)
hold on
for i=1:nexs
    il=idatasort(213-4*i);
    ist=indstart(il);
    ied=indend(il);
    ied=round(ee*(ied-ist))+ist;
    
    %confusing indexing
    il=dataIT(il);
    if ist>0 && ied>0  
         plot(meand1(ist:ied,il),nu1(ist:ied,il),'-','linewidth',3,'color',colors(i,:))
         plot(meand1(ist,il),nu1(ist,il),'o','linewidth',3,'color',colors(i,:))
         plot(meand1(ied,il),nu1(ied,il),'x','linewidth',3,'color',colors(i,:),'MarkerSize',8)
    end
end
xlabel('Mean Diam. (\mum)')
ylabel('Shape Parameter')
title('(c) Cloud')
set(gca,'xlim',[5 35],'ylim',[-2 15])

subplot(2,2,4)
hold on
for i=1:nexs
    il=idatasort(213-4*i);
    ist=indstart(il);
    ied=indend(il);
    ied=round(ee*(ied-ist))+ist;
    
    %confusing indexing
    il=dataIT(il);
    if ist>0 && ied>0
         plot(meand2(ist:ied,il),nu2(ist:ied,il),'-','linewidth',3,'color',colors(i,:))
         plot(meand2(ist,il),nu2(ist,il),'o','linewidth',3,'color',colors(i,:))
         plot(meand2(ied,il),nu2(ied,il),'x','linewidth',3,'color',colors(i,:),'MarkerSize',8)
    end
end
xlabel('Mean Diam. (\mum)')
ylabel('Shape Parameter')
title('(d) Rain')
set(gca,'xlim',[0 250],'ylim',[-2 15])

