% 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 all
%################## 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
            z=tt/(60*60*24)+datenum([1970,1,1,0,0,0]);
            jd=z-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_lb;
    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);
[nx,mx]=size(Ref);
dnzx=[];zzx=[];snrx=[];dopx=NaN*ones(4,mx);widx=NaN*ones(4,mx);
for i=1:mx
    if Nrangeh(i)==150
        rmx=max(Ref(:,i));
        if rmx<100 & rmx>0;
            ij=find(Ref(:,i)==rmx);;dnzx(i)=rmx;zzx(i)=ht(ij);snrx(i)=SNR(ij,i);dopx(:,i)=DV(ij-2:ij+1,i);;widx(:,i)=SW(ij-2:ij+1,i);
        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;

roll3 = interp1(P3_hh,P3_rollref,hh);alt3=interp1(P3_hh,P3_altref,hh);head3=interp1(P3_hh,P3_head,hh);
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,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)');

figure;plot(roll3,dopx(2:4,:),'.',[-20 20],[0 0 ],'--');ylim([-4 4]);grid; xlabel('Roll Angle (Deg)');ylabel('Doppler (m/s)');
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(head3(ii),(dopx(2,ii)-nanmean(dopx(2,ii)))./sin(roll3(ii)*pi/180),'.',head3(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');


