% Reads in W-Band nc moment data and plot
% range gates are 25 m
% timing is every 0.25 sec

%% initialize
clear;
close('all');
fclose('all');

run D:\DATA\BigSur_2018\Scientific_analysis\programs\read_parameters_BigSur_2018.m;
% colormap jet  %forces colormat to jet and be used full scale (has some glitches for some obscure reasons)
% cmap=colormap;
% cd([matlab_path, 'wband'])   

%################## Defines some options and parameters ##################################
Nrange=120;
RangeRes=0.031;
ht = (1:1:Nrange)*RangeRes; % range gates in km
zmxx=ht(Nrange);
ymaxx=2.5;
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;
graphformat='.png';  %select graphics format files
graphdevice='-dpng'; %select graphic device

motionread=0;%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)


% preallocate x-extended variables
% nxmax=175000;
nxmax=300000;
[Refx, SNRx, SWx, DVx]=deal(single(zeros(Nrange,nxmax))+NaN);
[jdx, hhx]=deal(single(zeros(1,nxmax))+NaN);

%%%%%%%%%%%%%%%%%%%%%%%%% Starts process %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~ automatically retrieve the previous day for process:
dd = clock;  % date vector today
edn = datenum(dd(1),dd(2),dd(3),0,0,0); % MATLAB day
bdn = edn-1;    % this is to get the prevuous days data.
ddv = datevec(bdn);
yyyy =  sprintf('%04i',ddv(1));
m = sprintf('%02i',ddv(2));
d = sprintf('%02i',ddv(3));
ddd = datenum(ddv(1),ddv(2),ddv(3))-datenum(ddv(1),1,0);
ddds= sprintf('%03i',ddd);
shr=0;ehr=23;
shrs= sprintf('%02i',shr);
ehrs= sprintf('%02i',ehr);
Vdate=ddv(1:3);

%~~~~ user inputs day to be processed:
% ddds = input('Input yearday  ','s');
% ddd = str2num(ddds);
% % [m,d]=yd2md(Vdate(1), ddd);
% % Vdate=[str2num(yyyy) m d];
% shrs = input('Input start hr  ','s');
% shr = str2num(shrs);
% ehrs = input('Input end hr  ','s');
% ehr = str2num(ehrs);

handl(1)=figure(1); % open 2 figures whether anything to plot
% handl(2)=figure(2);
handl(3)=figure(3);

% set(handl(:),'WindowStyle','docked');

% cldtop=[];DVx=[];Refx=[];hhx=[];SNRx=[];SWx=[];jdx=[];
nx=0; % count position in extended variables
for jam = shr:ehr
%     clf(handl(1:2));
    clear hh y z dt jd
    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])

        if exist(fn)  % Test to see if file exists
            %reads in variable of interest
            base_time=ncread(fn,'base_time');
            time_offset=ncread(fn,'time_offset');    
            z=ncread(fn,'Heights');
            Ref=ncread(fn,'Reflectivity'); %dBz
            SNR=ncread(fn,'SignalToNoiseRatio'); %dB
            SW=ncread(fn,'SpectralWidth'); %m/s
            DVa=ncread(fn,'MeanDopplerVelocity'); %m/s
            nrg=ncread(fn,'NumHeights'); %number of range gates
            %         Nrange=nrg(1); clear nrg

            % calculate time from time_offset and base_time
            tt = time_offset + double(base_time);
        
            %compute doy
            tt=time_offset+double(base_time);
            nta=length(tt); % time length of current file
            tti=tt/(60*60*24)+datenum([1970,1,1,0,0,0]);
            jd=tti-datenum(str2num(yyyy)-1,12,31); clear tti
            hh = (jd-(fix(jd(1))))*24; % subtract off JD and convert to an hour 
        
            %converts height to km and compute range resolution
            ht=z(:,1)/1e3;
            RangeRes=ht(2)-ht(1);
            %if RangeRes>0.04;ymaxx=ht(Nrange);end;

            %correct doppler for ship motion (if on ship)
            DV=single([]);
            if motionread%correct doppler for
                   read_Wband_motion_lb
                   try
    %                 DV_mot=single(interp1(hhk,vkons(9,:),hh));
                        fkactual=length(hhk)./3600;
                        hhk=vkons(4,:)+(0:1/fkactual:(length(hhk)-1)/fkactual)/3600;
                        DV_mot=single(interp1(hhk,vkons(9,:),sort(hh)));
                   catch
                    DV_mot=0*DVa(1,:)   ;
                    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;
            DV=single(DV);
        
            %insert data in extended daylong array(s)
            [m1,n1]=size(Ref);
            if m1>0 && n1>0
                try
                    jdx(:,nx+(1:nta))=jd;
                    hhx(:,nx+(1:nta))=hh;
                    Refx(1:m1,nx+(1:nta))=Ref;
                    DVx(1:m1,nx+(1:nta))=DV;
                    SNRx(1:m1,nx+(1:nta))=SNR;
                    SWx(1:m1,nx+(1:nta))=SW;
                    nx=nx+nta;
                catch
                    modeflag=1;  %true for dissimilar range modes in 1 day
                    [m0,n0] = size(Refx);[m1,n1] = size(Ref);
                    if m1<m0
                        Refx(1:m1,nx+(1:nta))=[Ref; NaN*ones(abs(m1-m0),n1)];               
                        DVx(1:m1,nx+(1:nta))=[DV; NaN*ones(abs(m1-m0),n1)];
                        SNRx(1:m1,nx+(1:nta))=[SNR; NaN*ones(abs(m1-m0),n1)];
                        SWx(1:m1,nx+(1:nta))=[SW; NaN*ones(abs(m1-m0),n1)];
                        nx=nx+nta;
    %                 elseif m1>m0  %case there is more gates than the 220
    %                 setup at beginning of code. Hopefuly doen's occur
    %                     Refx=single([Refx; NaN*ones(abs(m1-m0),n0)]);
    %                     Refx=[Refx single(Ref)];
                    end   
                end    
            end
             
            %% mask background noise from dpvl
            dbz0=-42;
            abovenoise = bsxfun(@minus, Ref, dbz0+20*log10(ht)) > 8.0;
            noisemask=zeros(size(abovenoise));
            noisemask(~abovenoise)=NaN;
%             noisemaskx(1:m1,nx+(1:nta))=noisemask;
            % compute profiles
            [rr,cc] = size(Ref);  % rr = # ht bins, cc = # timeseries samples
            z2 = NaN(cc,1); refl2 = z2; refl3 = z2; DVm = z2; nnt = z2;
            for k=1:cc
                kk = find(abovenoise(:,k) & ht>.1 & ht<zmxx & Ref(:,k)>-30);
                if ~isempty(kk)
                    z2(k) = ht(kk(end));
%                     refl2(k) = mean(10.^(Ref(kk,k)/10));
%                     refl3(k) = mean(10.^((Ref(kk,k)-dbz0-20*log10(ht(kk)))/10));
%                     DVm(k)=mean(DV(kk,k));
%                     nnt(k)=length(kk);
                else
                    z2(k) = NaN;
%                     refl2(k) = NaN;
%                     refl3(k) = NaN;
%                     DVm(k) = NaN;
%                     nnt(k) = 0;
                end
            end

            ncol = 32;
            crange = @(z1, z2, ncol) [z1-(z2-z1)/ncol, z2];
        
            hh = (jd-(fix(jd(1))))*24;
            
            %~~~~~~~~ do hourly plots
            %figure 1 with Reflectivity
            figure(1); clf
            colormap([[0 0 0]; b2rcolormap(ncol+1)]);
            ax(1)=subplot(4,1,1);imagesc(mod(hh,1)*60,ht, Ref,crange(-60, 20, ncol));  %
            axis xy;ylim([0 ymaxx]);
%             xlim([ddd+jam/24 ddd+(jam+1)/24]);
            ylabel('Ht (km)');
            title(sprintf('%s (%04i-%02i-%02i, DOY%03i, Hr-%02i).  W-Band',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,jam),'FontWeight','Bold','Interpreter','none')
            h1 = colorbar; set(h1,'ylim',[-60 20]); ylabel(h1,'Ref dBZ');
%             set(gca(gcf),'XTick',ddd+jam/24:10/60/24:ddd+(jam+1)/24); datetick('x','HH:MM','keepticks');
            hold on;

            ax(2)=subplot(4,1,2);imagesc(mod(hh,1)*60,ht, noisemask+DV,crange(-7, 7, ncol));  %
            axis xy;ylim([0 ymaxx]);%xlim([ddd+jam/24 ddd+(jam+1)/24]);
            ylabel('Ht (km)')
            h2 = colorbar; set(h2,'ylim',[-7 7]); ylabel(h2,'DopVel (m/s)')
%              set(gca(gcf),'XTick',ddd+jam/24:10/60/24:ddd+(jam+1)/24); datetick('x','HH:MM','keepticks');
            hold on;

            ax(3) = subplot(4,1,3);imagesc(mod(hh,1)*60,ht, noisemask+SW,crange(0, 5, ncol));  %
            axis xy;ylim([0 ymaxx]);%xlim([ddd+jam/24 ddd+(jam+1)/24]);
            ylabel('Ht (km)')
%             xlabel('Hrs (UTC)')
            h3 = colorbar; set(h3,'ylim',[0 5]); ylabel(h3, 'SpecWidth (m/s)');
%              set(gca(gcf),'XTick',ddd+jam/24:10/60/24:ddd+(jam+1)/24); datetick('x','HH:MM','keepticks');
            hold on;

            ax(4) = subplot(4,1,4);imagesc(mod(hh,1)*60,ht, SNR,crange(-25, 20, ncol));  %
            axis xy;ylim([0 ymaxx]);%xlim([ddd+jam/24 ddd+(jam+1)/24]);
            ylabel('Ht (km)');
            h4 = colorbar; set(h4,'ylim',[-25 20]); ylabel(h4,'SNR dBZ');
%              set(gca(gcf),'XTick',ddd+jam/24:10/60/24:ddd+(jam+1)/24); datetick('x','HH:MM','keepticks');
            hold on;
            xlabel('minutes')
            annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
                {'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',8,'FitBoxToText','off','LineStyle','none');      
            set(ax(:),'color',[0 0 0],'fontsize',16,'xlim',[-eps 60+eps],...
            'tickdir','out','ytick',0:6,'yticklabel',{'0','','2','','4','','6'});
            orient tall
            if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_Wband_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]);end

        end % test for file
    end % one or multiple files in hour
end % hr loop
% close all

% % truncate extended arrays
% Refx=Refx(:,1:nx);
% DVx=DVx(:,1:nx);
% SNRx=SNRx(:,1:nx);
% SWx=SWx(:,1:nx);
% jdx=jdx(:,1:nx);     
% hhx=hhx(:,1:nx);     
%05292018 - doesn't ehlp here as create out of memory issues...

%%%%%%%%%%%%%% PLOT full time series %%%%%%%%%%%%%% 
close('all');%clear temp DVx;
% ymaxx=7.0;
ii=1:2:length(jdx);
jd_ref=ddd:1/24:ddd+1;

%figure 4 with Ref, Doppler Vel and Spectral width
figure(4); 
colormap([[0 0 0]; b2rcolormap(ncol+1)]);
for ll=1:length(jd_ref)-1
    toto=find(jdx(ii)>=jd_ref(ll) & jdx(ii)<=jd_ref(ll+1));

    dbz0=-42;
    abovenoise = bsxfun(@minus, Refx(1:Nrange,ii(toto)), dbz0+20*log10(ht)) > 8.0;
    noisemask=zeros(size(abovenoise));
    noisemask(~abovenoise)=NaN;


    if length(toto)>2
        subplot(3,1,1),imagesc(jdx(ii(toto)),ht, Refx(1:Nrange,ii(toto)),crange(-60, 20,ncol));  %
        axis xy;ylim([0 ymaxx]);xlim([ddd ddd+1]);
        ylabel('Ht (km)');
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  W-Band',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        h = colorbar;        h.Label.String= 'Ref dBZ';
        set(gca(gcf),'XTick',ddd:1/24:ddd+1); datetick('x','HH','keepticks');
        hold on;

        subplot(3,1,2),imagesc(jdx(ii(toto)),ht, noisemask+DVx(1:Nrange,ii(toto)),crange(-7, 7,ncol));  %
        axis xy;ylim([0 ymaxx]);xlim([ddd ddd+1]);
        ylabel('Ht (km)')
        h = colorbar;        h.Label.String= 'DopVel (m/s)';
        set(gca(gcf),'XTick',ddd:1/24:ddd+1); datetick('x','HH','keepticks');
        hold on;

        subplot(3,1,3),imagesc(jdx(ii(toto)),ht, noisemask+SWx(1:Nrange,ii(toto)),crange(-5, 5,ncol));  %
        axis xy;ylim([0 ymaxx]);xlim([ddd ddd+1]);
        ylabel('Ht (km)')
        xlabel('Hrs (UTC)')
        h = colorbar;        h.Label.String= 'SpecWid (m/s)';
        set(gca(gcf),'XTick',ddd:1/24:ddd+1); datetick('x','HH','keepticks');
        hold on;
    end
end
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
% colormap(cmap)
if prtit; print(graphdevice,[way_proc_images_wband cruise '_Wband_ref_spec_summary_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end

clear temp DVx noisemask SWx SNRx
%figure 4 with Ref, Doppler Vel and Spectral width
%*******************  Load raw 15 s ceilo file
fceilo=[way_proc_data_ceilo sprintf('cloudbase_%04s_%02s_%02s_%03i.txt',yyyy,m,d,ddd)]; %[code height1 height2 height3 hour min sec yearday]
if exist(fceilo)
    ceilo=load(fceilo);
    doy_ceilo=ceilo(:,8)+(ceilo(:,5)+ceilo(:,6)/60+ceilo(:,7)/3600)/24;
else
    ceilo=single(zeros(24*60*60/15,8))+NaN;
    doy_ceilo=single(zeros(1,24*60*60/15))+NaN;
end;
clear fceilo 

%*******************  Load raw radiometer file
fmicro=fullfile(way_raw_data_mb,'ColInt_Obs.txt');
if exist(fmicro)
    mp = read_col_int(ddd, fmicro);
else
    mp=single(zeros(1000,4))+NaN;
end
micro_doy = mp(:,1);
micro_vap = mp(:,2); %cm
micro_liq = mp(:,3);  %mm
micro_cld_base = mp(:,4); %km
clear fmicro mp

%imagesc seems to have issues with dscoutinuos time series. It concatenates
%all and doesn't plot gap well. So do plot on an hourly basis for now.
% [b, m, n]=unique(jdx);
% jd_ref=ddd:.3/60/60/24:ddd+1;
% jd_ref=ddd:1/24:ddd+1;
figure(3);clf;
colormap([b2rcolormap(ncol+1)]);
h1=subplot(2,1,1);

for ll=1:length(jd_ref)-1
    toto=find(jdx(ii)>=jd_ref(ll) & jdx(ii)<=jd_ref(ll+1));
    if length(toto)>2
        hold on;hi=imagesc(jdx(ii(toto)),ht, Refx(1:Nrange,ii(toto)),crange(-60, 20,ncol));  %  
%         set(hi, 'AlphaData', ~isnan(Refx(1:Nrange,ii(toto)))); %set NaN to blank in case of any
    end    
end    
axis xy;ylim([0 ymaxx]);
xlim([ddd ddd+1]);
ylabel('Ht (km)');xlabel('Hour (UTC)');
title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  W-Band',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
h = colorbar;     
h.Label.String= 'Ref dBZ';
set(gca(gcf),'XTick',ddd:1/24:ddd+1); datetick('x','HH','keepticks');
% hold on;plot(micro_doy,micro_cld_base,'.',doy_ceilo,ceilo(:,2:4)/1000,'k+')
hold on;plot(doy_ceilo,ceilo(:,2:4)/1000,'k.')
hold off
% temp=get(h1,'Position')

pause(.1) %pause for 1/10 of a second as subplot(1) width doesn't seem to be detected right. 
% It still thinks it's the width without the colorbar...

h2=subplot(2,1,2);plot(micro_doy,micro_liq.*1e3,'.');box on;
h2.Position(3)=h1.Position(3)+0.035;%sets same width for figure
set(h2,'Position',h2.Position)
xlim([ddd ddd+1]);
% ylim([-20 300]);axis([ddd ddd+1 -20 300]);
set(gca(gcf),'XTick',ddd:1/24:ddd+1); datetick('x','HH','keepticks');
xlabel('Hour (UTC)');ylabel('LWC (g.m^{-2})');grid
% h3=get(gca,'position'); h3(3)=h1.Position(3); %sets same width for figure
% set(gca,'position',h3)
%       print('-dpng ',[way_raw_images_wband '\Wband_summary' num2str(ddd), '.png']);
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',8,'FitBoxToText','off','LineStyle','none'); 
% orient landscape
% colormap(cmap)
if prtit; print(graphdevice,[way_proc_images_wband cruise '_Wband_ref_cloud_IWC_summary_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end




exit