close all
clear all

%% Set values for isolume plotting as per Båtnes et al. (2013)

BatnesMin = log10(10^(-8));
BatnesMax = log10(10^(-6));
BatnesMean = log10(10^(-7));

%% Set up figure position
figure
set(gcf, 'Units', 'Centimeters', 'Pos', [1 1 20 24]);
ax1 = axes;
set(ax1,'Pos',[0.1 0.88 0.8 0.06]);
ax2 = axes;
set(ax2,'Pos',[0.1 0.68 0.8 0.18]);
ax3 = axes;
set(ax3,'Pos',[0.1 0.56 0.8 0.06]);
ax4 = axes;
set(ax4,'Pos',[0.1 0.36 0.8 0.18]);
ax5 = axes;
set(ax5,'Pos',[0.1 0.24 0.8 0.06]);
ax6 = axes;
set(ax6,'Pos',[0.1 0.04 0.8 0.18]);

%% Load data: ADCP (extracted midnight and midday), solar elevation (midnight and midday), time variable, depth variable, flu data
load DataFigure2

%% Load data

% Each year contains: 
% Sv_night, t, bins (Sv normalised between years as per Wallace et al., 2010)
% Flu, mtime


                                     %% 2013-14 %%

%% Load and filter (3 day) fluorometry data
mtime = Data2013.mtime;;
flu = Data2013.flu;

sid = 1/(mtime(2)-mtime(1));
cut_off_sid = 1/(3*24);
cut_off_Hz = cut_off_sid/24/60/60;
Hz = sid/24/60/60;
Wn = cut_off_Hz / Hz * 2;
[A,B] = butter(1,Wn,'low');
flu1 = filter(A,B,flu);

%% Scale fluorometry data by sample in Hesgeth et al. (2019)

FluJul = mean(flu1(find(mtime > datenum(2014,07,23) & mtime < datenum(2014,07,25))));
FluJulObs = 3.4; 
Scaling = FluJulObs/FluJul;
flu1 = flu1*Scaling;

%% Calculate Kd

Kd = 0.0313 * flu1 + 0.15;
Kd = interp1(mtime,Kd,Data2013.t);

%% Calculate solar elevation. Uses function SolarAzEl https://uk.mathworks.com/matlabcentral/fileexchange/23051-vectorized-solar-azimuth-and-elevation-estimation

[Az El_night] = SolarAzEl(Data2013.t', 78.9265, 11.8050, 0);
clear Az

%% Apply light model

a0    =  7.3168;
b0  =  3.2856;
El0  =   -4.9378;
y0   =  -5.5882;
logIrrad = y0+a0./(1+exp(-(El_night-El0)./b0)); 
E_night = 10.^logIrrad; 

[tt,zz] = meshgrid(Data2013.t,-Data2013.bins);
[kpar,zz] = meshgrid(Kd, -Data2013.bins);

[En,zz] = meshgrid(E_night,-Data2013.bins);
En = En .* exp(kpar .* zz);

%% Plotting

% ADCP + isolumes
axes(ax4)
hold on
contourf(tt,zz,Data2013.Sv_night,-95:5:45,'edgecolor','none');
ylim([-110 0]);
xlim([datenum(2013,10,07);datenum(2014,09,09)]);
datetick('x','mmm','keeplimits');
c = colorbar;
title(c,'S_v (dB)')
cpos = get(c,'Pos');
set(c,'Pos',[cpos(1)+0.08 cpos(2) cpos(3) cpos(4)]);
caxis([-95 -45]);
contour(tt,zz,log10(En),[BatnesMin BatnesMax],'k','LineStyle',':', 'LineWidth',0.6);
contour(tt,zz,log10(En),[BatnesMean BatnesMean],'k', 'LineWidth',1.2)
ylabel('Depth (m)');
colormap(warmrainbow)
text(datenum(2013,10,10),-10,'d)');

% Fluorometry
axes(ax3);
plot(mtime, flu1, 'k', 'LineWidth',1.2);
datetick('x','mmm','keeplimits');
xlim([datenum(2013,10,07);datenum(2014,09,09)]);
ylabel('Chl mg m^-^3');
hold on
box off
set(gca,'XTick',[],...
    'XColor','none');
text(datenum(2013,10,31),max(flu1)*0.8,strcat('Midnight  2013-14; Chl_m_a_x = ',num2str(round(max(flu1),3,'Significant')),' \mu gL^-^1'), 'Color', 'k', 'Fontweigh','bold', 'FontSize',10);
text(datenum(2013,10,10),max(flu1)*0.8,'c)', 'Color', 'k', 'FontSize',10);


                 %% 2007-08 %%

%% Load and filter (3 day) fluorometry data
mtime = Data2007.mtime;;
flu = Data2007.flu;

sid = 1/(mtime(2)-mtime(1));
cut_off_sid = 1/(3*24);
cut_off_Hz = cut_off_sid/24/60/60;
Hz = sid/24/60/60;
Wn = cut_off_Hz / Hz * 2;
[A,B] = butter(1,Wn,'low');
flu1 = filter(A,B,flu);

%% Scale fluorometry data by sample in Hesgeth et al. (2019)

FluJul = mean(flu1(find(mtime > datenum(2008,05,01) & mtime < datenum(2008,05,03))));
FluJulObs = 2; % From Hegseth paper
Scaling = FluJulObs/FluJul;
flu1 = flu1*Scaling;

%% Calculate Kd

Kd = 0.0313 * flu1 + 0.15;
Kd = interp1(mtime,Kd,Data2007.t);

%% Calculate solar elevation. Uses function SolarAzEl https://uk.mathworks.com/matlabcentral/fileexchange/23051-vectorized-solar-azimuth-and-elevation-estimation

[Az El_night] = SolarAzEl(Data2007.t', 78.9265, 11.8050, 0);
clear Az

%% Apply light model

a0    =  7.3168;
b0  =  3.2856;
El0  =   -4.9378;
y0   =  -5.5882;
logIrrad = y0+a0./(1+exp(-(El_night-El0)./b0)); 
E_night = 10.^logIrrad; 

[tt,zz] = meshgrid(Data2007.t,-Data2007.bins);
[kpar,zz] = meshgrid(Kd, -Data2007.bins);

[En,zz] = meshgrid(E_night,-Data2007.bins);
En = En .* exp(kpar .* zz);

%% Plotting

% ADCP + isolumes
axes(ax6)
hold on
contourf(tt,zz,Data2007.Sv_night,-95:5:45,'edgecolor','none');
ylim([-110 0]);
xlim([datenum(2007,10,07);datenum(2008,09,09)]);
datetick('x','mmm','keeplimits');
c = colorbar;
title(c,'S_v (dB)')
cpos = get(c,'Pos');
set(c,'Pos',[cpos(1)+0.08 cpos(2) cpos(3) cpos(4)]);
caxis([-95 -45]);
contour(tt,zz,log10(En),[BatnesMin BatnesMax],'k','LineStyle',':', 'LineWidth',0.6);
contour(tt,zz,log10(En),[BatnesMean BatnesMean],'k', 'LineWidth',1.2);
ylabel('Depth (m)');
colormap(warmrainbow)
text(datenum(2007,10,10),-10,'f)');

% Fluorometry
axes(ax5);
plot(mtime, flu1, 'k', 'LineWidth',1.2);
datetick('x','mmm','keeplimits');
xlim([datenum(2007,10,07) datenum(2008,09,09)]);
ylabel('Chl mg m^-^3');
hold on
box off
set(gca,'XTick',[],...
    'XColor','none');
text(datenum(2007,10,31),max(flu1)*0.8,strcat('Midnight  2007-08; Chl_m_a_x = ',num2str(round(max(flu1),3,'Significant')),' \mu gL^-^1'), 'Color', 'k', 'Fontweigh','bold', 'FontSize',10);
text(datenum(2007,10,10),max(flu1)*0.8,'e)', 'Color', 'k', 'FontSize',10);


                       %% 2008-09 %%

%% Load and filter (3 day) fluorometry data
mtime = Data2008.mtime;;
flu = Data2008.flu;

sid = 1/(mtime(2)-mtime(1));
cut_off_sid = 1/(3*24);
cut_off_Hz = cut_off_sid/24/60/60;
Hz = sid/24/60/60;
Wn = cut_off_Hz / Hz * 2;
[A,B] = butter(1,Wn,'low');
flu1 = filter(A,B,flu);

%% Scale fluorometry data by sample in Hesgeth et al. (2019)

FluJul = mean(flu1(find(mtime > datenum(2009,07,17) & mtime < datenum(2009,07,20))));
FluJulObs = 2.1; % From Hegseth paper
Scaling = FluJulObs/FluJul;
flu1 = flu1*Scaling;

%% Calculate Kd

Kd = 0.0313 * flu1 + 0.15;
Kd = interp1(mtime,Kd,Data2008.t);

%% Calculate solar elevation. Uses function SolarAzEl https://uk.mathworks.com/matlabcentral/fileexchange/23051-vectorized-solar-azimuth-and-elevation-estimation

[Az El_night] = SolarAzEl(Data2008.t', 78.9265, 11.8050, 0);
clear Az

%% Apply light model

a0    =  7.3168;
b0  =  3.2856;
El0  =   -4.9378;
y0   =  -5.5882;
logIrrad = y0+a0./(1+exp(-(El_night-El0)./b0)); 
E_night = 10.^logIrrad; 

[tt,zz] = meshgrid(Data2008.t,-Data2008.bins);
[kpar,zz] = meshgrid(Kd, -Data2008.bins);

[En,zz] = meshgrid(E_night,-Data2008.bins);
En = En .* exp(kpar .* zz);


%% Plotting

% ADCP + isolumes

axes(ax2)
hold on
contourf(tt,zz,Data2008.Sv_night,-95:5:45,'edgecolor','none');
ylim([-110 0]);
xlim([datenum(2008,10,07);datenum(2009,09,09)]);
datetick('x','mmm','keeplimits');
c = colorbar;
title(c,'S_v (dB)')
cpos = get(c,'Pos');
set(c,'Pos',[cpos(1)+0.08 cpos(2) cpos(3) cpos(4)]);
caxis([-95 -45]);
contour(tt,zz,log10(En),[BatnesMin BatnesMax],'k','LineStyle',':', 'LineWidth',0.6);
contour(tt,zz,log10(En),[BatnesMean BatnesMean],'k', 'LineWidth',1.2);
ylabel('Depth (m)');
colormap(warmrainbow)
text(datenum(2008,10,10),-10,'b)');

axes(ax1);
plot(mtime, flu1, 'k', 'LineWidth',1.2);
datetick('x','mmm','keeplimits');
xlim([datenum(2008,10,07);datenum(2009,09,09)]);
ylabel('Chl mg m^-^3');
hold on
box off
set(gca,'XTick',[],...
    'XColor','none');
text(datenum(2008,10,31),max(flu1)*0.8,strcat('Midnight 2008-09; Chl_m_a_x = ',num2str(round(max(flu1),3,'Significant')),' \mu gL^-^1'), 'Color', 'k', 'Fontweigh','bold', 'FontSize',10);
text(datenum(2008,10,10),max(flu1)*0.8,'a)', 'Color', 'k', 'FontSize',10);


%%
print Figure2 -dpng -r300

%% COLORMAP %%

function cm = warmrainbow(n,brightness,contrast);

% warmrainbow;
% warmrainbow(n);
% warmrainbow(n,brightness,contrast);
%
% better than jet. Still looks like a rainbow when n ~ 10; minimizes the
% luminance peak at cyan, and puts the luminance peak at yellow in the middle;
% and still looks like a rainbow when the value channel is made uniform, so
% shades of this colormap can be used to encode a second information dimension).
%
% neil banas, dec 2006
%
% default n=20
% default brightness (value of yellow) = .9
% default contrast (difference between yellow value and red-blue value) = .3

if nargin < 3, contrast = .3; end
if nargin < 2, brightness = .9; end
if nargin < 1, n = 20; end

resz = 0;
if ~mod(n,2) % if n is even, calculate an odd-length colormap and then shorten
	% 				the red side: this preserves the yellow color
	resz = 1;
	n=n+1;
end

Hbyr = [.55 1/6 0]; % blue yellow red hues
Vbyr = [brightness-contrast, brightness, brightness-contrast]; % blue yellow red values
W = 1.5; % width of the dip in saturation over green

x = 0:n-1;
nh = floor(n/2);
nh2 = ceil(n/2);
hue = nan.*ones(n,1); sat = hue; val = hue;

N=n-1;
b = (8*Hbyr(2)-7*Hbyr(1)-Hbyr(3))/N^2;
a = (Hbyr(3)-Hbyr(1))/N^3 - b/N;
hue = a.*x.^3 + b.*x.^2 + Hbyr(1);
hue(nh+1:end) = linspace(hue(nh+1),Hbyr(3),nh2);

x0 = x(find(abs(hue-1/3)==min(abs(hue-1/3))));
sat = 1-.5*exp(-(W.*(x./x0-1)).^2);

b = 4*(Vbyr(2)-Vbyr(3))/N;
c = Vbyr(3);
a = -b/N;
val = a.*x.^2 + b.*x + c;
if Vbyr(1)~=Vbyr(2)
	val(1:nh) = linspace(Vbyr(1),Vbyr(2),nh);
else 
	val(1:nh) = Vbyr(1);
end

cm = hsv2rgb([hue(:) sat(:) val(:)]);

if resz % subsampling red side because an even length was requested
	cm = interp1((1:n)',cm,[1:nh linspace(nh+1,n,nh)]');
end	
end