clear all
close all
addpath(genpath('/Users/Mau/Documents/MATLABhelp_new'));

%% Purpose: Plot current data from ADCP at M1

bid='M1';
% if downloading data for multiple years on disk, e.g. for climatology analysis
% path_data=['/Users/Mau/Documents/MBARIprojects/M1/'];
% files=dir([path_data,'*.nc4']);
% filename=files(1).name;
% NC=ncinfo(([path_data,filename]));
% nf=length(files);

% read latest data (for the current year) in ods
path_data=['http://dods.mbari.org:80/opendap/data/ssdsdata/deployments/m1/202106/'];
filename='OS_M1_20210614tenMinute_CMSTV.nc';
NC=ncinfo(([path_data,filename]));
nf=1;

% % read latest data (for 2020) in ods
% path_data=['http://dods.mbari.org:80/opendap/data/ssdsdata/deployments/m1/202008/'];
% filename='OS_M1_20200825tenMinute_CMSTV.nc';
% NC=ncinfo(([path_data,filename]));
% nf=1;

% select path to save plots
% % path_plots='/Users/Mau/Documents/MBARIprojects/M1/plots/';
% % % path_plots='/Users/Mau/Documents/MBARIprojects/Topo1/M1/plots/';
% path_plots='/Users/Mau/Documents/MBARIprojects/PIE/M1/plots/';
% dname='PIE'; % deployment names for plot titles and saving

path_plots='/Users/Mau/Documents/MBARIprojects/MR/Deployment_conditions/T1/';
dname='CPF-MR-T1'; % deployment names for plot titles and saving
mon=10; day=8;

path_plots='/Users/Mau/Documents/MBARIprojects/MR/Deployment_conditions/T2/';
dname='CPF-MR-T2'; % deployment names for plot titles and saving
mon=10; day=13;

path_plots='/Users/Mau/Documents/MBARIprojects/MR/Deployment_conditions/T3/';
dname='CPF-MR-T3'; % deployment names for plot titles and saving
mon=10; day=15; % changed PV to 0.15m/s

path_plots='/Users/Mau/Documents/MBARIprojects/MR/Deployment_conditions/T4/';
dname='CPF-MR-T4'; % deployment names for plot titles and saving
mon=10; day=18; % changed PV to 0.15m/s

%% DATE=cell(nf,1);
Sal=cell(nf,1);
Temp=cell(nf,1);
Depth=cell(nf,1);

LON=cell(nf,1);
LAT=cell(nf,1);
DATE=cell(nf,1);
MONTH=cell(nf,1);
DOY=cell(nf,1);

for f=1; %4; %nf % number of files
    
    % filename=files(f).name; % uncomment if processing multiple files
    
    % Load CTD data
    T=squeeze(ncread([path_data,filename],'SEA_WATER_SALINITY_TM')); % readtable(filename);
    S=squeeze(ncread([path_data,filename],'SEA_WATER_TEMPERATURE_TM'));
    D=squeeze(ncread([path_data,filename],'DEPTH'));
%     D=squeeze(ncread([path_data,filename],'DEPTH1'));
    
    % Load ADCP data
    U=squeeze(ncread([path_data,filename],'eastward_sea_water_velocity_TM')); % in cm/s
    V=squeeze(ncread([path_data,filename],'northward_sea_water_velocity_TM'));
    D_adcp=squeeze(ncread([path_data,filename],'TM_DEPTH_adcp'));
    time_adcp=squeeze(ncread([path_data,filename],'tm_time_adcp')); % gregorian, units: "seconds since 1970-01-01 00:00:00";
    
    % Load GPS
    lat=ncread([path_data,filename],'LATITUDE'); % 'GPS_LATITUDE_TM');
    lon=ncread([path_data,filename],'LONGITUDE');
    
    % Load time and convert to matlab time (for ts)
    time_ctd=ncread([path_data,filename],'tm_time_ts'); % gregorian, units: "seconds since 1970-01-01 00:00:00";
    to=datenum([1970 1 1 0 0 0]);
    time=time_ctd./86400; %in fracion of day
    mtime=time+to; % adjust start time
    
    % Load time and convert to matlab time (for adcp data)
    time_adcp=ncread([path_data,filename],'tm_time_adcp'); % gregorian, units: "seconds since 1970-01-01 00:00:00";
    to=datenum([1970 1 1 0 0 0]);
    time=time_adcp./86400; %in fracion of day
    mtime_adcp=time+to; % adjust start time
    
    % Convert to LOCAL TIME (-7 GMT)
    ltime=mtime-7/24;
    mDATE=datevec(ltime); % OK, data every ~10'
    month=mDATE(:,2);
    Date=datetime(mDATE);
%     doy=day(Date,'dayofyear'); % for climatology average
    
    ltime_adcp=mtime_adcp-7/24;
    mDATE_adcp=datevec(ltime_adcp); % OK, data every ~10'
    
    
    % Replace missing value with NaN
    missing_value = -1.0E34;
    T(T==missing_value)=NaN;
    S(S==missing_value)=NaN;
    
    U(U==missing_value)=NaN;
    V(V==missing_value)=NaN;
    
    % Save in structure array multiple years of data
    Temp{f}=T;
    Sal{f}=S;
    Depth{f}=D;
    
    LON{f}=lon;
    LAT{f}=lat;
    DATE{f}=mDATE;
    MONTH{f}=month;
%     DOY{f}=doy;
    
end

%% Select DATES coincident with a certain time period
% indt=find(mDATE(:,1)==2020 & mDATE(:,2)==12 & mDATE(:,3)<=16);
% imagesc(U(:,indt));

% Select DATES coincident with Topo1 Deploy 4
% indt=find(mDATE(:,1)==2021 & mDATE(:,2)==4 & mDATE(:,3)>=18);
% imagesc(U(:,indt));

% Select DAY to plot, for NRT ADCP data at M1
DATE=datevec(ltime_adcp);

pi=find(DATE(:,1)==2021 & DATE(:,2)==mon & DATE(:,3)==day,1,'first');
pf=find(DATE(:,1)==2021 & DATE(:,2)==mon & DATE(:,3)==day,1,'last');

%% Select depth range to plot velocity sections
[nZ,nP]=size(U);
zi=1; zf=nZ;

%% Select current component;
Z=U(:,pi:pf)/100; varname_plot='U';
cmin=-0.2; cmax=0.2;

Z=V(:,pi:pf)/100; varname_plot='V';
cmin=-0.2; cmax=0.2;

% Z=v/100; varname_plot='v_{adcp}';
% cmin=-0.1; cmax=0.1;

% Z=(sqrt(u.^2+v.^2)/100).*atan2(v,u); varname_plot='u_h_{adcp}';
% cmin=-1; cmax=1;

[TIME,~]=meshgrid(mtime_adcp(pi:pf),zi:zf);
[~,DEPTH]=meshgrid(mtime_adcp(pi:pf),D_adcp);

[nz,np]= size(Z);
close all

figure(1)
X=TIME;
Y=DEPTH;

h=pcolor(X,Y,Z);
h.EdgeColor = 'none';

% Colorbar
cmap = cmocean('balance',22);
colormap(cmap)
% shading flat; % works with pcolor

hc=colorbar; 
set(get(hc,'Title'),'String','[m/s]','fontsize',20);
caxis([cmin cmax]); % upcast

datetick('x','mm/dd'); %,'keeplimits')
dt=1/24; % 10 min spacing for adcp data
xticks([X(1,1):dt:X(1,end)]);
xlim([X(1,1),X(1,end)]);
xtickangle(30)
% axis tight
axis ij

datetick('x','mm/dd HH:MM','keeplimits','keepticks')
 
ylabel('Depth [m]');
xlabel([num2str(DATE(pi,1)),' Local Time']);

% ylim([nanmin(Y(:)) nanmax(Y)]);
title([varname_plot,' ADCP @ M1 during ',dname])

print_setup(50,30)
% % maximize;
print([path_plots,varname_plot,'_',dname,'_4-',num2str(day),'.eps'],'-depsc2');

%% Plot Vertical Section of Current Speed (magnitude)

uu=U(:,pi:pf); vv=V(:,pi:pf);
Z=sqrt(uu.^2+vv.^2)/100; varname_plot='W current';
cmin=0; cmax=0.2;


[TIME,~]=meshgrid(ltime_adcp(pi:pf),zi:zf);
[~,DEPTH]=meshgrid(ltime_adcp(pi:pf),D_adcp);

[nz,np]= size(Z);
close all

figure(1)
X=TIME;
Y=DEPTH;


subplot(3,1,1)
dq=1;
mu=nanmean(uu(1:dq:end,1:dq:end),1)/100;
mv=nanmean(vv(1:dq:end,1:dq:end),1)/100;
quiver(TIME(1,:),TIME(1,:).*0,mu,mv,0,'k');

dt=1/24; % 10 min spacing for adcp data
xticks([X(1,1):dt:X(1,end)]);
xlim([X(1,1),X(1,end)]);
% datetick('x','','keepticks')
% xtickangle(30)

title([varname_plot,' ADCP @ M1 during ',dname,' ',num2str(mon),'-',num2str(day)])


subplot(3,1,2:3)

h=pcolor(X,Y,Z);
h.EdgeColor = 'none';

% hold on
% dq=10;
% quiver(X(1:dq:end,1:dq:end),Y(1:dq:end,1:dq:end),uu(1:dq:end,1:dq:end)/100,vv(1:dq:end,1:dq:end)/100,0.001,'k');

% To omit space between subplots
ha=get(gcf,'children');
set(ha(1),'position',[0.1 0.1 0.8 0.6])
set(ha(2),'position',[0.1 0.7 0.8 0.2])

% Colorbar
cmap = flipud(cbrewer('div','Spectral',22));
colormap(cmap)
% shading flat; % works with pcolor

hc=colorbar('horizontal'); 
% set(get(hc,'Title'),'String','[m/s]','fontsize',20);
ylabel(hc,'[m/s]','fontsize',20);
caxis([cmin cmax]); % upcast

datetick('x','mm/dd'); %,'keeplimits')
dt=1/24; % 10 min spacing for adcp data
xticks([X(1,1):dt:X(1,end)]);
xlim([X(1,1),X(1,end)]);
% xtickangle(30)
% axis tight
axis ij

datetick('x','HH:MM','keeplimits','keepticks')
 
ylabel('Depth [m]');
xlabel([num2str(DATE(pi,2)),'/',num2str(DATE(pi,3)),'/',num2str(DATE(pi,1)),' Local Time']);
xtickangle(30)

% % print_setup(40,30)
% % % maximize;
% % print([path_plots,varname_plot,'_',dname,'_',num2str(mon),'-',num2str(day),'.eps'],'-depsc2');
% % % postcript conversion often breaks!
% printpdf([path_plots,varname_plot,'_',dname,'_',num2str(mon),'-',num2str(day)]);

print_setup(60,40)
print([path_plots,varname_plot,'_',dname,'_',num2str(mon),'-',num2str(day),'.jpg'],'-djpeg');

%% Wind Rose from M1 data 

% Select depth range for wind rose, 20 bins ~=100m
z1=1;
zf=20*1; drange='z1-100m';% top 100m
% zf=20*5; drange='z1-500m';% top 500m

% wind rose of currents from surface to zf
v=V(z1:zf,pi:pf);
u=U(zi:zf,pi:pf);

% % depth averaged currents
% v=nanmean(V(z1:zf,pi:pf),1);
% u=nanmean(U(zi:zf,pi:pf),1);
% drange='zint_1-500m';

wdir=rad2deg(atan2(v(:),u(:)));
wspd=sqrt(u(:).^2+v(:).^2)/100; %in m/s

% Define colormap
w_edges=[0:0.05:0.3];
nw=length(w_edges);
cmap=cbrewer('div','Spectral',nw);

WindRose(wdir,wspd,'vWinds',w_edges,'colors',cmap,'lablegend','Current Speeds in m/s',...
    'TitleString',{['M1 ADCP ',drange,' during ',dname,' ',num2str(mon),'-',num2str(day)];''}) %'LegendType', 1 % for colorbar
                                                % 'LegendType', 0 % no
                                                % legend
print_setup(50,40)                                              
% print([path_plots,'SurfaceCurrentsRose_',drange,'_adcp_M1_during_',dname,'_',num2str(mon),'-',num2str(day),'.eps'],'-depsc2');
print([path_plots,'SurfaceCurrentsRose_',drange,'_adcp_M1_during_',dname,'_',num2str(mon),'-',num2str(day),'.jpg'],'-djpeg');
