% This program plots E(w) and E(t) (it's fourrier transform)
% with : E(w)=sech(1.76*omega)exp(alpha*cos(omega*beta+Phi).
% and : E(t)=fftshift(fft(fftshift(E_w)));
% and : E2(w)=fftshift(ifft(fftshift(E(t)^2)));


function coherent_2_photon;

alpha=3;
beta=3;
Phi=0;

fac= 16;                        % for "hitkansut".
nop = 32*fac;                   % number of points/2.
N=2*nop;                        % number of points.
dw=0.3/fac;                          
w = [-nop:nop-1]*dw;  
omega=w;
dt=2*pi/(N*dw);           
t=[-nop:nop-1]*dt;
w0 = 1;
delta_w = 1;
  
  phi_w = alpha*cos(beta*omega +Phi );  
  A=sech( (1.76*omega)/delta_w);
  E_w=A.*exp(i*phi_w);
  
  E_t=fftshift(fft(fftshift(E_w)));
  E_t=E_t/sqrt(N);
 
  E2_w=fftshift(ifft(fftshift(E_t.*E_t)));
  
  figure;
 
  subplot(2,2,1);
  myplot(w,abs(real(E_w)),phase(E_w)); %for plot in colors as function of the phase
  hold on;
  plot(w,abs(E_w),'color','k'); 
  axis([-5 5 0 1]);
  xlabel('w');
  ylabel('|E(w)|   |Re(E(w))|');
  
  
  subplot(2,2,2);
  myplot(t,abs(real(E_t)),phase(E_t));
  xlabel('t');
  ylabel('Re(E(t))');
  axis([-60 60 -2 3]);
  hold on;
  
  
  subplot(2,2,4);
  myplot(t,real(E_t.^2),phase(E_t.^2));
  xlabel('t');
  ylabel('Re(E(t)^2)');
  axis([-60 60 -2 6]);
  hold on;
 
  
  subplot(2,2,3);
  myplot(w,abs(real(E2_w)),phase(E2_w));
  hold on;
  plot(w,abs(E2_w),'color','k');
  xlabel('w');
  ylabel('|E_{2}(w)|   |Re(E_{2}(w))|'); 
  axis([-5 5 0 0.06]);
  hold on;
 

return;