%% Relate the observed variability to various hydrodynamical and
%% topographic parameters

clear all 
close all
%%

load AlongshoreVariabilityVB_chi2.mat  % made with the chi2Vbeach.m
load Xo_and_volBeach.mat
load Barcrest.mat
load HighWaveDays.mat

Mpath='C:\Madeschipper\check_out\OpenEarthRawData\tudelft\vlugtenburg\';

load ([Mpath 'scripts\Matthieu\WaveHeightAnalysis\Pperiod.mat'])
load ([Mpath 'scripts\Matthieu\WaveHeightAnalysis\Pperiody.mat'])


load C:\Users\Matthieu\Dropbox\Matthieu\TU\PHD\Projects\AlongshoreVariability\AllVbeachsurveys.mat  % fix this!


%% 

for i_date=1:length(Vsurvey)

    datesmatlab(i_date)=Vsurvey(i_date).matlabdate;
    
end


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

%% 

VerticalAxisParam=1;  % 1= chi^2 , 2= d chi^2 / d t ,
HorizontalAxisParam=1;

for VerticalAxisParam=1%[1:3 ]
for HorizontalAxisParam=2 %[1 2 4 5 17 18]%1:19

if VerticalAxisParam==1
Y_resp=chi_squared
yl_string=('Alongshore topographic variability \chi^2');
end
if VerticalAxisParam==2
Y_resp=[NaN diff(chi_squared)] ;
yl_string=('Change in alongshore variability \Delta\chi^2');
end
if VerticalAxisParam==3
Y_resp=([NaN diff(chi_squared)./diff(datesmatlab)] );
yl_string=('Change in alongshore variability d\chi^2/dt');
end

% wave power
if HorizontalAxisParam==1
X_force=Pperiod(:,4);
xl_string=('Average wave power P_{mean}');
end
if HorizontalAxisParam==2
X_force=Pperiod(:,1);
xl_string=('Cumulative wave power P_{sum}');
end
if HorizontalAxisParam==3
X_force=Pperiod(:,5);
xl_string=('Max wave power P_{max}');
end
% alongshore wave power
if HorizontalAxisParam==4
X_force=Pperiody(:,4);
xl_string=('Average alongshore wave power P_{y,mean}');
end
if HorizontalAxisParam==5
X_force=Pperiody(:,1);
xl_string=('Cumulative alongshore wave power P_{y,sum}'); % check how -y and +y are added
end
if HorizontalAxisParam==6
X_force=Pperiody(:,7);
xl_string=('max alongshore power |Py_{max}|');
end
% ratio Along and total wave power
if HorizontalAxisParam==7
X_force=Pperiody(:,5);
xl_string=('Max Py');
end

% shoreline behavior
if HorizontalAxisParam==8
X_force=[NaN diff(mean(X0'))]';                            
xl_string=(' change in shoreline position');
end

if HorizontalAxisParam==9
X_force=mean(Volmin1toplus10')';                            
xl_string=('beach Vol');
end

if HorizontalAxisParam==10
X_force=[ NaN diff(mean(Volmin1toplus10'))]';                            
xl_string=('d Beach vol');
end

% bar crest
if HorizontalAxisParam==11
X_force=x_c';                               
xl_string=('bar crest distance');
end

if HorizontalAxisParam==12
X_force=(x_c'-mean(X0')');                               %%% remove X0
xl_string=('bar crest distance from shore');
end

if HorizontalAxisParam==13
X_force=[NaN diff(x_c)]';                               
xl_string=('change in bar crest distance');
end

if HorizontalAxisParam==14
X_force=d_c';
xl_string=('bar crest height');
end

if HorizontalAxisParam==15
X_force=n_high';
xl_string=('# high wave days');
end

if HorizontalAxisParam==16
X_force=n_high_oblique';
xl_string=('# high wave days oblique');
end

if HorizontalAxisParam==17
X_force=n_high_oblique'./n_high';
xl_string=('ratio high wave days / high wavedays oblique');
end

if HorizontalAxisParam==18
X_force=perc_high';
xl_string=('% days with wave height H_{s,0} above 2 m');
end

if HorizontalAxisParam==19
X_force=perc_high_oblique';
xl_string=('% days with wave height H_{s,0} above 2 m and wave angle wrt shore-normal \theta_0 larger than 30 deg');
end




cmap=colormap(jet(length(datesmatlab)));

figure(702)
plot([0 max(X_force)],[0 0],'k','linewidth',2)
hold on
for i_date=1:length(datesmatlab)
    scatter(X_force(i_date),Y_resp(i_date),20,'ok')
    scatter(X_force(i_date),Y_resp(i_date),20,cmap(i_date,:),'filled')
end
%plot(forcing,chi_squared,':k')
grid on
cl=colorbar

colorbarticks=2:4:length(datesmatlab)
for ic=1:length(colorbarticks)
colorbarlabels(ic,1:6)=xtickstring(colorbarticks(ic),1:6)
end
set(cl,'YTick',colorbarticks/length(datesmatlab))
set(cl,'YTickLabel',colorbarlabels)


xlabel(xl_string)
ylabel(yl_string)
hold off

%%% Correlation

param1=X_force(~isnan(X_force)&~isnan(Y_resp'));

param2=Y_resp(~isnan(X_force)&~isnan(Y_resp'))';

p = polyfit(param1,param2,1);

yfit = polyval(p,param1);
yresid = param2 - yfit;

SSresid = sum(yresid.^2);
SStotal = (length(param2)-1) * var(param2);
rsq = 1 - SSresid/SStotal
     
[R,P] = corrcoef(param1,param2);
    %R(1,2);
     rsquare=R(1,2)^2
     pvalue=P(1,2)

title(['R^2 = ' num2str(rsquare) ' , p value = ' num2str(pvalue)])     
     
     
     pause (0.1)


     if pvalue <= 0.05 
RsquaresAll(VerticalAxisParam,HorizontalAxisParam)=rsquare ;    
     else
RsquaresAll(VerticalAxisParam,HorizontalAxisParam)=NaN;
     end
% if savefigs==1
%     print ([ 'Chi_vsPymaxVbeach.png'],'-dpng')
% end
end
end




%% figure CD13






figure
set(gcf,'position',[  185         196        1293         505])
set(gcf,'PaperPosition',[0  0  40 21])
set(gcf,'color','w')

Y_resp=[chi_squared] ;
subplot(211)

plot([Pperiod(2,2); Pperiod(2:end,3)],Y_resp,'x-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(Pperiod)
    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({'Alongshore variability'; '\chi^2 [m2]'})
ylim([0 1.5E-2])

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(212)
X_force=Pperiod(:,4);

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

legend('Wave power P_0','Alongshore wave power P_{y0}')

%%% Plot survey dates
plot([Pperiod(2,2) Pperiod(2,2)],[0 4E4],':r') % Survey date
for iperiod=2:length(Pperiod)
    plot([Pperiod(iperiod,3) Pperiod(iperiod,3)],[0 4E4*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({'Mean wave power'; 'in between two surveys'})
ylim([0 4E4])





