clear all
close all

% load 
All=load('EURHSall.mat')
TAll=load('EURTall.mat')
Vb=load('EURHS.mat')
VbT=load('EURT.mat')
VbDir=load('EURdir.mat')

theta0=310; % angle of shor normal waves

% cut Vb data at first two years exactly
indVb=find(Vb.Hs(:,1)>=datenum([2009,08,01]) & Vb.Hs(:,1)<=datenum([2011,08,01]));
indVbT=find(VbT.T(:,1)>=datenum([2009,08,01]) & VbT.T(:,1)<=datenum([2011,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

%%

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))
%%


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 2011')

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))
    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
    Hi=Vb.Hs(ii,2);
    Ti=VbT.T(ind,2);
    Diri=VbDir.Dir(ind2,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;
    ipy=ipy+1;
    
    VbPy(ipy,1)=ti;
    VbPy(ipy,2)=Ei*cgi;
    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:7
    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(['10/04/2012'],'mm/dd/yyyy')])
xlabel('Date')
ylabel('Wave power')
ylim([0 13E8])


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

subplot(311)         
    
plot(dailyPy(:,1),dailyPy(:,2),'k')  % daily P signal
hold all
plot(dailyPy(:,1),abs(dailyPy(:,3)),'Color',[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('Wave power')
ylim([0 14E8])

% 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),'k')  %Total P
 hold all
 plot(dailyPy(:,1),PcumPeriod,'Color',[0.2 0.2 0.2]) 
 stem( Pperiody(2:end,3),abs( Pperiody(2:end,1)),'Color',[0 0.5 0],'Marker','d')  %Total P
plot(dailyPy(:,1),PycumPeriod,'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(['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),'ok','MarkerFaceColor','k') % max P
hold all

plot( Pperiody(2:end,6), Pperiody(2:end,7),'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)],'-k') % Mean P bar 
plot ([Pperiody(iperiod,2) Pperiody(iperiod,3)], [Pperiody(iperiod,4) Pperiody(iperiod,4)],'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];
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];
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(['10/04/2012'],'mm/dd/yyyy')])
xlabel('Date')
ylabel('Wave power')
ylim([0 14E8])








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


