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

LegNoStr = '4';

% set day-of-year
ddd = 240;
hhh = 21;

cruise = 'PISTON_MISOBOB_2018';     % Cruise name
yyyy = '2018';        % year string
ship = 'TGT';                       % 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 'byronb'
            data_drive = '/Users/byronb/Raw_Data_Files.noindex';
        case 'sdeszoek'
            data_drive = '/Users/sdeszoek/Data/cruises/';
    end
    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_',LegNoStr]);
path_plots = fullfile(data_drive,cruise,ship,'Wband',['NoiseLevelPlots_',LegNoStr]);
mkdir(path_plots);

dbz0 = -42;
zmax = 6.8;

%% load data
dddstr = sprintf('%03i',ddd);
hrstr = sprintf('%02i',hhh);
fdat = dir(fullfile(path_mom_data,[yyyy,dddstr,hrstr,'*MMCRMom.nc']));

base_time = ncread(fullfile(path_mom_data,fdat(1).name),'base_time');
time_offset = ncread(fullfile(path_mom_data,fdat(1).name),'time_offset');
if ~exist('z','var')
    z = ncread(fullfile(path_mom_data,fdat.name),'Heights');
    z = z(:,1); z = z/1000; % height bins in km
end
snr  = ncread(fullfile(path_mom_data,fdat.name),'SignalToNoiseRatio');
refl = ncread(fullfile(path_mom_data,fdat.name),'Reflectivity');
spw  = ncread(fullfile(path_mom_data,fdat.name),'SpectralWidth');
dpvl = ncread(fullfile(path_mom_data,fdat.name),'MeanDopplerVelocity');
noise = ncread(fullfile(path_mom_data,fdat.name),'NoiseLevel');
noise2 = noise + 20*log10(z);

base_time = double(base_time);
tt = base_time + time_offset;
ts = tt/86400 + datenum([1970,1,1,0,0,0]);
jd = ts - datenum(2017,12,31);
hh = (jd-(fix(jd(1))))*24;
tmin = mod(hh,1)*60;

ymaxx = 6.5;  % max altitude for plot y axes

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

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(2018,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, noise2,crange(50, 100, ncol));
axis xy; ylim([0 ymaxx]); ylabel('Ht (km)');
h2 = colorbar; set(h2,'ylim',[50 100]); ylabel(h2,'Noise Level');
ax(3)=subplot(4,1,3); imagesc(mod(hh,1)*60,z, 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
ppath = fullfile(path_plots,[cruise,'_',sprintf('%03i',ddd),'_',hrstr,'.png']);
print('-dpng',ppath);

