% Reads in W-Band nc moment data and plot
% range gates are variable
% Use jamn for multiple files in 1 hour.
% Preallocate arrays to combine full-day time series.
% Reads P3 flight level time series as P3nav 

close all
clear ddd ehr shr hh
%################## Cruise informations  ##################################################################
cruise='ATOMIC';             % Cruise name  
year='2020';               % year (deprecated--make year into a double) 
yyyy='2020';               % year (string)
ship='P3';               % Research vessel

%%%%%%%%%%%%%%%%%%%%%%%%% 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

big_dir='E:\DATA\ATOMIC\P3\wband\';%select address of main data files

disp('fockewe2')
% W-band radar files
way_raw_data_wband=[big_dir ddds '\mom\'];                  %location moment data
way_raw_images_wband=[big_dir ddds '\Raw_images\'];         %where put processed images
way_raw_motion_wband=[big_dir ddds '\motion\'];              %location P3 nav data


%##################%################# set control functions #%##################%##################%##################%##################

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;%use matlab net cdf software
srt=0;
%%%%%%%%%%%%%%%  Read aircraft flight level data
bigread=1;   
if bigread;
    %read_P3navdata_FLTLVL_atomic_nc;%nc file
    %read_P3navdata_FLTLVL_atomic_BB;%IWG dat file
    %%%%%%%%%%  use preprocessed mat file
    dr = dir([way_raw_motion_wband '*.mat']);
    filnx = dr.name;
    fname = [way_raw_motion_wband filnx];
    load(fname);
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    ii=find(P3_lat>10 & P3_lon>-60);
    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;
    %%%%%  compute mean profiles T,q,P, Un, Ue from P3 flight level data
    %%%%%  (could use dropsondes for this
    zmin=0;zmax=8000;dzp=100;nx=(zmax-zmin)/dzp+1;
   [Xp, Yp, M, SD, SE, SDM, SEM, N] = binave_4(P3_altref,P3_press,dzp,zmin,zmax);
   [Xt, Yt, M, SD, SE, SDM, SEM, N] = binave_4(P3_altref,P3_Ta,dzp,zmin,zmax);
   [Xq, Yq, M, SD, SE, SDM, SEM, N] = binave_4(P3_altref,P3_Q,dzp,zmin,zmax);
   [Xun, Yun, M, SD, SE, SDM, SEM, N] = binave_4(P3_altref,-P3_Windspeed.*cos(P3_WindDir*pi/180),dzp,zmin,zmax);%positive to n
   [Xue, Yue, M, SD, SE, SDM, SEM, N] = binave_4(P3_altref,-P3_Windspeed.*sin(P3_WindDir*pi/180),dzp,zmin,zmax);%postive to e
end;
cldtop=[];
% preallocate x-extended variables
nxmax=9*7200;%assumes no more than 9 hrs, 220 levels
Nrangep=238; % number ranges for re-mapped data top at 7.5 km
Nrange=220;  %number ranges of Wband file
Refxp=deal(single(zeros(Nrangep,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')
%define variables to append sucessive hours;
dnrxmm=[];dopxm=[];dopym=[];Vdopmm=[];swrxm=[];npm=[];ntm=[];snrxmm=[];swrxmm=[];refzxmm=[];nim=[];nmm=[];umm=[];npm=[];nrcsm=[];
Us3m=[];roll3m=[];attenm=[];reftzm=[];dnzxm=[];hhm=[];pitch3m=[];P3_latm=[];P3_lonm=[];P3_altrefm=[];P3_sstm=[];P3_hhm=[];dop0m=[];snr0m=[];
httm=[];reftm=[];Rs3m=[];Refpm=[];rainm=[];DVpmn=[];SWpm=[];SNRpm=[];DVpm=[];nrcfm=[];P3_pressm=[];P3_tasm=[];
Refm=[];DVm=[];SWm=[];SNRm=[];httbm=[];SSTm=[];snrxm=[];widxm=[];wvelm=[];sog3m=[];cog3m=[];
nx=0; % count position in extended variables

%%%%%%%%%%%%%%%%%%%%%%%%%%   Main loop reading hourly radar files, time hit
%%%%%%%%%%%%%%%%%%%%%%%%%%   cross section plots
for jam=shr:ehr
    isfile=false;
%    clf(handl(:));
    nh=0;
    clear y z dt jd
    hh=[]; Ref=[]; DVa=[]; SW=[]; SNR=[]; nrg=[];%zero 
    hr=sprintf('%02i',jam);
    dr=dir([way_raw_data_wband yyyy ddds hr '*MMCRMom.nc']);%find the hourly files jam
    for jamn=1:length(dr) % multiple files in one hour
        fn=[way_raw_data_wband dr(jamn).name];%file 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                    %Read each variable out of nc file                                  
                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
                NoiseLevel=netcdf.getVar(FID, 39)';
                %%%%%%%%%%   convert object to number (don't ask me why)
                Ref=single(Reflectivity)';
                SW=single(SpectralWidth)';
                DVa=single(MeanDopplerVelocity)';
                SNR=single(SignalToNoiseRatio)';
                Noise=single(NoiseLevel)';
                nrg=netcdf.getVar(FID,18);
                Nrange=nrg(1)
                netcdf.close(FID);  %close the file            
            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=-35;
    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;

     Konmotion=0;
    if motionread %correct doppler for ship motion
    %                 if exist([way_raw_motion_wband yyyy ddds hr '00Kongsberg.txt'],'file');
    
        if ~bigread  % reads hourly aircraft data files saved by Wband radar
            read_P3navdata_lb2_atomic_BB;
            figure;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');grid
            xlabel('Hour (UTC)');ylabel('Pitch');
            subplot(1,2,2);plot(P3_hh,P3_rollref,'r');grid
            xlabel('Hour (UTC)');ylabel('Roll');
            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
        end;
        %%%%%%%%   interpolated motion data onto radar time series, hh
        %%%%%%%%   radar hours
        ii=find(diff(P3_hh)<1e-4);P3_hh(ii)=P3_hh(ii)+1e-5;
        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);
        sst3=interp1(P3_hh,P3_TradD,hh);   press3=interp1(P3_hh,P3_press,hh);
        ii=find(abs(Us3)>30);Us3(ii)=NaN;

       if ~isempty(dir([way_raw_motion_wband '020\p3' ddds hr '_konesberg.txt'])); %no files at all
            %read_vertical velocity from PSD Konsberg sensor - only used
            %for ship system
            
            read_P3navdata_lb2_atomic_BB
            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 Reflectivity, Doppler, and Doppler width
    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)')
    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;
    
    
    if motionread
        zzx=(alt3/1e3);%altitude of the aircraft, converted to km
    else
        zzx=3.5;
    end;
    nnk=length(ht);dz=ht(2)-ht(1);htextr=ht(220):dz:7.6;htflp=[ht' htextr]'; %htflp=sortrows(htflp,'descend');
    figure(2); %clf
    ymaxx=ht(Nrange);
    zmxx=ymaxx;zzxx=zzx;
    %%%%%%%%%%%%%%%%%%%%%%%%  This code remaps range data onto a fixed altitude grid  
    Ref2=[];htx=[];DV2=[];SW2=[];SNR2=[];Refp=NaN*ones(length(htflp),length(Ref));htxx=[];DVp=Refp;SNRp=Refp;
    DVp=Refp;SWp=Refp;
    %%%%%%%%%%%%  restort data 
    for i=1:nnk
        %htx=(ht(nnk-i+1));
        Ref2(i,:)=Ref(nnk-i+1,:);
        DV2(i,:)=DV(nnk-i+1,:);
        SW2(i,:)=SW(nnk-i+1,:);
       SNR2(i,:)=SNR(nnk-i+1,:);
    end;
    %%%%%%%%%%%  Find the altitude of the closest range gate to the
    %%%%%%%%%%%  aircraft
    kk=floor((zzx-.15)/dz);
    ii=find(isnan(zzx));  if(~isempty(ii));kk(ii)=nanmedianx(kk);end;
    kk=max(kk,1);kk=min(kk,220);
    %%%%%%%%%%%  store the flipped data in the correct altitude
    for i=1:length(Ref2)
        Refp(1:kk(i),i)=Ref2(220-kk(i)+1:220,i);
        DVp(1:kk(i),i)=DV2(220-kk(i)+1:220,i);
        SWp(1:kk(i),i)=SW2(220-kk(i)+1:220,i);
        SNRp(1:kk(i),i)=SNR2(220-kk(i)+1:220,i);  
    end;
    
    
    gftp=7.5;%  Max altitude of graph
    figure; %clf
    subplot(3,1,1),
    imagesc((hh),htflp(1:Nrange,:), Refp(1:Nrange,:),[-50, 20]);
    axis xy; axis([jam jam+1 0 gftp]);
    ylabel('Ht (km)')
    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]);  %
    imagesc(hh,htflp, DVp(1:Nrange,:),[-7, 7]);  %
    axis xy;axis([jam jam+1 0 gftp]);
    ylabel('Ht (km)')
    h=colorbar;
    ylabel(h,'DopVel (m/s)')
    hold on
            
    subplot(3,1,3),imagesc(hh,htflp, SWp(1:Nrange,:),[0, 2]);  %
    axis xy;axis([jam jam+1 0 gftp]);
    ylabel('Ht (km)')
    xlabel('Hrs (UTC)')
    h=colorbar;
    axis xy;axis([jam jam+1 0 gftp]);
    ylabel(h,'SpecWid (m/s)')
 
    if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_Wband_Ref2_' 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(Refp);
     if m1>0 && n1>0
         hhx(:,nx+(1:n1))=single(hh);
         Refxp(1:m1,nx+(1:n1))=single(Refp);
         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 using altitude mapped data
        [htt reft reft2 wvel nnt nnk]=deal(single(zeros(1,size(Ref,2))));%put zeroes in the variables
        for k=1:size(Ref,2); % time dimension
           % logical for noisfor SNR exceeding a threshold and height in bounds
            kk=SNRp(:,k)>-10 & htflp>.5 & htflp<zmxx & Refp(:,k)>-30;
            if sum(kk)
                htt(k)=single(ht(find(kk,1,'last')));                               %Find cloud top height as the highest valid return
                reft(k)=single(nanmeanx(10.^(Refp(kk,k)/10)));                      %Z value at cloudtop
                reft2(k)=single(nanmeanx(10.^((Refp(kk,k)-dbz0-20*log10(ht(kk)))/10)));%Amount Z exceeds noise level
                wvel(k)=single(nanmeanx(DVp(kk,k)));                                 %Doppler
                nnt(k)=single(sum(kk));                                              %number of valid returns in the column
                nnk(k)=find(kk,1,'last');                                            %index of cloud top
            else
                htt(k)=NaN;
                reft(k)=NaN;
                reft2(k)=NaN;
                wvel(k)=NaN;
                nnt(k)=0;
                nnk(k)=238;
            end;
        end;
        
        y=prctilex(htt,[10 15 50 85 90]);                        %distribution of cloud top heights
        kk=htt>=y(1) & htt<=y(5);                                %times when cloud top in bounds
        httb=single(nanmedianx(htt(kk)));                        %mean cloud top height for hour
        reftb=single(10*log10(nanmedianx(reft(kk))));
        reft2b=single(10*log10(nanmedianx(reft2(kk))));
        wvelb=single(nanmeanx(wvel(kk)));                        %mean Doppler
        nntb(1,1:Nrange)=single(sum(nnt));                       %
        sigh=single(sqrt(nanmeanx(htt(kk).^2)-httb.^2));         %stnd dev cloudtop height
        disp(['Hour= ' num2str((hh(1))) '  Median cloud top ht= ' num2str(httb) '  std cloud top height=' num2str(sigh)])
        disp(['Mean dBZ= ' num2str(reftb) '   Median dBZ-noisedBZ= ' num2str(reft2b) '   Median 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
 

%Interpolated aircraft data again

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);
sst3=interp1(P3_hh,P3_TradD,hh);tas3=interp1(P3_hh,P3_TAS,hh);        ii=find(abs(Us3)>30);Us3(ii)=NaN;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Simple correction for pitch and aircraft speeed
ff=pi/180;
DVac=DVa-tas3.*sin(ff*pitch3);
DVpac=DVp-tas3.*sin(ff*pitch3);
Refv=Ref;
SNRv=SNR;
SWv=SW;
for i=1:220;
    ii=find(SNR(i,:)<-20);DVac(i,ii)=NaN;Ref(i,ii)=NaN;
    ii=find(SNRp(i,:)<-20);DVpac(i,ii)=NaN;Refp(i,ii)=NaN;
end;

%%%%%%%%%%%%%%%%%%%   Doppler veloc correctionss for aircraft pitch/roll
%%%%%%%%%%  not actually used here
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=NaN*ones(1,mx);dop0=dnzx;snr0=dop0;rrefx=dnzx;;zzx=[];snrx=[];dopx=NaN*ones(4,mx);rrefx=dopx;snrx=dopx;widx=NaN*ones(4,mx);dnrx1=[];reftz=[];dnrx2=[];dnrx3=[];zrx=[];dnrx1m=[];snrz=[];usfm=[];latm=[];lonm=[];swrx1=[];refzx1=[];reftz=[];snrx1=[];Vdopx1=[];
dij=3;
for i=1:mx
    if (Nrangeh(i)==150 | Nrangeh(i)==220)
        rng=alt3(i)./cos(roll3(i)*pi/180);
        kk=find(ht*1e3>rng-150 & ht*1e3<rng+200);qkk=size(kk);
        rmx=max(Ref(kk,i));ij=find(Ref(kk,i)==rmx);
        if (rmx<100 & rmx>0) & kk(ij)<217 ;
            ij=find(Ref(kk,i)==rmx);wij=ij-2:ij+1;
            if ij==2;wij=ij-1:ij+2;end;
            if ij==1;wij=ij:ij+3;end;  
            if ij==11;wij=ij-3:ij;end;
            if ij==12;wij=ij-3:ij;end;
            jk=kk(ij);dnzx(i)=rmx;zzx(i)=ht(kk(ij));dop0(i)=DV(kk(ij),i);snr0(i)=SNR(kk(ij),i);
            dnrx2(i)=Ref(jk-10-dij,i)+refcorr;dnrx3(i)=Ref(jk-5,i)+refcorr;snrz(i)=SNR(jk-3,i);
            dopx(:,i)=DV(kk(wij),i); snrx(:,i)=SNR(kk(wij),i);widx(:,i)=SW(kk(wij),i);rrefx(:,i)=Ref(kk(wij),i);
            %usfm(i)=Us3(jk);latm(i)=lat3(jk);lonm(i)=lon3(jk);
        else
             dnzx(i)=NaN;zzx(i)=NaN;snrx(i)=NaN;dnrx3(i)=NaN;snrz(i)=NaN;reftz(i)=NaN;%dopx(:,i)=NaN;dnrx1(:,i)=NaN
        end;
    else
        dnzx(i)=NaN;zzx(i)=NaN;snrx(i)=NaN;%dopx(:,i)=NaN;
    end;
end;
%Dvac=DVa-tas3.*sin(ff*pitch3);
snr_th=-12;
NRCS=dnzx-180+138.7+4.2;%+atten;
nrcf=9.94-.2*Us3-.004*Us3.^2;
rain=0.5*(nrcf-NRCS)*1.3/sqrt(1.2)./alt3*1e3;ii=find(abs(roll3)>1.5 | abs(pitch3)>2);rain(ii)=NaN;

% 
dopxm=[dopxm dopx];snrxm=[snrxm snrx];widxm=[widxm widx];%dopym=[dopym' dopyy']';Vdopmm=[Vdopmm' Vdopxx']';swrxmm=[swrxmm' swrxx']';refzxmm=[refzxmm' refzxx']';%

%npm=[npm' np']';ntm=[ntm' nt']';umm=[umm' Um']';
nrcsm=[nrcsm NRCS];Us3m=[Us3m Us3];roll3m=[roll3m roll3];pitch3m=[pitch3m pitch3];dnzxm=[dnzxm dnzx];hhm=[hhm hh];
P3_latm=[P3_latm lat3];P3_lonm=[P3_lonm lon3];P3_altrefm=[P3_altrefm alt3];P3_sstm=[P3_sstm sst3];dop0m=[dop0m dop0];snr0m=[snr0m snr0];
rainm=[rainm rain];Rs3m=[Rs3m Rs3];httm=[httm htt];reftm=[reftm reft];wvelm=[wvelm wvel];httbm=[httbm httb];P3_tasm=[P3_tasm tas3];
Refm=[Refm Ref];DVm=[DVm DVac];SWm=[SWm SW];SNRm=[SNRm SNR];P3_pressm=[P3_pressm press3];
Refpm=[Refpm Refp];DVpm=[DVpm DVpac];SWpm=[SWpm SWp];SNRpm=[SNRpm SNRp];nrcfm=[nrcfm nrcf];
 sog3m=[sog3m sog3]; cog3m=[cog3m cog3];
end; % hr loop
NRCSm=dnzxm-180+138.7+4.2;%+atten;
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)-nanmeanx(dopx(2,ii)))./sin(roll3(ii)*pi/180),'.',hdir3(ii),(dopx(3,ii)-nanmeanx(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;if ~isempty(ii);plot(mod(hdir3(ii)-90,360),NRCS(ii),'b.');end;
% hold on;hedx=-12.5;ii=find((roll3)>hedx-dhh & (roll3)<hedx+dhh & alt3<2.5e3);if ~isempty(ii);plot(mod(90+hdir3(ii),360),NRCS(ii),'.r');end;
% 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)-nanmeanx(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)-nanmeanx(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<5e3);
       dnh(i)=nanmedian((NRCS(jj)));
       rolh(i)=nanmedian(roll3(jj));
       dnhn(i)=nanmeanx(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<3.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)=nanmeanx(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)=nanmeanx(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');
% 
Refmx=[];DVmx=[];SWmx=[];SNRmx=[];Refr=[];DVr=[];SWr=[];SNRr=[];
if ddd==31;ik=find(hhm>19.55 & hhm<19.61);SNRm(:,ik)=-30;end;
jj=1:30;kk=find(P3_altrefm<3200 & P3_altrefm>1000);Refmx=Refm(jj,kk);DVmx=DVm(jj,kk);SWmx=SWm(jj,kk);SNRmx=SNRm(jj,kk);
[nn mm]=size(Refmx);
Refr=reshape(Refmx,1,nn*mm);
DVr=reshape(DVmx,1,nn*mm);
SWr=reshape(SWmx,1,nn*mm);
SNRr=reshape(SNRmx,1,nn*mm);
ii=find(SNRr>-10);[X, Y, M, SD, SE, SDM, SEM, N] = binave_3(Refr(ii),DVr(ii),5,-30,30);
ii=find(SNRr>-10);[Xr, Yr, Mr, SD, SE, SDM, SEM, N] = binave_3(Refr(ii),SNRr(ii),5,-30,30);
ii=find(SNRr>-10);[Xw, Yw, Mw, SD, SE, SDM, SEM, N] = binave_3(Refr(ii),SWr(ii),5,-30,30);


figure;ii=find(N>10);
subplot(2,1,1);plot(X(ii),M(ii),'-o');ylabel('Median Doppler Velocity (m/s)');grid;xlim([-30 30]);ylim([-5 1]);
title(sprintf('%s (%04i-%02i-%02i, DOY%03i ).  W-Band Stats down',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
subplot(2,1,2);semilogy(X(ii),N(ii)/nn/mm,'-o');ylabel('Fraction');grid
xlabel('Mean Reflectivity (dBZ)');
if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_Wband_RefStatsdown_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]); end;




%%  Compute dBZ vs DV stats
Refpmx=[];DVpmx=[];SWpmx=[];SNRpmx=[];Refpr=[];DVpr=[];SWpr=[];SNRpr=[];
if ddd==31;ik=find(hhm>19.55 & hhm<19.61);SNRpm(:,ik)=-30;end;
jj=40:100;kk=find(P3_altrefm<3200 & P3_altrefm>1000 & abs(roll3m)<2 & abs(pitch3m)<2);Refpmx=Refpm(jj,kk);DVpmx=DVpm(jj,kk);SWpmx=SWpm(jj,kk);SNRpmx=SNRpm(jj,kk);
[nn mm]=size(Refpmx);
Refpr=reshape(Refpmx,1,nn*mm);
DVpr=reshape(DVpmx,1,nn*mm);
SWpr=reshape(SWpmx,1,nn*mm);
SNRpr=reshape(SNRpmx,1,nn*mm); 
ii=find(SNRpr>-10);[Xpa, Ypa, Mpa, SD, SE, SDM, SEM, Npa] = binave_3(Refpr(ii),DVpr(ii),5,-30,30);

figure;ii=find(Npa>10);figure;
subplot(2,1,1);plot(Xpa(ii),Mpa(ii),'-o');ylabel('Median Doppler Velocity (m/s)');grid
title(sprintf('%s (%04i-%02i-%02i, DOY%03i ).  W-Band Stats Up',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
xlim([-30 30]);ylim([-5 1]);
subplot(2,1,2);semilogy(Xpa(ii),Npa(ii)/nn/mm,'-o');ylabel('Fraction');grid;
xlabel('Mean Reflectivity (dBZ)');xlim([-30 30]);ylim([1e-5 .1]);
if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_Wband_RefStats_up_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]); end;

%%%%%%  detect radar indication of clouds
cloudx=NaN*ones(1,length(hhm));cloudrf=cloudx;clouddp=cloudrf;
for i=1:length(hhm)-1
    [mx,ndx]=max(SNRpm(10:100,i));%% max SNR in column
    cloudx(i)=mx;
    if mx>-12;
        cloudrf(i)=Refpm(ndx,i);
        clouddp(i)=DVpm(ndx,i);
    end;%reflectivity if SNR>-12
end;

%%%%%%%%   use IR for cloud detection
ii=find(P3_sstm>30);P3_sstm(ii)=NaN;
ii=find(cloudx<-16); %locate no cloud conditions
kk=find(~isnan(P3_pressm(ii)));%pressure at aircraft height
PT=polyfit(P3_pressm(ii(kk))-1020,P3_sstm(ii(kk)),2);tir=polyval(PT,P3_pressm-1020);%fit no cloud ir signal from SST to aircraft altitude
%%%%%%%%%%%  fit to median to remove outliers
[Xg, Yg, Mg, SDg, SEg, SDMg, SEMg, Ng] = binave_4(P3_pressm(ii(kk)),P3_sstm(ii(kk)),20,390,1020);
nv=find(~isnan(Xg));PT2=polyfit(Xg(nv)-1020,Mg(nv),2);tir2=polyval(PT2,P3_pressm-1020);%fit no cloud ir signal from SST to aircraft altitude


prestm=[];tatm=[];uatm=[];
for i=1:length(hhm);
    iik=find(Xp>httm(i)*1e3);
    if ~isempty(iik);
        prestm(i)=Yp(iik(1));%estimate pressure at cloud top height
        tatm(i)=Yt(iik(1));
        uatm(i)=sqrt(Yue(iik(1)).^2+Yun(iik(1)).^2);
    else
        prestm(i)=NaN;tatm(i)=NaN;uatm(i)=NaN;
    end;
end;
tirct=P3_sstm+polyval(PT2,prestm-1020)-tir2;%compute ir temperature at cloud top.
tirsst=P3_sstm+polyval(PT2,0)-tir2;%compute ir sst , valid when no clouds below aircraft

ii=find(P3_altrefm<4000);figure;plot(10*log10(reftm(ii)),tatm(ii)-tirct(ii),'.');xlabel('Cloudtop dBZ');ylabel('Cloudtop Ta-Tir (C)');grid;

%ii=find(P3_altrefm>3500);tir(ii)=NaN;
cldindx=P3_sstm-tir2;
%%%%%%%  IR and Wband cloud fraction
ii=find(P3_altrefm<3500 & P3_altrefm>700);jj=find(cldindx(ii)<-1);%clouds when ac below 3.5 km and above cloud base
ii2=find(P3_altrefm<3500 & P3_altrefm>1000);kk=find(cloudx(ii2)>-12);%cloud when ac below 3.5 km and radar is on
[size(jj)/size(ii) size(kk)/size(ii2)]
ok=-20:1:-2;ok2=[ok -1.5 -1 -.5 0 .5 1 2 5];
[nx,X]=hist(cldindx(ii),ok2);i=21;[X(i) sum(nx(1:i))/sum(nx)]
figure;plot(X(1:i),nx(1:i)/sum(nx(1:i)),'-o');grid;xlabel('Cloudtop Temp, realtive to SST (C)');ylabel('Fraction');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i ).  Cloudtop Temperature-SST',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_WbandTir_cloudTemp_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]); end;


mj=5;mt=100;nkk=[];nkkj=[];
for i=mj:mt
    j=i-mj+1;
    kk=find(P3_altrefm>ht(i)*1e3 & P3_altrefm<3500 & P3_altrefm>1000);%find cases will radar on and below 3.5 km
    if isempty(kk)
        nkk(j)=0;
        nkkj(j)=1;
    else
       nkk(j)= length(kk);
       jjk=find(SNRpm(i,kk)>-15);nkkj(j)=length(jjk);%find valid returns
    end;
end;

figure;
subplot(3,1,1);plot(nanmeanx(Refpm(mj:mt,:)'),ht(mj:mt),'o');xlabel('Mean dBZ');ylabel('Height (m)');grid;xlim([-45 0]);
title(sprintf('%s (%04i-%02i-%02i, DOY%03i ).  W-Band median stats',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
subplot(3,1,2);plot(nanmedianx(DVpm(mj:mt,:)'),ht(mj:mt),'o');xlabel('Mean Doppler (m/s)');ylabel('Height (m)');grid;xlim([-5 .5]);
subplot(3,1,3);plot(nkkj./nkk,ht(mj:mt),'o');xlabel('Return Fraction');ylabel('Height (m)');grid;xlim([0 .1]);
if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_Wband_Prof_stats_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]); end;

[Xl, Yl, Ml, SD, SE, SDM, SEM, Nl] = binave_3(P3_lonm,cldindx,1,-60, -50);
ii=find(Nl>10);figure;plot(Xl(ii),-Yl(ii)+1,'-o',Xl(ii),SD(ii),'-o');xlabel('Longitude (deg)');ylabel('Mean IR cloud index (C)');grid;

figure;
subplot(2,1,1);plot(hhm,cloudx+16,'.',hhm,cldindx,'.');grid;ylabel('Cloud Index');legend('Wband','Tir');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i ).  W-Band IR Cloud Indices',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
subplot(2,1,2);plot(hhm, P3_altrefm);xlabel('Time (hr)');ylabel('Altitude (m)');
if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_Wband_CloudIndices_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]); end;

figure;plot(hhm,httm,'o');xlabel('Time (hr)');ylabel('Cloudtop Height (m)');grid
title(sprintf('%s (%04i-%02i-%02i, DOY%03i ).  W-Band Cloudtop Height (km)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_Wband_cloudTop_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]); end;
wvelma=wvelm-sin(pitch3m*pi/180).*P3_tasm;
reftmz=10*log10(reftm);cloudx2=cloudx+14;
%%%%%%%%%%   correct sfmr 10-m .04wind speed
u10c=5.7+.57*Us3m;
%%%%%%%%%%%%  correct nrcs 
Puc=10*(1-exp((P3_pressm-925)/150));nrcsmc=nrcsm-Puc+9.5;
ii=find(cloudx2>0 | P3_altrefm>3500 );jj=find(cldindx(ii)<-1);
tirsst(ii(jj))=NaN;

%save(['E:\DATA\ATOMIC\P3\wband\' ddds '\ATOMIC_P3_wband_' ddds '_cloudsum.mat'], 'hhm','P3_latm','P3_lonm','P3_altrefm','P3_tasm','pitch3m','roll3m','P3_pressm', 'P3_sstm','nrcsm','nrcsmc','Us3m','u10c','cloudx2','cldindx','httm','reftmz','wvelma','tatm','tirct','uatm','tirsst');
%                                                                                  jd      lat       lon       alt       airspd     pitch     roll       P         raw Tir    nrcs   U10     radar     IR       topht    topZ     topW    topT   topIR   utop   
%ddds='017';load(['E:\DATA\ATOMIC\P3\wband\' ddds '\ATOMIC_P3_wband_' ddds '_cloudsum.mat'], 'hhm','latm','lonm','P3_altrefm','P3_tasm','pitch3m','roll3m','P3_pressm', 'P3_sstm','nrcsm','cloudx2','cldindx','httm','reftmz','wvelma','tatm','tirct','uatm');