% plot_inwave_test_chan
%
% compute sima2 overf h2 for dx/L
%
% jcw 9/7/2021
%

cd E:\data\models\InWave\advect_test\output_testchan

%2 5 10 15 20 25

%here are teh dx grids for advection of wave from left
logo='c2';
logo='c4';
logo='u3';   %use this one
logo='a4'; 
logo='hs';
Hwave2=squeeze(ncread(['ocean_inwave_test_chan2v4_',logo,'_his.nc'],'Hwave',[1 5 1],[Inf 1 Inf]));
Hwave5=squeeze(ncread(['ocean_inwave_test_chan5v4_',logo,'_his.nc'],'Hwave',[1 5 1],[Inf 1 Inf]));
Hwave10=squeeze(ncread(['ocean_inwave_test_chan10v4_',logo,'_his.nc'],'Hwave',[1 5 1],[Inf 1 Inf]));
Hwave15=squeeze(ncread(['ocean_inwave_test_chan15v4_',logo,'_his.nc'],'Hwave',[1 5 1],[Inf 1 Inf]));
Hwave20=squeeze(ncread(['ocean_inwave_test_chan20v4_',logo,'_his.nc'],'Hwave',[1 5 1],[Inf 1 Inf]));
Hwave25=squeeze(ncread(['ocean_inwave_test_chan25v4_',logo,'_his.nc'],'Hwave',[1 5 1],[Inf 1 Inf]));
Hwave30=squeeze(ncread(['ocean_inwave_test_chan30v4_',logo,'_his.nc'],'Hwave',[1 5 1],[Inf 1 Inf]));

ot=ncread('ocean_inwave_test_chan5v4_c2_his.nc','ocean_time');
x_rho2=ncread('ocean_inwave_test_chan2v4_c2_his.nc','x_rho',[1 5],[Inf 1]);
x_rho5=ncread('ocean_inwave_test_chan5v4_c2_his.nc','x_rho',[1 5],[Inf 1]);
x_rho10=ncread('ocean_inwave_test_chan10v4_c2_his.nc','x_rho',[1 5],[Inf 1]);
x_rho15=ncread('ocean_inwave_test_chan15v4_c2_his.nc','x_rho',[1 5],[Inf 1]);
x_rho20=ncread('ocean_inwave_test_chan20v4_c2_his.nc','x_rho',[1 5],[Inf 1]);
x_rho25=ncread('ocean_inwave_test_chan25v4_c2_his.nc','x_rho',[1 5],[Inf 1]);
x_rho30=ncread('ocean_inwave_test_chan30v4_c2_his.nc','x_rho',[1 5],[Inf 1]);
%
if (0)
 figure
 tidx=2000;
 hold on
 grid on
 plot(x_rho2,squeeze(Hwave2(:,tidx)),'k')
 plot(x_rho5,squeeze(Hwave5(:,tidx)),'c')
 plot(x_rho10,squeeze(Hwave10(:,tidx)),'r')
 plot(x_rho15,squeeze(Hwave15(:,tidx)),'g')
 plot(x_rho20,squeeze(Hwave20(:,tidx)),'b')
 plot(x_rho25,squeeze(Hwave25(:,tidx)),'m')
 plot(x_rho30,squeeze(Hwave30(:,tidx)),'m')
 set(gca,'ylim',[0 3.5])
 set(gca,'xlim',[0 1000])  %= 6 wave forms repeated
end
%do x= 60 and 632  ??5590 = 5530+60
ie=630;
is2=find(x_rho2==60); ie2=find(x_rho2==ie);
is5=find(x_rho5==60); ie5=find(x_rho5==ie);
is10=find(x_rho10==60); ie10=find(x_rho10==ie);
is15=find(x_rho15==60); ie15=find(x_rho15==ie);
is20=find(x_rho20==60); ie20=find(x_rho20==ie+10);
is25=find(x_rho25==50); ie25=find(x_rho25==ie-5);
is30=find(x_rho30==60); ie30=find(x_rho30==ie);

%tidxs=2001;  tidxe=3024;
%figure
%plot(Hwave2(is2,tidxs:tidxe))
%hold on
%plot(Hwave2(ie2,tidxs:tidxe))

Lwave=squeeze(ncread('ocean_inwave_test_chan2v4long_his.nc','Lwave',[1 5  1],[Inf 1 Inf]));
cx2=squeeze(ncread('ocean_inwave_test_chan2v4long_his.nc','cx',[1 5 2 1],[Inf 1 1 Inf]));
cx=7.68;
Lwave=81.728;
cgr below = 7.68
so  7*Lwave/cx=74.51 sec delay
at a spot 7*81.728 m downthe channel =572.1 m
do x=60 and x = 632.25  ?
%isnce it havad +60 s in the cos term, it repeats every 120 sec.
%120*7.68m/s=921.6 m repeat cycle * 5 = 5530m
572.1m/7.68m/s = 74.5 sec

toff=74.5;   %720
figure
subplot(321)
hold off
plot(ot,Hwave2(is2,:))
hold on
plot(ot-toff,Hwave2(ie2,:))

subplot(322)
hold off
plot(ot,Hwave5(is5,:))
hold on
plot(ot-toff,Hwave5(ie5,:))

subplot(323)
hold off
plot(ot,Hwave10(is10,:))
hold on
plot(ot-toff,Hwave10(ie10,:))

subplot(324)
hold off
plot(ot,Hwave15(is15,:))
hold on
plot(ot-toff,Hwave15(ie15,:))

subplot(325)
hold off
plot(ot,Hwave20(is20,:))
hold on
plot(ot-toff,Hwave20(ie20,:))

subplot(326)
hold off
plot(ot,Hwave30(is30,:))
hold on
plot(ot-toff,Hwave30(ie30,:))



%%%%%now do calcs
toff=75;   %720
mean_2s=mean((Hwave2(is2,2576-toff:3600-toff)).^2);
mean_2e=mean((Hwave2(ie2,2576:3600)).^2);
std_2s=std((Hwave2(is2,2576-toff:3600-toff)).^2);
std_2e=std((Hwave2(ie2,2576:3600)).^2);
%
mean_5s=mean((Hwave5(is5,2576-toff:3600-toff)).^2);
mean_5e=mean((Hwave5(ie5,2576:3600)).^2);
std_5s=std((Hwave5(is5,2576-toff:3600-toff)).^2);
std_5e=std((Hwave5(ie5,2576:3600)).^2);
%
mean_10s=mean((Hwave10(is10,2576-toff:3600-toff)).^2);
mean_10e=mean((Hwave10(ie10,2576:3600)).^2);
std_10s=std((Hwave10(is10,2576-toff:3600-toff)).^2);
std_10e=std((Hwave10(ie10,2576:3600)).^2);
%
mean_15s=mean((Hwave15(is15,2576-toff:3600-toff)).^2);
mean_15e=mean((Hwave15(ie15,2576:3600)).^2);
std_15s=std((Hwave15(is15,2576-toff:3600-toff)).^2);
std_15e=std((Hwave15(ie15,2576:3600)).^2);
%
mean_20s=mean((Hwave20(is20,2576-toff:3600-toff)).^2);
mean_20e=mean((Hwave20(ie20,2576:3600)).^2);
std_20s=std((Hwave20(is20,2576-toff:3600-toff)).^2);
std_20e=std((Hwave20(ie20,2576:3600)).^2);
%
mean_25s=mean((Hwave25(is25,2576-toff:3600-toff)).^2);
mean_25e=mean((Hwave25(ie25,2576:3600)).^2);
std_25s=std((Hwave25(is25,2576-toff:3600-toff)).^2);
std_25e=std((Hwave25(ie25,2576:3600)).^2);
%
mean_30s=mean((Hwave30(is30,2576-toff:3600-toff)).^2);
mean_30e=mean((Hwave30(ie30,2576:3600)).^2);
std_30s=std((Hwave30(is30,2576-toff:3600-toff)).^2);
std_30e=std((Hwave30(ie30,2576:3600)).^2);

Ym=[mean_2e/mean_2s mean_5e/mean_5s mean_10e/mean_10s mean_15e/mean_15s mean_20e/mean_20s mean_25e/mean_25s mean_30e/mean_30s ];
YGF=[(std_2e/mean_2e)/(std_2s/mean_2s) (std_5e/mean_5e)/(std_5s/mean_5s) (std_10e/mean_10e)/(std_10s/mean_10s) ...
     (std_15e/mean_15e)/(std_15s/mean_15s) (std_20e/mean_20e)/(std_20s/mean_20s) (std_25e/mean_25e)/(std_25s/mean_25s) ...
      (std_30e/mean_30e)/(std_30s/mean_30s) ];
Ym_c2=Ym;
YGF_c2=YGF;
Ym_c4=Ym;
YGF_c4=YGF;
Ym_u3=Ym;
YGF_u3=YGF;
Ym_a4=Ym;
YGF_a4=YGF;
Ym_hs=Ym;
YGF_hs=YGF;
%  save inwave_test_chan_adv.mat      <-----------  i saved this work above, jcw 5.9.2022
%
figure
hold on
X=[2 5 10 15 20 25 30]/(81.75);
%plot(X,Ym_c2,'-+b')
plot(X,Ym_c4,'-+r','linewidth',2)
plot(X,Ym_a4,'-+k','linewidth',2)
plot(X,Ym_u3,'-+b','linewidth',2)
%plot(X,Ym_hs,'-+c')
hold on
%plot(X,YGF_c2,'-.b')
plot(X,YGF_c4,'-.r','linewidth',2)
plot(X,YGF_a4,'-.k','linewidth',2)
plot(X,YGF_u3,'-.b','linewidth',2)
%plot(X,YGF_hs,'-.c')
grid on
%legend('relative Hsig C2',...
%        'relative GF C2', ...
legend('relative Hsig C4',...
        'relative Hsig A4',...
        'relative Hsig U3',...
        'relative GF C4', ...
        'relative GF A4', ...
        'relative GF U3')
%        'relative GF HS')
%        'relative Hsig HS',...
grid on
set(gca,'xlim',[0 .4],'ylim',[0.5 1.1])
xlabel('dx/Lwave [-]')
ylabel('ratio [-]')
title([' Numerical error as function of grid size '])  %, ...
%      [' for the 3rd order upstream bias advection'])
set(gcf,'pos',[574         377        1285         302])
%print -dpng figx_numerror_vs_dx_test.png
%print -depsc figx_numerror_vs_dx_test.eps


figure
plot(ot(2576:3600),Hwave2(is5,2576:3600),'b--','linewidth',2)
hold on
toff=79% 74.2;%  75;   %720
plot(ot(2576:3600)-30,Hwave2(ie5,2576:3600),'k')
plot(ot(2576:3600)-toff,Hwave5(ie5,2576:3600),'r')
plot(ot(2576:3600)-toff,Hwave10(ie10,2576:3600),'g')
plot(ot(2576:3600)-toff,Hwave15(ie15,2576:3600),'c')
plot(ot(2576:3600)-toff,Hwave20(ie20,2576:3600),'m')
plot(ot(2576:3600)-toff,Hwave25(ie25,2576:3600),'b')
plot(ot(2576:3600)-toff,Hwave30(ie30,2576:3600),'k--')
grid on
legend('Initial shape','dx =  2m','dx =  5m','dx = 10m','dx = 15m','dx = 20m','dx = 25m','dx = 30m')
set(gca,'xlim',[2600 3080],'ylim',[0 3.5])
xlabel('time, s')
ylabel('Hmo, m')
title([' Time series of wave heigth after propagation 7 wave lengths'])
set(gcf,'pos',[574         377        1285         302])
%print -dpng figx_hwave_vs_time_test.png
%print -depsc figx_hwave_vs_time_test.eps


figure  %did not use this one
plot(x_rho2,Hwave2(:,3600),'b--','linewidth',2)
hold on
plot(x_rho2, Hwave2(:,3600),'k')
plot(x_rho5, Hwave5(:,3600),'r')
plot(x_rho10,Hwave10(:,3600),'g')
plot(x_rho15,Hwave15(:,3600),'c')
plot(x_rho20,Hwave20(:,3600),'m')
plot(x_rho25,Hwave25(:,3600),'b')
plot(x_rho30,Hwave30(:,3600),'k--')
grid on
legend('Initial shape','dx =  2m','dx =  5m','dx = 10m','dx = 15m','dx = 20m','dx = 25m','dx = 30m')
set(gca,'xlim',[0 1000],'ylim',[0 3.5])
xlabel('distance, m')
ylabel('Hsig, m')
title([' Wave shape after propagation 7 wave lengths'])












