% Università di Parma - Dipartimento di Ingegneria e Architettura
% Master Degree in Communications Engineering
% Course: Polarized Fiber Optic Transmission
% Prof. Armando Vannucci, 2018
%
% This script numerically simulates the (unit-determinant) Jones matrix of
% a randomly birefringent optical fiber, thorugh the Randow Waveplate Model.
% The Gisin-Pelleaux recursive rule [Opt.Commun.'92] is used to determine
% the local DGD of each waveplate, that is a standard PMF.
% If desired, the local DGD can be randomized (with Normal distribution) by
% setting the flag "randlocaldgd"

clear all;
% ***************** simulation parameters ***************** 
numFreq= 128;	% n. of optical frequencies (omega vector "om")
numplates= 100;		% n. of PMF waveplates to concatenate
AvgDGD= 10.0e-12;		% [s] expected (root-mean-square, rms) DGD value
linlocaleig=0;		% [0/1]: 1 the eigenstates of each waveplates must be linear polarizations
                    %        0 they can be elliptic polarizations as well
erroreDbeta0=0;		% [0/1]: 1 does not account for waveplates birefringence, at carrier freq. (omega=0)
                    %        0 includes a uniform [0,2*pi] birefringence, at omega=0
banda=10e9;           % [Hz] two-sided simulation bandwidth
semegen=2;         % generator seed for random sequences: rand() function
% semegen=input('generator seed? ');
if isempty(semegen)
else
 rand('seed',semegen);		% generator in Matlab 4.0; 
 				% use rand('state',semegen) in matlab 5.0
 randn('seed',semegen);		% generator in Matlab 4.0; 
end

% INITIALIZATION
% vector of optical frequencies
pulsfin= 2*pi*(banda/2);	% [rad/s]
deltaom= 2*pulsfin/numFreq;
% fft logic: do not start from -pulsfin (same value as final sample +pulsfin: periodic fft) 
% find frequency 0 after an even number of freq. samples
om=-pulsfin+deltaom:deltaom:pulsfin;    % [rad/s] opt. freq. vector ("numFreq" elements)
nomc= ceil(numFreq/2);               % index for reference freq. (om=0)


% ***************** generation of Jones matrix U *******************

% Unit-determinant Jones matrix: opt. freq. "om" in 3rd dimension
U=zeros(2,2,numFreq);
for i=1:numFreq, U(:,:,i)=eye(2); end;

% local mean dgd (Deltabeta1*z): concetto Gisin, <Dtau^2>=N <dtau^2>
localAvgDGD=AvgDGD/sqrt(numplates); % [s]

% cycle on waveplates: numbered input to output
for l=1:numplates,
  % azimuth and ellipticity
  the=rand*pi - pi/2;
  % uniform distribution on Poincaré sphere=> f_{2eps}(e)=(1/2)cos(e) ;  |e|<pi/2
  eps= (1/2)*asin(rand*2-1);
  % ellipticity set to 0 if linearly polarized local eigenstates are chosen
  if linlocaleig==1
    eps=0;
  end
  % Stokes parameters of local eigenstates
  e1= cos(2*the)*cos(2*eps);
  e2= sin(2*the)*cos(2*eps);
  e3= sin(2*eps);
  % local birefringence (Stokes rotation angle) at om=0 (Deltabeta0*z)
  % uniform in [-pi; pi]. Physically, equal to Deltabeta1*z*om0 (om0= carrier... THz)
  Dbeta0z=rand*2*pi - pi;
  % neglect borefringence at om=0, if flag "erroreDbeta0" is set to 1
  if erroreDbeta0==1
    Dbeta0z=0;
  end
  Dbeta1z=localAvgDGD;	% [s]	 Dbeta1z= ( d Dbeta(om)/d om )*z
  % For a standard PMF, the rotation angle due to birefringence is a linear function of frequency
  Dbetaz= Dbeta0z + Dbeta1z*om;	% [rad]
  % Jones matrix of n-th waveplate
  Un=zeros(2,2,numFreq);
  % Pauli matrices
  sig1=[1 0 ; 0 -1];
  sig2=[0 1 ; 1  0];
  sig3=complex(zeros(2),[0 -1 ; 1 0]);
  % Hermitian matrix: scalar product e.sigma
  eHerm= e1*sig1 + e2*sig2 + e3*sig3;
  for k=1:numFreq,
    % exponential form: exp(-i/2*Dbeta*z [e.sigma])
    Un(:,:,k)= cos(Dbetaz(k)/2)*eye(2)-complex(0,1)*sin(Dbetaz(k)/2)*eHerm;
    % update global Jones matrix
    U(:,:,k)= Un(:,:,k)*U(:,:,k); 
  end	% k (opt. freq. index)
end	% l (wveplate index)

% check
for k=1:numFreq,
  if (abs(det(U(:,:,k)) - 1) > 1e-9),
    sprintf('Beware: the Jones matrix is not unitary')
  end;
end	% k (opt. freq. index)
  
% ***************** Analysis of Jones matrix U *******************

% get rotation angle and axis (eigenmodes) from unitary matrix
delphi=zeros(1,numFreq);	% rotation angle (in Stokes)
b1=zeros(1,numFreq);		% rotation axis
b2=zeros(1,numFreq);
b3=zeros(1,numFreq);
for k=1:numFreq,
    
  % extract the parameters of interest
  tracciaU=trace(U(:,:,k));	             % the trace of U is 2cos(delphi/2)
  coshalf= real(tracciaU)/2;	     %  -1<coshalf<1
  delphi(k)= 2*acos(coshalf);	     %  0<delphi/2<pi
  senperbHerm= complex(0,1)*(U(:,:,k) - (1/2)*tracciaU*eye(2));	% sen(delphi/2)*[b.sigma]
  senqua= -1*det(senperbHerm);         % sen^2(delphi/2)
  senalfamez= sqrt(real(senqua));      % 0<delphi/2<pi
  if (abs(senqua) > 1e-9)
    % unit mag. vector b derived from Hermitian matrix [b.sigma]
    bHerm= senperbHerm / sqrt(senqua);
  else
    % we are here in the case alfa=0; fix vector a=[1 0 0] arbitrarily
    delphi=0;
    bHerm=[1 0; 0 -1];
    senqua=0;
  end;
  % use average, in case bHerm does not have the proper symmetries
  b1(k)= real(bHerm(1,1)-bHerm(2,2))/2;
  b2(k)= real(bHerm(1,2)+bHerm(2,1))/2;
  b3(k)= imag(bHerm(2,1)-bHerm(1,2))/2;
  
  % check and normalize rot. axis
  normab= sqrt(b1(k)^2+b2(k)^2+b3(k)^2);
  if (abs(normab-1) > 1e-9) | (abs(b1(k))>1) | (abs(b2(k))>1) | (abs(b3(k))>1),
    sprintf('Beware: The rotation axis b is not a unit magnitude Stokes vector')
  end;
  b1(k)= b1(k)/normab;
  b2(k)= b2(k)/normab;
  b3(k)= b3(k)/normab;
  
  % Each time delphi crosses zero, we invert the sign of previous (in freq.) angles
  % and axes, so as to get continuous plots. For the same reason, the
  % undetermined (DeltaPhi(k)==0) eigenstate b(k) is set to its previous value b(k-1)
  if (delphi(k)==0)
    % if entered, sen(DeltaPhi/2) is too close to 0 (and set to 0 by if statement above ): 
    b1(k)= b1(k-1);
    b2(k)= b2(k-1);
    b3(k)= b3(k-1);
    sprintf('sign inversion (angle and eigenstates) up to frequency # %d',k)
    for l=1:k,
      delphi(l)= -delphi(l);
      b1(l)= -b1(l);
      b2(l)= -b2(l);
      b3(l)= -b3(l);
    end;		% l (opt. frequencies)
  end;		% if (delphi(k)==0)

end;		% k (opt. frequencies)

% ***************** Plot of the Parameters of Jones matrix U *******************

% plot global angle and eigenstates of polarization
figure(1);
clf;
plot(om/(2*pi*1e9),delphi*(180/pi),'r.');
title('Rotation angle \Delta\phi (in Stokes space)');
xlabel('f [GHz]'); ylabel('\Delta\phi [deg.]');

% prepare sphere with illumination on angles the=0, eps=45 deg.
figure(2);
clf;
sphere;
axis equal;
axis ([-1 1 -1 1 -1 1]);
title(['Rotation axis b   seed=',num2str(semegen)]);
xlabel('s1'); ylabel('s2');zlabel('s3');
view(110,20);
lightangle(90,45);
colormap summer;
hold on;

% plot eigenstates
for i=1:numFreq,
  [azb,elb,mob]=cart2sph(b1(i),b2(i),b3(i));
  if abs(mob - 1) > 1e-9 
    sprintf('Beware: magnitude of b <> 1 at frequency # %4d',i) 
  end
  plot3(b1(i),b2(i),b3(i),'g.');
end	% i (opt. freq.)
% viewpoint=central SOP (+90 deg. due to odd difference between view and cart2sph)
[azb,elb,mob]=cart2sph(b1(nomc),b2(nomc),b3(nomc));
view(azb*180/pi + 90,elb*180/pi);


% ***************** Analysis of PMD vector Omega *******************

% evaluate Hermitian traceless matrix -i/2*[Omega(om).sigma] (so called "N-matrix")
Nmat= zeros(2,2,numFreq);
for k=1:numFreq-1, 
  % discrete approximation of derivative
  Uprime= (U(:,:,k+1)-U(:,:,k))/deltaom;
  Nmat(:,:,k)= Uprime*inv(U(:,:,k));
end;
Nmat(:,:,numFreq)= Nmat(:,:,numFreq-1);   % copy last value, for continuity

% etraction of PMD vector elements
Omega1=zeros(1,numFreq);
Omega2=zeros(1,numFreq);
Omega3=zeros(1,numFreq);
for k=1:numFreq,
    Omdotsigma= 2*complex(0,1)*Nmat(:,:,k);
    % average components... in case "Omdotsigma" is not numerically Hermitian
    Omega1(k)= real(Omdotsigma(1,1)-Omdotsigma(2,2))/2;
    Omega2(k)= real(Omdotsigma(1,2)+Omdotsigma(2,1))/2;
    Omega3(k)= imag(Omdotsigma(2,1)-Omdotsigma(1,2))/2;
end;		% k (opt. frequencies)

% calculate magnitude (DGD) and unit-magnitude vector (PSP)
deltau=zeros(1,numFreq);	% DGD Delta_tau(om)
q1=zeros(1,numFreq);		% (slow) PSP
q2=zeros(1,numFreq);
q3=zeros(1,numFreq);
for i=1:numFreq,
  % here magnitude in [rad/(2 pi Hz) == s]
  moduloOmega= sqrt(Omega1(i)^2 + Omega2(i)^2 + Omega3(i)^2);
  deltau(i)= moduloOmega;	% [s]
  q1(i)= Omega1(i)/moduloOmega;
  q2(i)= Omega2(i)/moduloOmega;
  q3(i)= Omega3(i)/moduloOmega;
end;	% i (pulsazioni)


% ***************** Plot of the Parameters of PMD vector *******************

% plot DGD
figure(3);
clf;
plot(om/(2*pi*1e9),deltau*1e12,'k.');
title(['DGD calculated from vector \Omega   \surd < \Delta\tau^2 >=',num2str(AvgDGD),' [ps]']);
xlabel('f [GHz]');
ylabel('\Delta\tau [ps] (black)');
axxis= axis; axis([axxis(1) axxis(2) 0 2*AvgDGD*1e12]);

% prepare sphere with light from the=0, eps=45 deg.
figure(4);
clf;
sphere;
axis equal;
axis ([-1 1 -1 1 -1 1]);
title(['PSP unit-mag. vector q (blue)   seed=',num2str(semegen)]);
xlabel('s1'); ylabel('s2');zlabel('s3');
view(110,20);
lightangle(90,45);
colormap summer;
hold on;

% plot PSP
for i=1:numFreq,
  [azq,elq,moq]=cart2sph(q1(i),q2(i),q3(i));
  if abs(moq - 1) > 1e-9 
    sprintf('q with magnitude <> 1 at freq. # %4d',i) 
  end
  plot3(q1(i),q2(i),q3(i),'b.');
end
% viewpoint=central SOP (+90 deg. due to odd difference between view and cart2sph)
[azq,elq,moq]=cart2sph(q1(nomc),q2(nomc),q3(nomc));
view(azq*180/pi + 90,elq*180/pi);