% 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 Polarization Maintaining Fiber (PMF), showing its retardation and
% (slow) axis of birefringence onto the Poincaré sphere

clear all;
% ***************** simulation parameters ***************** 
numFreq= 128;	% n. of optical frequencies (omega vector "om")
Lfiber= 100e3;        % [m] fiber length
banda=10e9;           % [Hz] two-sided simulation bandwidth
% semegen=2;         % generator seed for random sequences: rand() function
% % semegen=input('generator seed? ');
% rand('seed',semegen);		% generator in Matlab 4.0; 

% DGD is randomly fixed between two extreme values
dgd_low= 10;    dgd_high= 50;   %  [ps]
dgd= (dgd_low+rand*(dgd_high-dgd_low))*1e-12;          % [s]  Differential Group Delay (DGD) value
Dbeta1= dgd/Lfiber;	% [s/m] recall DGD: Dtau= Dbeta1*Lfiber , where Dbeta1= d Dbeta(om)/d om
 				
% 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);

  % azimuth and ellipticity
  the=rand*pi - pi/2;
  eps= rand*pi/2 - pi/4;
  % BEWARE: a uniform ellipticity does not imply a uniform statistical distribution on the Poincaré sphere
  % A uniform distribution is obtained when 2*epsilon is distributed according to 'the cosine law' => 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
  % Stokes parameters of local eigenstates
  e1= cos(2*the)*cos(2*eps);
  e2= sin(2*the)*cos(2*eps);
  e3= sin(2*eps);
  % Birefringence strength: for a PMF, the rotation angle is a linear function of frequency
  Dbetaz= Dbeta1*om*Lfiber;	% [rad]	
  % 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
  eDotSpin= e1*sig1 + e2*sig2 + e3*sig3;
  for k=1:numFreq,
    % exponential form: exp(-i/2*Dbeta*z [e.sigma])
    U(:,:,k)= cos(Dbetaz(k)/2)*eye(2)-complex(0,1)*sin(Dbetaz(k)/2)*eDotSpin;
  end	% k (opt. freq. index)

% check
for k=1:numFreq,
  if (abs(det(U(:,:,k)) - 1) > 1e-9),
    sprintf('The obtained Jones matrix in not unit-determinant')
  end;
end	% k (opt. freq. index)
  
% ***************** Analysis of Jones matrix U *******************

% get rotation angle and axis (Stokes eigenvector) 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
  traceU=trace(U(:,:,k));	     % the trace of U is 2cos(delphi/2)
  coshalf= real(traceU)/2;	     %  -1<coshalf<1
  delphi(k)= 2*acos(coshalf);	     %  0<delphi/2<pi
  sin_bDotSpin= complex(0,1)*(U(:,:,k) - (1/2)*traceU*eye(2));	% sen(delphi/2)*[b.sigma]
  sinSquare= -1*det(sin_bDotSpin);         % sen^2(delphi/2)
  if (abs(sinSquare) > 1e-9)
    % unit mag. vector b derived from Hermitian matrix [b.sigma]
    bDotSpin= sin_bDotSpin / sqrt(sinSquare);
    b1(k)=  real(bDotSpin(1,1));   % or  b1(k)= -real(bDotSpin(2,2));
    b2(k)=  real(bDotSpin(1,2));   % or  b2(k)=  real(bDotSpin(2,1));
    b3(k)= -imag(bDotSpin(1,2));   % or  b3(k)=  imag(bDotSpin(2,1));
  else
    % this is the case DeltaPhi=0 and the rotation axis is undefined  (here
    % arbitrarily fixed at the value of the previous frequency
    delphi(k)=0;
    b1(k)=  b1(k-1); b2(k)=  b2(k-1); b3(k)=  b3(k-1); 
  end;
  
  % check and normalize rot. axis
  norm_b= sqrt(b1(k)^2+b2(k)^2+b3(k)^2);
  if (abs(norm_b-1) > 1e-9) | (abs(b1(k))>1) | (abs(b2(k))>1) | (abs(b3(k))>1),
    sprintf('The rotation axis b is not a unit magnitude Stokes vector')
  end;
  b1(k)= b1(k)/norm_b;
  b2(k)= b2(k)/norm_b;
  b3(k)= b3(k)/norm_b;
  
%   % 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');
xlabel('s1'); ylabel('s2');zlabel('s3');
view(110,20);
lightangle(90,45);
colormap summer;
hold on;

% plot axis of birefringence
for i=1:numFreq,
  [azb,elb,mob]=cart2sph(b1(i),b2(i),b3(i));
  if abs(mob - 1) > 1e-9 
    sprintf('b a modulo <> 1 alla frequenza numero %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);