% The script file pwoprq03 calculates the Fourier series of a periodic
% soliton and discusses the problem of taking s points less compared
% with the perfect period. The approximation is not sensitive to such
% errors.

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(1:600));  
   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')
   pause
   close
   
   SA=sqrt(Ak.*Ak+Bk.*Bk); 
   plot(Nu,SA,'k o')
   title('amplitudes of components for x=0') 
   pause
   close
   PH=angle(Ak+i*Bk); PH=unwrap(PH);
   plot(Nu, PH,'k +')
   title('phases of components for x=0') 
   pause
   close
   
      RZT=RZ(1:600); WA1=zeros(600:20); WA2=zeros(600:20);
   a0=mean(RZT); RZT=RZT-a0*ones(size(RZT));   
   AT=0:dt:6*Td-dt; AT=AT';
   for r=1:20;
      WA1(:,r)=cos(2*pi/6/Td*r*AT);
   end
   for r=1:20;
      WA2(:,r)=sin(2*pi/6/Td*r*AT);
   end
   WA=[WA1,WA2];
   WAR=WA'*WA;
   IWAR=inv(WAR); WRZT=WA'*RZT;
   FCD0=IWAR*WRZT;
   AFD=FCD0(1:20); BFD=FCD0(21:40); Nu=1:20; Nu=Nu';
   plot(Nu,AFD,'r*',Nu,BFD,'k.')
   title('Fourier coefficients for x=0 by approximation') 
   pause 
   close
   ABFD=abs(AFD+i*BFD); PH=angle(AFD+i*BFD); PHU=unwrap(PH);
   plot(Nu,ABFD,'k o')
   title('Amplitudes of components for x=0 by approximation')
   pause 
   close
   plot(Nu,PHU,'k +')
   title('Phases of components for x=0 by approximation')
   pause 
   close

   [K,E]=ellipke(m);
   PER=2*K/beta/omegal/Td
   PEL=2*K/beta/2/kl/Ll;
   
for s=1:10
   RZC=RZ(1:600-s);   
   a0=mean(RZC); RZC=RZC-a0*ones(size(RZC)); 
   FS0=fft(RZC);  
   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 repeated for s=1:10, 3 periods')
   pause
   close
   
   RZT=RZC; WA1=zeros(600-s:20); WA2=zeros(600-s:20);
   a0=mean(RZT); RZT=RZT-a0*ones(size(RZT));   
   AT=0:dt:6*Td-(s+1)*dt; AT=AT';
   for r=1:20;
      WA1(:,r)=cos(2*pi/6/Td*r*AT);
   end
   for r=1:20;
      WA2(:,r)=sin(2*pi/6/Td*r*AT);
   end
   WA=[WA1,WA2];
   WAR=WA'*WA;
   IWAR=inv(WAR); WRZT=WA'*RZT;
   FCD0=IWAR*WRZT;
   AFD=FCD0(1:20); BFD=FCD0(21:40); Nu=1:20; Nu=Nu';
   plot(Nu,AFD,'r*',Nu,BFD,'k.')
   title('Coefficients by approximation repeated for s=1:10') 
   pause 
   close
end




