%  PISTON_MISOBOB_Wband_plot.m
%  reworked version of DYNAMO Wband script

%%  initialize run parameters
% clear;
% close('all');
fclose('all');
warning ('off','MATLAB:MKDIR:DirectoryExists');

% set leg number manually, 1-5, as defined below
% leg 2 155:174, MISOBOB1 Chennai to Colombo
% leg 3 181:201, MISOBOB2 Colombo to Colombo - no data for Wband
% leg 4 232:252, PISTON2 Taiwan to Palau
% leg 5 255:290, PISTON3 Palau to Taiwan

LegNo = 5;
LegNoStr = '5';

% set day-of-year
st_day = 267;  %%% no data before 253 hr 2
end_day = 267;
st_hour = 00;
end_hour = 23;

prtit = true;  % for saving plots--plotit must also be true
graphformat = '.png';  % select graphics format files
graphdevice = '-dpng'; % select graphic device
cruise = 'PISTON_2019';     % Cruise name
year = 2019;
yyyy = sprintf('%04d',year);        % year string
ship = 'Sally';                       % Research vessel

% system specific path defs
sysType = computer;
username=char(java.lang.System.getProperty('user.name'));
if strncmp(sysType,'MACI64',7)     % set Mac paths
%     switch username
%         case 'ludo'
%             data_drive = '/Users/eliz/DATA/';
%         case 'sdeszoek'
%             data_drive = '/Users/sdeszoek/Data/cruises/';
%     end
    data_drive = '/Users/eliz/DATA/';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs');
elseif strncmp(sysType,'PCWIN',7)  % set PSD DAC paths
    data_drive = 'D:\Data';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs');
end

% matlab script path
restoredefaultpath
cd(fullfile(path_prog,'Wband'));
addpath(genpath(fullfile(path_prog,'Wband')));
rehash toolboxcache;

% paths relative to data_drive
path_mom_data = fullfile(data_drive,cruise,ship,'Wband',['Mom']);
path_mru_data = fullfile(data_drive,cruise,ship,'Wband',['MRU']);
path_radar_plots = fullfile(data_drive,cruise,ship,'Wband',['HourlyImagePlots']);
mkdir(path_radar_plots);
path_motion_plots = fullfile(data_drive,cruise,ship,'Wband',['DailyMotionPlots']);
mkdir(path_motion_plots);
path_heave_plots = fullfile(data_drive,cruise,ship,'Wband',['HourlyHeavePlots']);
mkdir(path_heave_plots);

dbz0 = -42;
zmax = 6.8;

% load shift values for days prior to jd 251 on leg 4 (PISTON leg 2)
% lags = load(fullfile(data_drive,cruise,ship,'Wband','lag_all_leg5.txt'));
% ii = isfinite(lags(:,2));
% % default lag shift will be linear fit to time drift from prior runs
% coefs = polyfit(lags(ii,1),lags(ii,2),1);
%% loop through days

clear z
for ddd = st_day:end_day
    dddstr = sprintf('%03i',ddd);
    nn = 24*6; nn_idx=1;
    mot_stats = NaN(nn,10); % empty array for 10-min motion stats
    lag_hr = NaN(24,2); lag_hr(:,1) = 0:23; % array for saving lag shift values

    %% Loop through hours
    % cldtop=[]; DVx=[]; Refx=[]; hrx=[]; SNRx=[]; SWx=[]; jdx=[];
    for hhh = st_hour:end_hour

        %% read radar data
        jd_bin = (ddd+hhh/24):10/1440:(ddd+(hhh+1)/24); jd_bin = jd_bin';  % time bins for 10-min std dev
        hrstr = sprintf('%02i',hhh);
        fdat = dir(fullfile(path_mom_data,[yyyy,dddstr,hrstr,'*MMCRMom.nc']));
        if isempty(fdat)  % continue to next hour if no radar files for this hour
            nn_idx = nn_idx+6;
            continue
        end
        for filei = 1:length(fdat)  % handles multiple netCDF files per day...
            ncf = fullfile(path_mom_data,fdat(filei).name);
            if filei == 1
                base_time = double(ncread(ncf,'base_time'));
                time_offset = ncread(ncf,'time_offset');
                tt = base_time + time_offset;
                if ~exist('z','var')
                    z = ncread(ncf,'Heights');
                    z = z(:,1); z = z/1000; % height bins in km
                end
                snr  = ncread(ncf,'SignalToNoiseRatio');
                refl = ncread(ncf,'Reflectivity');
                spw  = ncread(ncf,'SpectralWidth');
                dpvl = ncread(ncf,'MeanDopplerVelocity');
            else
                base_time = double(ncread(ncf,'base_time'));
                time_offset = ncread(ncf,'time_offset');
                tt =   [tt  ; base_time + time_offset ];
                snr  = [snr , ncread(ncf,'SignalToNoiseRatio')];
                refl = [refl, ncread(ncf,'Reflectivity')];
                spw  = [spw , ncread(ncf,'SpectralWidth')];
                dpvl = [dpvl, ncread(ncf,'MeanDopplerVelocity')];
            end
        end
        ts = tt/86400 + datenum([1970,1,1,0,0,0]);
        [they, them, thed, theh, themin, thes] = datevec(ts(1));
        year_add = datenum(they-1,12,31,0,0,0);
        jd = ts - year_add;
        hh = (jd-(fix(jd(1))))*24;

        %% read motion file - here we assume only one Kongsberg file per hour.
        % In case of multiple files use text editor first to combine into single
        % .txt file for each hour.
        fname = dir(fullfile(path_mru_data,[yyyy,dddstr,hrstr,'*Kongsberg.txt']));
        dfl = fullfile(fname(1).folder, fname(1).name);  % assume only one mot file per day
%         disp(['reading ' dfl]);
        mot = read_kongsberg_PISTON_MISOBOB_2018(dfl,ddd,hhh);

        % min-max roll,pitch, & heave from Kongsberg
        minall = interval_min(mot(:,1),mot(:,2:4),jd_bin);
        maxall = interval_max(mot(:,1),mot(:,2:4),jd_bin);
        meanall = interval_avg(mot(:,1),mot(:,2:4),jd_bin);
        mot_stats(nn_idx:nn_idx+5,:) = [meanall(:,1:4),minall(:,2:4),maxall(:,2:4)];
    
        %%% example: 
%         mot_stats(1,:) = 
%   264.0000   -0.0259   -0.0052   -0.0031   -0.1690   -0.1096   -1.5900    1.0000    0.0897    1.5050
%        jd    r_mean     p_mean    h_mean    r_min          p_min         h_min         r_max        p_max         h_max
        nn_idx = nn_idx+6;

        % interp platform motion to radar timestamp
        % negate Kongsberg heave velocity since it's mounted upside down...
        % EJT is this correct? Is it really mounted upside down? 
        Vplat = interp1(mot(:,1),-mot(:,4), jd,'linear')';

        %% mask background noise from dpvl
        abovenoise = bsxfun(@minus, refl, dbz0+20*log10(z)) > 8.0;
        dpvl_mask = dpvl;
        dpvl_mask(~abovenoise) = NaN;
        dpvl_m = nanmean1(dpvl_mask,1);         % use only good dpvl returns . %% EJT 9/7 changed from nanmeanx to nanmean1
        dpvl_m(1,sum(abovenoise)<20) = NaN;
        dpvl_mx = dpvl_m;
        isf = isfinite(dpvl_mx);
        dpvl_mx(~isf) = mean(dpvl_m(isf));     % no NaNs can be in highpass
        hi = highpass(double(dpvl_mx),5,2/60); % to remove drift

        %% fix lag between motion and radar dpvl for PISTON 2018. ET commented out... ?? 
        
% %         if ddd<269     % default shift for jd < 251
% %             jj= find(lags(:,1)<269);
% %             ii = isfinite(lags(jj,2));
% %             % default lag shift will be linear fit to time drift from prior runs
% %             coefs = polyfit(lags(ii,1),lags(ii,2),1);
% %             dflt = round(coefs(1)*(ddd+hhh/24) + coefs(2));
% % %         else
% % %             dflt = 0;  % clock should be sync'ed after jd 251
% % %             dflt=12; %doen't seem to be the case for leg 5... LB - 09162018
% % %         end 
% %         else
% %             dflt = 1;  %LB reset manually on DOY 269
% %         end
        
           dflt = 1; 

        ii = (sum(abovenoise,1)>10 & isfinite(dpvl_m) & isfinite(Vplat));
        if sum(ii) > 1200   % sufficient cloud to compute xcov
            nlags = 150;
            [XYcorr,~] = xcov(hi(ii)',Vplat(ii)',nlags,'coeff');
            [minCorr,idx] = max(XYcorr);
            shift = idx - (nlags+1);  % shift is # data points @ 2 Hz
            if abs(shift-dflt)>4      % reject outlier shift values
                shift = dflt;                   % use default
                lag_hr(hhh+1,2) = 999+dflt;     % mark default by adding 999
            else
                lag_hr(hhh+1,2) = shift;        %  xcov shift is OK, so use
            end
        else    % insufficient cloud for xcov, so use default shift value
            shift = dflt;      % default for days when clock was drifting
            lag_hr(hhh+1,2) = 999+shift;  % mark as default by adding 999
        end

        Vplat2 = Vplat;
        if shift>0
            for ii=length(Vplat):-1:shift+1
                Vplat2(ii) = Vplat(ii-shift);
            end
        elseif shift<0
            for ii=1:length(Vplat)+shift
                Vplat2(ii) = Vplat(ii-shift);
            end
        end
        
        % plot vplat and mean dpvl as .fig to see if phase is correct
        %%% EJT: I'm not sure what is correct here.
        figure; 
        plot(mod(hh,1)*60,hi,'r-',mod(hh,1)*60,Vplat2,'b-'); grid;
        legend('dpvl','vplat'); xlabel('minutes'); ylabel('Vert Vel, m/s');
        title(sprintf('Heave Check: DOY %03i, Hr %s',ddd,hrstr),...
            'FontWeight','Bold','Interpreter','none');
        if prtit
            ppath = fullfile(path_heave_plots,['heave_check_',sprintf('%03i',ddd),'_',hrstr,graphformat]);            
            print(graphdevice,ppath,'-r300');
        end

        % remove heave from dpvl matrix
        mu = 1.0;  % assume unity calibration coefficient between velocities
        DV = dpvl - mu*Vplat2;

        %% mask the corrected Doppler velocity, and remove the residual mean
        DVmask = DV;
        DVmask(~abovenoise) = NaN;
        DVmean = nanmeanx(DVmask);  %% EJT 9/7 changed from nanmeanx to nanmean1
        DVfilt = bsxfun(@minus, DV, DVmean); % filters out residual column mean velocity
        %%% but this isn't used...? 

        % compute profiles
        [rr,cc] = size(refl);  % 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) & z>.1 & z<zmax & refl(:,k)>-30);
            if ~isempty(kk)
                z2(k) = z(kk(end));
                refl2(k) = mean(10.^(refl(kk,k)/10));
                refl3(k) = mean(10.^((refl(kk,k)-dbz0-20*log10(z(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

        %% some statistics on the radar data
        y = prctilex(z2,[10 15 50 85 90]);
        kk = find(z2>=y(1) & z2<=y(5));
        z2m = nanmeanx(z2(kk));
        refl2m = 10*log10(nanmeanx(refl2(kk)));
        refl3m = 10*log10(nanmeanx(refl3(kk)));
        DVmm = nanmeanx(DVm(kk));
        nntt = sum(nnt);
        sigh = sqrt(nanmeanx(z2(kk).^2)-z2m.^2);

        disp(['Hour = ',num2str(hhh),',  Mean cloud top ht = ',num2str(z2m),',  std cloud top height = ',num2str(sigh)])
        disp(['Mean dBZ = ',num2str(refl2m),',   Mean dBZ-noisedBZ = ',num2str(refl3m),',   Mean Vel = ',num2str(DVmm)])

        %     cldtop = [cldtop;[hhh,z2m,sigh,(y(4)-y(3))/2,refl2m,refl3m,DVmm,nntt]];
        %     jdx = [jdx,jd'];     % decimal day array for all hours
        %     hrx = [hrx,hh'];     % decimal hour array for all hours
        %     DVx = [DVx,DV];      % 220 x N doppler velocity array
        %     Refx = [Refx,refl];  % 220 x N reflectivity array
        %     SNRx = [SNRx,snr];   % 220 x N SNR array
        %     SWx = [SWx,spw];     % 220 x N spectral width array

        ymaxx = 6.5;  % max altitude for plot y axes

        noisemask=zeros(size(abovenoise));
        noisemask(~abovenoise)=NaN;

        ncol = 32;
        crange = @(z1, z2, ncol) [z1-(z2-z1)/ncol, z2];

        %% plot radar images
        figure;
        colormap([[0 0 0]; b2rcolormap(ncol+1)]);
        ax(1)=subplot(4,1,1); imagesc(mod(hh,1)*60,z, refl,crange(-60, 20, ncol));
        axis xy; ylim([0 ymaxx]); ylabel('Ht (km)');
        title(sprintf('PISTON %s DOY %03i, Hr %s.  W-Band',...
            datestr(datenum(they,0,0)+ddd),ddd,hrstr),...
            'FontWeight','Bold','Interpreter','none');
        h1 = colorbar; set(h1,'ylim',[-60 20]); ylabel(h1,'Ref dBZ');
        ax(2)=subplot(4,1,2); imagesc(mod(hh,1)*60,z, noisemask-DV,crange(-7, 7, ncol));
        axis xy; ylim([0 ymaxx]); ylabel('Ht (km)');
        h2 = colorbar; set(h2,'ylim',[-7 7]); ylabel(h2,'DopVel (m/s)');
        ax(3)=subplot(4,1,3); imagesc(mod(hh,1)*60,z, noisemask+spw,crange(0, 5, ncol));
        axis xy; ylim([0 ymaxx]); ylabel('Ht (km)');
        h3 = colorbar; set(h3,'ylim',[0 5]); ylabel(h3, 'SpecWidth (m/s)');
        ax(4)=subplot(4,1,4); imagesc(mod(hh,1)*60,z, snr,crange(-25, 20, ncol));
        axis xy; ylim([0 ymaxx]); ylabel('Ht (km)');
        h4 = colorbar; set(h4,'ylim',[-25 20]); ylabel(h4,'SNR dBZ');
        xlabel('minutes')
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
            {'NOAA/ESRL/PSD'},'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
            ppath = fullfile(path_radar_plots,[cruise,'_',sprintf('%03i',ddd),'_',hrstr,graphformat]);
            print(graphdevice,ppath,'-r300');
        end

        close all
    end

    % mot_stats = [jd,mean_roll,mean_pitch,mean_heave,min_roll,min_pitch,min_heave,max_roll,max_pitch,max_heave]
    %               1     2         3          4         5        6         7         8        9         10

    figure;
    set(gca,'fontsize',14);
    subplot(3,1,1); errorbar(mot_stats(:,1),mot_stats(:,2)*180/pi,mot_stats(:,5)*180/pi,mot_stats(:,8)*180/pi,'b.-'); grid;
    xlim([ddd,ddd+1]); ylabel('Roll, deg');ylim([-20 20]);
    set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
    xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
    title(sprintf('PISTON %s DOY %03i.  Radar Angles & Heave',datestr(datenum(they,0,0)+ddd),...
        ddd),'FontWeight','Bold','Interpreter','none','fontsize',16);
    subplot(3,1,2); errorbar(mot_stats(:,1),mot_stats(:,3)*180/pi,mot_stats(:,6)*180/pi,mot_stats(:,9)*180/pi,'b.-'); grid;
    xlim([ddd,ddd+1]);  ylabel('Pitch, deg');ylim([-20 20]);
    set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
    xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
    subplot(3,1,3); errorbar(mot_stats(:,1),mot_stats(:,4),mot_stats(:,7),mot_stats(:,10),'b.-'); grid; ylim([-2 2]);
    xlim([ddd,ddd+1]); ylabel('Heave, m/s'); xlabel('Hour of Day, UTC');
    set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
    xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
    if prtit
        ppath = fullfile(path_motion_plots,['motion_stats_',sprintf('%03i',ddd),graphformat]);
        print(graphdevice,ppath);
    end
    
    default = lag_hr; ii = default(:,2)>900;
    kk = lag_hr(:,2)<900;
    figure; 
    plot(lag_hr(kk,1),lag_hr(kk,2),'bd',default(ii,1),default(ii,2)-999,'ro',[0,24],[0,0],'k-'); grid;
    xlim([0,24]); xticks(0:2:24); ylim([-5,40]); yticks(-5:5:40); set(gca,'yminortick','on');
    xlabel('Hour'); ylabel('vplat & dpvl shift, points @ 2 Hz'); set(gca,'fontsize',14);
    legend('xcorr lag','default lag'); title('Platform : Radar W-Velocity Lag Correction');
    if prtit
        ppath = fullfile(path_motion_plots,['lag_correction_',sprintf('%03i',ddd),graphformat]);
        print(graphdevice,ppath);
    end
    close all
    
    %%% this is just the day... at 1 hr increments. 
    lag_hr(:,1) = ddd+lag_hr(:,1)/24;  % to jd
    lag_hr(ii,2) = NaN;  % save the xcorr lags for this day
    dfl = fullfile(path_motion_plots,['lag_',sprintf('%03i',ddd),'.txt']);
    save(dfl,'lag_hr','-ascii');
    
end

filecat([fullfile(path_motion_plots,['lag_*.txt'])],[fullfile(data_drive,cruise,ship,'Wband','lag_all_leg5.txt')])
