%%  Script to read in the wave data of Europlatform and compute the timeseries of wave power

clear all
close all

%% toggle to cut off the days with offshore waves.
Opt_maxangle=1;
iweg=[];

%% Load data
% load Long term data
All=load('EURHSall.mat')
TAll=load('EURTall.mat')

%%% Load old
 VbHs_old=load('EURHS.mat')
 VbT_old=load('EURT.mat')
 VbDir_old=load('EURdir.mat')

%%% Load new.
load('F:\check_outs\tudelft\vlugtenburg\MeteoHydro\Data Europlatform\Hmo.mat')
Vb.Hs(:,1)=date; Vb.Hs(:,2)=Hm0;
load('F:\check_outs\tudelft\vlugtenburg\MeteoHydro\Data Europlatform\Tm02.mat')
VbT.T(:,1)=date; VbT.T(:,2)=Tm02;
load('F:\check_outs\tudelft\vlugtenburg\MeteoHydro\Data Europlatform\wavedir.mat')
VbDir.Dir(:,1)=date; VbDir.Dir(:,2)=wavedir;

%% Fill gaps in data

%%% GAp in the data 25/7-25/8 2009 ,
% get the data out of other (Waterbase file) file.
ind_gat=find(VbT.T(:,1) >= datenum('2009_07_25', 'yyyy_mm_dd') & VbT.T(:,1) <= datenum('2009_08_25', 'yyyy_mm_dd'))
for ig=1:length(ind_gat)
i_fill=find (VbT_old.T(:,1) ==  VbT.T(ind_gat(ig),1)) 
   if ~isempty(i_fill)
       VbT.T(ind_gat(ig),2)=VbT_old.T(i_fill,2);
    VbT.T(ind_gat(ig),2)
   end
end
% get the data out of other file.
ind_gat=find(VbDir.Dir(:,1) >= datenum('2009_07_25', 'yyyy_mm_dd') & VbT.T(:,1) <= datenum('2009_08_25', 'yyyy_mm_dd'))
for ig=1:length(ind_gat)
i_fill=find (VbDir_old.Dir(:,1) ==  VbDir.Dir(ind_gat(ig),1)) 
   if ~isempty(i_fill)
       VbDir.Dir(ind_gat(ig),2)=VbDir_old.Dir(i_fill,2);
    VbDir.Dir(ind_gat(ig),2)
   end
end

%%% GAp in the data beginning 2011,
% get the data form betgin 2011 out of other (Waterbase file) file.
ind_gat=find(VbT.T(:,1) >= datenum('2011_02_09', 'yyyy_mm_dd') & VbT.T(:,1) <= datenum('2011_02_22', 'yyyy_mm_dd'))
for ig=1:length(ind_gat)
i_fill=find (VbT_old.T(:,1) ==  VbT.T(ind_gat(ig),1)) 
   if ~isempty(i_fill)
       VbT.T(ind_gat(ig),2)=VbT_old.T(i_fill,2);
    VbT.T(ind_gat(ig),2)
   end
end
% get the data form betgin 2011 out of other file.
ind_gat=find(VbDir.Dir(:,1) >= datenum('2011_02_09', 'yyyy_mm_dd') & VbDir.Dir(:,1) <= datenum('2011_02_22', 'yyyy_mm_dd'))
for ig=1:length(ind_gat)
i_fill=find (VbDir_old.Dir(:,1) ==  VbDir.Dir(ind_gat(ig),1)) 
   if ~isempty(i_fill)
       VbDir.Dir(ind_gat(ig),2)=VbDir_old.Dir(i_fill,2);
    VbDir.Dir(ind_gat(ig),2)
   end
end

%%% Gap in the data end 2012,
% % extend
ind_langer=find(VbHs_old.Hs(:,1) > max(Vb.Hs(:,1)))
 Vb.Hs=[ Vb.Hs ; VbHs_old.Hs(ind_langer,:)];
 VbT.T=[ VbT.T ; VbT_old.T(ind_langer,:)];
 VbDir.Dir=[ VbDir.Dir ; VbDir_old.Dir(ind_langer,:)];


%%
theta0=310; % angle of shore-normal waves

% cut Vb data at first three years exactly
indVb=find(Vb.Hs(:,1)>=datenum([2009,08,01]) & Vb.Hs(:,1)<=datenum([2012,08,01]));
indVbT=find(VbT.T(:,1)>=datenum([2009,08,01]) & VbT.T(:,1)<=datenum([2012,08,01]));
% cut all data to 20 years exactly
indAll=find(All.Hs(:,1)>=datenum([1991,09,01]) & All.Hs(:,1)<=datenum([2011,09,01]));
indAllT=find(TAll.T(:,1)>=datenum([1991,09,01]) & TAll.T(:,1)<=datenum([2011,09,01]));

%%  plot wave signals
figure
plot(All.Hs(:,1),All.Hs(:,2));
datetick('x')
title({ ' Sig. Waveheight at Europlatform' })
hold on
plot([All.Hs(indAll(1),1) All.Hs(indAll(1),1)],[0 700],':r')
plot([All.Hs(indAll(end),1) All.Hs(indAll(end),1)],[0 700],':r')
grid on
ylabel ('H_s [cm]')
print HsAllTimeserie.jpg -djpeg

%% plot wave period
figure
plot(TAll.T(:,1),TAll.T(:,2));
datetick('x')
title({ ' T_{m02} at Europlatform' })
hold on
plot([TAll.T(indAllT(1),1) TAll.T(indAllT(1),1)],[0 10],':r')
plot([TAll.T(indAllT(end),1) TAll.T(indAllT(end),1)],[0 10],':r')
grid on
ylabel ('T [s]')
print TAllTimeserie.jpg -djpeg
TmeanAllyears=nanmean(TAll.T(:,2))

%% plot
figure
plot(Vb.Hs(:,1),Vb.Hs(:,2));
datetick('x')
title({ ' Sig. Waveheight at Europlatform' })
hold on
plot([Vb.Hs(indVb(1),1) Vb.Hs(indVb(1),1)],[0 700],':r')
plot([Vb.Hs(indVb(end),1) Vb.Hs(indVb(end),1)],[0 700],':r')
grid on
ylabel ('H_s [cm]')
print HsVbTimeserie.jpg -djpeg

%%

figure

plot(VbT.T(:,1),VbT.T(:,2));
datetick('x')
title({ 'T_{m02} at Europlatform' })
hold on
plot([VbT.T(indVbT(1),1) VbT.T(indVbT(1),1)],[0 10],':r')
plot([VbT.T(indVbT(end),1) VbT.T(indVbT(end),1)],[0 10],':r')
grid on
ylabel ('T [s]')
print TVbTimeserie.jpg -djpeg

figure

plot(VbDir.Dir(:,1),VbDir.Dir(:,2));
datetick('x')
title({ 'wavedir at Europlatform' })
% hold on
% plot([VbT.T(indVbT(1),1) VbT.T(indVbT(1),1)],[0 10],':r')
% plot([VbT.T(indVbT(end),1) VbT.T(indVbT(end),1)],[0 10],':r')
grid on
ylabel ('Dir [deg]')
print DirVbTimeserie.jpg -djpeg



%% Plot both distributions

Hbins=0:15:700;
Tbins=0:0.3:10;
figure
[n,xout] =hist(All.Hs(indAll,2),Hbins);
[n2,xout2] =hist(Vb.Hs(indVb,2),Hbins);
[n3,xout3] =hist(TAll.T(indAllT,2),Tbins);
[n4,xout4] =hist(VbT.T(indVbT,2),Tbins);
close

figure
subplot(121)
plot(xout,n/sum(n),'.-')
hold all
plot(xout2,n2/sum(n2),'.-')
grid on
xlabel ('H_s [cm]')
ylabel ('P(H_s) [-]')
legend('1991 to 2011','2009 to 2012')

subplot(122)
plot(xout3,n3/sum(n3),'.-')
hold all
plot(xout4,n4/sum(n4),'.-')
grid on
xlabel ('T_{m02} [s]')
ylabel ('P(T_{m02}) [-]')

print ProbHs.jpg -djpeg

%% Wave Power, total
% P = E.cg
% P  =
rho=1030;
g=9.81;
d=32; % Depth at Eur
ip=0;
ipy=0;
disp('calculating wave power')
for ii=1:length(Vb.Hs(:,1))
    ii
    ti=Vb.Hs(ii,1);
    ind = find (VbT.T(:,1)==ti);
    ind2 = find (VbDir.Dir(:,1)==ti);
%     if length (ind) == 1
%     Hi=Vb.Hs(ii,2);
%     Ti=VbT.T(ind,2);
%     Ei=1/8*rho*g*Hi.^2;
%     wi=2*pi./Ti;
%     ki=disper (wi, d, g);
%     ci=wi./ki;
%     ni=0.5*(1+2*ki.*d./sinh(2*ki.*d));
%     cgi=ni.*ci;
%     ip=ip+1;
%     VbP(ip,1)=ti;
%     VbP(ip,2)=Ei*cgi;
%     end
%     
    if length (ind) == 1 & length (ind2) == 1  % so if there is wave height , period and direction on a time
    Hi=Vb.Hs(ii,2)/100;  % make to cm!
    Ti=VbT.T(ind,2);
    Diri=VbDir.Dir(ind2,2);
    Ei=1/16*rho*g*Hi.^2;  % use 1/16 here because it is HS!
    wi=2*pi./Ti;
    ki=disper (wi, d, g);
    ci=wi./ki;
    ni=0.5*(1+2*ki.*d./sinh(2*ki.*d));
    cgi=ni.*ci;
    ipy=ipy+1;
    
    VbPy(ipy,1)=ti;
    VbPy(ipy,2)=Ei*cgi;
    
    if Opt_maxangle==1
       if cosd(Diri-theta0) <=-0.5   % set wave energy to 0 if wave angle is more than 15 deg offshore!
       VbPy(ipy,2)=0;
       iweg=iweg+1
       end
       %disp(['data removed ' num2str(Hi)])
    end
    VbPy(ipy,3)=VbPy(ipy,2)*cosd(Diri-theta0)*sind(Diri-theta0);   % Py positive for Northern waves (southward transport)
      end
end
disp('done')


%% Filter signal in a daiy signal

opt=1;  %% choose 1 for mean daily value, 2 for maximum daily value
dailyPy(:,1)=(floor(VbPy(1,1)):1:VbPy(end,1))';  % t vector with timesteps of a day

dt=VbPy(2,1)-VbPy(1,1);
dt_new=dailyPy(2,1)-dailyPy(1,1);
for i=1:length(dailyPy)-1
    date_i=dailyPy(i,1);
    indh=find(VbPy(:,1) > date_i  & VbPy(:,1) <= date_i  + dt_new);
    if opt==1
    dailyPy(i,2)=nanmean(VbPy(indh,2));
    dailyPy(i,3)=nanmean(VbPy(indh,3));
    end
    if opt==2
    dailyPy(i,2)=max(VbPy(indh,2));
    dailyPy(i,3)=max(VbPy(indh,3));
    end

end
%save('Pdaily.mat','dailyP')

% format Pdaily
%column 1 = time
% column 2 = total Wave Power
% column 3 = wave power avaliable for alonsghore transport Py


%%
figure
plot(dailyPy(:,1),dailyPy(:,2),'k')
hold on
plot(dailyPy(:,1),dailyPy(:,3),'b')
datetick('x','keeplimits')
grid on

%% Store data per survey period
load F:\check_outs\tudelft\vlugtenburg\datesandpairs.mat

% Total wave power ()

for iperiod=1:length(pairs)-1
    
    periodstart=datenum(dates(pairs(iperiod,2),:),'mm/dd/yyyy');
    periodend=datenum(dates(pairs(iperiod+1,2),:),'mm/dd/yyyy');
    
    %Pb_i=VbP((VbP(:,1)>=periodstart & VbP(:,1)<=periodend),2);
    Pby_i=dailyPy((dailyPy(:,1)>=periodstart & dailyPy(:,1)<periodend),:);
        
    if length(Pby_i) >= 1
    Pperiod(iperiod+1,1)=sum(Pby_i(:,2));  % Ptotal
    Pperiod(iperiod+1,2)=periodstart;         % tstart
    Pperiod(iperiod+1,3)=periodend;           % tend
    if ~isnan(Pperiod(iperiod+1,1)) 
    [Pperiod(iperiod+1,5),ind]=max(Pby_i(:,2));   % Pmax
    Pperiod(iperiod+1,6)=Pby_i(ind,1);              % t_Pmax
    else
     Pperiod(iperiod+1,5)  =NaN ;
     Pperiod(iperiod+1,6)   =NaN;
    end
    
    else
    Pperiod(iperiod+1,1)=NaN;
    Pperiod(iperiod+1,2)=periodstart;
    Pperiod(iperiod+1,3)=periodend;
    Pperiod(iperiod+1,5)=NaN;
    Pperiod(iperiod+1,6)=NaN;
    end
    
    PcumPeriod((dailyPy(:,1)>=periodstart & dailyPy(:,1)<periodend))=cumsum(abs(Pby_i(:,2))) % Timeseries of wave power building up from the start of the survey period
    
  end

Pperiod(:,4)=Pperiod(:,1)./(Pperiod(:,3)-Pperiod(:,2)); % Pmean 
Pperiod (1,:)=NaN;


save('Pperiod.mat','Pperiod')
%%

ind1=find (dailyPy(:,3)<0)
ind2=find (dailyPy(:,3)>0)
dailyPyN=dailyPy(:,3);
dailyPyN (ind1)=NaN;
dailyPyS=-dailyPy(:,3);
dailyPyS (ind2)=NaN;

% Bulk values alongshore wave power available for transport.

for iperiod=1:length(pairs)-1
    
    periodstart=datenum(dates(pairs(iperiod,2),:),'mm/dd/yyyy');
    periodend=datenum(dates(pairs(iperiod+1,2),:),'mm/dd/yyyy');
    
    %Pb_i=VbP((VbP(:,1)>=periodstart & VbP(:,1)<=periodend),2);
    Pby_i=dailyPy((dailyPy(:,1)>=periodstart & dailyPy(:,1)<periodend),:); 
        
    if length(Pby_i) >= 1
    Pperiody(iperiod+1,1)=sum(abs(Pby_i(:,3))); %Py total
    Pperiody(iperiod+1,2)=periodstart;  
    Pperiody(iperiod+1,3)=periodend;
  
    if ~isnan(Pperiody(iperiod+1,1)) 
    [Pperiody(iperiod+1,7),ind]=max(abs(Pby_i(:,3))); % max (|Py|)
    Pperiody(iperiod+1,6)=Pby_i(ind,1); % t (|Py| max)
    Pperiody(iperiod+1,5)=Pby_i(ind,3); % (Py max)
    else
    Pperiody(iperiod+1,7)=NaN;    
    Pperiody(iperiod+1,6)=NaN;
    Pperiody(iperiod+1,5)=NaN;
    end
    [Pperiody(iperiod+1,8),ind]=max(Pby_i(:,3));      % max Southward Wave power Py
    [Pperiody(iperiod+1,9),ind]=min(Pby_i(:,3));      % max Northward Wave power Py
    
    else
    Pperiody(iperiod+1,1)=NaN;
    Pperiody(iperiod+1,2)=periodstart;
    Pperiody(iperiod+1,3)=periodend;
    Pperiody(iperiod+1,5)=NaN;
    Pperiody(iperiod+1,6)=NaN;
    Pperiody(iperiod+1,7)=NaN;
    Pperiody(iperiod+1,8)=NaN;
    Pperiody(iperiod+1,9)=NaN;
    end
    PycumPeriod((dailyPy(:,1)>=periodstart & dailyPy(:,1)<periodend))=cumsum(abs(Pby_i(:,3)))  % Timeseries of abs alongshore wave power available for transport building up from the start of the survey period
 end

Pperiody(:,4)=Pperiody(:,1)./(Pperiody(:,3)-Pperiody(:,2)); %Py mean
Pperiody(1,:)=NaN;

save('Pperiody.mat','Pperiody')

%%  % figure wave power timeseries
%xticks
monthv=[];
for i_mm=7:12
    iv=datenum([num2str(i_mm) '/01/2009'],'mm/dd/yyyy');
monthv=[ monthv iv];
end
for i_mm=1:12
    iv=datenum([num2str(i_mm) '/01/2010'],'mm/dd/yyyy');
monthv=[ monthv iv];
end
for i_mm=1:12
    iv=datenum([num2str(i_mm) '/01/2011'],'mm/dd/yyyy');
monthv=[ monthv iv];
end
for i_mm=1:12
    iv=datenum([num2str(i_mm) '/01/2012'],'mm/dd/yyyy');
monthv=[ monthv iv];
end

xtickstring=datestr(monthv,'mmm/yy')
xtickstring(1:2:end,:)=NaN;
figure
set(gcf,'position',[185         553        1064         388])
         set(gcf,'PaperPosition',[0  0  60 21])

hold all
plot(dailyPy(:,1),dailyPy(:,2),'k')  % daily P signal

stem( Pperiod(2:end,3), Pperiod(2:end,1)/10,'k')  %Total P
   stem( Pperiod(2:end,3), Pperiod(2:end,4),'b') % mean P
  plot( Pperiod(2:end,6), Pperiod(2:end,5),'om') % max P

plot([Pperiod(2,2) Pperiod(2,2)],[0 max(Pperiod(:,5))*1.1],':r') % Survey date
for iperiod=2:length(pairs)
plot([Pperiod(iperiod,3) Pperiod(iperiod,3)],[0 max(Pperiod(:,5))*1.1],':r') % Survey date
end
for iperiod=2:length(pairs)
    plot ([Pperiod(iperiod,2) Pperiod(iperiod,3)], [Pperiod(iperiod,4) Pperiod(iperiod,4)],'-b') % Mean P bar 
end


plot(dailyPy(:,1),dailyPyN,'.-c')  % daily P signal
plot(dailyPy(:,1),dailyPyS,'.-g')  % daily P signal


set(gca,'Xtick',monthv)

set(gca,'Xticklabel',xtickstring)
legend('Wave power P_d','0.1 * Total wave power between surveys P_{total}','Mean wave power between surveys P_{mean}','Max wave power P_{max} between surveys','Survey dates')
set(legend,'position',[0.6573    0.6838    0.2932    0.3144])
xlim([datenum(['06/01/2009'],'mm/dd/yyyy') datenum(['01/01/2013'],'mm/dd/yyyy')])
xlabel('Date')
ylabel('Wave power')
ylim([0 13E8])


%% Wave Power timeseries images
figure
set(gcf,'position',[  185         0     1293         745])
set(gcf,'PaperPosition',[0  0  40 21])
set(gcf,'color','w')

subplot(311)         
    
plot(dailyPy(:,1),dailyPy(:,2)/10^3,'k')  % daily P signal
hold all
plot(dailyPy(:,1),abs(dailyPy(:,3))/10^3,'Color',[0 0.5 0])  % daily P signal

legend('Wave power P','Alongshore wave power for sed. Transp. | P_{y} |')
set(legend,'position',[ 0.7455    0.8765    0.1780    0.0667])


plot( Pperiod(2:end,6), Pperiod(2:end,5)/10^3,'ok','MarkerFaceColor','k') % max P
plot( Pperiody(2:end,6), Pperiody(2:end,7)/10^3,'d','Color',[0 0.5 0],'MarkerFaceColor',[0 0.5 0]) % max (|Py|)

%%% Plot survey dates
plot([Pperiod(2,2) Pperiod(2,2)],[0 max(Pperiod(:,5))*1.1],':r') % Survey date
for iperiod=2:length(pairs)
plot([Pperiod(iperiod,3) Pperiod(iperiod,3)],[0 max(Pperiod(:,5))*1.1],':r') % Survey date
end

% get xticks ok
set(gca,'Xtick',monthv)
set(gca,'Xticklabel',xtickstring)

xlim([datenum(['06/01/2009'],'mm/dd/yyyy') datenum(['01/01/2013'],'mm/dd/yyyy')])
xlabel('Date')
ylabel('Wave power [kW/m]')
ylim([0 7E4/10^3])

% hatch years
x = [datenum(['01/01/2010'],'mm/dd/yyyy'); datenum(['12/31/2010'],'mm/dd/yyyy');datenum(['12/31/2010'],'mm/dd/yyyy') ; datenum(['01/01/2010'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

x = [datenum(['01/01/2012'],'mm/dd/yyyy'); datenum(['12/31/2012'],'mm/dd/yyyy');datenum(['12/31/2012'],'mm/dd/yyyy') ; datenum(['01/01/2012'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')


subplot(312)         
     
%plot(dailyPy(:,1),dailyPy(:,2),'k')  % daily P signal
stem( Pperiod(2:end,3), Pperiod(2:end,1)/10^3,'k')  %Total P
 hold all
 plot(dailyPy(1:length(PcumPeriod),1),PcumPeriod/10^3,'Color',[0.2 0.2 0.2]) 
 stem( Pperiody(2:end,3),abs( Pperiody(2:end,1))/10^3,'Color',[0 0.5 0],'Marker','d')  %Total P
plot(dailyPy(1:length(PycumPeriod),1),PycumPeriod/10^3,'Color',[0 0.5 0]) 

%%% Plot survey dates
plot([Pperiod(2,2) Pperiod(2,2)],[0 max(PcumPeriod)*1.1],':r') % Survey date
for iperiod=2:length(pairs)
plot([Pperiod(iperiod,3) Pperiod(iperiod,3)],[0 max(PcumPeriod)*1.1],':r') % Survey date
end

% get xticks ok
set(gca,'Xtick',monthv)
set(gca,'Xticklabel',xtickstring)


x = [datenum(['01/01/2010'],'mm/dd/yyyy'); datenum(['12/31/2010'],'mm/dd/yyyy');datenum(['12/31/2010'],'mm/dd/yyyy') ; datenum(['01/01/2010'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

x = [datenum(['01/01/2012'],'mm/dd/yyyy'); datenum(['12/31/2012'],'mm/dd/yyyy');datenum(['12/31/2012'],'mm/dd/yyyy') ; datenum(['01/01/2012'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')


xlim([datenum(['06/01/2009'],'mm/dd/yyyy') datenum(['01/01/2013'],'mm/dd/yyyy')])
xlabel('Date')
ylabel({'Cumulative wave power'; 'between surveys [kW/m]'})
ylim([0 5E5/10^3])

subplot(313)         
% plot( Pperiod(2:end,6), Pperiod(2:end,5)/10^3,'ok','MarkerFaceColor','k') % max P
hold all

% plot( Pperiody(2:end,6), Pperiody(2:end,7)/10^3,'d','Color',[0 0.5 0],'MarkerFaceColor',[0 0.5 0]) % max (|Py|)

x = [datenum(['01/01/2010'],'mm/dd/yyyy'); datenum(['12/31/2010'],'mm/dd/yyyy');datenum(['12/31/2010'],'mm/dd/yyyy') ; datenum(['01/01/2010'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

x = [datenum(['01/01/2012'],'mm/dd/yyyy'); datenum(['12/31/2012'],'mm/dd/yyyy');datenum(['12/31/2012'],'mm/dd/yyyy') ; datenum(['01/01/2012'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

for iperiod=2:length(pairs)
plot ([Pperiod(iperiod,2) Pperiod(iperiod,3)], [Pperiod(iperiod,4) Pperiod(iperiod,4)]/10^3,'-k') % Mean P bar 
plot ([Pperiody(iperiod,2) Pperiody(iperiod,3)], [Pperiody(iperiod,4) Pperiody(iperiod,4)]/10^3,'Color',[0 0.5 0]) % Mean Py bar 


x = [Pperiod(iperiod,2) ; Pperiod(iperiod,3) ;Pperiod(iperiod,3)  ; Pperiod(iperiod,2) ];
y = [Pperiod(iperiod,4);Pperiod(iperiod,4);0; 0]/10^3;
tcolor(1,1,1:3) = 0.7*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

x = [Pperiod(iperiod,2) ; Pperiod(iperiod,3) ;Pperiod(iperiod,3)  ; Pperiod(iperiod,2) ];
y = [Pperiody(iperiod,4);Pperiody(iperiod,4);0; 0]/10^3;
tcolor(1,1,1:3) = 0.7*[0 1 0];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

end

%%% Plot survey dates
plot([Pperiod(2,2) Pperiod(2,2)],[0 max(Pperiod(:,5))*1.1],':r') % Survey date
for iperiod=2:length(pairs)
plot([Pperiod(iperiod,3) Pperiod(iperiod,3)],[0 max(PcumPeriod)*1.1],':r') % Survey date
end

% get xticks ok
set(gca,'Xtick',monthv)
set(gca,'Xticklabel',xtickstring)

xlim([datenum(['06/01/2009'],'mm/dd/yyyy') datenum(['01/01/2013'],'mm/dd/yyyy')])
xlabel('Date')
ylabel({'Mean wave power'; 'between surveys [kW/m]'})
ylim([0 2E4/10^3])

%%% Add a,b,c
annotation('textbox',[0.08 0.95 0.08 0.0407],...
    'String',{'a )'},...
    'LineStyle','none');
% Create textbox
annotation('textbox',[0.08 0.65 0.08 0.0407],...
    'String',{'b )'},...
    'LineStyle','none');
% Create textbox
annotation('textbox',[0.08 0.35 0.08 0.0407],...
    'String',{'c )'},...
    'LineStyle','none');



set (gcf,'renderer','zbuffer')
print WavePowerTimeseries.jpg -djpeg


%% Wave Power timeseries only maxima per survey period image
figure
set(gcf,'position',[  185         196        1293         745])
set(gcf,'PaperPosition',[0  0  40 21])
set(gcf,'color','w')

subplot(311)         
    
plot(Pperiod(:,3),Pperiod(:,4),'o-k','MarkerFaceColor','k')  % daily P signal
hold all
plot(Pperiody(:,3),Pperiody(:,4),'d-','Color',[0 0.5 0],'MarkerFaceColor',[0 0.5 0])  % daily P signal

legend('Wave power P','Alongshore wave power P_y')
set(legend,'position',[ 0.7455    0.8765    0.1780    0.0667])

%%% Plot survey dates
plot([Pperiod(2,2) Pperiod(2,2)],[0 max(Pperiod(:,5))*1.1],':r') % Survey date
for iperiod=2:length(pairs)
plot([Pperiod(iperiod,3) Pperiod(iperiod,3)],[0 max(Pperiod(:,5))*1.1],':r') % Survey date
end

% get xticks ok
set(gca,'Xtick',monthv)
set(gca,'Xticklabel',xtickstring)

xlim([datenum(['06/01/2009'],'mm/dd/yyyy') datenum(['10/04/2012'],'mm/dd/yyyy')])
xlabel('Date')
ylabel('Mean Wave power')
ylim([0 4E8])

x = [datenum(['01/01/2010'],'mm/dd/yyyy'); datenum(['12/31/2010'],'mm/dd/yyyy');datenum(['12/31/2010'],'mm/dd/yyyy') ; datenum(['01/01/2010'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

x = [datenum(['01/01/2012'],'mm/dd/yyyy'); datenum(['12/31/2012'],'mm/dd/yyyy');datenum(['12/31/2012'],'mm/dd/yyyy') ; datenum(['01/01/2012'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')


subplot(312)         
plot(Pperiod(:,3),Pperiod(:,1),'o-k','MarkerFaceColor','k')  % daily P signal
hold all
plot(Pperiody(:,3),Pperiody(:,1),'d-','Color',[0 0.5 0],'MarkerFaceColor',[0 0.5 0])  % daily P signal

%%% Plot survey dates
plot([Pperiod(2,2) Pperiod(2,2)],[0 max(PcumPeriod)*1.1],':r') % Survey date
for iperiod=2:length(pairs)
plot([Pperiod(iperiod,3) Pperiod(iperiod,3)],[0 max(PcumPeriod)*1.1],':r') % Survey date
end

% get xticks ok
set(gca,'Xtick',monthv)
set(gca,'Xticklabel',xtickstring)


x = [datenum(['01/01/2010'],'mm/dd/yyyy'); datenum(['12/31/2010'],'mm/dd/yyyy');datenum(['12/31/2010'],'mm/dd/yyyy') ; datenum(['01/01/2010'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

x = [datenum(['01/01/2012'],'mm/dd/yyyy'); datenum(['12/31/2012'],'mm/dd/yyyy');datenum(['12/31/2012'],'mm/dd/yyyy') ; datenum(['01/01/2012'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')


xlim([datenum(['06/01/2009'],'mm/dd/yyyy') datenum(['10/04/2012'],'mm/dd/yyyy')])
xlabel('Date')
ylabel({'Cumulative wave power'; 'of the previous period'})
ylim([0 10E9])

subplot(313)         
        

plot( Pperiod(2:end,6), Pperiod(2:end,5),'o-k','MarkerFaceColor','k') % max P
hold all
plot( Pperiody(2:end,6), Pperiody(2:end,5),'-d','Color',[0 0.5 0],'MarkerFaceColor',[0 0.5 0]) % max P

plot( Pperiody(2:end,6), Pperiody(2:end,8),'.','Color',[0 0.5 0],'MarkerFaceColor',[0 0.5 0]) % max P
plot( Pperiody(2:end,6), Pperiody(2:end,9),'.','Color',[0 0.5 0],'MarkerFaceColor',[0 0.5 0]) % max P

x = [datenum(['01/01/2010'],'mm/dd/yyyy'); datenum(['12/31/2010'],'mm/dd/yyyy');datenum(['12/31/2010'],'mm/dd/yyyy') ; datenum(['01/01/2010'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

x = [datenum(['01/01/2012'],'mm/dd/yyyy'); datenum(['12/31/2012'],'mm/dd/yyyy');datenum(['12/31/2012'],'mm/dd/yyyy') ; datenum(['01/01/2012'],'mm/dd/yyyy')];
y = [ 10E9; 10E9;0; 0];
tcolor(1,1,1:3) = 0.9*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')



%%% Plot survey dates
plot([Pperiod(2,2) Pperiod(2,2)],[0 max(Pperiod(:,5))*1.1],':r') % Survey date
for iperiod=2:length(pairs)
plot([Pperiod(iperiod,3) Pperiod(iperiod,3)],[0 max(PcumPeriod)*1.1],':r') % Survey date
end

% get xticks ok
set(gca,'Xtick',monthv)
set(gca,'Xticklabel',xtickstring)

xlim([datenum(['06/01/2009'],'mm/dd/yyyy') datenum(['10/04/2012'],'mm/dd/yyyy')])
xlabel('Date')
ylabel('Wave power')
ylim([-14E8 14E8])








set (gcf,'renderer','zbuffer')
% print WavePowerTimeseries.jpg -djpeg

%% Store wave height per survey period

% Total wave power ()

for iperiod=1:length(pairs)-1
    
    periodstart=datenum(dates(pairs(iperiod,2),:),'mm/dd/yyyy');
    periodend=datenum(dates(pairs(iperiod+1,2),:),'mm/dd/yyyy');
    
    %Pb_i=VbP((VbP(:,1)>=periodstart & VbP(:,1)<=periodend),2);
    Hs_i=Vb.Hs((Vb.Hs(:,1)>=periodstart & Vb.Hs(:,1)<periodend),2);
        
    if length(Hs_i) >= 1
    Hsperiod(iperiod+1,1)=nanmean(Hs_i);  % Ptotal
    Hsperiod(iperiod+1,2)=periodstart;         % tstart
    Hsperiod(iperiod+1,3)=periodend;           % tend
    if ~isnan(Hsperiod(iperiod+1,1)) 
    [Hsperiod(iperiod+1,5),ind]=max(Hs_i);   % Pmax
%     Hsperiod(iperiod+1,6)=Hs_i(ind,1);              % t_Pmax
    else
     Hsperiod(iperiod+1,5)  =NaN ;
%      Hsperiod(iperiod+1,6)   =NaN;
    end
    
    else
    Hsperiod(iperiod+1,1)=NaN;
    Hsperiod(iperiod+1,2)=periodstart;
    Hsperiod(iperiod+1,3)=periodend;
    Hsperiod(iperiod+1,5)=NaN;
%     Hsperiod(iperiod+1,6)=NaN;
    end
    
%     HscumPeriod((Vb.Hs(:,1)>=periodstart & Vb.Hs(:,1)<periodend))=cumsum(abs(Hs_i)); % Timeseries of wave power building up from the start of the survey period
    
  end

 Hsperiod(:,4)= Hsperiod(:,1)./( Hsperiod(:,3)- Hsperiod(:,2)); % hsmean 
 Hsperiod (1,:)=NaN;

 figure
  plot(Hsperiod(:,3),Hsperiod(:,1))
  datetick('x')