% DEW KM read in W-Band nc moment data and plot
% range gates are 25 m
% timing is every 0.25 sec
%Version_3 reads kongsberg motion data and corrects vertical velocity
%   cd c:\data\cwf\matlabstf\cwf\VOCALS2008_programs_cwf\wband
%The read_parameters_VOCALSz.m file sets the paths for the radar data
% 96  W-band radar files
% 97  way_raw_data_wband=[big_dir '\RHB\radar\wband\Raw\'];
% 98  way_raw_images_wband=[big_dir '\RHB\radar\wband\Raw_images\'];
% You must go to lines 97-98 and set the path for your data archive
%close all
clear
%read_parameters_1;%this program just calls read_parameters_VOCALSz
%cd F:\CALNEX_2010\Atlantis\Scientific_Analysis\programs\Wband\
% read_parameters;
% disp('flocmoi');

% %******************   set paths
% cruise='CALNEX_2010';     % Cruise name, for example 'STRATUS_2009', 'TOTO_1969', etc...
% year='2010';               % year (string)
% ship='Atlantis';               % Research vessel
% bigdisc='F:\';
% %******************* W-band radar files  *********************************
% disc_radar='F:';
% way_raw_data_wband=[disc_radar '\CALNEX_2010\Atlantis\radar\wband\Raw\'];
% way_raw_images_wband=[disc_radar  '\CALNEX_2010\Atlantis\radar\wband\Raw_images\'];
% way_cloud_data_wband=[disc_radar '\CALNEX_2010\Atlantis\radar\wband\cloud\'];
% way_raw_motion_wband=[disc_radar '\CALNEX_2010\Atlantis\radar\wband\kongs\'];
% %##################%  ceilo files %####################################
% way_images_ceilo=[bigdisc cruise  '\' ship '\ceilometer\Raw_Images\']; % daily saved impates
% way_raw_data_ceilo=[bigdisc cruise  '\' ship '\ceilometer\Raw\'];  % Raw Data 
% way_proc_data_ceilo=[bigdisc cruise  '\' ship '\ceilometer\Processed\']; % Processed Data

%################## Cruise informations  ##################################################################
cruise='Hurricane';             % Cruise name  
year='2013';               % year (deprecated--make year into a double) 
yyyy='2013';               % year (string)
ship='P3';               % Research vessel
big_dir='H:\';
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\'];

% %  ceilo files
% way_images_ceilo=[big_dir '\RHB\ceilometer\Raw_Images\']; % daily saved impates
% way_raw_data_ceilo=[big_dir '\RHB\ceilometer\Raw\'];  % Raw Data 
% way_proc_data_ceilo=[big_dir '\RHB\ceilometer\Processed\']; % Processed Data
% 
% way_proc_data_micro=[big_dir '\RHB\radiometer\Mailbox-UM\Processed\'];  % Raw Data 
% 
% 



motionread=0;%set to zero to skip motion correction

ddds = input('Input yearday  ','s');
ddd = str2num(ddds);
shrs = input('Input start hr  ','s');
shr = str2num(shrs);
ehrs = input('Input end hr  ','s');
ehr = str2num(ehrs);
%eval(['cd ', way_raw_data_wband]);%
eval(['cd ', 'C:\data\cwf\dropbox\data\VOCALS_2008\']);
save dew1 ddds ddd shrs ehrs shr ehr
%eval(['cd ', 'F:\VOCALS_2008\RHB\radar\wband\']);

%way_raw_images_wband= 'F:\VOCALS_2008\RHB\radar\wband\Raw_images\';
%way_raw_images_wband= 'c:\data\cwf\data\vocals\';

%ncload('C:\Data\VOCALS\Data\wband\NetCDFRadarmfiles\WbandMom.nc')  % First rain shower

%20082821346MMCRMom.nc  20082821400MMCRMom.nc% 2nd rain shower

ht = (1:1:120)*.025; % range gates in km
fn = [way_raw_data_wband '\20132071600MMCRMom.nc'];
fn(41:43) = ddds;

dr=dir([way_raw_data_wband '\2013' ddds shrs '*.nc']);
film=dr.name;
jamn=1;
cldtop=[];DVx=[];Refx=[];hhx=[];SNRx=[];SWx=[];jdx=[];
for jam = shr:ehr
    save dew2 jam
    %close all
    clear hh y z dt jd
    %read_parameters;%read_parameters_1
    load dew1
    load dew2


    %way_raw_images_wband= 'F:\VOCALS_2008\RHB\radar\wband\Raw_images\';
    %way_raw_images_wband= 'C:\Data\cwf\DATA\VOCALS_2008\cband_images';

    %ncload('C:\Data\VOCALS\Data\wband\NetCDFRadarmfiles\WbandMom.nc')  % First rain shower

    %20082821346MMCRMom.nc  20082821400MMCRMom.nc% 2nd rain shower

    ht = (1:1:150)*.025; % range gates in km
    zmxx=ht(150);
    dbz0=-38;
    %fn(28:30) = ddds;

    if jam<10,
        hr=['0' num2str(jam)];
    else
        hr=num2str(jam);
    end;
    jamn=jamn+1;
    film=dr(jamn).name
    fn = [way_raw_data_wband '\2013' ddds hr film(10:11) 'MMCRMom.nc'];
    FID = netcdf.open(fn,'NC_NOWRITE');%fopen(fn,'r');
    if FID > 0  % Test to see if file exists
       % ncload(fn)
         base_time=netcdf.getVar(FID,0); % bug for multiple files in 1 hour!   0
        % read 1D vectors
        time_offset=netcdf.getVar(FID,1)';
        tt = time_offset + double(base_time);
        % calculate time form time_offset and base_time
        y=[];z=[];scs=[];minx=[];sec0=[];
        for i = 1:length(tt)
            %X = base_time;
            y(i) = tt(i)/(60*60*24);
            z(i) = y(i)+datenum([1970,1,1,0,0,0]);
            dt = datestr(z(i), 31);    %2008-10-07 16:35:53
            yr = str2num(dt(1:4));
            month = str2num(dt(6:7));
            day = str2num(dt(9:10));
            hour = str2num(dt(12:13));
            minutes = str2num(dt(15:16));
            if i==1;sec0=str2num(dt(18:19));end;
            secs = mod(time_offset(i)+sec0,60);
            jd(i) = datenum(yr,month,day,hour,minutes,secs)-datenum(yr-1,12,31);
            minx(i)=minutes;
            scs(i)=secs;
            %correct for
        end
        hh = (jd-(fix(jd(1))))*24; % subtract off JD and convert to an hour
        NNX=ncread(fn,'NumHeights');nnx=NNX(1);
        z=netcdf.getVar(FID,27);% 27
        ht=z(:,1)'/1e3;
        Ref=netcdf.getVar(FID,40);  %40
        SNR=netcdf.getVar(FID,36); %36
        SW=netcdf.getVar(FID,38);  % 38
        DVa=netcdf.getVar(FID,35); %35
        NoiseLevel=netcdf.getVar(FID,39); %35
        netcdf.close(FID);
%         Ref = Reflectivity';
%         SW = SpectralWidth';
%         DVa = MeanDopplerVelocity';
%         SNR = SignalToNoiseRatio';
        DV=[];
        if motionread%correct doppler for
               read_Wband_motion_P3
               DV_mot=interp1(hhk,vkons(9,:),hh);
               ii=find(isnan(DV_mot));DV_mot(ii)=0;
               for im=1:length(ht);
                    DV(im,:)=DVa(im,:)+DV_mot;
               end;
        else
            DV=DVa;
        end;
        
        clear Reflectivity SpectralWidth MeanDopplerVelocity SignalToNoiseRatio
        [nx,mx]=size(Ref);
        for k=1:mx;
            kk=find(Ref(:,k)-dbz0-20*log10(ht')>8 & ht'>.5 & ht'<zmxx & Ref(:,k)>-30);
            if ~isempty(kk)
                htt(k)=ht(kk(length(kk)));
                reft(k)=mean(10.^(Ref(kk,k)/10));
                reft2(k)=mean(10.^((Ref(kk,k)-dbz0-20*log10(ht(kk)'))/10));
                wvel(k)=mean(DV(kk,k));
                nnt(k)=length(kk);
            else
                htt(k)=NaN;
                reft(k)=NaN;
                reft2(k)=NaN;
                wvel(k)=NaN;
                nnt(k)=0;
            end;
        end;
        y = prctilex(htt,[10 15 50 85 90]);
        kk=find(htt>=y(1) & htt<=y(5));
        httb=nanmeanx(htt(kk));
        reftb=10*log10(nanmeanx(reft(kk)));
        reft2b=10*log10(nanmeanx(reft2(kk)));
        wvelb=nanmeanx(wvel(kk));
        nntb=sum(nnt);
        sigh=sqrt(nanmeanx(htt(kk).^2)-httb.^2);
        ['Hour= ' num2str((hh(1))) '  Mean cloud top ht= ' num2str(httb) '  std cloud top height = ' num2str(sigh)]
        ['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=[hhx hh];
        DVx=[DVx DV];
        Refx=[Refx Ref];
        SNRx=[SNRx SNR];
        SWx=[SWx SW];
        jdx=[jdx jd];
        ymaxx=max(ht);
        figure(4)
        clf
        subplot(3,1,1),imagesc(hh,ht, Ref(1:nnx,:),[-60, 20]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        tit = {[cruise  ' ' year ' W-Band'' YD ' ddds ' Hr ' hr]};
        title(tit)
        h = colorbar;
        axes(h)
        ylabel('Ref dBZ')

        subplot(3,1,2),imagesc(hh,ht, DV(1:nnx,:),[-7, 7]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        h = colorbar;
        axes(h)
        ylabel('DopVel (m/s)')


        subplot(3,1,3),imagesc(hh,ht, SW(1:nnx,:),[0, 5]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        xlabel('Hrs (UTC)')
        h = colorbar;
        axes(h)
        ylabel('SpecWid (m/s)')
        print('-djpeg90 ',[way_raw_images_wband '\Wband_ref_' num2str(ddd), '_' hr '.jpg']);

        figure(5)
        clf
        subplot(3,1,1),imagesc(hh,ht, SNR(1:nnx,:),[-25, 20]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        %axis([15 15.25 0 3])
        tit = {[cruise ' ' year ' W-Band'' YD ' ddds ' Hr ' hr]};
        title(tit)
        h = colorbar;
        axes(h)
        ylabel('SNR dBZ')

        subplot(3,1,2),imagesc(hh,ht, DV(1:nnx,:),[-3, 5]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        %axis([15 15.25 0 3])
        h = colorbar;
        axes(h)
        ylabel('DopVel (m/s)')
%         
%         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;
%         axes(h)
%         ylabel('DopVel (m/s)')

        subplot(3,1,3),imagesc(hh,ht, SW(1:nnx,:),[0, 5]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        xlabel('Hrs (UTC)')
        h = colorbar;
        axes(h)
        ylabel('SpecWid (m/s)')
        print('-djpeg90 ',[way_raw_images_wband '\Wband_snr_' num2str(ddd), '_' hr '.jpg']);

        %fclose(FID);

    end % test for file

end % hr loop

%       
