% The script file pwoprp03 calculates the Fourier series of a periodic
% soliton as a function of time for different shifts in space x and
% compares with solutions of approximations by trigonometric functions.
% The solutions when the orthogonal properties of the set of functions
% are not used are better.

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;
   
   x1=Ll/6; X1=2*kl*x1*ones(size(X))+X;
   Y1=kl*x1*ones(size(X))+Y;
   [sn,cn,dn1]=ellipj(beta*X1,m);   
   Z1=100*sqrt(2)/kl*beta*dn1.*exp(i*Y1);
   AZ1=abs(Z1); RZ1=real(Z1); IZ1=imag(Z1);
   plot(t,RZ1,'k',t,AZ1,'r',t,-AZ1,'r')
   title( 'Surface elevation as function of t for x=Ll/6')   
   pause
   close
   a01=mean(RZ1(1:601)); RZ1=RZ1-a01*ones(size(RZ1));
   FS1=fft(RZ1(1:601));  
   Ak1=2*real(FS1(2:21))/600; Bk1=-2*imag(FS1(2:21))/600; 
   Nu=1:20;Nu=Nu';
   plot(Nu,Ak1,'r*',Nu,Bk1,'k.')
   title('Coefficients of Fourier series x=Ll/6')
   pause
   close
   SA1=sqrt(Ak1.*Ak1+Bk1.*Bk1);
   plot(Nu,SA1,'k o')
   title('Amplitudes of components for x=Ll/6') 
   pause
   close
   PH1=angle(Ak1+i*Bk1); PHU1=unwrap(PH1);
   plot(Nu,PHU1,'k +')
   title('Phases of components x=Ll/6') 
   pause
   close
   
   RZ1T=RZ1(1:600);
   a01=mean(RZ1T); RZ1T=RZ1T-a01*ones(size(RZ1T));   
   WRZ1T=WA'*RZ1T;
   FCD1=IWAR*WRZ1T;
   AFD1=FCD1(1:20); BFD1=FCD1(21:40); Nu=1:20; Nu=Nu';
   plot(Nu,AFD1,'r*',Nu,BFD1,'k.')
   title('Coefficients of Fourier series x=Ll/6 by approximation')   
   pause 
   close
   ABFD1=abs(AFD1+i*BFD1); PH1=angle(AFD1+i*BFD1); PHU1=unwrap(PH1);
   plot(Nu,ABFD1,'k o')
   title('Amplitudes of components x=Ll/6 by approximation')
   pause 
   close
   plot(Nu,PHU1,'k +')
   title('Phases of components x=Ll/6 by approximation')
   pause 
   close

   
   x2=2*Ll/6; X2=2*kl*x2*ones(size(X))+X;
   Y2=kl*x2*ones(size(X))+Y;
   [sn,cn,dn2]=ellipj(beta*X2,m);   
   Z2=100*sqrt(2)/kl*beta*dn2.*exp(i*Y2);
   AZ2=abs(Z2); RZ2=real(Z2); 
   plot(t,[RZ2,AZ2,-AZ2])
   title( 'surface elevation as function of t for x=2*Ll/6')      
   pause
   close
   a02=mean(RZ2(1:601)); RZ2=RZ2-a02*ones(size(RZ2));
   FS2=fft(RZ2(1:601));  
   Ak2=2*real(FS2(2:21))/6/N; Bk2=-2*imag(FS2(2:21))/6/N; 
   Nu=1:20;Nu=Nu';
   plot(Nu,Ak2,'r*',Nu,Bk2,'k.')
   title('Coefficients of Fourier series x=2*Ll/2')
   pause
   close
   SA2=sqrt(Ak2.*Ak2+Bk2.*Bk2);
   plot(Nu,SA2,'k o')
   title('amplitudes of components for x=2*Ll/6') 
   pause
   close
   PH2=angle(Ak2+i*Bk2); PHU2=unwrap(PH2);
   plot(Nu,PHU2,'k +')
   title('phases of components x=2*Ll/6') 
   pause
   close
   
   RZ2T=RZ2(1:600);
   a02=mean(RZ2T); RZ2T=RZ2T-a02*ones(size(RZ2T));   
   WRZ2T=WA'*RZ2T;
   FCD2=IWAR*WRZ2T;
   AFD2=FCD2(1:20); BFD2=FCD2(21:40); Nu=1:20; Nu=Nu';
   plot(Nu,AFD2,'r*',Nu,BFD2,'k.')
   title('Coefficients of Fourier series x=2*Ll/6 by approximation')   
   pause 
   close
   ABFD2=abs(AFD2+i*BFD2); PH2=angle(AFD2+i*BFD2); PHU2=unwrap(PH2);
   plot(Nu,ABFD2,'k o')
   title('Amplitudes of components x=2*Ll/6 by approximation')
   pause 
   close
   plot(Nu,PHU2,'k +')
   title('Phases of components x=2*Ll/6 by approximation')
   pause 
   close

   

   
   x3=3*Ll/6; X3=2*kl*x3*ones(size(X))+X;
   Y3=kl*x3*ones(size(X))+Y;
   [sn,cn,dn3]=ellipj(beta*X3,m);   
   Z3=100*sqrt(2)/kl*beta*dn3.*exp(i*Y3);
   AZ3=abs(Z3); RZ3=real(Z3);
   plot(t,[RZ3,AZ3,-AZ3])
   title( 'surface elevation as function of t for x=3*Ll/6')      
   pause
   close
   a03=mean(RZ3(1:601)); RZ3=RZ3-a03*ones(size(RZ3));
   FS3=fft(RZ3(1:601));  
   Ak3=2*real(FS3(2:21))/6/N; Bk3=-2*imag(FS3(2:21))/6/N; 
   Nu=1:20;Nu=Nu';
   plot(Nu,Ak3,'r*',Nu,Bk3,'k.')
   title('Coefficients of Fourier series x=3*Ll/2')
   pause
   close
   SA3=sqrt(Ak3.*Ak3+Bk3.*Bk3);
   plot(Nu,SA3,'k o')
   title('amplitudes of components for x=3*Ll/6') 
   pause
   close
   PH3=angle(Ak3+i*Bk3); PHU3=unwrap(PH3);
   plot(Nu,PHU3,'k +')
   title('phases of components x=3*Ll/6') 
   pause
   close
   
   RZ3T=RZ3(1:600);
   a03=mean(RZ3T); RZ3T=RZ3T-a03*ones(size(RZ3T));   
   WRZ3T=WA'*RZ3T;
   FCD3=IWAR*WRZ3T;
   AFD3=FCD3(1:20); BFD3=FCD3(21:40); Nu=1:20; Nu=Nu';
   plot(Nu,AFD3,'r*',Nu,BFD3,'k.')
   title('Coefficients of Fourier series x=3*Ll/6 by approximation')   
   pause 
   close
   ABFD3=abs(AFD3+i*BFD3); PH3=angle(AFD3+i*BFD3); PHU3=unwrap(PH3);
   plot(Nu,ABFD3,'k o')
   title('Amplitudes of components x=3*Ll/6 by approximation')
   pause 
   close
   plot(Nu,PHU3,'k +')
   title('Phases of components x=3*Ll/6 by approximation')
   pause 
   close




   
   x4=4*Ll/6; X4=2*kl*x4*ones(size(X))+X;
   Y4=kl*x4*ones(size(X))+Y;
   [sn,cn,dn4]=ellipj(beta*X4,m);   
   Z4=100*sqrt(2)/kl*beta*dn4.*exp(i*Y4);
   AZ4=abs(Z4); RZ4=real(Z4);
   plot(t,[RZ4,AZ4,-AZ4]) 
   title( 'surface elevation as function of t for x=4*Ll/6')      
   pause
   close
   a04=mean(RZ4(1:601)); RZ4=RZ4-a04*ones(size(RZ4));
   FS4=fft(RZ4(1:601));  
   Ak4=2*real(FS4(2:21))/6/N; Bk4=-2*imag(FS4(2:21))/6/N; 
   Nu=1:20;Nu=Nu';
   plot(Nu,Ak4,'r*',Nu,Bk4,'k.')
   title('Coefficients of Fourier series x=4*Ll/6')
   pause
   close
   SA4=sqrt(Ak4.*Ak4+Bk4.*Bk4);
   plot(Nu,SA4,'k o')
   title('amplitudes of components for x=4*Ll/6') 
   pause
   close
   PH4=angle(Ak4+i*Bk4); PHU4=unwrap(PH4);
   plot(Nu,PHU4,'k +')
   title('phases of components x=4*Ll/6') 
   pause
   close
   
   RZ4T=RZ4(1:600);
   a04=mean(RZ4T); RZ4T=RZ4T-a04*ones(size(RZ4T));   
   WRZ4T=WA'*RZ4T;
   FCD4=IWAR*WRZ4T;
   AFD4=FCD4(1:20); BFD4=FCD4(21:40); Nu=1:20; Nu=Nu';
   plot(Nu,AFD4,'r*',Nu,BFD4,'k.')
   title('Coefficients of Fourier series x=4*Ll/6 by approximation')   
   pause 
   close
   ABFD4=abs(AFD4+i*BFD4); PH4=angle(AFD4+i*BFD4); PHU4=unwrap(PH4);
   plot(Nu,ABFD4,'k o')
   title('Amplitudes of components x=4*Ll/6 by approximation')
   pause 
   close
   plot(Nu,PHU4,'k +')
   title('Phases of components x=4*Ll/6 by approximation')
   pause 
   close
   
   x5=5*Ll/6; X5=2*kl*x5*ones(size(X))+X;
   Y5=kl*x5*ones(size(X))+Y;
   [sn,cn,dn5]=ellipj(beta*X5,m);   
   Z5=100*sqrt(2)/kl*beta*dn5.*exp(i*Y5);
   AZ5=abs(Z5); RZ5=real(Z5);
   plot(t,[RZ5,AZ5,-AZ5])
   title( 'surface elevation as function of t for x=5*Ll/6')      
   pause
   close
   a05=mean(RZ5(1:601)); RZ5=RZ5-a05*ones(size(RZ5));
   FS5=fft(RZ5(1:601));  
   Ak5=2*real(FS5(2:21))/6/N; Bk5=-2*imag(FS5(2:21))/6/N; 
   Nu=1:20;Nu=Nu';
   plot(Nu,Ak5,'r*',Nu,Bk5,'k.')
   title('Coefficients of Fourier series x=5*Ll/2')
   pause
   close
   SA5=sqrt(Ak5.*Ak5+Bk5.*Bk5);
   plot(Nu,SA5,'k o')
   title('amplitudes of components for x=5*Ll/6') 
   pause
   close
   PH5=angle(Ak5+i*Bk5); PHU5=unwrap(PH5);
   plot(Nu,PHU5,'k +')
   title('phases of components x=5*Ll/6') 
   pause
   close
   
   RZ5T=RZ5(1:600);
   a05=mean(RZ5T); RZ5T=RZ5T-a05*ones(size(RZ5T));   
   WRZ5T=WA'*RZ5T;
   FCD5=IWAR*WRZ5T;
   AFD5=FCD5(1:20); BFD5=FCD5(21:40); Nu=1:20; Nu=Nu';
   plot(Nu,AFD5,'r*',Nu,BFD5,'k.')
   title('Coefficients of Fourier series x=5*Ll/6 by approximation')   
   pause 
   close
   ABFD5=abs(AFD5+i*BFD5); PH5=angle(AFD5+i*BFD5); PHU5=unwrap(PH5);
   plot(Nu,ABFD5,'k o')
   title('Amplitudes of components x=5*Ll/6 by approximation')
   pause 
   close
   plot(Nu,PHU5,'k +')
   title('Phases of components x=5*Ll/6 by approximation')
   pause 
   close
      
   x6=6*Ll/6; X6=2*kl*x6*ones(size(X))+X;
   Y6=kl*x6*ones(size(X))+Y;
   [sn,cn,dn6]=ellipj(beta*X6,m);   
   Z6=100*sqrt(2)/kl*beta*dn6.*exp(i*Y6);
   AZ6=abs(Z6); RZ6=real(Z6);
   plot(t,[RZ6,AZ6,-AZ6])
   title( 'surface elevation as function of t for x=6*Ll/6')      
   pause
   close
   a06=mean(RZ6(1:601)); RZ6=RZ6-a06*ones(size(RZ6));
   FS6=fft(RZ6(1:601));  
   Ak6=2*real(FS6(2:21))/6/N; Bk6=-2*imag(FS6(2:21))/6/N; 
   Nu=1:20;Nu=Nu';
   plot(Nu,Ak6,'r*',Nu,Bk6,'k.')
   title('Coefficients of Fourier series x=6*Ll/2')
   pause
   close
   SA6=sqrt(Ak6.*Ak6+Bk6.*Bk6);
   plot(Nu,SA6,'k o')
   title('amplitudes of components for x=6*Ll/6') 
   pause
   close
   PH6=angle(Ak6+i*Bk6); PHU6=unwrap(PH6);
   plot(Nu,PHU6,'k +')
   title('phases of components x=6*Ll/6') 
   pause
   close
      
   RZ6T=RZ6(1:600);
   a06=mean(RZ6T); RZ6T=RZ5T-a06*ones(size(RZ6T));   
   WRZ6T=WA'*RZ6T;
   FCD6=IWAR*WRZ6T;
   AFD6=FCD6(1:20); BFD6=FCD6(21:40); Nu=1:20; Nu=Nu';
   plot(Nu,AFD6,'r*',Nu,BFD6,'k.')
   title('Coefficients of Fourier series x=6*Ll/6 by approximation')   
   pause 
   close
   ABFD6=abs(AFD6+i*BFD6); PH6=angle(AFD6+i*BFD6); PHU6=unwrap(PH6);
   plot(Nu,ABFD6,'k o')
   title('Amplitudes of components x=6*Ll/6 by approximation')
   pause 
   close
   plot(Nu,PHU6,'k +')
   title('Phases of components x=6*Ll/6 by approximation')
   pause 
   close
   
   
   x7=PER*omegal*Td/kl/4; X7=2*kl*x7*ones(size(X))+X;
   Y7=kl*x7*ones(size(X))+Y;
   [sn,cn,dn1]=ellipj(beta*X7,m);   
   Z7=100*sqrt(2)/kl*beta*dn1.*exp(i*Y7);
   AZ7=abs(Z7); RZ7=real(Z7); IZ7=imag(Z7);
   plot(t,RZ7,'k',t,AZ7,'r',t,-AZ7,'r')
   title( 'Surface elevation as function of t for half period')   
   pause
   close
   
   x8=PER*omegal*Td/kl/2; X8=2*kl*x8*ones(size(X))+X;
   Y8=kl*x8*ones(size(X))+Y;
   [sn,cn,dn1]=ellipj(beta*X8,m);   
   Z8=100*sqrt(2)/kl*beta*dn1.*exp(i*Y8);
   AZ8=abs(Z8); RZ8=real(Z8); IZ8=imag(Z8);
   plot(t,RZ8,'k',t,AZ8,'r',t,-AZ8,'r')
   title( 'Surface elevation as function of t for full period')   
   pause
   close
   
   
   

   
   
   


   
   
   
   


  
  

