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

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

way_raw_images_wband= 'H:\RHB\radar\wband\Raw_images\';
%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 = 'C:\Data\VOCALS\data\wband\raw\20082821400MMCRMom.nc';
fn(35:37) = ddds;

for jam = shr:ehr
    save dew2 jam
    close all
    clear
    read_parameters_VOCALS
    load dew1
    load dew2


    way_raw_images_wband= 'H:\RHB\radar\wband\Raw_images\';
    %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 = 'C:\Data\VOCALS\data\wband\raw\20082821400MMCRMom.nc';
    fn(35:37) = ddds;

    if jam<10,
        hr=['0' num2str(jam)];
    else
        hr=num2str(jam);
    end;
    fn(38:39) = hr
    FID = fopen(fn,'r');
    if FID > 0  % Test to see if file exists
        ncload(fn)

        tt = time_offset + base_time;
        % calculate time form time_offset and base_time
        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';
        DV = MeanDopplerVelocity';
        SNR = SignalToNoiseRatio';

        clear Reflectivity SpectralWidth MeanDopplerVelocity SignalToNoiseRatio

        figure(4)
        subplot(3,1,1),imagesc(hh,ht, Ref(1:120,:),[-60, 20]);  %
        axis xy
        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:120,:),[-7, 7]);  %
        axis xy
        ylabel('Ht (km)')
        h = colorbar;
        axes(h)
        ylabel('DopVel (m/s)')


        subplot(3,1,3),imagesc(hh,ht, SW(1:120,:),[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(5)
        clf
        subplot(3,1,1),imagesc(hh,ht, SNR(1:120,:),[-25, 20]);  %
        axis xy
        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: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:120,:),[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);

    end % test for file

end % hr loop
