% This Filename:  stoke.mp   [MetaPost source]
% Creation time:  Thu Jan 23 20:19:34 2014
%
% Copyright (C) 1997-2005, Fredrik Jonsson <fj@optics.kth.se>
%
% Input Filename [Stokes parameters]:  
% This MetaPost source code was automatically generated by poincare
% Full set of command line options that generated this code:
%     --verbose --normalize --axislengths 0.3 1.6 0.3
%     2.7 0.3 1.5 --axislabels S_1 bot
%     S_2 bot S_3 rt --rotatephi 15.0
%     --rotatepsi -70.0 --shading 0.75 0.99 --rhodivisor
%     50 --phidivisor 80 --scalefactor 25.0 --outputfile
%     stoke.mp --auxsource polstate.mp
%
% Description:  Map of Stokes parameters, visualized as trajectories
%               onto the Poincare sphere. This file contains MetaPost
%               source code, to be compiled with John Hobby's MetaPost
%               compiler or used with anything that understands MetaPost
%               source code.
%
% If you want to create PostScript output, or include the resulting
% output in a TeX document, this example illustrates the procedure,
% assuming 'poincaremap.mp' to be the name of the file containing the
% MetaPost code to be visualized: (commands run on command-line)
%
%       mp poincaremap.mp;
%       echo "\input epsf\centerline{\epsfbox{poincaremap.1}}\bye" > tmp.tex;
%       tex tmp.tex;
%       dvips tmp.dvi -o poincaremap.ps;
%
% Here, the first command compiles the MetaPost source code, and leaves
% an Encapsulated PostScript file named 'poincaremap.1', containing TeX
% control codes for characters, etc. This file does not contain any
% definitions for characters or TeX-specific items, and it cannot be
% viewed or printed simply as is stands; it must rather be included into
% TeX code in order to provide something useful.
%     The second command creates a temporary minimal TeX-file 'tmp.tex',
% that only includes the previously generated Encapsulated PostScript
% code.
%     The third command compiles the TeX-code into device-independent,
% or DVI, output, stored in the file 'tmp.dvi'.
%     Finally, the last command converts the DVI output into a free-
% standing PostScript file 'poincaremap.ps', to be printed or viewed
% with some PostScript viewer, such as GhostView.
%
scalefactor := 25.000000 mm;
rot_psi := -70.000000;  % Rotation angle round z-axis (first rotation)
rot_phi := 15.000000;  % Rotation angle round y-axis (second rotation)
alpha := -35.416613;    % == arctan(sin(rot_phi)*tan(rot_psi))
beta  := -5.381520;    % == arctan(sin(rot_phi)/tan(rot_psi))

%
% Parameters specifying the location of the light source; for Phong
% shading of the sphere.
%
%    phi_source:  Angle (in deg.) to light source counterclockwise
%                 'from three o'clock', viewed from the observer.
%
%  theta_source:  Angle (in deg.) between light source and observer,
%                 seen from the centre of the sphere.
%
% Parameters specifying the shading 'intensity' in terms of maximum
% (for the highlighs) and minimum (for the deep shadowed regions)
% values for the Phong shading.  '0.0' <=> 'black'; '1.0' <=> 'white'
%
%   upper_value:  Maximum value of whiteness.
%   lower_value:  Minimum value of whiteness.
%
phi_source := 30.000000;
theta_source := 30.000000;
upper_value := 0.990000;
lower_value := 0.750000;
radius := scalefactor;
delta_rho := radius/50.000000;
delta_phi := 360.0/80.000000;
beginfig(1);
  path p;
  path equator;
  transform T;
  c1:=lower_value;
  c2:=upper_value-lower_value;
  nx_source := sind(theta_source)*cosd(phi_source);
  ny_source := sind(theta_source)*sind(phi_source);
  nz_source := cosd(theta_source);
  phistop := 360.0;
  rhostop := radius - delta_rho/2.0;
%
% Draw the shaded Poincare sphere projected on 2D screen coordinates
%
  for rho=0.0cm step delta_rho until rhostop:
    for phi=0.0 step delta_phi until phistop:
      rhomid := rho + delta_rho/2.0;
      phimid := phi + delta_phi/2.0;
      x1 := rho*cosd(phi);
      y1 := rho*sind(phi);
      x2 := (rho+delta_rho)*cosd(phi);
      y2 := (rho+delta_rho)*sind(phi);
      x3 := (rho+delta_rho)*cosd(phi+delta_phi);
      y3 := (rho+delta_rho)*sind(phi+delta_phi);
      x4 := rho*cosd(phi+delta_phi);
      y4 := rho*sind(phi+delta_phi);
      p:=makepath makepen ((x1,y1)--(x2,y2)--(x3,y3)--(x4,y4)--(x1,y1));
      quot := (rhomid/radius);
      nx_object := quot*cosd(phimid);
      ny_object := quot*sind(phimid);
      nz_object := sqrt(1-quot*quot);
      prod:=nx_object*nx_source+ny_object*ny_source
            +nz_object*nz_source;
      if prod < 0.0:
         value := c1;
      else:
         value := c1 + c2*prod*prod;
      fi
      fill p withcolor value[black,white];
    endfor
  endfor

%
% Draw the 'equators' of the Poincare sphere
%
   equator := halfcircle scaled (2.0*radius);
   eqcolval := .45;    % '0.0' <=> 'white';  '1.0' <=> 'black'

   pickup pencircle scaled 0.600000 pt;
%
% Draw equator $S_3=0$...
%
   T := identity yscaled sind(rot_phi) rotated 180.0;
   draw equator transformed T withcolor eqcolval [white,black];

%
% ... then equator $S_2=0$...
%
   T := identity yscaled (cosd(rot_phi)*sind(rot_psi))
                 rotated (270.0 + alpha);
   draw equator transformed T withcolor eqcolval [white,black];

%
% ... and finally equator $S_1=0$.
%
   T := identity yscaled (cosd(rot_phi)*cosd(rot_psi))
                 rotated (270.0 - beta);
   draw equator transformed T withcolor eqcolval [white,black];

%
% Draw the $S_1$-, $S_2$- and $S_3$-axis of the Poincare sphere.
% First of all, calculate the transformations of the intersections
% for the unity sphere.
%
% Used variables:
%
%    behind_distance : Specifies the relative distance of the coordi-
%                      axes to be plotted behind origo (in negative di-
%                      rection of respective axis.
%
%   outside_distance_s1 : The relative distance from origo to the point
%                         of the arrow head of the coordinate axis S1.
%                         If this is set to 1.0, the arrow head will
%                         point directly at the Poincare sphere.
%
%   outside_distance_s2 : Same as above, except that this one controls
%                         the S2 coordinate axis instead.
%
%   outside_distance_s3 : Same as above, except that this one controls
%                         the S3 coordinate axis instead.
%
%    insidecolval :    Specifies the shade of gray to use for the parts
%                      of the coordinate axes that are inside the Poin-
%                      care sphere. Values must be between 0 and 1,
%                      where:  '0.0' <=> 'white';  '1.0' <=> 'black'
%
   behind_distance_s1  := -0.300000;
   behind_distance_s2  := -0.300000;
   behind_distance_s3  := -0.300000;
   outside_distance_s1 :=  1.600000;
   outside_distance_s2 :=  2.700000;
   outside_distance_s3 :=  1.500000;
   insidecolval := .85;    % '0.0' <=> 'white';  '1.0' <=> 'black'

   pickup pencircle scaled 0.600000 pt;
%
% Start with drawing the x-axis...
%
   x_bis_start :=  radius*behind_distance_s1*cosd(rot_psi)*cosd(rot_phi);
   y_bis_start :=  radius*behind_distance_s1*sind(rot_psi);
   z_bis_start := -radius*behind_distance_s1*cosd(rot_psi)*sind(rot_phi);
   x_bis_intersect :=  radius*cosd(rot_psi)*cosd(rot_phi);
   y_bis_intersect :=  radius*sind(rot_psi);
   z_bis_intersect := -radius*cosd(rot_psi)*sind(rot_phi);
   p := makepath makepen (y_bis_intersect,z_bis_intersect)--
             (outside_distance_s1*y_bis_intersect,
              outside_distance_s1*z_bis_intersect);
   drawarrow p;
   label.bot(btex $S_1$ etex,
             (outside_distance_s1*y_bis_intersect,
              outside_distance_s1*z_bis_intersect));

%
% ... then draw the y-axis ...
%
   x_bis_start := -radius*behind_distance_s2*sind(rot_psi)*cosd(rot_phi);
   y_bis_start :=  radius*behind_distance_s2*cosd(rot_psi);
   z_bis_start :=  radius*behind_distance_s2*sind(rot_psi)*sind(rot_phi);
   x_bis_intersect := -radius*sind(rot_psi)*cosd(rot_phi);
   y_bis_intersect :=  radius*cosd(rot_psi);
   z_bis_intersect :=  radius*sind(rot_psi)*sind(rot_phi);
   p := makepath makepen (y_bis_intersect,z_bis_intersect)--
             (outside_distance_s2*y_bis_intersect,
              outside_distance_s2*z_bis_intersect);
   drawarrow p;
   label.bot(btex $S_2$ etex,
             (outside_distance_s2*y_bis_intersect,
              outside_distance_s2*z_bis_intersect));

%
% ... then, finally, draw the z-axis.
%
   x_bis_start := radius*behind_distance_s3*sind(rot_phi);
   y_bis_start := 0.0;
   z_bis_start := radius*behind_distance_s3*cosd(rot_phi);
   x_bis_intersect := radius*sind(rot_phi);
   y_bis_intersect := 0.0;
   z_bis_intersect := radius*cosd(rot_phi);
   p := makepath makepen (y_bis_intersect,z_bis_intersect)--
             (outside_distance_s3*y_bis_intersect,
              outside_distance_s3*z_bis_intersect);
   drawarrow p;
   label.rt(btex $S_3$ etex,
             (outside_distance_s3*y_bis_intersect,
              outside_distance_s3*z_bis_intersect));

%
% The following external file is included (using the --auxsource option): polstate.mp  [MetaPost source]
%
   input polstate.mp
   endfig;
end
