clear data
colorsb=parula(6);
colorsb=colorsb([1,3,5],:);
colorsl=brighten(colorsb,0.5);
colorsd=brighten(colorsb,-0.5);
field='c';
nvec=.999:-.001:.01;
proc='cc';
titles={'Cloud M3','Cloud M0','Cloud M6','Rain M0','Rain M6','Total M0','Total M6'};
letts='a':'f';

gtime=50; %tercile group evaluation time

if length(imomset)==3
    ncol=5;
else
    ncol=1;
end

for ifield=1:2
    if ifield==1; field='c';elseif ifield==2; field='r';else; field='f';end
    for imom=imomset%[3,0,6]
        if (strcmp(field,'f') || strcmp(field,'r')) && imom==3
            continue
        end
        jp=jp+1;
        
        nonparametric_error_evolution
        
        if imom==3
            
            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); %sort by error at the normtime=.75
            %     [datasort,idatasort]=sort((data(gtime,:,1)-data(gtime,:,2))); %sort by error at the normtime=.75
            
            %             inds1=idatasort(1:71);
            %             inds2=idatasort(72:144);
            %             inds3=idatasort(145:214);
            
            i1=round(lastreal/3);
            inds1=idatasort(1:i1);
            inds2=idatasort(i1+1:2*i1);
            inds3=idatasort(2*i1+1:lastreal);
        end
        
        subplot(ns,ncol,jp)
        hold on
        
        for ipc=1:3
            %plot(normtime,data(:,ip(ipc),1),'-','linewidth',3,'color',colors((ipc-1)*2+1,:))
            %plot(normtime,data(:,ip(ipc),2),'--','linewidth',3,'color',colors((ipc-1)*2+2,:))
            
        end
        %Should do this with a loop instead
        
        plot(normtime,nanmean(data(:,inds1,1),2),'-','linewidth',3,'color',colorsd(1,:))
        plot(normtime,nanmean(data(:,inds1,2),2),'--','linewidth',3,'color',colorsl(1,:))
        plot(normtime,nanmean(data(:,inds1,1)-data(:,inds1,2),2),':','color',colorsb(1,:),'linewidth',3)
        
        plot(normtime,nanmean(data(:,inds2,1),2),'-','linewidth',3,'color',colorsd(2,:))
        plot(normtime,nanmean(data(:,inds2,2),2),'--','linewidth',3,'color',colorsl(2,:))
        plot(normtime,nanmean(data(:,inds2,1)-data(:,inds2,2),2),':','color',colorsb(2,:),'linewidth',3)
        
        plot(normtime,nanmean(data(:,inds3,1),2),'-','linewidth',3,'color',colorsd(3,:))
        plot(normtime,nanmean(data(:,inds3,2),2),'--','linewidth',3,'color',colorsl(3,:))
        plot(normtime,nanmean(data(:,inds3,1)-data(:,inds3,2),2),':','color',colorsb(3,:),'linewidth',3)
        
        if jp>(ns-1)*ncol
            xlabel('Normalized Time')
        end
        if jp==1            
            legend('AMP','BIN','diff','Location','south','Fontsize',12,'Orientation','horizontal','NumColumns',2,'Color','none')
            legend('boxoff')
        end
        if mod(jp,ncol)==1 || ncol==1
            ylabel(['\bf',simname,10,'\rm Normalized Moment'])
        end
        if jp<=ncol;title(titles{jp});end
        set(gca,'xlim',[0 1])
        set(gca,'ylim',[-.4 1.2],'ytick',-.4:.4:1.2)
        grid on
        plotlabel(['(',letts(ins),num2str(mod(jp-1,ncol)+1),')'],'FontSize',12)
        
    end
end

%% Evaporation - error evolution and normalized evolution
% clear data
% colors=parula(9);%brighten(linspecer(9,'sequential'),-.4);
% field='c';
% nvec=.999:-.001:.01;
% proc='cc';
% %figure('position',[1 1 1125 675])
% figure('position',[1 1 930 300])
% letts='a':'f';
% j=0;
% % for imom=[0,3,6]
% %     j=j+1;
% %        
% %     nonparametric_error_evolution
% %     
% %     dummy=(data(:,:,1:end-1)-data(:,:,end))./data(:,:,end)*100;
% %     dummy=abs(dummy);
% %     everrs=squeeze(prctile(dummy,50,2));
% %     everrs25=squeeze(prctile(dummy,25,2));
% %     everrs75=squeeze(prctile(dummy,75,2));
% %     
% %     subplot(2,3,j)
% %     hold on
% %     
% %     inds=[1,25:25:900];
% %     tdat=1:length(everrs);
% %     plot(normtime,everrs,'r-','linewidth',3)
% %     hold on
% %     plot(normtime,everrs25,'r-','linewidth',3)
% %     plot(normtime,everrs75,'r-','linewidth',3)
% %         
% %     %set(gca,'xscale','linear','xdir','reverse','xlim',[0.1 1],'yscale','log');
% % 
% %     grid on
% %     if j==1
% %         title('0th Moment')
% %         ylabel(['Percent Error',10,'25,50,75^{th} Percentiles'])
% %     elseif j==2
% %         %
% %         title('3rd Moment')
% %     elseif j==3
% %         title('6th Moment')
% %     end
% %     %plotlabel(['(',letts(j),')'],'upperleft')
% %     xlabel('Normalized Time')
% %     set(gca,'yscale','log')
% %     set(gca,'ylim',[.1 1e4])
% % end
% 
% clear data
% j=0;
% for imom=[0,3,6]
%     j=j+1;
%     
%     nonparametric_error_evolution
%     
%     dummy=data;
%     everrs=squeeze(prctile(dummy,50,2));
%     everrs25=squeeze(prctile(dummy,25,2));
%     everrs75=squeeze(prctile(dummy,75,2));
%     
%     subplot(1,3,j)
%     hold on
% 
%     tdat=1:length(everrs);
%     
%     for ipc=[25,50,75]
%        % if imom==3
%        %     plot(normtime,1-squeeze(prctile(dummy(:,:,1),ipc,2)),'b-','linewidth',3)
%        %     plot(normtime,1-squeeze(prctile(dummy(:,:,2),ipc,2)),'k-','linewidth',3)
%        % else
%             plot(normtime,squeeze(prctile(dummy(:,:,1),ipc,2)),'b-','linewidth',3)
%             plot(normtime,squeeze(prctile(dummy(:,:,2),ipc,2)),'k-','linewidth',3)
%        % end
%     end
%     
%     %set(gca,'xscale','linear','xdir','reverse','xlim',[0 1],'xtick',0:.2:1);
%     xlabel('Normalized Time')
%     if j==1
%         ylabel(['Normalized Moment Value',10,'25,50,75^{th}',' Percentiles'])
%         legend('AMP','BIN')
%     end
%     set(gca,'xlim',[0 1])
%     grid on
%      if imom==0 || imom==3
%          set(gca,'ylim',[0 1])
%      end
%     %plotlabel(['(',letts(j),')'],'upperleft')
% end
% 
% %minspace(gcf,2,3)
% %%
% figure
% % subplot(1,3,1)
% % plot(mom(1).c(:,:,2)./mom(1).c(1,:,2),mom(1).c(:,:,1)./mom(1).c(1,:,2),'-','linewidth',2)
% % set(gca,'xdir','reverse')
% % hold on
% % plot([1e-3 1.2],[1e-3 1.2],'k-')
% % 
% % subplot(1,3,2)
% plot(mom(4).c(:,:,2)./mom(4).c(1,:,2),mom(4).c(:,:,1)./mom(4).c(1,:,2),'-','linewidth',1)
% set(gca,'xdir','reverse')
% hold on
% plot([1e-3 1],[1e-3 1],'k-')
% set(gca,'xlim',[0 1],'ylim',[0 1])
% xlabel('Normalized BIN Mass')
% ylabel('Normalized AMP Mass')
% 
% % subplot(1,3,3)
% % plot(mom(7).c(:,:,2)./mom(7).c(1,:,2),mom(7).c(:,:,1)./mom(7).c(1,:,2),'-','linewidth',2)
% % set(gca,'xdir','reverse')
% % hold on
% % plot([1e-3 1.2],[1e-3 1.2],'k-')
% 
% 
