% 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);
        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
             
%         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')



        fclose(FID);
        close all
    end % test for file

end % hr loop


%Read Ceilo daily
ceilo_data=load([way_proc_data_ceilo 'cloudbase20090713_194.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==10);
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*');
print('-djpeg90 ',[way_raw_images_wband '\Wband_ceilo_' num2str(ddd), '_' hr '.jpg']);
   
