% The script file pwoprr03 calculates the Fourier series of a periodic
% soliton and discusses the problem when additionally two and three
% periods are considered. Such a change does not change a solution for
% a perfect period.

g=9.81;Ll=1.20; kl=2*pi/Ll; omegal=sqrt(g*kl); Tw=2*pi/omegal;  
m=0.995445; amax=0.06;

   ar=amax*sqrt(1-m*m/2);
   gamma=kl*ar/sqrt(2);
   omegad=omegal*(1+gamma*gamma);Tl=2*pi/omegal;
   Td=2*pi/omegad;
   beta=gamma*sqrt(2/(2-m*m));
   
   N=100; dt=Td/N; t=0:dt:3*6*Td; t=t'; X=-omegal*t;
   [sn,cn,dn]=ellipj(beta*X,m);
   Y=-omegad*t;
   Z=100*sqrt(2)/kl*beta*dn.*exp(i*Y);
   AZ=abs(Z); RZ=real(Z); IZ=imag(Z); 
   plot(t,RZ,'k',t,AZ,'r',t,-AZ,'r')
   title( 'surface elevation as function of t for x=0')  
   pause   
   close
   
   a0=mean(RZ); RZ=RZ-a0*ones(size(RZ)); 
   FS0=fft(RZ);  
   Ak=2*real(FS0(2:61))/3/6/N; Bk=-2*imag(FS0(2:61))/3/6/N;
   Nu=1:60; Nu=Nu';
   plot(Nu,Ak,'r*',Nu,Bk,'k.')
   title('Coefficients of Fourier series x=0, 3 periods')
   pause
   close
   
   RZA=RZ(1:1200);
   a0=mean(RZA); RZA=RZA-a0*ones(size(RZA)); 
   FS0=fft(RZA);  
   Ak=2*real(FS0(2:41))/2/6/N; Bk=-2*imag(FS0(2:41))/2/6/N;
   Nu=1:40; Nu=Nu';
   plot(Nu,Ak,'r*',Nu,Bk,'k.')
   title('Coefficients of Fourier series x=0, 2 periods')
   pause
   close
   
   RZB=RZ(1:600);
   a0=mean(RZB); RZB=RZB-a0*ones(size(RZB)); 
   FS0=fft(RZB);  
   Ak=2*real(FS0(2:21))/6/N; Bk=-2*imag(FS0(2:21))/6/N;
   Nu=1:20; Nu=Nu';
   plot(Nu,Ak,'r*',Nu,Bk,'k.')
   title('Coefficients of Fourier series x=0, 1 period')
   pause
   close
   
for s=1:10
   RZC=RZ(1:1800-s);   
   a0=mean(RZC); RZ=RZC-a0*ones(size(RZC)); 
   FS0=fft(RZC);  
   Ak=2*real(FS0(2:61))/3/6/N; Bk=-2*imag(FS0(2:61))/3/6/N;
   Nu=1:60; Nu=Nu';
   plot(Nu,Ak,'r*',Nu,Bk,'k.')
   title('Coefficients of Fourier series repeated for s=1:10, 3 periods')
   pause
   close
end
