% KM read in W-Band nc moment data and plot
% range gates are 25 m
% timing is every 0.25 sec --> check??
close all
clear

%%%%%%%%%%%  Read parameters  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% addpath('G:\WHOTS_2009\Kilo_Moana\WHOTS09_programs\') %add path for the read_parameters function
read_parameters_WHOTS09

motionread=1;%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 ', 'H:\RHB\radar\wband\']);
% save dew1 ddds ddd shrs ehrs shr ehr

ngates=198; %define number of gates to plot.

ht = (1:1:ngates)*.025; % range gates in km

for jam = shr:ehr
    hr=sprintf('%02i',jam);   
    fn=[way_raw_data_wband yyyy ddds hr '00MMCRMom.nc'];
    disp(['Reading: ' fn])
    FID = fopen(fn,'r');
    if FID > 0  % Test to see if file exists
%         ncload(fn,'Reflectivity', 'SpectralWidth','MeanDopplerVelocity','SignalToNoiseRatio')
        ncload(fn)
        
        tt = time_offset + base_time;
        % calculate time from time_offset and base_time
        dt=datestr(tt/(60*60*24)+datenum([1970,1,1,0,0,0]), 31);
        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));
        secs = str2num(dt(:,18:19));
        jd = datenum(yr,month,day,hour,minutes,secs)-datenum(yr-1,12,31);
 % Unnecessary looped removed! --> too slow        
% %         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));
% %             secs = str2num(dt(18:19));
% %             jd(i) = datenum(yr,month,day,hour,minutes,secs)-datenum(yr-1,12,31);
% %             %correct for
% %         end
        hh = (jd-(fix(jd(1))))*24; % subtract off JD and convert to an hour

        Ref = Reflectivity';
        SW = SpectralWidth';
        DVa = MeanDopplerVelocity';
        SNR = SignalToNoiseRatio';

        DV=[];
        if motionread%correct doppler for
               read_Wband_motion
               try
                   DV_mot=interp1(hhk,vkons(9,:),hh);
               catch
                   %fix 'he values of X should be distinct.' error message 071109 - LB
                   disp('error fixed')
                   [B,I,J] = UNIQUE(hhk);
                   DV_mot=interp1(hhk(I),vkons(9,I),hh);
               end;
               
               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
     
%         for k=1:length(Ref);
%             kk=find(Ref(:,k)+39-20*log10(ht')>8 & ht'>.5 & ht'<3 & 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)+39-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 = prctile(htt,[10 15 50 85 90]);
%         kk=find(htt>=y(1) & htt<=y(5));
%         httb=nanmean(htt(kk));
%         reftb=10*log10(nanmean(reft(kk)));
%         reft2b=10*log10(nanmean(reft2(kk)));
%         wvelb=nanmean(wvel(kk));
%         nntb=sum(nnt);
%         sigh=sqrt(nanmean(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]']';
        
        
        
        figure
        subplot(3,1,1),imagesc(hh,ht, Ref(1:ngates,:),[-60, 20]);  %
        axis xy
        ylabel('Ht (km)')
        tit = {[cruise '-' year ' W-Band / DOY ' ddds ' Hr ' hr]};
        title(tit)
        h = colorbar;
        axes(h)
        ylabel('Ref dBZ')

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

        subplot(3,1,3),imagesc(hh,ht, SW(1:ngates,:),[0, 5]);  %
        axis xy
        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
        clf
        subplot(3,1,1),imagesc(hh,ht, SNR(1:ngates,:),[-25, 20]);  %
        axis xy
        ylabel('Ht (km)')
        %axis([15 15.25 0 3])
        tit = {[cruise '-' year ' W-Band / DOY ' ddds ' Hr ' hr]};
        title(tit)
        h = colorbar;
        axes(h)
        ylabel('SNR dBZ')

        subplot(3,1,2),imagesc(hh,ht, DV(1:ngates,:),[-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:ngates,:),[0, 5]);  %
        axis xy
        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);
        close all
    end % test for file

end % hr loop




%Read Ceilo daily
% ceilo_data=load([way_proc_data_ceilo 'cloudbase20090710_191.txt']);
% cloudcode=ceilo_data(:,1);
% cloud1=ceilo_data(:,2);
% cloud2=ceilo_data(:,3);
% cloud3=ceilo_data(:,4);
% ceil_hr=ceilo_data(:,5);
% ceil_mn=ceilo_data(:,6);
% ceil_sc=ceilo_data(:,7);
% ceilo_doy=ceilo_data(:,8);
% 
% ii=find(ceil_hr==9);
% time_ceilo=ceil_hr+ceil_mn/60+ceil_sc/60/60;
% 
% figure
% imagesc(hh,ht, Ref(1:ngates,:),[-60, 20]);  %
% axis xy
% ylabel('Ht (km)')
% tit = {[cruise '-' year ' W-Band / DOY ' ddds ' Hr ' hr]};
% title(tit)
% h = colorbar;
% axes(h)
% ylabel('Ref dBZ')
% 
% hold on;plot(time_ceilo(ii),cloud1(ii)/1000,'k*')
%