% Reads in W-Band nc moment data and plot
% range gates are variable
% Version_3 and above read kongsberg motion data and corrects vertical velocity
% Use jamn for multiple files in 1 hour.
% Preallocate arrays to combine full-day time series.

close all
clear ddd ehr shr hh
%################## Cruise informations  ##################################################################
cruise='HFIP_13';             % Cruise name  
year='2013';               % year (deprecated--make year into a double) 
yyyy='2013';               % year (string)
ship='P3';               % Research vessel
big_dir='H:\';
big_dir='C:\Data\cwf\Dropbox\DATA\HFIP_13\';

disp('fockewe2')
% W-band radar files
way_raw_data_wband=[big_dir 'Raw\'];
way_raw_images_wband=[big_dir 'Raw_images\'];
way_raw_motion_wband=[big_dir 'Raw\Motion\'];
%way_cloud_data_wband=[big_dir '\RHB\reports\'];
%##################%##################%##################%##################%##################%##################

graphformat='.png';  %select graphics format files
graphdevice='-dpng'; %select graphic device
prtit=true;  % for printing plots--ineffective unless plotit is also true going to a temp directory
prtit_hr=1;  %to save individual hourly graphs
modeflag=0;
motionread=1;%set to zero to skip motion correction
% ceilo=1;  %to add ceilo data to the Wband plots (note you need to run the ceilo codes first)
% micro=1;  %to add mailbox data to the Wband plots (note you need to run the mailbox codes first)
Matlab_netcdf=1;


%%%%%%%%%%%%%%%%%%%%%%%%% Starts process %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% use known day/hour; prompt for it if unknown
if ~exist('ddd','var')
    ddds=input('Input yearday  ','s');
    ddd=str2double(ddds);
else
    ddds=sprintf('%03i',ddd);
end
[m,d]=yd2md(str2double(yyyy), ddd);
Vdate=[str2double(yyyy) m d];
if ~exist('shr','var')
%     if 0
        shrs=input('Input start hr  ','s');
        shr=str2double(shrs);
    else
        shr=0;
    end
% end
if ~exist('ehr','var')
%     if 0
        ehrs=input('Input end hr  ','s');
        ehr=str2double(ehrs);
    else % use hours 0-23 by default
        ehr=23;
    end
% end

cldtop=[];
% preallocate x-extended variables
nxmax=175000;
Nrange=220; % may only use 190 of the ranges or less or more
Refx=deal(single(zeros(Nrange,nxmax))+NaN);
%[SNRx, SWx, DVx]=deal(Refx); % not used
[hhx, rangemode]=deal(single(zeros(1,nxmax))+NaN);

Refx=deal(single(zeros(Nrange,nxmax))+NaN);
[Refh,SNRh, SWh, DVah,]=deal(single(zeros(Nrange,15000))+NaN);
[hh]=deal(single(zeros(1,15000))+NaN);


handl(1)=figure(1); % open 2 figures whether anything to plot
handl(2)=figure(2);
% set(handl(:),'WindowStyle','docked');

nx=0; % count position in extended variables
for jam=shr:ehr
    isfile=false;
    clf(handl(:));
    nh=0;
    clear y z dt jd
    hh=[]; Ref=[]; DVa=[]; SW=[]; SNR=[]; nrg=[];
    hr=sprintf('%02i',jam);
    dr=dir([way_raw_data_wband yyyy ddds hr '*MMCRMom.nc']);
    for jamn=1:length(dr) % multiple files in one hour
        fn=[way_raw_data_wband dr(jamn).name];
        disp(['Reading file: ' fn])
        try
            if Matlab_netcdf
                FID=netcdf.open(fn,'NC_NOWRITE');
            else
                FID=fopen(fn,'r');
            end
        catch
            FID=-1;
        end
        
        if FID>0  % Test to see if file exists
            isfile=true;
            if Matlab_netcdf
                base_time=netcdf.getVar(FID,0);
                time_offset=netcdf.getVar(FID,1)';
            else
                ncload(fn);
            end
            % calculate time from time_offset and base_time
            tt=time_offset+double(base_time);
            nta=length(tt); % time length of current file
            zc=tt/(60*60*24)+datenum([1970,1,1,0,0,0]);
            jd=zc-datenum(str2double(yyyy)-1,12,31);
            %         dt = datestr(z, 31);    %2008-10-07 16:35:53
            %         [yr, month, day, hour, minutes, secs] = datevec(dt);
           
            hh0=(jd-(fix(jd(1))))*24; % subtract off JD and convert to an hour
            
            if Matlab_netcdf
                z=netcdf.getVar(FID,27);% 27
                ht=z(:,1)/1e3;
                Reflectivity=netcdf.getVar(FID,40)';  %40
                SignalToNoiseRatio=netcdf.getVar(FID,36)'; %36
                SpectralWidth=netcdf.getVar(FID,38)';  % 38
                MeanDopplerVelocity=netcdf.getVar(FID,35)'; %35
                
                Ref=single(Reflectivity)';
                SW=single(SpectralWidth)';
                DVa=single(MeanDopplerVelocity)';
                SNR=single(SignalToNoiseRatio)';
                nrg=netcdf.getVar(FID,18);
                    Nrange=nrg(1)
                netcdf.close(FID);              
            else
                z=single(Heights);
                ht=z(1,:)/1e3;
                Nrange=single(NumHeights(1));
                Ref=single(Reflectivity)';
                SW=single(SpectralWidth)';
                DVa=single(MeanDopplerVelocity)';
                SNR=single(SignalToNoiseRatio)';
                fclose(FID);
            end
        end % test for file
       
        %concatenate multiples files into one hour    
        [m0,n0]=size(Ref);
        hh(nh+(1:nta))=hh0;
%         hth(1:m0)=ht;
        Refh(1:m0,nh+(1:nta))=Ref;
        DVah(1:m0,nh+(1:nta))=DVa;
        SWh(1:m0,nh+(1:nta))=SW;
        SNRh(1:m0,nh+(1:nta))=SNR;
        Nrangeh(nh+(1:nta))=Nrange;
        nh=nh+nta;    
        
    end % one or multiple files in hour
     
    % truncate hourly arrays
    clear hh0 Ref DVa SW SNR 
    Ref=Refh(1:Nrange,1:nh);
    DVa=DVah(1:Nrange,1:nh);
    SW=SWh(1:Nrange,1:nh);   
    SNR=SNRh(1:Nrange,1:nh); 
%     ht=hth(1:nh);
  %  clear Refh DVah SWh SNRh hth
    [nx,mx]=size(Ref);jj=1:nx;
    ht=(z(1)/1e3+0.025*jj)';
    dbz0=-43;
    RangeRes=ht(2)-ht(1);
    if RangeRes>0.04;
        dbz0=-51;
    end

    % store 2 height variables for the 2 modes if necessary
    if length(ht)==220
        ht1=ht;
    elseif length(ht)==190
        ht0=ht;
    end


    if motionread %correct doppler for ship motion
    %                 if exist([way_raw_motion_wband yyyy ddds hr '00Kongsberg.txt'],'file');
       if ~isempty(dir([way_raw_motion_wband '2013' ddds hr '*Kongsberg.txt'])); %no files at all
            read_Wband_motion_P3_lb
            pitchk=vkonsa(7,:);
            rollk=vkonsa(8,:);
            try
                DV_mot=single(interp1(hhka,vkonsa(9,:),hh));
            catch
                DV_mot=0*DVa(1,:);% interpolation failed
            end
            DV_mot(isnan(DV_mot))=0;
            DV=DVa-repmat(DV_mot,Nrange,1);
            Konmotion=1;
        else
            DV=DVa;% no motion file
            Konmotion=0;
            [hhka,pitchk, rollk]=deal(single(zeros(1,36000))+NaN);
        end
    else
        DV=DVa; %motion correction turned off
    end
    DV=single(DV);
    
  
    
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %do hourly plots
             ymaxx=ht(Nrange);
            zmxx=ymaxx;
            
            figure(1); %clf
            subplot(3,1,1),
           imagesc(hh,ht, Ref(1:Nrange,:),[-60, 20]);
            axis xy; axis([jam jam+1 0 ymaxx]);
            ylabel('Ht (km)')
            %         tit={[cruise ' W-Band'' YD ' ddds ' Hr ' hr]};
            %         title(tit,'Interpreter','none')
            title(sprintf('%s (%04i-%02i-%02i, DOY%03i, Hr-%s).  W-Band (motionread=%i, Kongsberg=%i)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr,motionread,Konmotion),'FontWeight','Bold','Interpreter','none')
            h=colorbar;
            ylabel(h,'Ref dBZ') 
            hold on
            
            subplot(3,1,2)
            imagesc(hh,ht, DV(1:Nrange,:),[-7, 7]);  %
            axis xy;axis([jam jam+1 0 ymaxx]);
            ylabel('Ht (km)')
            h=colorbar;
            ylabel(h,'-DopVel (m/s)')
            hold on
            
            subplot(3,1,3),imagesc(hh,ht, SW(1:Nrange,:),[0, 5]);  %
            axis xy;axis([jam jam+1 0 ymaxx]);
            ylabel('Ht (km)')
            xlabel('Hrs (UTC)')
            h=colorbar;
            ylabel(h,'SpecWid (m/s)')
            % print('-djpeg90 ',[way_raw_images_wband '\Wband_ref_' num2str(ddd), '_' hr '.jpg']);
            annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
            hold on
            if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_Wband_ref_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]);end
            
            figure(2); %clf
            subplot(3,1,1),imagesc(hh,ht, SNR(1:Nrange,:),[-25, 20]);  %
            axis xy; axis([jam jam+1 0 ymaxx]);
            ylabel('Ht (km)')
            %axis([15 15.25 0 3])
            %         tit={[cruise ' W-Band'' YD ' ddds ' Hr ' hr]};
            %         title(tit,'Interpreter','none')
            title(sprintf('%s (%04i-%02i-%02i, DOY%03i, Hr-%s).  W-Band (motionread=%i, Kongsberg=%i)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr,motionread,Konmotion),'FontWeight','Bold','Interpreter','none')
            h=colorbar;
            ylabel(h,'SNR dBZ')
            hold on

            subplot(3,1,2),imagesc(hh,ht, DV(1:Nrange,:),[-3, 5]);  %
            axis xy;axis([jam jam+1 0 ymaxx]);
            ylabel('Ht (km)')
            %axis([15 15.25 0 3])
            h=colorbar;
            ylabel(h,'DopVel (m/s)')
            hold on
            %         subplot(3,1,3),imagesc(hh,ht, DVa(1:120,:),[-3, 5]);  %
            %         axis xy
            %         ylabel('Ht (km)')
            %         %axis([15 15.25 0 3])
            %         h=colorbar;
            %         ylabel(h,'DopVel (m/s)')

            subplot(3,1,3),imagesc(hh,ht, SW(1:Nrange,:),[0, 5]);  %
            axis xy;axis([jam jam+1 0 ymaxx]);
            ylabel('Ht (km)')
            xlabel('Hrs (UTC)')
            h=colorbar;
            ylabel(h,'SpecWid (m/s)')
            %print('-djpeg90 ',[way_raw_images_wband '\Wband_snr_' num2str(ddd), '_' hr '.jpg']);
            annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
            hold on
            if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_Wband_snr_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]);                 end
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
     
    %reads in the P3 data
    read_P3navdata_lb2;
    figure(4);clf;
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Pitch/Roll',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    subplot(1,2,1);plot(P3_hh,P3_pitchref,'r',hhka,pitchk);grid
    xlabel('Hour (UTC)');ylabel('Pitch');legend('P3','kongsberg')
    subplot(1,2,2);plot(P3_hh,P3_rollref,'r',hhka,rollk);grid
    xlabel('Hour (UTC)');ylabel('Roll');legend('P3','kongsberg')
    annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
    if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_pitch_roll_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]);end

    
    
    

     clear Reflectivity SpectralWidth MeanDopplerVelocity SignalToNoiseRatio Power NoiseLevel SkyNoiseLevel vkons CircularDepolarizationRatio
    clear RangeCorrectedPower AvgNoiseLevel Azimuth
   

     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %insert data in extended daylong array(s)
            [m1,n1]=size(Ref);
            if m1>0 && n1>0
                hhx(:,nx+(1:n1))=single(hh);
                Refx(1:m1,nx+(1:n1))=single(Ref);
                rangemode(1,nx+(1:n1))=m1==220; % rangemode=1 for 220 gates, 0 for 190 gates
                nx=nx+nta;
            end
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                   
    if isfile
        % hourly height statistics
        [htt reft reft2 wvel nnt]=deal(single(zeros(1,size(Ref,2))));
        for k=1:size(Ref,2); % time dimension
            kk=Ref(:,k)-dbz0-20*log10(ht)>8 & ht>.1 & ht<zmxx & Ref(:,k)>-30; % logical for noise&height in bounds
            if sum(kk)
                htt(k)=single(ht(find(kk,1,'last')));
                reft(k)=single(nanmean(10.^(Ref(kk,k)/10)));
                reft2(k)=single(nanmean(10.^((Ref(kk,k)-dbz0-20*log10(ht(kk)))/10)));
                wvel(k)=single(nanmean(DV(kk,k)));
                nnt(k)=single(sum(kk));
            else
                htt(k)=NaN;
                reft(k)=NaN;
                reft2(k)=NaN;
                wvel(k)=NaN;
                nnt(k)=0;
            end
        end
        
        y=prctile(htt,[10 15 50 85 90]);
        kk=htt>=y(1) & htt<=y(5); % times when cloud top in bounds
        httb=single(nanmean(htt(kk)));
        reftb=single(10*log10(nanmean(reft(kk))));
        reft2b=single(10*log10(nanmean(reft2(kk))));
        wvelb=single(nanmean(wvel(kk)));
        nntb(1,1:Nrange)=single(sum(nnt));
        sigh=single(sqrt(nanmean(htt(kk).^2)-httb.^2));
        disp(['Hour= ' num2str((hh(1))) '  Mean cloud top ht= ' num2str(httb) '  std cloud top height=' num2str(sigh)])
        disp(['Mean dBZ= ' num2str(reftb) '   Mean dBZ-noisedBZ= ' num2str(reft2b) '   Mean Vel= ' num2str(wvelb)])
        cldtop=[cldtop' [hh(1) httb sigh (y(4)-y(3))/2 reftb reft2b wvelb nntb]']';
        %hhx=[single(hhx) single(hh)];
    end % if there is a file
    
%      clear hh Ref DV SW SNR  vkonsa
%      clear hhka pitchk rollk
     
end % hr loop

% % truncate extended arrays
% Refx=Refx(:,1:nx);
% hhx=hhx(:,1:nx);
% rangemode=rangemode(:,1:nx);

roll3 = interp1(P3_hh,P3_rollref,hh);alt3=interp1(P3_hh,P3_altref,hh);pitch3 = interp1(P3_hh,P3_pitchref,hh);
Us3 = interp1(P3_hh,P3_sfmrU,hh);Rs3=interp1(P3_hh,P3_sfmrR,hh);cog3=interp1(P3_hh,P3_cog,hh);sog3=interp1(P3_hh,P3_GS,hh);
lon3=interp1(P3_hh,P3_lon,hh);lat3=interp1(P3_hh,P3_lat,hh);wdir3=interp1(P3_hh,P3_WindDir,hh);hdir3=interp1(P3_hh,P3_Hed,hh);
ff=pi/180;
Gn=sog3.*cos(cog3*ff);Ge=sog3.*sin(cog3*ff);
Uwn=Us3.*cos((wdir3+180)*ff);Uwe=Us3.*sin((wdir3+180)*ff);
Dopx=(sin(hdir3*ff).*sin(roll3*ff)+cos(hdir3*ff).*sin(pitch3*ff).*cos(roll3*ff)).*(Gn-Uwn);
Dopy=(sin(hdir3*ff).*sin(pitch3*ff).*cos(roll3*ff)-cos(hdir3*ff).*sin(roll3*ff)).*(Ge-Uwe);

refcorr=4.2;%+1.6
%  1.6 is loss of window, Attenuation is est as 4.2 dB; ken says 1.6 is
%  already in the files
[nx,mx]=size(Ref);
dnzx=[];zzx=[];snrx=[];dopx=NaN*ones(4,mx);widx=NaN*ones(4,mx);dnrx1=[];reftz=[];dnrx2=[];dnrx3=[];zrx=[];dnrx1m=[];snrz=[];usfm=[];latm=[];lonm=[];swrx1=[];refzx1=[];reftz=[];
dij=3;
for i=1:mx
    if Nrangeh(i)==150
        rng=alt3(i)./cos(roll3(i)*pi/180);
        kk=find(ht*1e3>rng-150 & ht*1e3<rng+150);
        rmx=max(Ref(kk,i));
        if rmx<100 & rmx>0;
            ij=find(Ref(kk,i)==rmx);jk=kk(ij);dnzx(i)=rmx;zzx(i)=ht(kk(ij));snrx1(:,i)=SNR(jk-10-dij:jk-dij,i);swrx1(:,i)=SW(jk-10-dij:jk-dij,i);swrx1(:,i)=SW(jk-10-dij:jk-dij,i);zrx=1e3*(ht(jk)-ht(jk-10-dij:jk-dij));dnrx1(:,i)=Ref(jk-10-dij:jk-dij,i)+refcorr;Vdopx1(:,i)=DVa(jk-10-dij:jk-dij,i);dnrx2(i)=Ref(jk-10,i)+refcorr;dnrx3(i)=Ref(jk-15,i)+refcorr;refzx1(:,i)=10.^((Ref(jk-10-dij:jk-dij,i)+refcorr)/10);reftz(:,i)=nanmean(10.^((Ref(jk-10-dij:jk-dij,i)+refcorr)/10));snrz(i)=SNR(jk-3,i);dopx(:,i)=DV(kk(ij-2:ij+1),i);widx(:,i)=SW(kk(ij-2:ij+1),i);
            %usfm(i)=Us3(jk);latm(i)=lat3(jk);lonm(i)=lon3(jk);
        else
             dnzx(i)=NaN;zzx(i)=NaN;snrx(i)=NaN;%dopx(:,i)=NaN;
        end;
    else
        dnzx(i)=NaN;zzx(i)=NaN;snrx(i)=NaN;%dopx(:,i)=NaN;
    end;
end;
ii=find(abs(roll3)<5 & abs(pitch3)<3);P = polyfit(10*log10(reftz(ii))-min(10*log10(reftz)),max(dnzx)-dnzx(ii)-4,1);atten=polyval(P,10*log10(reftz)-min(10*log10(reftz)));
P2 = polyfit((reftz(ii))-min((reftz)),max(dnzx)-dnzx(ii)-4,1);atten2=polyval(P2,(reftz)-min((reftz)));
atten2=max(atten2-1.5,0)/.434;
ii=find(P3_lat<50 & P3_lon<-90);
figure;
subplot(2,1,1);scatter(P3_lon(ii),P3_lat(ii),20,P3_sfmrU(ii),'o');
colorbar;
subplot(2,1,2);scatter(P3_lon(ii),P3_lat(ii),20,P3_WindDir(ii),'o');
colorbar;

snr_th=-12;

ii=find(dnrx3<-18 & snrz>snr_th & Us3>15 & abs(roll3)<3 & atten<5);dnrx1m=mean(dnrx1(:,ii)');snrx1m=mean(snrx1(:,ii)');Um=mean(Us3(ii));lonm=mean(lon3(ii));latm=mean(lat3(ii));np=length(ii);ij=find(dnrx3<-18  & Us3>15);nt=length(ij);
dirm=mean(wdir3(ii));hedm=mean(hdir3(ii));swrx1m=mean(widx(:,ii)');refzx1m=mean(refzx1(:,ii)');
figure;plot(mean(dnrx1(:,ii)'),zrx','-o',10*log10(refzx1m),zrx','x');xlabel('dBZ');ylabel('z (m)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  dBZ vs Altitude',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');
figure;plot(-mean(Vdopx1(:,ii)')+mean(Dopx(ii))+mean(Dopy(ii)),zrx','-o', mean(swrx1(:,ii)'),zrx','-x');xlabel('Fall Velocity (m/s)');ylabel('z (m)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  dBZ vs Altitude',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');
NRCS=dnzx-180+137.9+4.2;%+atten;

%
%dnrxmm=[];dopxm=[];dopym=[];Vdopmm=[];swrxm=[];npm=[];ntm=[];snrxmm=[];swrxmm=[];refzxmm=[];nim=[];nmm=[];umm=[];npm=[];nrcsm=[];Us3m=[];roll3m=[];attenm=[];reftzm=[];dnzxm=[];hhm=[];pitch3m=[];
dnrxmm=[dnrxmm' dnrx1m']';snrxmm=[snrxmm' snrx1m']';dopxm=[dopxm' mean(Dopx(ii))']';dopym=[dopym' mean(Dopy(ii))']';Vdopmm=[Vdopmm' mean(Vdopx1(:,ii)')']';swrxmm=[swrxmm' mean(swrx1(:,ii)')']';refzxmm=[refzxmm' mean(refzx1(:,ii)')']';%Us3m=[Us3m' Us3']';
npm=[npm' np']';ntm=[ntm' nt']';umm=[umm' Um']';nrcsm=[nrcsm NRCS];Us3m=[Us3m Us3];roll3m=[roll3m roll3];pitch3m=[pitch3m pitch3];attenm=[attenm atten];reftzm=[reftzm reftz];dnzxm=[dnzxm dnzx];hhm=[hhm hh];


figure;
%subplot(2,1,1);plot(mean(dnrxmm),zrx','-o',dnrxmm(1,:),zrx','-',dnrxmm(2,:),zrx','--',dnrxmm(3,:),zrx',':');xlabel('dBZ');ylabel('z (m)');
%subplot(2,1,1);plot(10*log10(nanmean(refzxmm)),zrx','-o',10*log10(refzxmm(1,:)),zrx','-',10*log10(refzxmm(2,:)),zrx','--',10*log10(refzxmm(3,:)),zrx',':',nanmean(snrxmm),zrx','-x');xlabel('dBZ');ylabel('z (m)');
subplot(2,1,1);plot(10*log10(nanmean(refzxmm)),zrx','-o',10*log10(refzxmm(1,:)),zrx','-',10*log10(refzxmm(3,:)),zrx',':',nanmean(snrxmm),zrx','-x');xlabel('dBZ');ylabel('z (m)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  dBZ and Vg vs Altitude',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,'7-10'),'FontWeight','Bold','Interpreter','none');
legend('Mean dBZ','hr 07','hr 09','hr 10','Mean SNR');
legend('Mean dBZ','hr 07','hr 10','Mean SNR');

subplot(2,1,2);plot(-nanmean(Vdopmm)+nanmean(dopxm)+nanmean(dopym),zrx','-o',-Vdopmm(1,:)+dopxm(1)+dopym(1),zrx','-',-Vdopmm(2,:)+dopxm(2)+dopym(2),zrx','--',-Vdopmm(3,:)+dopxm(3)+dopym(3),zrx',':',mean(swrxmm),zrx','-x');xlabel('Fall Velocity (m/s)');ylabel('z (m)');xlim([0 2])
subplot(2,1,2);plot(-nanmean(Vdopmm)+nanmean(dopxm)+nanmean(dopym),zrx','-o',-Vdopmm(1,:)+dopxm(1)+dopym(1),zrx','-',-Vdopmm(3,:)+dopxm(3)+dopym(3),zrx',':',nanmean(swrxmm),zrx','-x');xlabel('Fall Velocity (m/s)');ylabel('z (m)');xlim([0 2])
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  Fall Velocity vs Altitude',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,'7-10'),'FontWeight','Bold','Interpreter','none');
legend('Mean Vg','hr 07','hr 09','hr 10', 'Mean Width');
legend('Mean Vg','hr 07','hr 10', 'Mean Width');

vb=[];for i=1:3;vb(i,:)=-Vdopmm(i,:)+dopxm(i)+dopym(i);end;
figure;
subplot(2,1,1);plot(nanmean(dnrxmm),zrx','-o',nanstd(dnrxmm),zrx','-x',nanmean(snrxmm),zrx','-d');xlabel('dBZ');ylabel('z (m)');legend('dBZ mean','\sigma_{dBZ}','Mean SNR');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  dBZ and Vg vs Altitude',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,'7-10'),'FontWeight','Bold','Interpreter','none');
subplot(2,1,2);plot(-nanmean(Vdopmm)+nanmean(dopxm)+nanmean(dopym),zrx','-o',nanstd(vb),zrx,'-x');xlabel('Velocity (m/s)');ylabel('z (m)');legend('Vg mean','\sigma_{Vg}')

%%%%%%%%%%%%%%  frisch moment method
r=2:2:1000;%particle radius, micron
D=2*r;%diameter
for i=1:length(r)
   vf(i)=wf2(r(i)/1e6);%droplet fall velocity at radius r, Prupaccher and Klett
   vfp(i)=wf2(r(i)*2/1e6);
end;

a=1.2e-4;b=1e-5;
vgm=-nanmean(Vdopmm)+nanmean(dopxm)+nanmean(dopym);zzm=nanmean(dnrxmm);%mean vert doppler
svm=nanmean(swrxmm);%mean width
dbzm=nanmean(dnrxmm);
zzm=nanmean(refzxmm);
sgx=sqrt(log(1+svm.^2./(vgm+b/a).^2));%gaussian droplet width
sigx=max(.4,sgx);
r0=(a*vgm+b)./exp(13*sgx(length(sgx)).^2/2);%mode radius
n0=(zzm*1e-18)./((2.*r0).^6.*exp(18*sigx.^2));%number drops/m^3
dfrac=mean(npm)/mean(ntm);%fraction of acceptable conditions with valid returns (nsr>-12)
ql=zzm./(a*vgm).^3.*10.^((-132+26*sgx.^2)/10)*dfrac;
ij=find(zrx<300);PP=polyfit(zrx(ij)',log10(ql(ij)),1);
zrz=[zrx' 50 25 0]';qlf=polyval(PP,zrz);qlf=10.^(qlf);
%ql2=zzm./(a*vgm).^3.*10.^((-132+26*sgx.^2)/10)*dfrac;
mq=length(ql);
figure;semilogx(ql,zrx,'-o',qlf,zrz,'--',[.45e-3 .48e-3 .54e-3 .62e-3 .73e-3 .97e-3],[300 200 100 50 25 10],':', 1.2*[1e-3 3e-3 5e-3 .02 .06 .08]*(16/20)^3,[250 200 150 100 50 25] ,'-x');xlabel('Liquid Water (g/m^3)');ylabel('z (m)');xlim([1e-4 1.e-1]);
hold on;semilogx(sqrt(.4*4)*[1e-3 1e-3],[70 200],'-k',[.4*1e-3 4*1e-3],[150 150],'-k',sqrt(.4*4)*1e-3,150,'dk')
legend('Frisch et al. 1995 Retrieval','Least Square','Jones & Andreas 2013','Shpund et al. 2011','CBLAST CIP 2003');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  Droplet Liquid Water vs Altitude',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,'7-10'),'FontWeight','Bold','Interpreter','none');


figure;plot(roll3,dnzx,'.');grid; xlabel('Roll Angle (Deg)');ylabel('Peak Reflectivity (dBZ)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  dBZ vs Roll',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');

figure;plot(roll3,NRCS,'.');grid; xlabel('Roll Angle (Deg)');ylabel('NRCS (dB)');
%%%%%%%%%%%%%   137.9 is 10*log10(dR*pi^5*K^2/lambda^4),  1.6 is loss of
%%%%%%%%%%%%%   window, Attenuation is est as 4.2 dB
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  NRCS vs Roll',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');

figure;plot(roll3,dopx(2:4,:),'.',[-20 20],[0 0 ],'--');ylim([-4 4]);grid; xlabel('Roll Angle (Deg)');ylabel('Doppler (m/s)');
figure;plot(roll3,zzx*1e3-alt3,'.',roll3,+25+alt3.*(1./cos(roll3*pi/180)-1));grid; xlabel('Roll Angle (Deg)');ylabel('Range-Altitude (m)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  Alt-Range vs Roll',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');
legend('ReturnRange-Altitude','25m+Altitude*(1/cos(roll)-1)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  Doppler vs Roll',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');
legend('-1 Gate','Max Return Gate','+1 Gate');

dhh=2;hedx=10;ii=find((roll3)>hedx-dhh & alt3<2.5e3);figure;plot(hdir3(ii),(dopx(2,ii)-nanmean(dopx(2,ii)))./sin(roll3(ii)*pi/180),'.',hdir3(ii),(dopx(3,ii)-nanmean(dopx(3,ii)))./sin(roll3(ii)*pi/180),'.',[0 360],[0 0],'--');ylim([-7.5 7.5]);
grid; xlabel('Heading Angle (Deg)');ylabel('Doppler/sin(roll) (m/s)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  Doppler vs Heading',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');
legend('Max Return Gate','+1 Gate');
         
dhh=3;hedx=12.5;ii=find((roll3)>hedx-dhh & (roll3)<hedx+dhh & alt3<2.5e3);figure;plot(mod(hdir3(ii)-90,360),NRCS(ii),'b.');
hold on;hedx=-12.5;ii=find((roll3)>hedx-dhh & (roll3)<hedx+dhh & alt3<2.5e3);plot(mod(90+hdir3(ii),360),NRCS(ii),'.r');
grid; xlabel('Direction (Deg)');ylabel('NRCS (dB)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  NRCS vs Radial Direction ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');
legend('Roll Positive','Roll Negative');

dhh=2;hedx=10;ii=find((roll3)>hedx-dhh  & alt3<2.5e3);figure;plot(mod(hdir3(ii)-90,360),(dopx(2,ii)-nanmean(dopx(2,ii)))./sin(roll3(ii)*pi/180),'b.');ylim([-7.5 7.5]);
hold on;hedx=-10.;ii=find((roll3)<hedx+dhh & alt3<2.5e3);plot(mod(90+hdir3(ii),360),(dopx(2,ii)-nanmean(dopx(2,ii)))./sin(roll3(ii)*pi/180),'.r');ylim([-7.5 7.5]);
grid; xlabel('Direction (Deg)');ylabel('Doppler/sin(roll) (m/s)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s). Doppler vs Radial Direction ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');
legend('Roll Positive','Roll Negative');


dhh=3;hedx=12.5;ii=find((roll3)>hedx-dhh & (roll3)<hedx+dhh & alt3<2.5e3);figure;plot(mod(hdir3(ii)-90,360),10.^((NRCS(ii))/10),'b.');
hold on;hedx=-12.5;ii=find((roll3)>hedx-dhh & (roll3)<hedx+dhh & alt3<2.5e3);plot(mod(90+hdir3(ii),360),10.^((NRCS(ii))/10),'.r');
grid; xlabel('Direction (Deg)');ylabel('NRCS');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  NRCS vs Radial Direction ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');
legend('Roll Positive','Roll Negative');

ddh=[0 2 4 6 8 11 14 18  26 30];
mx=length(ddh);dnh=[];rolh=[];dnhn=[];nnd=[];
for i=1:mx-1
       jj=find(roll3>ddh(i) & roll3<ddh(i+1) & alt3<2.5e3);
       dnh(i)=nanmedian((NRCS(jj)));
       rolh(i)=nanmedian(roll3(jj));
       dnhn(i)=nanmean(10.^((NRCS(jj))/10));
       nnd(i)=length(jj);
end;    
kk=find(nnd>25);figure;plot(rolh(kk),dnh(kk),'-ob',rolh(kk),10*log10(dnhn(kk)),'-xr');grid; xlabel('Roll Angle (Deg)');ylabel('NRCS (dB)');
ij=find(alt3<2.5e3);figure;plot(abs(roll3(ij)),NRCS(ij),'.',rolh,dnh,'-ob',rolh,10*log10(dnhn),'-xr');grid; xlabel('Roll Angle (Deg)');ylabel('NRCS (dB)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  NRCS vs Roll Angle ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');
legend('Median','Mean Linear');

kkf=find(nnd>25 & (tan(rolh*pi/180)).^2<.1);[P,S] = polyfit((tan(rolh(kkf)*pi/180)).^2,10*log10(dnhn(kkf)),1);
kk=find(nnd>25);dnhfit=polyval(P,(tan(rolh(kk)*pi/180)).^2);figure;plot(tan(rolh(kk)*pi/180).^2,dnh(kk),'-ob',tan(rolh(kk)*pi/180).^2,10*log10(dnhn(kk)),'-xr',tan(rolh(kk)*pi/180).^2,dnhfit,'--');grid; xlabel('(Tan(Roll))^2');ylabel('NRCS (dB)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  NRCS vs [tan(roll)]^2 ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');
legend('Median','Mean Linear');

dirdopp=NaN*ones(size(roll3));
jj=find(roll3>2);dirdopp(jj)=mod(hdir3(jj)-90,360);
jj=find(roll3<-2);dirdopp(jj)=mod(hdir3(jj)+90,360);

ddh1=0:10:360;
mx=length(ddh1);
dnh1=[];dirh1=[];dnhn1=[];
for i=1:mx-1
       jj=find(dirdopp>ddh1(i) & dirdopp<ddh1(i+1) & alt3<2.5e3 & abs(roll3)>8 & abs(roll3)<12);
       dnh1(i)=nanmedian((NRCS(jj)));
       dirh1(i)=nanmedian(dirdopp(jj));
       dnhn1(i)=nanmean(10.^((NRCS(jj))/10));
end;
dnh2=[];dirh2=[];dnhn2=[];
for i=1:mx-1
       jj=find(dirdopp>ddh1(i) & dirdopp<ddh1(i+1) & alt3<2.5e3 & abs(roll3)>12 & abs(roll3)<18);
       dnh2(i)=nanmedian((NRCS(jj)));
       dirh2(i)=nanmedian(dirdopp(jj));
       dnhn2(i)=nanmean(10.^((NRCS(jj))/10));
end;

figure;plot(dirh1,10*log10(dnhn1),'-ob',dirh2,10*log10(dnhn2),'-or');grid; xlabel('Doppler Direction (Deg)');ylabel('NRCS (dB)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i hr-%s).  NRCS vs Direction, 10/15 deg roll ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none');
legend('roll=10 deg','roll=15 deg');


figure;
subplot(2,1,1);plot(hh-10,10.^(dnrx1(1,:)/10),'o',hh-10,10.^(dnrx1(6,:)/10),'d',hh-10,10.^(dnrx1(10,:)/10),'x');xlim([.310 0.3125]);
subplot(2,1,2);plot(hh-10,dnrx3);xlim([.3110 0.3125]);
            Ref2=[];htx=[];DV2=[];SW2=[];
            for i=1:150
             htx(i)=-(ht(150-i+1)-mean(zzx));
             Ref2(i,:)=Ref(150-i+1,:);
             DV2(i,:)=DV(150-i+1,:);
             SW2(i,:)=SW(150-i+1,:);
             end;
            figure;
           imagesc((hh-hh(1))*3600,htx(79:Nrange), Ref2(79:Nrange,:),[-40, 20]);
            axis xy; %axis([jam jam+1 0 ymaxx]);
            ylabel('Ht (km)');xlabel('Time (s)');
            %         tit={[cruise ' W-Band'' YD ' ddds ' Hr ' hr]};
            %         title(tit,'Interpreter','none')
            title(sprintf('%s (%04i-%02i-%02i, DOY%03i, Hr-%s).  W-Band (motionread=%i, Kongsberg=%i)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr,motionread,Konmotion),'FontWeight','Bold','Interpreter','none')
            h=colorbar;
            ylabel(h,'Ref dBZ') 

        ii=find (abs(roll3m)<2 & attenm<2 );     
        nrcf=9.94-.2*Us3m-.004*Us3m.^2;
%         sgnrc(1)=0.75;%signrc=1+.04*Us3m;
%         sgnrc(2)=0.85;
%         sgnrc(3)=1.1;%
%         sgnrc(4)=1.1;
%         sgnrc(5)=1.1;
        signrc=.8+.03*Us3m;
        
        dk=2;nrcsn=[];frc=[];frc2=[];uk=[];
 for k=1:dk:20
     jj=find (Us3m(ii)>k & Us3m(ii)<=k+dk);
     if ~isempty(ii)
          nrcsn(k)=nanmedian(nrcsm(ii(jj)));
          frc(k)=length(jj);
          kk=find((nrcf(ii(jj))-nrcsm(ii(jj)))>2.5*signrc(ii(jj)));
          frc2(k)=length(kk);
          uk(k)=nanmedian(Us3m(ii(jj)));
     else
         nrcsn(k)=NaN;
         frc(k)=NaN;
         frc2(k)=NaN;
         uk(k)=NaN;
     end;
 end;
 ij=find(frc>300);figure;plot(uk(ij),frc2(ij)./frc(ij),'o');
 ij=find(frc>100);figure;plot(uk(ij),nrcsn(ij),'o',Us3m,nrcf,'.',[5.5 8 11.5],[11 9 8],'-x');
 ii=find(10*log10(reftzm)<-30);reftzm(ii)=10.^((-28)/10);ii=find(isnan(dnzxm));dnzxm(ii)=45;
 ii=find(abs(roll3m)<3 );P = polyfit(10*log10(reftzm(ii))-min(10*log10(reftzm)),max(dnzxm)-dnzxm(ii)-4,1);attenx=polyval(P,10*log10(reftzm)-min(10*log10(reftzm)))-1.5;
 attenz=1.25*2*2.3/50*reftzm;Rainz=(reftzm/50).^1.25;

 xk=-5:.25:10;
 k=0;jj=find (Us3m(ii)>k & Us3m(ii)<=k+4 & attenx(ii)<2 );figure;subplot(5,1,1);hist((nrcf(ii(jj))-nrcsm(ii(jj))),xk);xlim([-5 10]); xlabel('<NRCS>-NRCS');ylabel('Number');legend('D-NRCS 0<U_{10}<4 ms^{-1}');
 k=4;jj=find (Us3m(ii)>k & Us3m(ii)<=k+8 & attenx(ii)<2 );subplot(5,1,2);hist((nrcf(ii(jj))-nrcsm(ii(jj))),xk);xlim([-5 10]); xlabel('<NRCS>-NRCS');ylabel('Number');legend('D-NRCS 4<U_{10}<8 ms^{-1}');
 k=10;jj=find (Us3m(ii)>k & Us3m(ii)<=k+3 & attenx(ii)<2 );subplot(5,1,3);hist((nrcf(ii(jj))-nrcsm(ii(jj))),xk);xlim([-5 10]); xlabel('<NRCS>-NRCS');ylabel('Number');legend('D-NRCS 10<U_{10}<13 ms^{-1}');
 k=13;jj=find (Us3m(ii)>k & Us3m(ii)<=k+3 & attenx(ii)<2 );subplot(5,1,4);hist((nrcf(ii(jj))-nrcsm(ii(jj))),xk);xlim([-5 10]); xlabel('<NRCS>-NRCS');ylabel('Number');legend('D-NRCS 13<U_{10}<16 ms^{-1}');  
 k=16;jj=find (Us3m(ii)>k & Us3m(ii)<=k+5 & attenx(ii)<2 );subplot(5,1,5);hist((nrcf(ii(jj))-nrcsm(ii(jj))),xk);xlim([-5 10]); xlabel('<NRCS>-NRCS');ylabel('Number');legend('D-NRCS 16<U_{10}<20 ms^{-1}');  

 xk=-5:.25:15;
 k=0;jj=find (Us3m(ii)>k & Us3m(ii)<=k+4 & attenx(ii)<2 );figure;subplot(5,1,1);hist((nrcsm(ii(jj))),xk);xlim([-5 15]); xlabel('NRCS');ylabel('Number');legend('NRCS 0<U_{10}<4 ms^{-1}');
 k=4;jj=find (Us3m(ii)>k & Us3m(ii)<=k+4 & attenx(ii)<2 );subplot(5,1,2);hist((nrcsm(ii(jj))),xk);xlim([-5 15]); xlabel('NRCS');ylabel('Number');legend('NRCS 4<U_{10}<8 ms^{-1}');
 k=10;jj=find (Us3m(ii)>k & Us3m(ii)<=k+3 & attenx(ii)<2 );subplot(5,1,3);hist((nrcsm(ii(jj))),xk);xlim([-5 15]); xlabel('NRCS');ylabel('Number');legend('NRCS 10<U_{10}<13 ms^{-1}');
 k=13;jj=find (Us3m(ii)>k & Us3m(ii)<=k+3 & attenx(ii)<2 );subplot(5,1,4);hist((nrcsm(ii(jj))),xk);xlim([-5 15]); xlabel('NRCS');ylabel('Number');legend('NRCS 13<U_{10}<16 ms^{-1}');  
 k=16;jj=find (Us3m(ii)>k & Us3m(ii)<=k+5 & attenx(ii)<2 );subplot(5,1,5);hist((nrcsm(ii(jj))),xk);xlim([-5 15]); xlabel('NRCS');ylabel('Number');legend('NRCS 16<U_{10}<20 ms^{-1}');  

ii=find(attenx<2 & abs(roll3m)<2);figure;plot(Us3m(ii),nrcsm(ii),'.',Us3m(ii),nrcf(ii),'.');ylim([-5 15]);xlabel('U_{10} (ms^{-1})');ylabel('NRCS (dB)');



            ymaxx=ht(Nrange);
            zmxx=ymaxx;
             Ref2=[];htx=[];DV2=[];SW2=[];
            for i=1:150
             htx(i)=-(ht(150-i+1)-mean(zzx));
             Ref2(i,:)=Ref(150-i+1,:);
             DV2(i,:)=DV(150-i+1,:);
             SW2(i,:)=SW(150-i+1,:);
             end;
            figure;
            figure; %clf
            subplot(3,1,1),
           %imagesc(hh,ht, Ref(1:Nrange,:),[-60, 20]);
            imagesc((hh),htx(1:Nrange), Ref2(1:Nrange,:),[-30, 20]);
            axis xy; axis([jam jam+1 0 2.3]);
            ylabel('Ht (km)')
            %         tit={[cruise ' W-Band'' YD ' ddds ' Hr ' hr]};
            %         title(tit,'Interpreter','none')
            title(sprintf('%s (%04i-%02i-%02i, DOY%03i, Hr-%s).  W-Band (motionread=%i, Kongsberg=%i)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr,motionread,Konmotion),'FontWeight','Bold','Interpreter','none')
            h=colorbar;
            ylabel(h,'Ref dBZ') 
            hold on
            
            subplot(3,1,2)
            imagesc(hh,htx, DV2(1:Nrange,:),[-2, 7]);  %
            axis xy;axis([jam jam+1 0 2.3]);
            ylabel('Ht (km)')
            h=colorbar;
            ylabel(h,'-DopVel (m/s)')
            hold on
            
            subplot(3,1,3),imagesc(hh,htx, SW2(1:Nrange,:),[0, 2]);  %
            axis xy;axis([jam jam+1 0 2.3]);
            ylabel('Ht (km)')
            xlabel('Hrs (UTC)')
            h=colorbar;
            axis xy;axis([jam jam+1 0 2.3]);
            ylabel(h,'SpecWid (m/s)')
            % print('-djpeg90 ',[way_raw_images_wband '\Wband_ref_' num2str(ddd), '_' hr '.jpg']);
%            annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
            