% Here is my matlab source code ... 
% this scripts computes the wave elevation and bottom pressure for a
% shoaling wave train, and estimates the k=0 spectral density 
% which comes into the primary microseisms signal. 

%clear; close all;
addpath /export/home/ardhuin/TOOLS/MATLAB
nx=2586;
dx1=65.757904031020715;

nz2=2268;
X2=linspace(0,nz2-1,nz2)*dx1./10;

nT=20
fall=linspace(0.005,0.1,nT);
   
wnum(1)=dispNewtonTH(fall(17),100);
   
Z2=100+20.*cos(wnum(1)*1.05.*X2);
for it=17:17 %1:nT
% Parameters for wind wave
   a1=1;surfig=0;gam=0.5;gi=9.81;f=fall(it);imu=sqrt(-1);h1=4800;h2=0.001;om=2*pi*f;
 %  a1=1;surfig=0;gam=0.5;gi=9.81;f=fall(it);imu=sqrt(-1);h1=2400;h2=0.001;om=2*pi*f;
%   a1=1;surfig=1;gam=0.5;gi=9.81;f=fall(it);imu=sqrt(-1);h1=3000;h2=0.001;om=2*pi*f;

   nn=64;
   Nx=21000*nn;
   x=linspace(X2(1),X2(end),Nx);
   ZZ=interp1(X2,Z2,x)';
   dx=x(2)-x(1);
   L=dx*Nx;

   a=zeros(Nx,1);
   iz2=find(ZZ < h2);
   iz1=find(ZZ > h1);
   h3=100;
   ZZ(iz1)=h1;
   if(length(iz2) > 0) 
   ZZ(iz2(1):end)=h2;
   end

   %ZZs = smoothn(ZZ,1E16);
   %I=find(x >2E5 & x < 3.15E5);
   %ZZ(I)=ZZs(I);

   ZZ(iz2)=h2;
   Hp=0.01;

   % solves dispersion relation for estimating Airy wave solutions
   wnum=zeros(Nx,1);
   kx=zeros(Nx,1);
   wnum(1)=dispNewtonTH(f,ZZ(1));a(1)=a1;
   C=om/wnum(1);kh=wnum(1)*ZZ(1);pp=0.5*(1+(2*kh)/sinh(2*kh));Cginf=pp*C;
   kinf=wnum(1);
   for i=2:Nx
      wnum(i)=dispNewtonTH(f,ZZ(i));
      if i==2
          kx(i)=kx(i-1)+wnum(i)*dx;
      else
          kx(i)=kx(i-1)+0.5*(wnum(i)+wnum(i-1))*dx;
      end
      if (mod(kx(i),80) < mod(kx(i-1),80)) 
          %[i kx(i)]
          %kx(i) = kx(i)+0.2*(rand-0.5);
      end
      %kh=wnum(i)*ZZ(i);
      %pp=0.5*(1+(2*kh)/sinh(2*kh));
      %Cg=pp*om/wnum(i);
      %%[ZZ(i) Cg pp]
      %dkm=0.05*(om/Cg)*cos(kx(i)/(3*pi^2));
      %kx(i)=kx(i)+dkm*dx;
      %%[wnum(i) dkm cos(kx(i)/pi^2)]
%%% THE AMPLITUDE IS CALCULATED from   
%%% conservation of action: d/dx(Cg*E/sigma)=0--> Cg*A^2=const over the x-axis
%%% 
      C=om/wnum(i);kh=wnum(i)*ZZ(i);pp=0.5*(1+(2*kh)/sinh(2*kh));Cg=pp*C;
%
% limitation by breaking for wind waves ... 
      a(i)=a1.*min(sqrt(Cginf/Cg),ZZ(i)*gam);
   end
   Ldeep=2*pi/wnum(1);

%% This is for IG waves: constant amplitude in surf zone
   indsurf=find (a > a1.*ZZ*gam*0.99);
   if (surfig == 1 & length(indsurf) > 0)
      a(indsurf)=a(indsurf(1));
   end

inddry=find(ZZ < 0.0015);
a(inddry)=0.;

%kx=kx-kx(iz2); %-om/(sqrt(gi*Hp)).*2.*sqrt(abs(h2/Hp));
z=a.*(cos(kx));
zquad=a.*(sin(kx));

nphase=120;
Kphase=zeros(nphase,1);
Kqphase=zeros(nphase,1);

pb=a.*(cos(kx))./cosh(wnum.*ZZ);
pbquad=a.*(sin(kx))./cosh(wnum.*ZZ);

%%% Creates an animation of the surface elevation and bottom pressure
nfft=Nx;

figure(10)
fig=figure('Position',[1 1 600 200])
set(fig,'Color',[ 1 1 1]);
nt=36;
prefix='bottom2';
for ii=1:nt
    ii
    clf
    phi=(ii-1)*2*pi/nt;
    zphi=z.*cos(phi)+zquad.*sin(phi);
    pphi=pb.*cos(phi)+pbquad.*sin(phi);
   hold on
    plot(x./1000,-ZZ./60-2.,'k-','LineWidth',1);
    plot(x./1000,zphi./4,'b-','LineWidth',2);
    plot(x./1000,pphi.*4-1.5,'r-','LineWidth',2);
   hold off

 xlabel('x (km)')
 ht=text(0.5,0.8,'sea surface elevation / 4 (m)','Fontsize',16,'Color',[0 0 1])

ht=text(0.5,-2.9,'bottom pressure \times 4 - 1.5 (m)','Fontsize',16,'Color',[1 0 0])
ht=text(10,-2.9,'depth/60 -2 (m)','Fontsize',16,'Color',[0 0 0])
   number=sprintf('%.4d',ii);
   oname=[prefix '_' number '.png' ];
   set(gcf, 'PaperPositionMode', 'auto')
   set(gca,'Xlim',[0 15],'Ylim',[-4 0.9])
   set(gca,'position',[0.05 0.2 0.92 0.74])
   saveas(gcf, oname, 'png')

end