Search:

Return to previous page

Contents of file 'poincare/poincare.c':




Deprecated: Function split() is deprecated in /storage/content/45/2011745/jonsson.eu/public_html/php/htmlicise.php on line 46
    1   /*-----------------------------------------------------------------------------
    2   | File: poincare.c [ANSI-C conforming source code]                            |
    3   | Created:       November 17, 1997, Fredrik Jonsson <fj@optics.kth.se>        |
    4   | Last modified: January 24, 2014, Fredrik Jonsson <http://jonsson.eu>        |
    5   | Description:                                                                |
    6   |       This program creates maps of Stokes parameters, visualized as         |
    7   |       trajectories on the Poincare sphere. The program generates MetaPost   |
    8   |       source code as output, to be compiled into PostScript or Encapsulated |
    9   |       PostScript, using John Hobby's MetaPost compiler, or to be used with  |
   10   |       anything that understands MetaPost source code.                       |
   11   |                                                                             |
   12   | Copyright (C) 1997-2014 Fredrik Jonsson <fj@optics.kth.se>                  |
   13   |                                                                             |
   14   | Compile with:                                                               |
   15   |                                                                             |
   16   |        gcc -O2 -g -Wall -pedantic -ansi ./poincare.c -o ./poincare -lm      |
   17   |                                                                             |
   18   | or use the enclosed Makefile, designed for use with GCC:                    |
   19   |      ----------------- CUT HERE ------------------                          |
   20   |      CC     = gcc                                                           |
   21   |      CCOPTS = -O2 -Wall -pedantic -ansi                                     |
   22   |      LNOPTS = -lm                                                           |
   23   |                                                                             |
   24   |      all:                                                                   |
   25   |              make poincare                                                  |
   26   |                                                                             |
   27   |      poincare: ./poincare.o                                                 |
   28   |              $(CC) $(CCOPTS) -o ./poincare ./poincare.o $(LNOPTS)           |
   29   |                                                                             |
   30   |      poincare.o: ./poincare.c                                               |
   31   |              $(CC) $(CCOPTS) -c ./poincare.c                                |
   32   |      ----------------- CUT HERE ------------------                          |
   33   |                                                                             |
   34   |                              Code by:  Fredrik Jonsson <fj@optics.kth.se>   |
   35   |                                        Department of Physics II (Optics)    |
   36   |                                        The Royal Institute of Technology    |
   37   |                                        S-100 44  Stockholm,  Sweden         |
   38   |                                                                             |
   39   | Revision history:                                                           |
   40   |                                                                             |
   41   |  971117:  After a busy weekend I have now finished a program 'poincare.c'   |
   42   |  [v.1.0]  that generates the MetaPost-code for visualizing the output       |
   43   |           Stokes-parameters -- from traject.c -- as trajectory maps         |
   44   |           onto a Phong shaded Poincare sphere.                              |
   45   |                                                                             |
   46   |  971216:  Visiting Ecole Polytechnique for a week.  Added the '-n' option   |
   47   |  [v.1.1]  to 'poincare', which makes it capable of plotting normalized      |
   48   |           Stokes-parameters (S1/S0,S2/S0,S3/S0) without preparing initial   |
   49   |           data.  Also added the feature that 'poincare' writes comments     |
   50   |           into the resulting output MetaPost source file, making it a       |
   51   |           bit easier for the masochist who enjoy reading MetaPost source    |
   52   |           [but -- of course -- Real Men read PostScript directly instead].  |
   53   |                                                                             |
   54   |  980216:  Implemented support for drawing arrows on the poincare sphere     |
   55   |  [v.1.2]  directly by specifications in terms of command-line options.      |
   56   |           In this first version of the arrow drawing command, only          |
   57   |           solid-line- and dashed-line-arrows are supported, though with     |
   58   |           arbitrary shade of gray from white to black.                      |
   59   |                                                                             |
   60   |  980608:  Changed the Stokes-parameter-path-drawing parts, so that hidden   |
   61   |  [v.1.3]  parts of the paths are properly taken into account, and drawn     |
   62   |           with dashed lines (instead of the regular in-front solid lines.   |
   63   |           I also changed the format specifications of the input files, so   |
   64   |           that an arbitrary number of separated paths may be present.       |
   65   |                                                                             |
   66   |  980610:  Added the '--bezier' option, to give the possibility of changing  |
   67   |  [v.1.4]  between regular piecewise straight-line type paths and Bezier     |
   68   |           interpolated paths with continuous first-derivatives.             |
   69   |           Added the '--paththickness' option.                               |
   70   |                                                                             |
   71   |           Added the '--draw_hidden_dashed' option, which gives the user     |
   72   |           a choice between drawing hidden parts of the paths (specified     |
   73   |           with the '--input' option) with dashed or solid lines (of diffe-  |
   74   |           rent shade than the visible, 'in-front' parts.  This required     |
   75   |           some restructuring of the program, since I up until now have      |
   76   |           used only black, dashed lines for hidden parts of paths, in which |
   77   |           case the order in which the paths are written doesn't matter.     |
   78   |           Now the path-drawing statements instead are called upon twice,    |
   79   |           the first time for drawing hidden parts, and the second time for  |
   80   |           drawing visible, 'in-front' parts.  This ensures that no hidden   |
   81   |           part will be written over a visible part, giving a proper         |
   82   |           appearance of the mapped trajectory of the specified path.        |
   83   |                                                                             |
   84   |           Also added the '--axislength' option which gives the possibility  |
   85   |           tuning the axis lengths, depending on the viewpoint for the       |
   86   |           Poincare sphere.                                                  |
   87   |                                                                             |
   88   |  981113:  (Ecole Polytechnique) Added the possibility to visualize an       |
   89   |  [v.1.5]  additional coordinate system, rotated to the original one.        |
   90   |           Useful for visualizing solutions that are best expressed in       |
   91   |           a rotated system, but without loosing the reference to the        |
   92   |           original one.  This feature is accessed via the somewhat          |
   93   |           messy '--xtracoordsys' option, and its relatives                  |
   94   |           '--xtra_axislabels', '--xtra_axislengths' et. al.                 |
   95   |                                                                             |
   96   |           Also added the '--axislabels' option, enabling the user to        |
   97   |           specify his/her own axis labels instead of the default 'S_1',     |
   98   |           'S_2' and 'S_3'.                                                  |
   99   |                                                                             |
  100   |  000419:  (Stockholm, busy writing my PhD thesis) Added the '--auxsource'   |
  101   |  [v.1.6]  option, causing an auxiliary source file to be included at the    |
  102   |           end of the generated MetaPost source. Useful for including        |
  103   |           additional comments, labels etc.                                  |
  104   |                                                                             |
  105   |  030717:  (Stockholm) When compiling the source code of the poincare the    |
  106   |  [v.1.7]  program, the GCC compiler that I now use complains about that     |
  107   |           the string lengths of my arguments to fprintf() are greater than  |
  108   |           the length of 509 characters that ISO C89 compilers are required  |
  109   |           to support. Well, it is no secret that I consider the source to   |
  110   |           be quite messy (since it produces quite a lot of text output, and |
  111   |           in addition produces commented and readable MetaPost source).     |
  112   |               In order to fix this potential future source of trouble, I    |
  113   |           have chopped up all charcter strings greater than 500 characters  |
  114   |           into smaller pieces. This, on the other hand, has made the code   |
  115   |           look even worse; this is something that I plan to take care of    |
  116   |           sometime in the future, whenever I have the time for it.          |
  117   |                                                                             |
  118   |  030721:  (Stockholm) Corrected a bug in the block that append an external  |
  119   |  [v.1.8]  MetaPost file to the generated one (via the --auxsource option).  |
  120   |           (The bug caused strange characters to be included into the gene-  |
  121   |           rated MetaPost even when the --auxsource option was entirely      |
  122   |           omitted at the command line.) Also changed the --scalefactor      |
  123   |           argument to yield the radius of the Poincare sphere in milli-     |
  124   |           metres, instead of centimetres as previously used.                |
  125   |                                                                             |
  126   |  030722:  (Stockholm) Wrote a program 'stokeproc' for preprocessing of      |
  127   |  [v.1.9]  files containing Stokes parameters into the format as accepted    |
  128   |           by the Poincare program. The 'stokproc' program essentially       |
  129   |           three ASCII files, containing lists of the S1, S2, and S3 Stokes  |
  130   |           parameters, respectively, and generates an output file (out of    |
  131   |           these parameters), that is valid as direct input to the poincare  |
  132   |           program.                                                          |
  133   |                                                                             |
  134   |  030724:  (Stockholm) Finished with the blocks for parsing the input file   |
  135   | [v.1.10]  for possible labels at the end points of a trajectory on the      |
  136   |           Poincare sphere. The syntax of the input file is now:             |
  137   |               p [b <pos> <TeX labelstring at begin point of trajectory>]    |
  138   |               <S1[1]>  <S2[1]>  <S3[1]>   (1st Stokes param. triplet)       |
  139   |               <S1[2]>  <S2[2]>  <S3[2]>   (2nd Stokes param. triplet)       |
  140   |                 .        .        .                                         |
  141   |                 .        .        .                                         |
  142   |               <S1[N]>  <S2[N]>  <S3[N]>   (Nth Stokes param. triplet)       |
  143   |               q [e <pos> <TeX labelstring at end point of trajectory>]      |
  144   |           where <pos> is any of                                             |
  145   |               top   (puts label in top relative current point)              |
  146   |               ulft  (puts label in upper left relative current point)       |
  147   |               lft   (puts label in left relative current point)             |
  148   |               llft  (puts label in lower left relative current point)       |
  149   |               bot   (puts label in bottom relative current point)           |
  150   |               lrt   (puts label in lower right relative current point)      |
  151   |               rt    (puts label in right relative current point)            |
  152   |               urt   (puts label in upper right relative current point)      |
  153   |           This syntax is now also incorporated into the data preprocessor   |
  154   |           'stokeproc', which now accepts the command line options           |
  155   |           '--beginlabel' and '--endlabel' for simplifying the end point     |
  156   |           labeling. See comments in the C source code of 'stokeproc' for    |
  157   |           further information.                                              |
  158   |                                                                             |
  159   |  031217:  (Wuppertal) Added a few blocks to ensure that any spaces,         |
  160   | [v.1.11]  linefeeds, or other auxiliary characters are properly read        |
  161   |           away at the end of every Stokes-parameter path that is supplied   |
  162   |           as input to the program. The breaking of a path into subpaths     |
  163   |           works properly whenever the path is passing from the visible to   |
  164   |           the hidden hemisphere of the Poincare sphere; however, there      |
  165   |           still seem to remain some bug that prohibits "makapath makepen"   |
  166   |           to be written to the output file whenever a new path (beginning   |
  167   |           with the character 'p') starts after the very first path is       |
  168   |           finished. In addition, endlabel statements (e.g. 'q e $t=10$')    |
  169   |           seem to fail in the positioning, so there still seem to be plenty |
  170   |           to fix before considering all bugs to be removed.                 |
  171   |                                                                             |
  172   |  040318:  (Wuppertal) Added the '--draw_paths_as_arrows option', enabling   |
  173   | [v.1.12]  the poincare program to draw multiple paths as arrows, rather     |
  174   |           than just as solid curves. Nice option to use if one wish to show |
  175   |           a general direction of polarization state evolution on the sphere |
  176   |           irregardless of the orientation of the polarization ellipse.      |
  177   |           The source code was also transferred to my new Mac OS X machine   |
  178   |           and compiled with the Mac port of GCC without any complaints even |
  179   |           when fed with the '-Wall -pedantic' option.                       |
  180   |                                                                             |
  181   |  050330:  (Cork) Today started the catharsis of the poincare program by     |
  182   | [v.1.13]  collecting all variables related to the Poincare map generation   |
  183   |           into one single C struct, and using this struct for initializa-   |
  184   |           tion and command-line parsing of variables. This is the first     |
  185   |           step towards obtaining a more readable program, and doing so I    |
  186   |           separated the code into the following, more self-explanatory      |
  187   |           routines:                                                         |
  188   |               parse_command_line(argc,argv);                                |
  189   |               display_arrow_specs(map);                                     |
  190   |               open_outfile(map);                                            |
  191   |               write_header(outfileptr,map,argc,argv);                       |
  192   |               write_euler_angle_specs(outfileptr,map);                      |
  193   |               write_sphere_shading_specs(outfileptr,map);                   |
  194   |               write_shaded_sphere(outfileptr,map);                          |
  195   |               write_equators(outfileptr,map);                               |
  196   |               write_scanned_trajectories(outfileptr,map);                   |
  197   |               write_additional_arrows(outfileptr,map);                      |
  198   |               write_coordinate_axes(outfileptr,map);                        |
  199   |               write_additional_coordinate_axes(outfileptr,map);             |
  200   |               write_included_auxiliary_source(outfileptr,map);              |
  201   |               generate_eps_image(map);                                      |
  202   |           There still remains to be something done about the huge           |
  203   |           write_scanned_trajectories() routine, which deals with flushing   |
  204   |           the Stokes trajectories to the MetaPost file. This routine        |
  205   |           involves the splitting of the trajectories into visible and       |
  206   |           hidden parts, as well as all labeling and tick marking, and is    |
  207   |           at the time of writing quite hard to read. This is the first      |
  208   |           thing that needs attention in the progress of cleaning up the     |
  209   |           poincare program.                                                 |
  210   |           In the verification of todays version I ran the 'examples' block  |
  211   |           of the Makefile enclosed with the source code, without any        |
  212   |           deviations from the expected result. However, there still seem    |
  213   |           to be lurking a bug around which causes unpredictable generation  |
  214   |           of multiple Stokes trajectories in the same Poincare map,         |
  215   |           something that needs to be attended in the future.                |
  216   |                                                                             |
  217   |  050402:  (Cork) Today I took the giant leap forward and entirely rewrote   |
  218   | [v.1.14]  the algorithms for the generation of trajectory maps. In the form |
  219   |           as now employed, the entire work is performed with the following  |
  220   |           newly written routines:                                           |
  221   |               short new_trajectory()                                        |
  222   |               short end_of_trajectory()                                     |
  223   |               short beginlabel()                                            |
  224   |               short endlabel()                                              |
  225   |               void readaway_comments_and_blanks()                           |
  226   |               void scan_label()                                             |
  227   |               void scan_beginlabel()                                        |
  228   |               void scan_endlabel()                                          |
  229   |               void scan_for_stokes_triplet()                                |
  230   |               short visible()                                               |
  231   |               short point_just_became_hidden()                              |
  232   |               short point_just_became_visible()                             |
  233   |               void get_screen_coordinates()                                 |
  234   |               void add_subtrajectory()                                      |
  235   |               void sort_out_visible_and_hidden()                            |
  236   |               void add_hidden_subtrajectories()                             |
  237   |               void add_visible_subtrajectories()                            |
  238   |               void add_scanned_trajectory()                                 |
  239   |               void add_scanned_labels()                                     |
  240   |               void write_scanned_trajectories()                             |
  241   |           Checked multiple path specifications within one single map with   |
  242   |           success, using the {\tt example-c} block of the Makefile enclosed |
  243   |           with the poincare source. Label positioning and multiple labels   |
  244   |           functionality works as expected.                                  |
  245   |                                                                             |
  246   |  050403:  (Cork) Added the --hiddengraytone option, and also added the      |
  247   | [v.1.15]  functionality that one in the input trajectory file anywhere      |
  248   |           after a Stokes parameter triplet can specify that this point      |
  249   |           should get a tick mark (simply using a 't' letter after the       |
  250   |           triplet (s1,s2,s3). This functionality is highly useful whenever  |
  251   |           one wish to put the scale of evolution of the Stokes vector       |
  252   |           explicitly in the map. (For example with the Stokes vector as     |
  253   |           function of wavelength or strain, as I currently am in the stage  |
  254   |           of investigating for the sinusoidal gratings with phase           |
  255   |           discontinuities.)                                                 |
  256   |                                                                             |
  257   |  050404:  (Cork) Added the option of including tick mark labeling in the    |
  258   | [v.1.16]  input file syntax, so that tick marks may be given labels if      |
  259   |           wanted. This adds to yesterdays implemented tick marks, so that   |
  260   |           one now can have the absolute scale of measure explicitly present |
  261   |           in the figure, with numbers and all. The current input file       |
  262   |           syntax after this inclusion yields (assuming M points forming the |
  263   |           trajectory):                                                      |
  264   |               p [b <pos> "<TeX label string at begin point>"]               |
  265   |               <s1(1)> <s2(1)> <s3(1)> [options]                             |
  266   |               <s1(2)> <s2(2)> <s3(2)> [options]                             |
  267   |                  .       .       .       .                                  |
  268   |               <s1(M)> <s2(M)> <s3(M)> [options]                             |
  269   |               q [e <pos> "<TeX label string at end point>"]                 |
  270   |           where each of the rows can have individually specified options of |
  271   |           the form                                                          |
  272   |               options = [t [l <pos> "<text>"]] [% <comments>]               |
  273   |                                                                             |
  274   |  050419:  (Southampton) Today fixed a bug in the program which in the case  |
  275   | [v.1.17]  of multiple sets of Stokes trajectories caused previously written |
  276   |           visible parts of earlier trajectories to be crossed by hidden     |
  277   |           parts of trajectories added later on. The issue was simply solved |
  278   |           by adding the |viewtype| flag as input variable to the            |
  279   |           |add_scanned_trajectory()|, |add_scanned_tickmarks()|, and        |
  280   |           |write_scanned_trajectories()| routines, so as to be able to call |
  281   |           |write_scanned_trajectories()| twice, causing the routine to      |
  282   |           perform the parsing of all trajectories twice and adding the      |
  283   |           visible parts only after all hidden parts have been written to    |
  284   |           file. Also fixed a minor bug which under some circumstances       |
  285   |           caused the program to write any additional Stokes trajectories    |
  286   |           with a thinner pen stroke than otherwise specified by the         |
  287   |           numerical value of |map.paththickness|.                           |
  288   |                                                                             |
  289   |  050420:  (Southampton) Fixed a bug which potentially would cause the       |
  290   | [v.1.18]  program to crash for a very large number of tickmarks in the      |
  291   |           Stokes trajectory.                                                |
  292   |                                                                             |
  293   |  050514:  (Southampton) Added a routine which in verbose mode displays      |
  294   | [v.1.19]  the command-line options as they are parsed. This is useful in    |
  295   |           debugging batch-mode calls from, for example, Makefiles.          |
  296   |                                                                             |
  297   |  050515:  (Southampton) Added the --reverse_arrow_paths option, which       |
  298   | [v.1.20]  is used to change the direction of traversal in drawing paths     |
  299   |           as arrows using the --draw_paths_as_arrows option. Useful when    |
  300   |           sampled data have an order which do not coincide with the         |
  301   |           natural order of traversal in the Poincare map.                   |
  302   |                                                                             |
  303   |  050810:  (Southampton) Wrote the strip_away_path() routine in order to     |
  304   | [v.1.21]  finally solve the problem with long path strings that appear in   |
  305   |           the program name string whenever poincare is called with an       |
  306   |           explicit path specified at the command line. The call to the      |
  307   |           strip_away_path() routine is located in the beginning of the      |
  308   |           block for parsing the command line.                               |
  309   |                                                                             |
  310   |  111215:  (Stockholm) Fixed a few rather clumsy blocks in the hidden path   |
  311   | [v.1.22]  generation. Updated the Makefile.                                 |
  312   |                                                                             |
  313   |  140124:  (Stockholm) This morning I received an email from a person in     |
  314   | [v.1.23]  Germany who tried out the program but found that for some reason  |
  315   |           all S1 coordinates seemed to have their signs dropped, i.e. as    |
  316   |           if its absolute value had been extracted. It turned out that the  |
  317   |           sign-dropping was due to a change in the way the standard         |
  318   |           isalnum() function interprets signs. Previously--at least in the  |
  319   |           GNU development environment distributed with Linux (1996-2000)    |
  320   |           and OSX 10.3 (2003)--I had no problems whatsoever to parse the    |
  321   |           signs of signed numbers as belonging to the alphanumeric set;     |
  322   |           however, with newer compilers this behaviour of isalnum() seems   |
  323   |           to have changed, despite that we here follow a strict ISO C90     |
  324   |           standard. Anyway, thanks for pointing this out, and the code is   |
  325   |           now fixed and operational again.                                  |
  326   |                                                                             |
  327   | Example of usage (the figure on the front page of my PhD thesis):           |
  328   |                                                                             |
  329   |    poincare --normalize --verbose  --bezier --draw_hidden_dashed \          |
  330   |               --axislengths 0.3 1.5 0.3 2.7 0.3 1.5 \                       |
  331   |               --axislabels  "s_1"  "s_2"  "s_3,w_3" \                       |
  332   |               --rotatephi 15.0  --rotatepsi -70.0 \                         |
  333   |               --xtracoordsys  7.0181217  0.0 \                              |
  334   |               --xtracoordsys_axislabel_x "w_1" \                            |
  335   |               --xtracoordsys_axislabel_y "w_2" \                            |
  336   |               --xtracoordsys_axislengths 0.3 1.8 0.3 1.9 0.3 1.8 \          |
  337   |               --shading 0.75 0.99 \                                         |
  338   |               --rhodivisor 50  --phidivisor 80  --scalefactor 2.5 \         |
  339   |               --inputfile tfig2.dat  --outputfile tfig2.mp \                |
  340   |               --paththickness 0.8 --arrowthickness 0.4                      |
  341   |                                                                             |
  342   =============================================================================*/
  343   #include <math.h>
  344   #include <stdio.h>
  345   #include <string.h>
  346   #include <stdlib.h>
  347   #include <time.h>
  348   #include <ctype.h>
  349   
  350   #define VERSION_NUMBER "1.23"
  351   
  352   /*-----------------------------------------------------------------------------
  353   | PI should have been defined from <math.h>, but let's be on the safe side
  354   | (since I have found strange behaviour on some SPARC:s).
  355   -----------------------------------------------------------------------------*/
  356   #ifndef PI
  357   #define PI (3.14159265358979323846)
  358   #endif
  359   
  360   /*-----------------------------------------------------------------------------
  361   | Definitions of maximum number of allowed coordinates and labels per
  362   | trajectory, as used in the allocation of memory. These fixed parameters
  363   | determine the following:
  364   |    MAX_NUM_STOKE_COORDS    Maximum number of coordinates per trajectory
  365   |    MAX_NUM_TICKMARKS       Maximum number of tick marks per trajectory
  366   |    MAX_NUM_LABELS          Maximum number of text labels per trajectory
  367   |    MAX_LABEL_TEXTLENGTH    Maximum number of characters per text label
  368   -----------------------------------------------------------------------------*/
  369   #define MAX_NUM_STOKE_COORDS (5000)
  370   #define MAX_NUM_TICKMARKS ((int)(MAX_NUM_STOKE_COORDS/10))
  371   #define MAX_NUM_LABELS ((int)(MAX_NUM_TICKMARKS/10))
  372   #define MAX_LABEL_TEXTLENGTH (256)
  373   #define MAX_FILENAME_TEXTLENGTH (256)
  374   
  375   #define SUCCESS 0 /* Return code for successful program termination */
  376   #define FAILURE 1 /* Return code for unsuccessful program termination */
  377   
  378   #define DEFAULT_OUTFILENAME "aout.mp"
  379   #define DEFAULT_EPSJOBNAME "aout"
  380   #define DEFAULT_AXISLABEL_S1 "S_1"
  381   #define DEFAULT_AXISLABEL_S2 "S_2"
  382   #define DEFAULT_AXISLABEL_S3 "S_3"
  383   #define DEFAULT_AXISLABELPOSITION_S1 "urgt"
  384   #define DEFAULT_AXISLABELPOSITION_S2 "urgt"
  385   #define DEFAULT_AXISLABELPOSITION_S3 "urgt"
  386   
  387   #define DEFAULT_ROT_PSI   (-40.0*(PI/180))      /* Angle in radians */
  388   #define DEFAULT_ROT_PHI   (15.0*(PI/180))       /* Angle in radians */
  389   
  390   #define DEFAULT_PHI_SOURCE   (30.0*(PI/180))    /* Angle in radians */
  391   #define DEFAULT_THETA_SOURCE (30.0*(PI/180))    /* Angle in radians */
  392   
  393   #define DEFAULT_MAX_WHITENESS (0.99)  /*  '0.0' <=> black; '1.0' <=> white */
  394   #define DEFAULT_MIN_WHITENESS (0.75)  /*  '0.0' <=> black; '1.0' <=> white */
  395   #define DEFAULT_HIDDEN_GRAYTONE (0.65)
  396   
  397   #define DEFAULT_RHO_DIVISOR  (50.0)
  398   #define DEFAULT_PHI_DIVISOR  (80.0)
  399   
  400   #define DEFAULT_POSITIVE_AXIS_LENGTH (1.5)
  401   #define DEFAULT_NEGATIVE_AXIS_LENGTH (0.1)
  402   
  403   #define DEFAULT_PATH_THICKNESS  (1.0)
  404   #define DEFAULT_ARROW_THICKNESS (0.6)
  405   #define DEFAULT_ARROW_HEADANGLE (30.0)
  406   #define DEFAULT_TICKSIZE (4*DEFAULT_PATH_THICKNESS)
  407   
  408   /* Number of coordinates per line in the generated MetaPost code for the map */
  409   #define NUM_COORDS_PER_METAPOST_LINE (3)
  410   
  411   /* Definitions of flags used for identifying label positions */
  412   #define NOLABEL         (0)
  413   #define TOPLABEL        (1)
  414   #define UPPERLEFTLABEL  (2)
  415   #define LEFTLABEL       (3)
  416   #define LOWERLEFTLABEL  (4)
  417   #define BOTTOMLABEL     (5)
  418   #define LOWERRIGHTLABEL (6)
  419   #define RIGHTLABEL      (7)
  420   #define UPPERRIGHTLABEL (8)
  421   
  422   #define MAKE_VERBOSE_REALLY_VERBOSE (1)  /* Avoid TOO much of ASCII */
  423   
  424   /* Definition of values taken by flags determining whether hidden or visible
  425    * parts of trajectories should be flushed to file.
  426    */
  427   #define HIDDEN (0)
  428   #define VISIBLE (1)
  429   
  430   /*---------------------------------------------------------------------
  431   | The only global variables allowed in my programs are `optarg`, which is
  432   | the string of characters that specified the call from the command line,
  433   | and `progname`, which simply is the string containing the name of the
  434   | program, as it was invoked from the command line.
  435   ---------------------------------------------------------------------*/
  436   extern char *optarg;
  437   char *progname;
  438   
  439   typedef struct {
  440      double **arrows;
  441      int numarrows;
  442      short verbose;
  443      short save_memory;
  444      short use_normalized_stokes_params;
  445      short use_bezier_curves;
  446      short user_specified_inputfile;
  447      short user_specified_auxfile;
  448      short user_specified_axislabels;
  449      short user_specified_additional_coordinate_system;
  450      short user_specified_xtra_axislabel_x;
  451      short user_specified_xtra_axislabel_y;
  452      short user_specified_xtra_axislabel_z;
  453      short draw_hidden_dashed;
  454      short draw_paths_as_arrows;
  455      short reverse_arrow_paths;
  456      short last_point_infront;
  457      short current_point_is_a_beginlabelpoint,
  458            current_point_is_an_endlabelpoint;
  459      short draw_axes_inside_sphere;
  460      short currently_drawing_path;
  461      short generate_eps_output;
  462      char infilename[MAX_FILENAME_TEXTLENGTH];
  463      char outfilename[MAX_FILENAME_TEXTLENGTH];
  464      char auxfilename[MAX_FILENAME_TEXTLENGTH];
  465      char epsjobname[MAX_FILENAME_TEXTLENGTH];
  466      char axislabel_s1[MAX_LABEL_TEXTLENGTH];
  467      char axislabel_s2[MAX_LABEL_TEXTLENGTH];
  468      char axislabel_s3[MAX_LABEL_TEXTLENGTH];
  469      char axislabelposition_s1[8],
  470           axislabelposition_s2[8],
  471           axislabelposition_s3[8];
  472      char xtra_axislabel_x[MAX_LABEL_TEXTLENGTH];
  473      char xtra_axislabel_y[MAX_LABEL_TEXTLENGTH];
  474      char xtra_axislabel_z[MAX_LABEL_TEXTLENGTH];
  475      double xtra_neg_axis_length_x;
  476      double xtra_neg_axis_length_y;
  477      double xtra_neg_axis_length_z;
  478      double xtra_pos_axis_length_x;
  479      double xtra_pos_axis_length_y;
  480      double xtra_pos_axis_length_z;
  481      char labelstr_beginpoint[MAX_LABEL_TEXTLENGTH];
  482      char labelstr_endpoint[MAX_LABEL_TEXTLENGTH];
  483      double scalefactor;
  484      double rot_psi, rot_phi;
  485      double delta_rot_psi, delta_rot_phi;
  486      double phi_source, theta_source;
  487      double upper_whiteness_value;
  488      double lower_whiteness_value;
  489      double hiddengraytone;
  490      double rho_divisor;
  491      double phi_divisor;
  492      double xpos_beginpoint;
  493      double ypos_beginpoint;
  494      double xpos_endpoint;
  495      double ypos_endpoint;
  496      double neg_axis_length_s1;
  497      double neg_axis_length_s2;
  498      double neg_axis_length_s3;
  499      double pos_axis_length_s1;
  500      double pos_axis_length_s2;
  501      double pos_axis_length_s3;
  502      double paththickness;
  503      double arrowthickness;
  504      double arrowheadangle;
  505      double coordaxisthickness;
  506      double ticksize;
  507   } pmap;
  508   
  509   typedef struct {
  510      long numcoords;
  511      double *s1,*s2,*s3;
  512      short *visible;
  513      int numtickmarks;
  514      long *tickmark;
  515      int numlabels;
  516      long *label;
  517      char **labeltext;
  518      int *labellength;
  519      short *labelpos;
  520   } stoketraject;
  521   
  522   /*----------------------------------------------------------------------------
  523   | The |svector| routine allocates a vector of short integer precision,
  524   | with vector index ranging from |nl| to |nh|.
  525   ----------------------------------------------------------------------------*/
  526   short *svector(long nl, long nh) {
  527      short *v;
  528      v=(short *)malloc((size_t) ((nh-nl+2)*sizeof(short)));
  529      if (!v) {
  530         fprintf(stderr,"Error: Allocation failure in svector()\n");
  531         exit(1);
  532      }
  533      return v-nl+1;
  534   }
  535   
  536   /*----------------------------------------------------------------------------
  537   | The |ivector| routine allocates a vector of integer precision,
  538   | with vector index ranging from |nl| to |nh|.
  539   ----------------------------------------------------------------------------*/
  540   int *ivector(long nl, long nh) {
  541      int *v;
  542      v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int)));
  543      if (!v) {
  544         fprintf(stderr,"Error: Allocation failure in ivector()\n");
  545         exit(1);
  546      }
  547      return v-nl+1;
  548   }
  549   
  550   /*----------------------------------------------------------------------------
  551   | The |lvector| routine allocates a vector of long integer precision,
  552   | with vector index ranging from |nl| to |nh|.
  553   ----------------------------------------------------------------------------*/
  554   long *lvector(long nl, long nh) {
  555      long *v;
  556      v=(long *)malloc((size_t) ((nh-nl+2)*sizeof(long)));
  557      if (!v) {
  558         fprintf(stderr,"Error: Allocation failure in lvector()\n");
  559         exit(1);
  560      }
  561      return v-nl+1;
  562   }
  563   
  564   /*----------------------------------------------------------------------------
  565   | The |free_svector| routine release the memory occupied by the
  566   | short integer vector |v[nl..nh]|.
  567   ----------------------------------------------------------------------------*/
  568   void free_svector(short *v, long nl, long nh) {
  569      free((char*) (v+nl-1));
  570   }
  571   
  572   /*----------------------------------------------------------------------------
  573   | The |free_ivector| routine release the memory occupied by the
  574   | integer vector |v[nl..nh]|.
  575   ----------------------------------------------------------------------------*/
  576   void free_ivector(long *v, long nl, long nh) {
  577      free((char*) (v+nl-1));
  578   }
  579   
  580   /*----------------------------------------------------------------------------
  581   | The |free_lvector| routine release the memory occupied by the
  582   | long integer vector |v[nl..nh]|.
  583   ----------------------------------------------------------------------------*/
  584   void free_lvector(long *v, long nl, long nh) {
  585      free((char*) (v+nl-1));
  586   }
  587   
  588   /*----------------------------------------------------------------------------
  589   | The |free_dvector| routine release the memory occupied by the
  590   | real-valued vector |v[nl..nh]|.
  591   ----------------------------------------------------------------------------*/
  592   void free_dvector(double *v, long nl, long nh) {
  593      free((char*) (v+nl-1));
  594   }
  595   
  596   /*----------------------------------------------------------------------------
  597   | The |dvector| routine allocates a real-valued vector of double precision,
  598   | with vector index ranging from |nl| to |nh|.
  599   ----------------------------------------------------------------------------*/
  600   double *dvector(long nl, long nh) {
  601      double *v;
  602      v=(double *)malloc((size_t) ((nh-nl+2)*sizeof(double)));
  603      if (!v) {
  604         fprintf(stderr,"Error: Allocation failure in dvector()\n");
  605         exit(1);
  606      }
  607      return v-nl+1;
  608   }
  609   
  610   /*-------------------------------------------------------------------------
  611   | The scan_for_boundingbox(infilename,llx,lly,urx,ury) routine takes the
  612   | name of a regular ASCII text file (infilename) containing Encapsulated
  613   | PostScript (EPS) as input, and returns the bounding box of the figure,
  614   | in terms of the corner coordinates given by llx (lower left x), lly
  615   | (lower left y), urx (upper right x), and ury (upper right y).
  616   | Typical usage (as a sample program for extraction of bounding box, similar
  617   | to the UNIX standard program psbb):
  618   |
  619   |     #include <stdlib.h>
  620   |     #include <stdio.h>
  621   |     #include <string.h>
  622   |
  623   |     extern char *optarg;
  624   |     char *progname;
  625   |
  626   |     void scan_for_boundingbox(char infilename[256],
  627   |        long int *llx,long int *lly,long int *urx,long int *ury) {
  628   |        ...
  629   |     }
  630   |
  631   |     int main(int argc,char *argv[]) {
  632   |        char infilename[256];
  633   |        int no_arg;
  634   |        long int llx,lly,urx,ury;
  635   |
  636   |        progname=argv[0];
  637   |        no_arg=argc;
  638   |        if (argc<2) {
  639   |           fprintf(stderr,"%s: Error! Please specify input (EPS) filename.\n",
  640   |              progname);
  641   |           fprintf(stderr,"Usage: %s <filename>\n",progname);
  642   |           exit(FAILURE);
  643   |        }
  644   |        --argc;
  645   |        strcpy(infilename,argv[no_arg-argc]);
  646   |        scan_for_boundingbox(infilename,&llx,&lly,&urx,&ury);
  647   |        fprintf(stdout,"%ld %ld %ld %ld\n",llx,lly,urx,ury);
  648   |        exit(0);
  649   |     }
  650   |
  651   -------------------------------------------------------------------------*/
  652   void scan_for_boundingbox(char infilename[256],
  653      long int *llx,long int *lly,long int *urx,long int *ury) {
  654      short int done;
  655      char tmpstr[256];
  656      int tmpch;
  657      FILE *infile;
  658   
  659      if ((infile=fopen(infilename,"r"))==NULL) {
  660         fprintf(stderr,"%s: Error! Could not open file %s for reading.\n",
  661            progname,infilename);
  662         exit(FAILURE);
  663      }
  664      fseek(infile,0L,SEEK_SET);
  665      done=0;
  666      while (((tmpch=getc(infile))!=EOF)&&(!done)) {
  667         ungetc(tmpch,infile);
  668         fscanf(infile,"%s",tmpstr);
  669         if (!strcmp(tmpstr,"%%BoundingBox:")) {
  670            done=1;
  671            fscanf(infile,"%ld",llx);
  672            fscanf(infile,"%ld",lly);
  673            fscanf(infile,"%ld",urx);
  674            fscanf(infile,"%ld",ury);
  675         }
  676      }
  677      if (tmpch==EOF) {
  678         fprintf(stderr,"%s: Error! End of file reached without ",progname);
  679         fprintf(stderr,"finding any %%%%BoundingBox statement!\n");
  680         fprintf(stderr,"%s: (Does %s really contain ",progname,infilename);
  681         fprintf(stderr,"Encapsulated PostScript?)\n");
  682         exit(FAILURE);
  683      }
  684   } /* end of scan_for_boundingbox() */
  685   
  686   double **matrix(long nrl, long nrh, long ncl, long nch) {
  687      long i, j, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  688      double **m;
  689      m=(double **) malloc((size_t)((nrow+1)*sizeof(double*)));
  690      if (!m) {
  691         printf("%s: Allocation failure 1 in matrix()\n",progname);
  692         exit(FAILURE);
  693      }
  694      m += 1;
  695      m -= nrl;
  696      m[nrl]=(double *) malloc((size_t)((nrow*ncol+1)*sizeof(double)));
  697      if (!m[nrl]) {
  698         printf("%s: Allocation failure 2 in matrix()\n",progname);
  699         exit(FAILURE);
  700      }
  701      m[nrl] += 1;
  702      m[nrl] -= ncl;
  703      for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
  704      for(i=nrl;i<=nrh;i++) for(j=ncl;j<=nch;j++) m[i][j]=0.0;
  705      return m;
  706   }
  707   
  708   void free_matrix(double **m, long nrl, long nrh, long ncl, long nch) {
  709      free((char*) (m[nrl]+ncl-1));
  710      free((char*) (m+nrl-1));
  711   }
  712   
  713   char **cmatrix(long nrl, long nrh, long ncl, long nch) {
  714      long i, j, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  715      char **m;
  716      m=(char **) malloc((size_t)((nrow+1)*sizeof(char*)));
  717      if (!m) {
  718         printf("%s: Allocation failure 1 in cmatrix()\n",progname);
  719         exit(1);
  720      }
  721      m += 1;
  722      m -= nrl;
  723      m[nrl]=(char *) malloc((size_t)((nrow*ncol+1)*sizeof(char)));
  724      if (!m[nrl]) {
  725         printf("%s: Allocation failure 2 in cmatrix()\n",progname);
  726         exit(1);
  727      }
  728      m[nrl] += 1;
  729      m[nrl] -= ncl;
  730      for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
  731      for(i=nrl;i<=nrh;i++) for(j=ncl;j<=nch;j++) m[i][j]=0.0;
  732      return m;
  733   }
  734   
  735   void free_cmatrix(char **m, long nrl, long nrh, long ncl, long nch) {
  736      free((char*) (m[nrl]+ncl-1));
  737      free((char*) (m+nrl-1));
  738   }
  739   
  740   void show_banner(void) {
  741      fprintf(stdout,"This is %s v.%s.  ",progname,VERSION_NUMBER);
  742      fprintf(stdout,"Copyright (C) 1997-2005, Fredrik Jonsson\n");
  743   }
  744   
  745   void showsomehelp(void) {
  746      fprintf(stdout,"Usage: %s [options]\n",progname);
  747      fprintf(stdout,
  748    "where [options] include:\n"
  749    " -h, --help              Show this help-message and exit clean.\n"
  750    " -v, --verbose           Toggle verbose mode; show beautiful ASCII.\n"
  751    "                         Default: off.\n"
  752    " -s, --save_memory       Toggle memory save mode; spare some RAM.\n"
  753    "                         Default: off.\n"
  754    " -V, --version           Show version and exit clean.\n\n");
  755      fprintf(stdout,
  756    " -f, --inputfile <name>  Read input Stokes-parameters from file <name>.\n"
  757    "                         The input file may contain an arbitrary number of\n"
  758    "                         separate paths, and the format of the file is\n");
  759      fprintf(stdout,
  760    "                               p\n"
  761    "                               <s1> <s2> <s3>       [1:st triplet, path 1]\n"
  762    "                               <s1> <s2> <s3>       [2:nd triplet, path 1]\n"
  763    "                                  .  .  .\n"
  764    "                               <s1> <s2> <s3>       [M:th triplet, path 1]\n"
  765    "                               q\n");
  766      fprintf(stdout,
  767    "                               p\n"
  768    "                               <s1> <s2> <s3>       [1:st triplet, path 2]\n"
  769    "                                  .  .  .\n"
  770    "                               <s1> <s2> <s3>       [N:th triplet, path 2]\n"
  771    "                               q\n"
  772    "                                  .  .  .                   [etc.]\n");
  773      fprintf(stdout,
  774    "                         Thus, each separate path should be defined by an\n"
  775    "                         initial 'p', after which each following row con-\n"
  776    "                         tains a triplet of Stokes parameters.  After the\n"
  777    "                         Stokes parameter triplet, comments and additional\n"
  778    "                         information may be written (ignored by program)\n"
  779    "                         until linefeed.  Finally, each separate path is\n"
  780    "                         is ended with a 'q' on a separate line.\n"
  781    "\n");
  782      fprintf(stdout,
  783    " --paththickness <val>   Specifies the thickness in PostScript points (pt)\n"
  784    "                         of the path to draw.  Default: <val> = 1.0 [pt].\n"
  785    "                         [1 pt == 1/72 inch]\n"
  786    "\n");
  787      fprintf(stdout,
  788    " --draw_hidden_dashed    Toggles between drawing of hidden parts of the\n"
  789    "                         specified path with dashed and solid lines.\n"
  790    "                         Default: off. (Solid lines)\n"
  791    "\n");
  792      fprintf(stdout,
  793    " --draw_paths_as_arrows  Draw all specified trajectories as arrowed\n"
  794    "                         curves, with arrowheads at the and point.\n"
  795    "                         This option is useful whenever one wish to,\n"
  796    "                         for example, show on the direction of evolution\n"
  797    "                         of a certain trajectory, or the direction of\n"
  798    "                         rotation of the Stokes vector in a circular\n");
  799      fprintf(stdout,
  800    "                         path. With this option it is often useful to\n"
  801    "                         chop up the trajectory of the original input\n"
  802    "                         file into subtrajectories, so as to create\n"
  803    "                         multiple arrow heads in the same trajectory.\n"
  804    "                         See also the --reverse_arrow_paths option.\n");
  805      fprintf(stdout,
  806    " --reverse_arrow_paths   Reverse the direction of all arrows drawn using\n"
  807    "                         the --draw_paths_as_arrows option. This is useful\n"
  808    "                         if the sampled trajectory data are not ordered in\n"
  809    "                         the natural direction of trajectory traversal.\n");
  810      fprintf(stdout,
  811    " --auxsource <name>      Causes the auxiliary file <name> to be included\n"
  812    "                         at the end of the generated MetaPost source.\n"
  813    "                         Useful for including additional comments, labels\n"
  814    "                         etc. in the figure.\n"
  815    "\n");
  816      fprintf(stdout,
  817    " --arrowthickness <val>  Analogous to the '--paththickness' option, but\n"
  818    "                         with the difference that this one applies to\n"
  819    "                         (eventually occuring) the thickness of additional\n"
  820    "                         arrows to be drawn with the '--arrow' option.\n"
  821    "                         Default: <val> = 0.6 [pt].\n"
  822    "\n");
  823      fprintf(stdout,
  824    " --arrowheadangle <deg>  Specifies the head angle of any arrows used in\n"
  825    "                         the mapping of Stokes parameters on the Poincare\n"
  826    "                         sphere. Notice that this does not affect the\n"
  827    "                         head angles of the arrows of the coordinate\n"
  828    "                         axes. Default value: 30 degrees.\n");
  829      fprintf(stdout,
  830    " -b, --bezier            Toggle Bezier mode, in which Bezier interpolation\n"
  831    "                         is used in order to obtain smooth paths for the\n"
  832    "                         input trajectory(-ies), specified with the '-f'\n"
  833    "                         option.  Otherwise regular piecewise stright-line\n"
  834    "                         type lines are used.   Default: off.\n"
  835    "\n");
  836      fprintf(stdout,
  837    " -o, --outputfile <name> Write output MetaPost-code [1] to file <name>.\n"
  838    "\n");
  839      fprintf(stdout,
  840    " -e, --epsoutput <name>  In addition to just generating MetaPost-code for\n"
  841    "                         the figure, also try to generate a complete EPS\n"
  842    "                         (Encapsulated PostScript) figure, using <name>\n"
  843    "                         as the base name for the job. This option relies\n"
  844    "                         on system calls for TeX, MetaPost, and DVIPS, and\n"
  845    "                         relies on that they are properly installed in the\n"
  846    "                         system environment.\n");
  847      fprintf(stdout,
  848    "                         The EPS output and the intermediate TeX, DVI, and\n"
  849    "                         and log files will from the base name be named\n"
  850    "                         <name>.eps,<name>.tex,<name>.dvi, and <name>.log,\n"
  851    "                         respectively.\n"
  852    "\n");
  853      fprintf(stdout,
  854    "--psi, --rotatepsi <val> When mapping Poincare-sphere and corresponding\n"
  855    "                         coordinate-system (S_1,S_2,S_3), first rotate\n"
  856    "                         angle psi == <val> around the 'z'-axis (S_3).\n"
  857    "                         Default: -40.0 (Degrees)\n"
  858    "\n");
  859      fprintf(stdout,
  860    "--phi, --rotatephi <val> When mapping Poincare-sphere and corresponding\n"
  861    "                         coordinate-system (S_1,S_2,S_3), after the first\n"
  862    "                         rotation (psi above), rotate angle phi == <val>\n"
  863    "                         around the 'y'-axis (S_2).\n"
  864    "                         Default: 15.0 (Degrees)\n"
  865    "\n");
  866      fprintf(stdout,
  867    " --rhodivisor  <val>     Number of segments in radial direction of the 2D-\n"
  868    "                         mapped Poincare sphere.  Default: 50.\n"
  869    "\n");
  870      fprintf(stdout,
  871    " --phidivisor  <val>     Number of segments in tangential direction of the\n"
  872    "                         2D-mapped Poincare sphere.  Default: 80.\n"
  873    "\n");
  874      fprintf(stdout,
  875    " --scalefactor <val>     Specifies the radius of the printed Poincare\n"
  876    "                         sphere (Encapsulated PostScript) in millimetres.\n"
  877    "\n");
  878      fprintf(stdout,
  879    " --shading <w1> <w2>     Specifies the minimum (<w1>) and maximum (<w2>)\n"
  880    "                         whiteness values of the Poincare sphere to draw\n"
  881    "                         (using the Phong shading algorithm).\n"
  882    "                         Here:\n"
  883    "                             <wx> == 0.0  corresponds to 'white'\n"
  884    "                             <wx> == 1.0  corresponds to 'white'\n"
  885    "                         Default values:  <w1> == 0.65,  <w2> == 0.99\n"
  886    "\n");
  887      fprintf(stdout,
  888    " --hiddengraytone <w>    Specifies the whiteness to be used in drawing\n"
  889    "                         trajectory parts that are hidden behind the\n"
  890    "                         Poincare sphere.\n"
  891    "                             <w> == 0.0  corresponds to black,\n"
  892    "                             <w> == 1.0  corresponds to white,\n");
  893      fprintf(stdout,
  894    " --axislengths <v>       Specifies the lengths of negative and positive\n"
  895    "                         parts of the coordinate axes, on the form\n"
  896    "                           <v> = <xmin> <xmax> <ymin> <ymax> <zmin> <zmax>\n"
  897    "                         with 'x' as the s1-axis, 'y' as the s2-axis, and\n"
  898    "                         'z' as the s3-axis. All values are taken relative\n"
  899    "                         to the radius of the Poincare sphere; thus <v>=1\n");
  900      fprintf(stdout,
  901    "                         correspond to the radius, while <v>=1.5 corre-\n"
  902    "                         spond to an axis length such that 50 percent of\n"
  903    "                         the axis is showed outside the Poincare sphere.\n"
  904    "                         Default:  <xmin> = <ymin> = <zmin> = 0.3 (30 %%)\n"
  905    "                                   <xmax> = <ymax> = <zmax> = 1.5 (150 %%)\n"
  906    "\n");
  907      fprintf(stdout,
  908    " --axislabels <s>        Specifies the labels of the coordinate axes, on\n"
  909    "                         the form\n"
  910    "                           <s> = <s1> <p1> <s2> <p2> <s3> <p3>\n"
  911    "                         where <s1>, <s2>, and <s3> are strings to use for\n"
  912    "                         the s1-, s2-, and s3-labels, respectively, and\n");
  913      fprintf(stdout,
  914    "                         where the strings <p1>, <p2>, <p3> determine the\n"
  915    "                         position of respective label, relative the end\n"
  916    "                         point of the arrow of respective axis. The label\n"
  917    "                         position is determined by the following syntax:\n");
  918      fprintf(stdout,
  919    "                              lft    Left\n"
  920    "                              rt     Right\n"
  921    "                              top    Top\n"
  922    "                              bot    Bottom\n"
  923    "                              ulft   Upper left\n"
  924    "                              urt    Upper Right\n"
  925    "                              llft   Lower left\n"
  926    "                              lrt    Lower right\n");
  927      fprintf(stdout,
  928    "                         The label strings should be expressed in plain\n"
  929    "                         TeX [2] mathmode syntax.\n");
  930      fprintf(stdout,
  931    "                         Default: <s1> = $S_1$, <s2> = $S_2$, <s3> = $S_3$\n"
  932    "                         Imortant note:  No blank spaces are allowed in\n"
  933    "                         the strings.\n"
  934    "\n");
  935      fprintf(stdout,
  936    " --draw_axes_inside      Toggles drawing (with dashed lines) of coordinate\n"
  937    "                         axes inside Poincare sphere.  Default: off.\n"
  938    "\n");
  939      fprintf(stdout,
  940    " -n, --normalize         Instead of making a trajectory plot of the para-\n"
  941    "                         meters (s1,s2,s3), contained in the file speci-\n"
  942    "                         fied by the '-f' option, instead use the norma-\n");
  943      fprintf(stdout,
  944    "                         lized parameter (s1/s0,s2/s0,s3/s0), which for\n"
  945    "                         completely polarized light corresponds to a tra-\n"
  946    "                         jectory mapped directly on the Poincare sphere,\n"
  947    "                         without any deviations fromthe spheres surface.\n");
  948      fprintf(stdout,
  949    "                         This option is particularly useful when only the\n"
  950    "                         state of polarization (and not the intensity) of\n"
  951    "                         the light is of interest.\n"
  952    "\n");
  953      fprintf(stdout,
  954    " --arrow <pa> <pb> <v>   Display an arrow, in Stokes parameter space, from\n"
  955    "                         point <pa>, at the command-line specified as the\n"
  956    "                         triple of floats <s1a> <s2a> <s3a>, to the point\n"
  957    "                         <pb>, similarly specified as <s1b> <s2b> <s3b>.\n");
  958      fprintf(stdout,
  959    "                         Useful for pointing out certain operation cycles\n"
  960    "                         in polarization domain, or just as an easy direct\n"
  961    "                         way of creating paths on the Poincare sphere\n"
  962    "                         without having to use external input files.\n");
  963      fprintf(stdout,
  964    "                         The arrow is drawn as a circular arc onto the\n"
  965    "                         Poincare sphere, through the closest path between\n"
  966    "                         the points.  The '--arrow' statement may appear\n"
  967    "                         repeated times,for producing multiple arrows.\n"
  968    "                         Currently there is a limit of 24 arrows in one\n"
  969    "                         single Poincare map (which should do for most\n"
  970    "                         people).\n");
  971      fprintf(stdout,
  972    "                           The last argument <v> is a pair of float values\n"
  973    "                         which determines the style of the drawn arrow.\n"
  974    "                         The pair <v> should be specified as <v1> <v2> on\n"
  975    "                         the command-line.\n");
  976      fprintf(stdout,
  977    "                         The first parameter <v1> determines the line-type\n"
  978    "                         of the arrow to draw.  The rules are:\n"
  979    "                            -0.5 <= <v1> < 0.5   -   Solid line\n"
  980    "                             0.5 <= <v1> < 1.5   -   Dashed line\n"
  981    "                         The second parameter, <v2>, determines the black-\n"
  982    "                         ness of the arrow to draw, where <v2> == 0 corre-\n"
  983    "                         sponds to white and <v2> == 1 to black.\n"
  984    "\n");
  985      fprintf(stdout,
  986    "Suffix conventions of the files:\n\n"
  987    "    .mp    - MetaPost source code (ASCII) [1]\n"
  988    "    .tex   - TeX source code (ASCII) [2]\n"
  989    "    .dvi   - Device independent output file from TeX [2]\n"
  990    "    .ps    - PostScript [3]\n"
  991    "    .eps   - Encapsulated PostScript [3]\n\n");
  992      fprintf(stdout,"References\n\n");
  993      fprintf(stdout,
  994    " [1] For information on the MetaPost program for typesetting figures,\n"
  995    "     see for example John Hobbys page, at\n"
  996    "     http://cm.bell-labs.com/who/hobby/MetaPost.html.\n\n");
  997      fprintf(stdout,
  998    " [2] For information on the TeX typesetting system, as well as references\n"
  999    "     to the dvips program, see for example the homepage of the TeX Users\n"
 1000    "     Group, at http://www.tug.org.\n\n");
 1001      fprintf(stdout,
 1002    " [3] For information on the PostScript programming language, see for\n"
 1003    "     example the homepage of Adobe Systems Inc., at\n"
 1004    "     http://www.adobe.com/products/postscript/main.html,\n"
 1005    "     or 'PostScript Language - Tutorial and Cookbook' (Addison-Wesley,\n"
 1006    "     Reading, Massachusetts, 1985), ISBN 0-201-10179-3.\n\n");
 1007      fprintf(stdout,
 1008    "Please report bugs to Fredrik Jonsson <fredrik@physics.kth.se>\n"
 1009    "Copyright (C) 1997-2005, Fredrik Jonsson <fredrik@physics.kth.se>\n");
 1010   } /* end of showsomehelp() */
 1011   
 1012   void initialize_variables(pmap *map) {
 1013      (*map).arrows=NULL;
 1014      (*map).numarrows=0;
 1015      (*map).save_memory=0;
 1016      (*map).use_normalized_stokes_params=0;
 1017      (*map).use_bezier_curves=0;
 1018      (*map).user_specified_inputfile=0;
 1019      (*map).user_specified_auxfile=0;
 1020      (*map).user_specified_axislabels=0;
 1021      (*map).user_specified_additional_coordinate_system=0;
 1022      (*map).user_specified_xtra_axislabel_x=0;
 1023      (*map).user_specified_xtra_axislabel_y=0;
 1024      (*map).user_specified_xtra_axislabel_z=0;
 1025      (*map).draw_hidden_dashed=0;
 1026      (*map).draw_paths_as_arrows=0;
 1027      (*map).reverse_arrow_paths=0;
 1028      (*map).last_point_infront=1;
 1029      (*map).current_point_is_a_beginlabelpoint=0;
 1030      (*map).current_point_is_an_endlabelpoint=0;
 1031      (*map).draw_axes_inside_sphere=0;
 1032      (*map).currently_drawing_path=0;
 1033      (*map).generate_eps_output=0;
 1034      (*map).xtra_neg_axis_length_x=DEFAULT_NEGATIVE_AXIS_LENGTH;
 1035      (*map).xtra_neg_axis_length_y=DEFAULT_NEGATIVE_AXIS_LENGTH;
 1036      (*map).xtra_neg_axis_length_z=DEFAULT_NEGATIVE_AXIS_LENGTH;
 1037      (*map).xtra_pos_axis_length_x=DEFAULT_POSITIVE_AXIS_LENGTH;
 1038      (*map).xtra_pos_axis_length_y=DEFAULT_POSITIVE_AXIS_LENGTH;
 1039      (*map).xtra_pos_axis_length_z=DEFAULT_POSITIVE_AXIS_LENGTH;
 1040      (*map).scalefactor=6.0;
 1041      (*map).rot_psi=DEFAULT_ROT_PSI;
 1042      (*map).rot_phi=DEFAULT_ROT_PHI;
 1043      (*map).delta_rot_psi=0.0;
 1044      (*map).delta_rot_phi=0;
 1045      (*map).phi_source=DEFAULT_PHI_SOURCE;
 1046      (*map).theta_source=DEFAULT_THETA_SOURCE;
 1047      (*map).upper_whiteness_value=DEFAULT_MAX_WHITENESS;
 1048      (*map).lower_whiteness_value=DEFAULT_MIN_WHITENESS;
 1049      (*map).hiddengraytone=DEFAULT_HIDDEN_GRAYTONE;
 1050      (*map).rho_divisor=DEFAULT_RHO_DIVISOR;
 1051      (*map).phi_divisor=DEFAULT_PHI_DIVISOR;
 1052      (*map).xpos_beginpoint=0.0;
 1053      (*map).ypos_beginpoint=0.0;
 1054      (*map).xpos_endpoint=0.0;
 1055      (*map).ypos_endpoint=0.0;
 1056      (*map).neg_axis_length_s1=DEFAULT_NEGATIVE_AXIS_LENGTH;
 1057      (*map).neg_axis_length_s2=DEFAULT_NEGATIVE_AXIS_LENGTH;
 1058      (*map).neg_axis_length_s3=DEFAULT_NEGATIVE_AXIS_LENGTH;
 1059      (*map).pos_axis_length_s1=DEFAULT_POSITIVE_AXIS_LENGTH;
 1060      (*map).pos_axis_length_s2=DEFAULT_POSITIVE_AXIS_LENGTH;
 1061      (*map).pos_axis_length_s3=DEFAULT_POSITIVE_AXIS_LENGTH;
 1062      (*map).paththickness=DEFAULT_PATH_THICKNESS;
 1063      (*map).arrowthickness=DEFAULT_ARROW_THICKNESS;
 1064      (*map).arrowheadangle=DEFAULT_ARROW_HEADANGLE;
 1065      (*map).coordaxisthickness=DEFAULT_ARROW_THICKNESS;
 1066      (*map).ticksize=DEFAULT_TICKSIZE;
 1067      strcpy((*map).outfilename,DEFAULT_OUTFILENAME);
 1068      strcpy((*map).epsjobname,DEFAULT_EPSJOBNAME);
 1069      strcpy((*map).axislabel_s1,DEFAULT_AXISLABEL_S1);
 1070      strcpy((*map).axislabel_s2,DEFAULT_AXISLABEL_S2);
 1071      strcpy((*map).axislabel_s3,DEFAULT_AXISLABEL_S3);
 1072      strcpy((*map).axislabelposition_s1,DEFAULT_AXISLABELPOSITION_S1);
 1073      strcpy((*map).axislabelposition_s2,DEFAULT_AXISLABELPOSITION_S2);
 1074      strcpy((*map).axislabelposition_s3,DEFAULT_AXISLABELPOSITION_S3);
 1075      strcpy((*map).xtra_axislabel_x,"");
 1076      strcpy((*map).xtra_axislabel_y,"");
 1077      strcpy((*map).xtra_axislabel_z,"");
 1078      strcpy((*map).labelstr_beginpoint,"");
 1079      strcpy((*map).labelstr_endpoint,"");
 1080   } /* end of initialize_variables() */
 1081   
 1082   /*
 1083    * In the initialization of the |stoketraject| struct, MAX_NUM_LABELS is the
 1084    * maximum number of allowed labels along each trajectory.
 1085    * However, we reserve two elements of the arrays containing the label data
 1086    * for possibly present labels at begin and end points of the trajectory,
 1087    * hence the total number of elements in the label-related arrays are
 1088    * MAX_NUM_LABELS+2.
 1089    */
 1090   void initialize_stoke_trajectory(stoketraject *tr) {
 1091      (*tr).numcoords=0;
 1092      (*tr).s1=dvector(1,MAX_NUM_STOKE_COORDS);
 1093      (*tr).s2=dvector(1,MAX_NUM_STOKE_COORDS);
 1094      (*tr).s3=dvector(1,MAX_NUM_STOKE_COORDS);
 1095      (*tr).visible=svector(1,MAX_NUM_STOKE_COORDS);
 1096      (*tr).numtickmarks=0;
 1097      (*tr).tickmark=lvector(1,MAX_NUM_TICKMARKS);
 1098      (*tr).numlabels=0;
 1099      (*tr).label=lvector(1,(MAX_NUM_LABELS+2));
 1100      (*tr).labeltext=cmatrix(1,(MAX_NUM_LABELS+2),1,MAX_LABEL_TEXTLENGTH);
 1101      (*tr).labellength=ivector(1,(MAX_NUM_LABELS+2));
 1102      (*tr).labelpos=svector(1,(MAX_NUM_LABELS+2));
 1103   } /* end of initialize_stoke_trajectory() */
 1104   
 1105   void reset_stokes_trajectory_struct(stoketraject *st) {
 1106      long int j,k;
 1107      (*st).numcoords=0;
 1108      for (k=1;k<=MAX_NUM_STOKE_COORDS;k++) {
 1109         (*st).s1[k]=0.0;
 1110         (*st).s2[k]=0.0;
 1111         (*st).s3[k]=0.0;
 1112         (*st).visible[k]=0;
 1113      }
 1114      (*st).numtickmarks=0;
 1115      for (k=1;k<=MAX_NUM_TICKMARKS;k++) (*st).tickmark[k]=0;
 1116      (*st).numlabels=0;
 1117      for (k=1;k<=(MAX_NUM_LABELS+2);k++) {
 1118         (*st).label[k]=0;
 1119         for (j=1;j<=MAX_LABEL_TEXTLENGTH;j++) (*st).labeltext[k][j]=' ';
 1120         (*st).labellength[k]=0;
 1121         (*st).labelpos[k]=0;
 1122      }
 1123   }
 1124   
 1125   /*---------------------------------------------------------------------
 1126   | The pathcharacter() routine takes one character `ch` as argument,
 1127   | and returns 1 ('true') if the character is valid character of a path
 1128   | string, otherwise 0 ('false') is returned.
 1129   ---------------------------------------------------------------------*/
 1130   short pathcharacter(int ch) {
 1131      return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
 1132         (ch=='-')||(ch=='+'));
 1133   }
 1134   
 1135   /*---------------------------------------------------------------------
 1136   | The strip_away_path() routine takes a character string `filename` as
 1137   | argument, and returns a pointer to the same string but without any
 1138   | preceding path segments. This routine is, for example, useful for
 1139   | removing paths from program names parsed from the command line.
 1140   ---------------------------------------------------------------------*/
 1141   char *strip_away_path(char filename[]) {
 1142      int j,k=0;
 1143      while (pathcharacter(filename[k])) k++;
 1144      j=(--k); /* this is the uppermost index of the full path+file string */
 1145      while (isalnum((int)(filename[j]))) j--;
 1146      j++; /* this is the lowermost index of the stripped file name */
 1147      return(&filename[j]);
 1148   }
 1149   
 1150   void display_parsed_command_line_option(pmap *map,char *optstr) {
 1151      if ((*map).verbose)
 1152         fprintf(stdout,"%s: Parsing '%s' option.\n",progname,optstr);
 1153   }
 1154   
 1155   /*-----------------------------------------------------------------------------
 1156   | Parsing for command-line arguments. This routine also takes care of default
 1157   | initialization of the variables of the pmap struct.
 1158   | For any additionally specified arrows, we allocate a [7x24] matrix for
 1159   | (potential) arrow coordinates. The organization is as follows:
 1160   |    Each column c=1,2,...,24, of the matrix corresponds to one arrow. The
 1161   |    first three elements are the (s1,s2,s3)-coordinates of the first point
 1162   |    of the arrow (the arrow base), the next three elements are the
 1163   |    (s1,s2,s3)-coordinates of the endpoint of the arrow. The seventh element
 1164   |    determines the type of arrow to draw: if  0 < arrows[...][7] < 1, the
 1165   |    arrow will be of solid-line type,  if  1 < arrows[...][7] < 2, the arrow
 1166   |    will be dashed.
 1167   |    Finally, the eighth element determines the blackness of the arrow, where
 1168   |    arrows[...][8] == 0 corresponds to white, and arrows[...][8] = 1 to
 1169   |    black.
 1170   -----------------------------------------------------------------------------*/
 1171   pmap parse_command_line(int argc, char *argv[]) {
 1172      int no_arg;
 1173      pmap map;
 1174   
 1175      initialize_variables(&map);
 1176      progname=strip_away_path(argv[0]); /* Strip from preceeding path string */
 1177      no_arg=argc;
 1178      map.arrows=matrix(1,8,1,24);
 1179      while(--argc) { /* while any arguments still are left */
 1180         if (!strcmp(argv[no_arg-argc],"-v") ||
 1181            !strcmp(argv[no_arg-argc],"--verbose")) {
 1182            map.verbose=(map.verbose?0:1);
 1183            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1184         } else if (!strcmp(argv[no_arg-argc],"--save_memory")) {
 1185            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1186            map.save_memory=(map.save_memory?0:1);
 1187         } else if (!strcmp(argv[no_arg-argc],"-n") ||
 1188                 !strcmp(argv[no_arg-argc],"--normalize")) {
 1189            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1190            map.use_normalized_stokes_params
 1191               =(map.use_normalized_stokes_params?0:1);
 1192         } else if (!strcmp(argv[no_arg-argc],"--paththickness")) {
 1193            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1194            --argc;
 1195            if (!sscanf(argv[no_arg-argc],"%lf",&map.paththickness)) {
 1196               fprintf(stderr,"%s: Couldn't get path thickness!\n",progname);
 1197               exit(FAILURE);
 1198            }
 1199         } else if (!strcmp(argv[no_arg-argc],"--draw_hidden_dashed")) {
 1200            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1201            map.draw_hidden_dashed=(map.draw_hidden_dashed?0:1);
 1202         } else if (!strcmp(argv[no_arg-argc],"--draw_paths_as_arrows")) {
 1203            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1204            map.draw_paths_as_arrows=(map.draw_paths_as_arrows?0:1);
 1205         } else if (!strcmp(argv[no_arg-argc],"--reverse_arrow_paths")) {
 1206            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1207            map.reverse_arrow_paths=(map.reverse_arrow_paths?0:1);
 1208         } else if (!strcmp(argv[no_arg-argc],"--arrowthickness")) {
 1209            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1210            --argc;
 1211            if (!sscanf(argv[no_arg-argc],"%lf",&map.arrowthickness)) {
 1212               fprintf(stderr,\
 1213                 "%s: Couldn't get arrow thickness in [pt]!\n",progname);
 1214               exit(FAILURE);
 1215            }
 1216         } else if (!strcmp(argv[no_arg-argc],"--arrowheadangle")) {
 1217            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1218            --argc;
 1219            if (!sscanf(argv[no_arg-argc],"%lf",&map.arrowheadangle)) {
 1220               fprintf(stderr,\
 1221                 "%s: Couldn't get arrow head angle in [deg]!\n",progname);
 1222               exit(FAILURE);
 1223            }
 1224         } else if (!strcmp(argv[no_arg-argc],"-b") ||
 1225                 !strcmp(argv[no_arg-argc],"--bezier")) {
 1226            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1227            map.use_bezier_curves=(map.use_bezier_curves?0:1);
 1228         } else if (!strcmp(argv[no_arg-argc],"-h") ||
 1229                 !strcmp(argv[no_arg-argc],"--help")) {
 1230            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1231            showsomehelp();
 1232            exit(FAILURE);
 1233         } else if (!strcmp(argv[no_arg-argc],"-V") ||
 1234                 !strcmp(argv[no_arg-argc],"--version")) {
 1235            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1236            show_banner();
 1237            exit(0);
 1238         } else if (strcmp(argv[no_arg-argc],"-f")==0 ||
 1239                strcmp(argv[no_arg-argc],"--inputfile")==0) {
 1240            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1241            --argc;
 1242            strcpy(map.infilename,argv[no_arg-argc]);
 1243            map.user_specified_inputfile=1;
 1244         } else if (strcmp(argv[no_arg-argc],"-e")==0 ||
 1245                strcmp(argv[no_arg-argc],"--epsoutput")==0) {
 1246            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1247            --argc;
 1248            strcpy(map.epsjobname,argv[no_arg-argc]);
 1249            map.generate_eps_output=1;
 1250         } else if (strcmp(argv[no_arg-argc],"-o")==0 ||
 1251                strcmp(argv[no_arg-argc],"--outputfile")==0) {
 1252            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1253            --argc;
 1254            strcpy(map.outfilename,argv[no_arg-argc]);
 1255         } else if (strcmp(argv[no_arg-argc],"--auxsource")==0) {
 1256            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1257            --argc;
 1258            strcpy(map.auxfilename,argv[no_arg-argc]);
 1259            map.user_specified_auxfile=1;
 1260         } else if (!strcmp(argv[no_arg-argc],"--psi") ||
 1261                 !strcmp(argv[no_arg-argc],"--rotatepsi")) {
 1262            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1263            --argc;
 1264            if (!sscanf(argv[no_arg-argc],"%lf",&map.rot_psi)) {
 1265               fprintf(stderr,\
 1266                 "%s: Couldn't get value for psi (rotation round z)!\n",progname);
 1267               exit(FAILURE);
 1268            }
 1269            map.rot_psi *= (PI/180);
 1270         } else if (!strcmp(argv[no_arg-argc],"--phi") ||
 1271                 !strcmp(argv[no_arg-argc],"--rotatephi")) {
 1272            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1273            --argc;
 1274            if (!sscanf(argv[no_arg-argc],"%lf",&map.rot_phi)) {
 1275               fprintf(stderr,\
 1276                 "%s: Couldn't get value for phi (rotation round y)!\n",progname);
 1277               exit(FAILURE);
 1278            }
 1279            map.rot_phi *= (PI/180);
 1280         } else if (!strcmp(argv[no_arg-argc],"--rhodivisor")) {
 1281            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1282            --argc;
 1283            if (!sscanf(argv[no_arg-argc],"%lf",&map.rho_divisor)) {
 1284               fprintf(stderr,\
 1285                 "%s: Couldn't get value for rho divisor!\n",progname);
 1286               exit(FAILURE);
 1287            }
 1288         } else if (!strcmp(argv[no_arg-argc],"--phidivisor")) {
 1289            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1290            --argc;
 1291            if (!sscanf(argv[no_arg-argc],"%lf",&map.phi_divisor)) {
 1292               fprintf(stderr,\
 1293                 "%s: Couldn't get value for phi divisor!\n",progname);
 1294               exit(FAILURE);
 1295            }
 1296         } else if (!strcmp(argv[no_arg-argc],"--scalefactor")) {
 1297            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1298            --argc;
 1299            if (!sscanf(argv[no_arg-argc],"%lf",&map.scalefactor)) {
 1300               fprintf(stderr,\
 1301                 "%s: Couldn't get value for scalefactor!\n",progname);
 1302               exit(FAILURE);
 1303            }
 1304         } else if (!strcmp(argv[no_arg-argc],"--hiddengraytone")) {
 1305            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1306            --argc;
 1307            if (!sscanf(argv[no_arg-argc],"%lf",&map.hiddengraytone)) {
 1308               fprintf(stderr,\
 1309                 "%s: Couldn't get whiteness value of hidden parts!\n",progname);
 1310               exit(FAILURE);
 1311            }
 1312         } else if (!strcmp(argv[no_arg-argc],"--shading")) {
 1313            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1314            --argc;
 1315            if (!sscanf(argv[no_arg-argc],"%lf",&map.lower_whiteness_value)) {
 1316               fprintf(stderr,\
 1317                 "%s: Couldn't get lower value of sphere whiteness!\n",progname);
 1318               exit(FAILURE);
 1319            }
 1320            --argc;
 1321            if (!sscanf(argv[no_arg-argc],"%lf",&map.upper_whiteness_value)) {
 1322               fprintf(stderr,\
 1323                 "%s: Couldn't get upper value of sphere whiteness!\n",progname);
 1324               exit(FAILURE);
 1325            }
 1326         } else if (!strcmp(argv[no_arg-argc],"--axislengths")) {
 1327            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1328            --argc;
 1329            if (!sscanf(argv[no_arg-argc],"%lf",&map.neg_axis_length_s1)) {
 1330               fprintf(stderr,"%s: Couldn't get minimum for s1 axis!\n",progname);
 1331               fprintf(stderr,"%s: Check the '--axislength' option\n", progname);
 1332               exit(FAILURE);
 1333            }
 1334            --argc;
 1335            if (!sscanf(argv[no_arg-argc],"%lf",&map.pos_axis_length_s1)) {
 1336               fprintf(stderr,"%s: Couldn't get maximum for s1 axis!\n",progname);
 1337               fprintf(stderr,"%s: Check the '--axislength' option\n", progname);
 1338               exit(FAILURE);
 1339            }
 1340            --argc;
 1341            if (!sscanf(argv[no_arg-argc],"%lf",&map.neg_axis_length_s2)) {
 1342               fprintf(stderr,"%s: Couldn't get minimum for s2 axis!\n",progname);
 1343               fprintf(stderr,"%s: Check the '--axislength' option\n", progname);
 1344               exit(FAILURE);
 1345            }
 1346            --argc;
 1347            if (!sscanf(argv[no_arg-argc],"%lf",&map.pos_axis_length_s2)) {
 1348               fprintf(stderr,"%s: Couldn't get maximum for s2 axis!\n",progname);
 1349               fprintf(stderr,"%s: Check the '--axislength' option\n", progname);
 1350               exit(FAILURE);
 1351            }
 1352            --argc;
 1353            if (!sscanf(argv[no_arg-argc],"%lf",&map.neg_axis_length_s3)) {
 1354               fprintf(stderr,"%s: Couldn't get minimum for s3 axis!\n",progname);
 1355               fprintf(stderr,"%s: Check the '--axislength' option\n", progname);
 1356               exit(FAILURE);
 1357            }
 1358            --argc;
 1359            if (!sscanf(argv[no_arg-argc],"%lf",&map.pos_axis_length_s3)) {
 1360               fprintf(stderr,"%s: Couldn't get maximum for s3 axis!\n",progname);
 1361               fprintf(stderr,"%s: Check the '--axislength' option\n", progname);
 1362               exit(FAILURE);
 1363            }
 1364         } else if (!strcmp(argv[no_arg-argc],"--axislabels")) {
 1365            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1366            --argc;
 1367            strcpy(map.axislabel_s1,argv[no_arg-argc]);
 1368            --argc;
 1369            strcpy(map.axislabelposition_s1,argv[no_arg-argc]);
 1370            --argc;
 1371            strcpy(map.axislabel_s2,argv[no_arg-argc]);
 1372            --argc;
 1373            strcpy(map.axislabelposition_s2,argv[no_arg-argc]);
 1374            --argc;
 1375            strcpy(map.axislabel_s3,argv[no_arg-argc]);
 1376            --argc;
 1377            strcpy(map.axislabelposition_s3,argv[no_arg-argc]);
 1378            map.user_specified_axislabels=1;
 1379         } else if (!strcmp(argv[no_arg-argc],"--draw_axes_inside")) {
 1380            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1381            map.draw_axes_inside_sphere=(map.draw_axes_inside_sphere?0:1);
 1382         } else if (!strcmp(argv[no_arg-argc],"--xtracoordsys")) {
 1383            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1384            --argc;
 1385            if (!sscanf(argv[no_arg-argc],"%lf",&map.delta_rot_psi)) {
 1386               fprintf(stderr,"%s: Couldn't get delta_rot_psi!\n",progname);
 1387               fprintf(stderr,"%s: Check the '--xtracoordsys' option\n",progname);
 1388               fprintf(stderr,"%s: (first argument)\n",progname);
 1389               exit(FAILURE);
 1390            }
 1391            map.delta_rot_psi *= (PI/180);
 1392            --argc;
 1393            if (!sscanf(argv[no_arg-argc],"%lf",&map.delta_rot_phi)) {
 1394               fprintf(stderr,"%s: Couldn't get delta_rot_phi!\n",progname);
 1395               fprintf(stderr,"%s: Check the '--xtracoordsys' option\n",progname);
 1396               fprintf(stderr,"%s: (second argument)\n",progname);
 1397               exit(FAILURE);
 1398            }
 1399            map.delta_rot_phi *= (PI/180);
 1400            map.user_specified_additional_coordinate_system=1;
 1401         } else if (!strcmp(argv[no_arg-argc],"--xtracoordsys_axislabel_x")) {
 1402            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1403            --argc;
 1404            strcpy(map.xtra_axislabel_x,argv[no_arg-argc]);
 1405            map.user_specified_xtra_axislabel_x=1;
 1406         } else if (!strcmp(argv[no_arg-argc],"--xtracoordsys_axislabel_y")) {
 1407            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1408            --argc;
 1409            strcpy(map.xtra_axislabel_y,argv[no_arg-argc]);
 1410            map.user_specified_xtra_axislabel_y=1;
 1411         } else if (!strcmp(argv[no_arg-argc],"--xtracoordsys_axislabel_z")) {
 1412            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1413            --argc;
 1414            strcpy(map.xtra_axislabel_z,argv[no_arg-argc]);
 1415            map.user_specified_xtra_axislabel_z=1;
 1416         } else if (!strcmp(argv[no_arg-argc],"--xtracoordsys_axislengths")) {
 1417            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1418            --argc;
 1419            if (!sscanf(argv[no_arg-argc],"%lf",&map.xtra_neg_axis_length_x)) {
 1420               fprintf(stderr,"%s: Couldn't get minimum for x axis!\n",progname);
 1421               fprintf(stderr,
 1422                  "%s: Check the '--xtracoordsys_axislengths' option", progname);
 1423               exit(FAILURE);
 1424            }
 1425            --argc;
 1426            if (!sscanf(argv[no_arg-argc],"%lf",&map.xtra_pos_axis_length_x)) {
 1427               fprintf(stderr,"%s: Couldn't get maximum for x axis!\n",progname);
 1428               fprintf(stderr,
 1429                  "%s: Check the '--xtracoordsys_axislengths' option", progname);
 1430               exit(FAILURE);
 1431            }
 1432            --argc;
 1433            if (!sscanf(argv[no_arg-argc],"%lf",&map.xtra_neg_axis_length_y)) {
 1434               fprintf(stderr,"%s: Couldn't get minimum for y axis!\n",progname);
 1435               fprintf(stderr,
 1436                  "%s: Check the '--xtracoordsys_axislengths' option", progname);
 1437               exit(FAILURE);
 1438            }
 1439            --argc;
 1440            if (!sscanf(argv[no_arg-argc],"%lf",&map.xtra_pos_axis_length_y)) {
 1441               fprintf(stderr,"%s: Couldn't get maximum for y axis!\n",progname);
 1442               fprintf(stderr,
 1443                  "%s: Check the '--xtracoordsys_axislengths' option", progname);
 1444               exit(FAILURE);
 1445            }
 1446            --argc;
 1447            if (!sscanf(argv[no_arg-argc],"%lf",&map.xtra_neg_axis_length_z)) {
 1448               fprintf(stderr,"%s: Couldn't get minimum for z axis!\n",progname);
 1449               fprintf(stderr,
 1450                  "%s: Check the '--xtracoordsys_axislengths' option", progname);
 1451               exit(FAILURE);
 1452            }
 1453            --argc;
 1454            if (!sscanf(argv[no_arg-argc],"%lf",&map.xtra_pos_axis_length_z)) {
 1455               fprintf(stderr,"%s: Couldn't get maximum for z axis!\n",progname);
 1456               fprintf(stderr,
 1457                  "%s: Check the '--xtracoordsys_axislengths' option", progname);
 1458               exit(FAILURE);
 1459            }
 1460         } else if (!strcmp(argv[no_arg-argc],"--arrow")) {
 1461            display_parsed_command_line_option(&map,argv[no_arg-argc]);
 1462            map.numarrows++;
 1463            --argc;
 1464            if (!sscanf(argv[no_arg-argc],"%lf",&map.arrows[1][map.numarrows])) {
 1465               fprintf(stderr,\
 1466                 "%s: Couldn't get S1 coordinate for starting point of arrow"
 1467                 " No. %d!\n", progname, map.numarrows);
 1468               exit(FAILURE);
 1469            }
 1470            --argc;
 1471            if (!sscanf(argv[no_arg-argc],"%lf",&map.arrows[2][map.numarrows])) {
 1472               fprintf(stderr,\
 1473                 "%s: Couldn't get S2 coordinate for starting point of arrow"
 1474                 " No. %d!\n", progname, map.numarrows);
 1475               exit(FAILURE);
 1476            }
 1477            --argc;
 1478            if (!sscanf(argv[no_arg-argc],"%lf",&map.arrows[3][map.numarrows])) {
 1479               fprintf(stderr,\
 1480                 "%s: Couldn't get S3 coordinate for starting point of arrow"
 1481                 " No. %d!\n", progname, map.numarrows);
 1482               exit(FAILURE);
 1483            }
 1484            --argc;
 1485            if (!sscanf(argv[no_arg-argc],"%lf",&map.arrows[4][map.numarrows])) {
 1486               fprintf(stderr,\
 1487                 "%s: Couldn't get S1 coordinate for ending point of arrow"
 1488                 " No. %d!\n", progname, map.numarrows);
 1489               exit(FAILURE);
 1490            }
 1491            --argc;
 1492            if (!sscanf(argv[no_arg-argc],"%lf",&map.arrows[5][map.numarrows])) {
 1493               fprintf(stderr,\
 1494                 "%s: Couldn't get S2 coordinate for ending point of arrow"
 1495                 " No. %d!\n", progname, map.numarrows);
 1496               exit(FAILURE);
 1497            }
 1498            --argc;
 1499            if (!sscanf(argv[no_arg-argc],"%lf",&map.arrows[6][map.numarrows])) {
 1500               fprintf(stderr,\
 1501                 "%s: Couldn't get S3 coordinate for ending point of arrow"
 1502                 " No. %d!\n", progname, map.numarrows);
 1503               exit(FAILURE);
 1504            }
 1505            --argc;
 1506            if (!sscanf(argv[no_arg-argc],"%lf",&map.arrows[7][map.numarrows])) {
 1507               fprintf(stderr,\
 1508                 "%s: Couldn't get line type of arrow No. %d!\n",\
 1509                                                        progname, map.numarrows);
 1510               exit(FAILURE);
 1511            }
 1512            --argc;
 1513            if (!sscanf(argv[no_arg-argc],"%lf",&map.arrows[8][map.numarrows])) {
 1514               fprintf(stderr,\
 1515                 "%s: Couldn't get blackness of arrow No. %d!\n",\
 1516                                                        progname, map.numarrows);
 1517               exit(FAILURE);
 1518            }
 1519         } else { /* if no valid option is to be found... */
 1520            fprintf(stderr,"%s: Error: Specified option '%s' invalid!\n",
 1521               progname,argv[no_arg-argc]);
 1522            showsomehelp();
 1523            exit(FAILURE);
 1524         }
 1525      }
 1526      return map; /* return all parameter values as a struct of type |pmap| */
 1527   } /* end of parse_command_line() */
 1528   
 1529   void display_arrow_specs(pmap map) {
 1530      int i;
 1531      if (map.numarrows > 0) {
 1532         if (map.verbose) {
 1533            printf("%s: You specified the following arrows to draw:\n",progname);
 1534            for (i=1;i<=map.numarrows;i++) {
 1535               printf("%s:    Arrow No. %d: (%1.2f,%1.2f,%1.2f) --> "
 1536                      "(%1.2f,%1.2f,%1.2f) [%1.2f,%1.2f]\n", progname, i,
 1537                            map.arrows[1][i], map.arrows[2][i], map.arrows[3][i],
 1538                            map.arrows[4][i], map.arrows[5][i], map.arrows[6][i],
 1539                            map.arrows[7][i], map.arrows[8][i]);
 1540            }
 1541         }
 1542      }
 1543   } /* end of display_arrow_specs() */
 1544   
 1545   FILE *open_outfile(pmap map) {
 1546      FILE *outfileptr;
 1547      if ((outfileptr=fopen(map.outfilename,"w")) == NULL) {
 1548         fprintf(stderr,"Couldn't open file %s for output!\n",map.outfilename);
 1549         exit(FAILURE);
 1550      } else {
 1551         if (map.verbose)
 1552            printf("%s: Writing MetaPost code to %s\n",progname,map.outfilename);
 1553      }
 1554      return outfileptr;
 1555   } /* end of open_outfile() */
 1556   
 1557   FILE *open_infile(pmap map) {
 1558      FILE *infileptr;
 1559      if (!map.user_specified_inputfile) {
 1560         fprintf(stderr,"%s: No input trajectory file specified.\n",progname);
 1561         infileptr=NULL;
 1562      } else {
 1563      if ((infileptr=fopen(map.infilename,"r")) == NULL) {
 1564         fprintf(stderr,"%s: Couldn't open trajectory file %s for reading\n",
 1565            progname,map.infilename);
 1566         fprintf(stderr,"%s: Please check -f or --inputfile option arguments\n",
 1567                            progname);
 1568         exit(FAILURE);
 1569      } else {
 1570         if (map.verbose) {
 1571            printf("%s: Reading Stokes parameters from %s\n",\
 1572                                              progname, map.infilename);
 1573         }
 1574      }
 1575      }
 1576      return infileptr;
 1577   } /* end of open_infile() */
 1578   
 1579   /*
 1580    * Write out the heading comments of the MetaPost file that is about
 1581    * to be generated.
 1582    */
 1583   void write_header(FILE *outfileptr,pmap map,int argc, char *argv[]) {
 1584      time_t now=time(NULL);
 1585      int i=0,no_arg=argc;
 1586      fprintf(outfileptr,
 1587    "%% This Filename:  %s   [MetaPost source]\n"
 1588    "%% Creation time:  %s"                   /* ctime() takes care of '\n'. */
 1589    "%%\n"
 1590    "%% Copyright (C) 1997-2005, Fredrik Jonsson <fj@optics.kth.se>\n"
 1591    "%%\n"
 1592    "%% Input Filename [Stokes parameters]:  %s\n"
 1593    "%% This MetaPost source code was automatically generated by %s\n",
 1594      map.outfilename, ctime(&now), map.infilename, progname);
 1595   
 1596      /*
 1597       * In many cases, it is quite useful to keep the full command-line
 1598       * that generated a certain figure, for future reference.
 1599       */
 1600      argc=no_arg-1;
 1601      fprintf(outfileptr,
 1602         "%% Full set of command line options that generated this code:\n");
 1603      while (argc) {
 1604         if (i==0) fprintf(outfileptr,"%%    ");
 1605         fprintf(outfileptr," %s",argv[no_arg-argc]);
 1606         argc--;
 1607         i++;
 1608         if (i==6) {
 1609            fprintf(outfileptr,"\n");
 1610            i=0;
 1611         }
 1612      }
 1613      fprintf(outfileptr,"\n%%\n");
 1614      fprintf(outfileptr,
 1615    "%% Description:  Map of Stokes parameters, visualized as trajectories\n"
 1616    "%%               onto the Poincare sphere. This file contains MetaPost\n"
 1617    "%%               source code, to be compiled with John Hobby's MetaPost\n"
 1618    "%%               compiler or used with anything that understands MetaPost\n"
 1619    "%%               source code.\n"
 1620    "%%\n");
 1621      fprintf(outfileptr,
 1622    "%% If you want to create PostScript output, or include the resulting\n"
 1623    "%% output in a TeX document, this example illustrates the procedure,\n"
 1624    "%% assuming 'poincaremap.mp' to be the name of the file containing the\n"
 1625    "%% MetaPost code to be visualized: (commands run on command-line)\n"
 1626    "%%\n");
 1627      fprintf(outfileptr,
 1628    "%%       mp poincaremap.mp;\n"
 1629    "%%       echo \042\\input epsf\\centerline{\\epsfbox{poincaremap.1}}\\bye\042 > tmp.tex;\n"
 1630    "%%       tex tmp.tex;\n"
 1631    "%%       dvips tmp.dvi -o poincaremap.ps;\n"
 1632    "%%\n");
 1633      fprintf(outfileptr,
 1634    "%% Here, the first command compiles the MetaPost source code, and leaves\n"
 1635    "%% an Encapsulated PostScript file named 'poincaremap.1', containing TeX\n"
 1636    "%% control codes for characters, etc. This file does not contain any\n"
 1637    "%% definitions for characters or TeX-specific items, and it cannot be\n"
 1638    "%% viewed or printed simply as is stands; it must rather be included into\n"
 1639    "%% TeX code in order to provide something useful.\n");
 1640      fprintf(outfileptr,
 1641    "%%     The second command creates a temporary minimal TeX-file 'tmp.tex',\n"
 1642    "%% that only includes the previously generated Encapsulated PostScript\n"
 1643    "%% code.\n");
 1644      fprintf(outfileptr,
 1645    "%%     The third command compiles the TeX-code into device-independent,\n"
 1646    "%% or DVI, output, stored in the file 'tmp.dvi'.\n"
 1647    "%%     Finally, the last command converts the DVI output into a free-\n"
 1648    "%% standing PostScript file 'poincaremap.ps', to be printed or viewed\n"
 1649    "%% with some PostScript viewer, such as GhostView.\n"
 1650    "%%\n");
 1651   } /* end of write_header() */
 1652   
 1653   void write_euler_angle_specs(FILE *outfileptr,pmap map) {
 1654      fprintf(outfileptr,
 1655    "scalefactor := %f mm;\n"
 1656    "rot_psi := %f;  %% Rotation angle round z-axis (first rotation)\n"
 1657    "rot_phi := %f;  %% Rotation angle round y-axis (second rotation)\n"
 1658    "alpha := %f;    %% == arctan(sin(rot_phi)*tan(rot_psi))\n"
 1659    "beta  := %f;    %% == arctan(sin(rot_phi)/tan(rot_psi))\n\n",
 1660      map.scalefactor, (180/PI)*map.rot_psi, (180/PI)*map.rot_phi,
 1661      (180/PI)*atan(sin(map.rot_phi)*tan(map.rot_psi)),
 1662      (180/PI)*atan(sin(map.rot_phi)/tan(map.rot_psi)));
 1663   } /* end of write_euler_angle_specs() */
 1664   
 1665   void write_sphere_shading_specs(FILE *outfileptr,pmap map) {
 1666      /*
 1667       * Parameters specifying the location of the light source.
 1668       */
 1669      fprintf(outfileptr,
 1670        "%%\n"
 1671        "%% Parameters specifying the location of the light source; for Phong\n"
 1672        "%% shading of the sphere.\n"
 1673        "%%\n"
 1674        "%%    phi_source:  Angle (in deg.) to light source counterclockwise\n"
 1675        "%%                 'from three o'clock', viewed from the observer.\n"
 1676        "%%\n"
 1677        "%%  theta_source:  Angle (in deg.) between light source and observer,\n"
 1678        "%%                 seen from the centre of the sphere.\n"
 1679        "%%\n");
 1680      fprintf(outfileptr,
 1681        "%% Parameters specifying the shading 'intensity' in terms of maximum\n"
 1682        "%% (for the highlighs) and minimum (for the deep shadowed regions)\n"
 1683        "%% values for the Phong shading.  '0.0' <=> 'black'; '1.0' <=> 'white'\n"
 1684        "%%\n"
 1685        "%%   upper_value:  Maximum value of whiteness.\n"
 1686        "%%   lower_value:  Minimum value of whiteness.\n"
 1687        "%%\n");
 1688      fprintf(outfileptr,
 1689        "phi_source := %f;\n"
 1690        "theta_source := %f;\n"
 1691        "upper_value := %f;\n"
 1692        "lower_value := %f;\n",
 1693          (180/PI)*map.phi_source,(180/PI)*map.theta_source,
 1694          map.upper_whiteness_value,map.lower_whiteness_value);
 1695      fprintf(outfileptr,
 1696        "radius := scalefactor;\n"
 1697        "delta_rho := radius/%f;\n"
 1698        "delta_phi := 360.0/%f;\n"
 1699        "beginfig(1);\n"
 1700        "  path p;\n"
 1701        "  path equator;\n"
 1702        "  transform T;\n"
 1703        "  c1:=lower_value;\n"
 1704        "  c2:=upper_value-lower_value;\n",map.rho_divisor, map.phi_divisor);
 1705   
 1706   /*-----------------------------------------------------------------------------
 1707   | Here follows the x-, y- and z-components of the unit normal vector pointing
 1708   | from the center of the sphere to the position of the (point-like) light
 1709   | source.
 1710   -----------------------------------------------------------------------------*/
 1711      fprintf(outfileptr,
 1712        "  nx_source := sind(theta_source)*cosd(phi_source);\n"
 1713        "  ny_source := sind(theta_source)*sind(phi_source);\n"
 1714        "  nz_source := cosd(theta_source);\n"
 1715        "  phistop := 360.0;\n"
 1716        "  rhostop := radius - delta_rho/2.0;\n");
 1717   } /* end of write_sphere_shading_specs() */
 1718   
 1719   /*-----------------------------------------------------------------------------
 1720   | Now generate the Phong-shaded poincare sphere as it is projected down
 1721   | in 2D.  This is done by creating a number of >>parallelltrapetser<<
 1722   | mapped centro-symmmetrically around origo.
 1723   | The steps involved in this are as follows:
 1724   | 1. Take the coordinate of the center of each >>parallelltrapets<< to be
 1725   |    the reference for calculating the coordinates of the corners.
 1726   | 2. Calculate the coordinates of the corners of >>parallelltrapetserna<<.
 1727   | 3. Create the path p of the >>parallelltrapets<<:
 1728   |
 1729   |          (x3,y3)                       (x2,y2)
 1730   |               +------------<------------+
 1731   |                \                       /
 1732   |                 \                     /
 1733   |                  v                   ^
 1734   |                   \                 /
 1735   |                    \               /
 1736   |                     +------>------+
 1737   |                 (x4,y4)          (x1,y1)
 1738   |
 1739   | 4. Check if the >>parallelltrapets<< just created is shaded by the sphere or
 1740   |    not, i.e. if the position of the >>parallelltrapets<< is situated on the
 1741   |    opposite side of the Poincare sphere, seen from the direction of the light
 1742   |    source. If so, make the >>parallelltrapets<< shaded with the darkest tone,
 1743   |    specified by c1; otherwise make it shaded with the tone of gray specified
 1744   |    by prod, the scalar product of the normal vector to the object with the
 1745   |    normal vector to the lightsource.
 1746   -----------------------------------------------------------------------------*/
 1747   void write_shaded_sphere(FILE *outfileptr,pmap map) {
 1748      fprintf(outfileptr,
 1749        "%%\n"
 1750        "%% Draw the shaded Poincare sphere projected on 2D screen coordinates\n"
 1751        "%%\n");
 1752      fprintf(outfileptr,
 1753        "  for rho=0.0cm step delta_rho until rhostop:\n"
 1754        "    for phi=0.0 step delta_phi until phistop:\n");
 1755      fprintf(outfileptr,
 1756        "      rhomid := rho + delta_rho/2.0;\n"
 1757        "      phimid := phi + delta_phi/2.0;\n");
 1758      fprintf(outfileptr,
 1759        "      x1 := rho*cosd(phi);\n"
 1760        "      y1 := rho*sind(phi);\n"
 1761        "      x2 := (rho+delta_rho)*cosd(phi);\n"
 1762        "      y2 := (rho+delta_rho)*sind(phi);\n"
 1763        "      x3 := (rho+delta_rho)*cosd(phi+delta_phi);\n"
 1764        "      y3 := (rho+delta_rho)*sind(phi+delta_phi);\n"
 1765        "      x4 := rho*cosd(phi+delta_phi);\n"
 1766        "      y4 := rho*sind(phi+delta_phi);\n");
 1767      fprintf(outfileptr,
 1768        "      p:=makepath makepen ((x1,y1)--(x2,y2)--(x3,y3)--(x4,y4)--(x1,y1));\n"
 1769        "      quot := (rhomid/radius);\n"
 1770        "      nx_object := quot*cosd(phimid);\n"
 1771        "      ny_object := quot*sind(phimid);\n"
 1772        "      nz_object := sqrt(1-quot*quot);\n");
 1773      fprintf(outfileptr,
 1774        "      prod:=nx_object*nx_source+ny_object*ny_source\n"
 1775        "            +nz_object*nz_source;\n");
 1776      fprintf(outfileptr,
 1777        "      if prod < 0.0:\n"
 1778        "         value := c1;\n"
 1779        "      else:\n"
 1780        "         value := c1 + c2*prod*prod;\n"
 1781        "      fi\n"
 1782        "      fill p withcolor value[black,white];\n"
 1783        "    endfor\n"
 1784        "  endfor\n\n");
 1785   } /*end of write_shaded_sphere() */
 1786   
 1787   /*-----------------------------------------------------------------------------
 1788   | Draw the 'equators' $S_1=0$, $S_2=0$ and $S_3=0$ on the sphere.
 1789   | If specified by the user, also draw the additional 'equators' (corresponding
 1790   | to a coordinate system rotated around the original one by the Euler angles
 1791   | delta_psi, delta_phi) on the sphere.
 1792   -----------------------------------------------------------------------------*/
 1793   void write_equators(FILE *outfileptr,pmap map) {
 1794      fprintf(outfileptr,
 1795        "%%\n"
 1796        "%% Draw the 'equators' of the Poincare sphere\n"
 1797        "%%\n"
 1798        "   equator := halfcircle scaled (2.0*radius);\n"
 1799        "   eqcolval := .45;    %% '0.0' <=> 'white';  '1.0' <=> 'black'\n"
 1800        "\n");
 1801      fprintf(outfileptr,
 1802        "   pickup pencircle scaled %f pt;\n",map.coordaxisthickness);
 1803      fprintf(outfileptr,
 1804        "%%\n"
 1805        "%% Draw equator $S_3=0$...\n"
 1806        "%%\n"
 1807        "   T := identity yscaled sind(rot_phi) rotated 180.0;\n"
 1808        "   draw equator transformed T withcolor eqcolval [white,black];\n"
 1809        "\n");
 1810      fprintf(outfileptr,
 1811        "%%\n"
 1812        "%% ... then equator $S_2=0$...\n"
 1813        "%%\n"
 1814        "   T := identity yscaled (cosd(rot_phi)*sind(rot_psi))\n"
 1815        "                 rotated (270.0 + alpha);\n"
 1816        "   draw equator transformed T withcolor eqcolval [white,black];\n"
 1817        "\n");
 1818      fprintf(outfileptr,
 1819        "%%\n"
 1820        "%% ... and finally equator $S_1=0$.\n"
 1821        "%%\n"
 1822        "   T := identity yscaled (cosd(rot_phi)*cosd(rot_psi))\n"
 1823        "                 rotated (270.0 - beta);\n"
 1824        "   draw equator transformed T withcolor eqcolval [white,black];\n"
 1825        "\n");
 1826      if (map.user_specified_additional_coordinate_system) { /* whoaaaou ... */
 1827         fprintf(outfileptr,
 1828           "%%\n"
 1829           "%% Some handy parameters used in calculations below.\n"
 1830           "%%\n");
 1831         fprintf(outfileptr,
 1832           "delta_rot_psi := %f; %% Additional 1st rotation angle round z-axis\n"
 1833           "delta_rot_phi := %f;  %% Additional 2nd rotation angle round y-axis\n"
 1834           "delta_alpha := %f;    %% == arctan(sin(rot_phi)*tan(rot_psi))\n"
 1835           "delta_beta  := %f;    %% == arctan(sin(rot_phi)/tan(rot_psi))\n\n",
 1836             (180/PI)*(map.delta_rot_psi),(180/PI)*(map.delta_rot_phi),
 1837             (180/PI)*atan(sin(map.rot_phi + map.delta_rot_phi)
 1838               *tan(map.rot_psi + map.delta_rot_psi)),
 1839             (180/PI)*atan(sin(map.rot_phi + map.delta_rot_phi)
 1840               /tan(map.rot_psi + map.delta_rot_psi)));
 1841         fprintf(outfileptr,
 1842           "%%\n"
 1843           "%% Draw the additional 'equators' of the Poincare sphere,\n"
 1844           "%% corresponding to a system rotated by the Euler-angles\n"
 1845           "%%   delta_psi=, delta_phi=\n"
 1846           "%%\n"
 1847           "   equator := halfcircle scaled (2.0*radius);\n"
 1848           "   eqcolval := .45;    %% '0.0' <=> 'white';  '1.0' <=> 'black'\n"
 1849           "\n");
 1850         fprintf(outfileptr,
 1851           "%%\n"
 1852           "%% Draw equator $W_3=0$...\n"
 1853           "%%\n"
 1854           "   T := identity yscaled sind(rot_phi+delta_rot_phi) rotated 180.0;\n"
 1855           "   draw equator transformed T withcolor eqcolval [white,black];\n"
 1856           "\n");
 1857         fprintf(outfileptr,
 1858           "%%\n"
 1859           "%% ... then equator $W_2=0$...\n"
 1860           "%%\n"
 1861           "   T := identity yscaled (cosd(rot_phi + delta_rot_phi)\n"
 1862           "             *sind(rot_psi + delta_rot_psi))\n"
 1863           "             rotated (270.0 + delta_alpha);\n"
 1864           "   draw equator transformed T withcolor eqcolval [white,black];\n"
 1865           "\n");
 1866         fprintf(outfileptr,
 1867           "%%\n"
 1868           "%% ... and finally equator $W_1=0$.\n"
 1869           "%%\n"
 1870           "   T := identity yscaled (cosd(rot_phi + delta_rot_phi)\n"
 1871           "             *cosd(rot_psi + delta_rot_psi))\n"
 1872           "             rotated (270.0 - delta_beta);\n"
 1873           "   draw equator transformed T withcolor eqcolval [white,black];\n"
 1874           "\n");
 1875      }
 1876   } /* end of write_equators() */
 1877   
 1878   void display_stokes_trajectory(stoketraject st) {
 1879      long k;
 1880   
 1881      fprintf(stdout,"Stokes trajectory [%ld coordinates]:\n",st.numcoords);
 1882      for (k=1;k<=st.numcoords;k++) {
 1883         fprintf(stdout," S[%ld]=(%f, %f ,%f), ",k,st.s1[k],st.s2[k],st.s3[k]);
 1884         fprintf(stdout,"[%s]\n",(st.visible[k]?"visible":"hidden"));
 1885      }
 1886   }
 1887   
 1888   short new_trajectory(FILE *infileptr) {
 1889      int ch;
 1890      if ((ch=getc(infileptr))=='p') {
 1891         return 1;
 1892      } else {
 1893         ungetc(ch,infileptr);
 1894         return 0;
 1895      }
 1896   }
 1897   
 1898   short end_of_trajectory(FILE *infileptr) {
 1899      int ch;
 1900      if ((ch=getc(infileptr))=='q') {
 1901         return 1;
 1902      } else {
 1903         ungetc(ch,infileptr);
 1904         return 0;
 1905      }
 1906   }
 1907   
 1908   short beginlabel(FILE *infileptr) {
 1909      int ch;
 1910      if ((ch=getc(infileptr))=='b') {
 1911         return 1;
 1912      } else {
 1913         ungetc(ch,infileptr);
 1914         return 0;
 1915      }
 1916   }
 1917   
 1918   short tickmark(FILE *infileptr) {
 1919      int ch;
 1920      if ((ch=getc(infileptr))=='t') {
 1921         return 1;
 1922      } else {
 1923         ungetc(ch,infileptr);
 1924         return 0;
 1925      }
 1926   }
 1927   
 1928   short tickmarklabel(FILE *infileptr) {
 1929      int ch;
 1930      if ((ch=getc(infileptr))=='l') {
 1931         return 1;
 1932      } else {
 1933         ungetc(ch,infileptr);
 1934         return 0;
 1935      }
 1936   }
 1937   
 1938   short endlabel(FILE *infileptr) {
 1939      int ch;
 1940      if ((ch=getc(infileptr))=='e') {
 1941         return 1;
 1942      } else {
 1943         ungetc(ch,infileptr);
 1944         return 0;
 1945      }
 1946   }
 1947   
 1948   /*
 1949    * The readaway_comments_and_blanks(FILE *infileptr,long *linenum) routine,
 1950    * well... , reads away comments and blanks from the input file (character
 1951    * stream) until next alphanumeric character, including any leading sign, or
 1952    * until end-of-file.
 1953    *
 1954    * Important: The "while ((!(isalnum(ch)||(ch=='-')||(ch=='+')))&&(ch!=EOF))"
 1955    * construction is absolutely necessary, since any preceding signs are not
 1956    * covered by isalnum(). Previously, isalnum() used to cover this as well,
 1957    * at least in the GNU development environment distributed with Linux (1996-
 1958    * 2000) and OSX 10.3 (2003), but somehow this seems to have changed, despite
 1959    * that we here follow a strict ISO C90 standard.
 1960    */
 1961   void readaway_comments_and_blanks(FILE *infileptr,long *linenum) {
 1962      int ch;
 1963      ch=getc(infileptr);
 1964      while ((!(isalnum(ch)||(ch=='-')||(ch=='+')))&&(ch!=EOF)) {
 1965         if (ch=='\n') {
 1966            (*linenum)++;
 1967         } else if (ch=='%') { /* if comment found, then read away the row */
 1968            while (((ch=getc(infileptr))!='\n')&&(ch!=EOF));
 1969            (*linenum)++;
 1970         }
 1971         ch=getc(infileptr);
 1972      }
 1973      if (ch!=EOF) ungetc(ch,infileptr);
 1974   }
 1975   
 1976   void scan_label(FILE *infileptr,stoketraject *st,int lnum,long *linenum,
 1977         pmap *map,long coordnum) {
 1978      char tmpstr[256];
 1979      int ch,k;
 1980   
 1981      /* Scan for positioning of label */
 1982      readaway_comments_and_blanks(infileptr,linenum);
 1983      if ((*map).verbose) {
 1984         fprintf(stdout,
 1985            "%s: Scanning label text starting at line %ld of trajectory file\n",
 1986               progname,*linenum);
 1987      }
 1988      /* tells at which coordinate index the label position is to be found */
 1989      (*st).label[lnum]=coordnum;
 1990      fscanf(infileptr,"%s",tmpstr);
 1991      if (!strcmp(tmpstr,"top")) {
 1992         (*st).labelpos[lnum]=TOPLABEL;
 1993      } else if (!strcmp(tmpstr,"ulft")) {
 1994         (*st).labelpos[lnum]=UPPERLEFTLABEL;
 1995      } else if (!strcmp(tmpstr,"lft")) {
 1996         (*st).labelpos[lnum]=LEFTLABEL;
 1997      } else if (!strcmp(tmpstr,"llft")) {
 1998         (*st).labelpos[lnum]=LOWERLEFTLABEL;
 1999      } else if (!strcmp(tmpstr,"bot")) {
 2000         (*st).labelpos[lnum]=BOTTOMLABEL;
 2001      } else if (!strcmp(tmpstr,"lrgt")) {
 2002         (*st).labelpos[lnum]=LOWERRIGHTLABEL;
 2003      } else if (!strcmp(tmpstr,"rgt")) {
 2004         (*st).labelpos[lnum]=RIGHTLABEL;
 2005      } else if (!strcmp(tmpstr,"urgt")) {
 2006         (*st).labelpos[lnum]=UPPERRIGHTLABEL;
 2007      } else {
 2008         fprintf(stderr,
 2009            "%s: Invalid string '%s' found at line %ld of trajectory file.\n",
 2010               progname,tmpstr,*linenum);
 2011         exit(1);
 2012      }
 2013      if ((*map).verbose) {
 2014         fprintf(stdout,
 2015            "%s: Scanned label positioning '%s' at line %ld of trajectory file\n",
 2016               progname,tmpstr,*linenum);
 2017      }
 2018   
 2019      /* Scan for label string to be printed */
 2020      while ((ch=getc(infileptr))==' ');
 2021      ungetc(ch,infileptr);
 2022      if ((ch=getc(infileptr))!='\"') {
 2023         fprintf(stderr,"%s: Error in line %ld of trajectory file. [ch=%c]\n",
 2024            progname,*linenum,ch);
 2025         fprintf(stderr,"%s: Use enclosing quote marks (\") around label text.\n",
 2026            progname);
 2027         exit(1);
 2028      }
 2029      k=0;
 2030      while ((ch=getc(infileptr))!='\"') {
 2031         k++;
 2032         if (ch!='\n') {
 2033            tmpstr[k]=ch;
 2034         } else {
 2035            fprintf(stderr,
 2036               "%s: Error: Reached end of line %ld without closing quote mark.\n",
 2037               progname,*linenum);
 2038            fprintf(stderr,"%s: (Check this label statement.)\n",progname);
 2039            exit(1);
 2040         }
 2041      }
 2042      (*st).labellength[lnum]=k;
 2043      while (0<k) {
 2044         (*st).labeltext[lnum][k]=tmpstr[k];
 2045         k--;
 2046      }
 2047   }
 2048   
 2049   /*
 2050   The scan_beginlabel() scans for a label string immediatelty after a
 2051   statement for a new trajectory, and places the string text, string length,
 2052   and label position relative the first point in thye data structure st.
 2053   The routine uses the more general scan_label() routine for the implementation,
 2054   in order to keep a comapct and consistent behaviour of the algorithm.
 2055   */
 2056   void scan_beginlabel(FILE *infile,stoketraject *st,long *linenum,pmap *map,
 2057         long coordnum) {
 2058      scan_label(infile,st,1,linenum,map,coordnum);
 2059   }
 2060   
 2061   void scan_ticklabel(FILE *infile,stoketraject *st,long *linenum, pmap *map,
 2062         long coordnum) {
 2063      if ((*map).verbose)
 2064         fprintf(stdout,"%s: Scanning label No %d\n",progname,(*st).numlabels);
 2065      if ((1<(*st).numlabels)&&((*st).numlabels<MAX_NUM_LABELS+2)) {
 2066         scan_label(infile,st,(*st).numlabels,linenum,map,coordnum);
 2067      } else {
 2068         fprintf(stderr,"%s: Error in scan_ticklabel routine:\n",progname);
 2069         fprintf(stderr,"%s: Index number %d out of range.\n",progname,
 2070            (*st).numlabels);
 2071      }
 2072   }
 2073   
 2074   void scan_endlabel(FILE *infile,stoketraject *st,long *linenum,pmap *map,
 2075         long coordnum) {
 2076      scan_label(infile,st,MAX_NUM_LABELS+2,linenum,map,coordnum);
 2077   }
 2078   
 2079   /*------------------------------------------------------------------------
 2080   | In this routine we do not only like to access the numerical values of the
 2081   | input data structure <st> which contains all information of the Stokes
 2082   | vector trajectory, but we also wish to update the values.
 2083   | Hence this routine takes a pointer to the structure as input, and internally
 2084   | dereferences this pointer whenever a value is to be updated. Typical
 2085   | usage of the routine, for example as when called from within the
 2086   | write_scanned_trajectories routine, is:
 2087   |    stoketraject st;
 2088   |     FILE *infileptr;
 2089   |    long int linenum=0;
 2090   |    st=initialize_stoke_trajectory();
 2091   |    infileptr=open_infile(map);
 2092   |    scan_for_stokes_triplet(infileptr,&st,&linenum);
 2093   |        . . .
 2094   ------------------------------------------------------------------------*/
 2095   void scan_for_stokes_triplet(FILE *infileptr,stoketraject *st,long *linenum) {
 2096      float s1,s2,s3;
 2097      if (!fscanf(infileptr,"%f",&s1)) {
 2098         fprintf(stderr,"%s: Error: Faulty S1 in line %ld of trajectory file.\n",
 2099            progname,*linenum);
 2100         exit(1);
 2101      }
 2102      if (!fscanf(infileptr,"%f",&s2)) {
 2103         fprintf(stderr,"%s: Error: Faulty S2 in line %ld of trajectory file.\n",
 2104            progname,*linenum);
 2105         exit(1);
 2106      }
 2107      if (!fscanf(infileptr,"%f",&s3)) {
 2108         fprintf(stderr,"%s: Error: Faulty S3 in line %ld of trajectory file.\n",
 2109            progname,*linenum);
 2110         exit(1);
 2111      }
 2112      ((*st).numcoords)++;
 2113      (*st).s1[(*st).numcoords]=s1;
 2114      (*st).s2[(*st).numcoords]=s2;
 2115      (*st).s3[(*st).numcoords]=s3;
 2116   }
 2117   
 2118   short visible(double s1,double s2,double s3,pmap *map) {
 2119      double vprod;
 2120      vprod=s1*cos((*map).rot_psi)*cos((*map).rot_phi)
 2121              -s2*sin((*map).rot_psi)*cos((*map).rot_phi)
 2122              +s3*sin((*map).rot_phi);
 2123      if (vprod>=0.0) { /* point (s1,s2,s3) located in visible region */
 2124         return 1;
 2125      } else { /* point (s1,s2,s3) located in hidden region */
 2126         return 0;
 2127      }
 2128   }
 2129   
 2130   /*
 2131   The point_just_became_hidden(&st,k) routine takes an index k as input
 2132   and checks if the corresponding Stokes vector (st.s1[k],st.s2[k],st.s3[k])
 2133   just became hidden.
 2134   The routine returns 1 (true) if (st.s1[k],st.s2[k],st.s3[k]) is hidden
 2135   while (st.s1[k-1],st.s2[k-1],st.s3[k-1]) is visible, or, for k=1, if
 2136   (st.s1[k],st.s2[k],st.s3[k]) is hidden; otherwise 0 (false) is returned.
 2137   This routine is useful for finding the start and end points of sub-trajectories
 2138   of hidden and visible parts of a Stokes trajectory scanned from file.
 2139   */
 2140   short point_just_became_hidden(stoketraject *st,long int k) {
 2141      if (k==1) {
 2142         if ((*st).visible[k]==0) {
 2143            return 1;
 2144         } else {
 2145            return 0;
 2146         }
 2147      } else if ((2<=k)&&(k<=(*st).numcoords)) {
 2148         if (((*st).visible[k]==0)&&((*st).visible[k-1]==1)) {
 2149            return 1;
 2150         } else {
 2151            return 0;
 2152         }
 2153      } else {
 2154         fprintf(stderr,"%s: Error in routine point_just_became_hidden()!\n",
 2155            progname);
 2156         fprintf(stderr,"%s: Index %ld is out of range of current trajectory.\n",
 2157            progname,k);
 2158         fprintf(stderr,"%s: (Maximum possible index is %ld.)\n",
 2159            progname,(*st).numcoords);
 2160         exit(1);
 2161      }
 2162   }
 2163   
 2164   short point_just_became_visible(stoketraject *st,long int k) {
 2165      if (k==1) {
 2166         if ((*st).visible[k]==1) {
 2167            return 1;
 2168         } else {
 2169            return 0;
 2170         }
 2171      } else if ((2<=k)&&(k<=(*st).numcoords)) {
 2172         if (((*st).visible[k]==1)&&((*st).visible[k-1]==0)) {
 2173            return 1;
 2174         } else {
 2175            return 0;
 2176         }
 2177      } else {
 2178         fprintf(stderr,"%s: Error in routine point_just_became_visible()!\n",
 2179            progname);
 2180         fprintf(stderr,"%s: Index %ld is out of range of current trajectory.\n",
 2181            progname,k);
 2182         fprintf(stderr,"%s: (Maximum possible index is %ld.)\n",
 2183            progname,(*st).numcoords);
 2184         exit(1);
 2185      }
 2186   }
 2187   
 2188   /*
 2189    * The get_screen_coordinates() routine calculates the projected screen
 2190    * coordinates (x,y) in the image from a Stokes triplet (s1,s2,s3).
 2191    * In order to update the output parameters we use these as pointers,
 2192    * to give a behaviour as when using reference variables in languages
 2193    * as Pascal.
 2194    */
 2195   void get_screen_coordinates(double *x,double *y,
 2196         double s1,double s2,double s3,pmap *map) {
 2197      double snorm;
 2198      (*x)= s1*sin((*map).rot_psi)+s2*cos((*map).rot_psi);
 2199      (*y)=-s1*cos((*map).rot_psi)*sin((*map).rot_phi)
 2200             +s2*sin((*map).rot_psi)*sin((*map).rot_phi)+s3*cos((*map).rot_phi);
 2201      if ((*map).use_normalized_stokes_params) {
 2202         snorm=sqrt(s1*s1+s2*s2+s3*s3);
 2203         (*x)=(*x)/snorm;
 2204         (*y)=(*y)/snorm;
 2205      }
 2206   }
 2207   
 2208   void add_subtrajectory(FILE *outfileptr,stoketraject *st,
 2209         long int ka,long int kb,pmap *map,char type[256]) {
 2210      long int k;
 2211      short j;
 2212      double x,y; /* trash variables for screen coordinates */
 2213      j=1;
 2214      fprintf(outfileptr,"   pickup pencircle scaled %f pt;\n",
 2215         (*map).paththickness);
 2216      if (ka<kb) { /* only draw paths of two points or more */
 2217         for (k=ka;k<=kb;k++) {
 2218            j++;
 2219            get_screen_coordinates(&x,&y,(*st).s1[k],(*st).s2[k],(*st).s3[k],map);
 2220            if (k==ka) {
 2221      /*
 2222               fprintf(outfileptr,"%%\n%% Drawing path No %d\n%%\n",pathnum);
 2223      */
 2224               fprintf(outfileptr,"   p := makepath makepen ");
 2225            }
 2226            if (j==(NUM_COORDS_PER_METAPOST_LINE+1)) {
 2227               fprintf(outfileptr,"\n    ");
 2228               j=1;
 2229            }
 2230            if (k>ka)
 2231               fprintf(outfileptr,"%s",((*map).use_bezier_curves)?"..":"--");
 2232            fprintf(outfileptr,"(%1.4f,%1.4f)",x,y);
 2233            if (k==kb) {
 2234               fprintf(outfileptr,";\n");
 2235               if ((kb==(*st).numcoords)&&((*map).draw_paths_as_arrows)) {
 2236                  if ((*map).reverse_arrow_paths) {
 2237                     fprintf(outfileptr,"   drawarrow reverse p scaled radius");
 2238                  } else {
 2239                     fprintf(outfileptr,"   drawarrow p scaled radius");
 2240                  }
 2241               } else {
 2242                  fprintf(outfileptr,"   draw p scaled radius");
 2243               }
 2244               if (!strcmp(type,"hidden")) { /* end of a hidden segment */
 2245                  if ((*map).draw_hidden_dashed) {
 2246                     fprintf(outfileptr," dashed evenly withcolor black;\n");
 2247                  } else {
 2248                     fprintf(outfileptr," withcolor %f [black,white];\n",
 2249                        (*map).hiddengraytone);
 2250                  }
 2251               } else { /* end of a visible segment */
 2252                  fprintf(outfileptr," withcolor black;\n");
 2253               }
 2254            }
 2255         }
 2256      }
 2257   }
 2258   
 2259   /*
 2260    * Sort out visible from hidden parts of the trajectory.
 2261    */
 2262   void sort_out_visible_and_hidden(FILE *outfile,stoketraject *st,pmap *map) {
 2263      long int k;
 2264      for (k=1;k<=(*st).numcoords;k++) {
 2265         if (visible((*st).s1[k],(*st).s2[k],(*st).s3[k],map)) {
 2266            (*st).visible[k]=1;
 2267         } else {
 2268            (*st).visible[k]=0;
 2269         }
 2270      }
 2271   }
 2272   
 2273   /*
 2274    * Write hidden parts of the trajectory to file.
 2275    */
 2276   void add_hidden_subtrajectories(FILE *outfile,stoketraject *st,pmap *map) {
 2277      long int k,ka,kb;
 2278      k=ka=kb=1;
 2279      while (k<=(*st).numcoords) {
 2280         if (point_just_became_hidden(st,k)) {
 2281            ka=k; /* first index of hidden sub-trajectory */
 2282            kb=ka-1; /* initialize kb<ka as check for the following while loop */
 2283            while ((k<=(*st).numcoords)&&(kb<ka)) {
 2284               if (point_just_became_visible(st,k)) {
 2285                  kb=k-1;
 2286               } else {
 2287                  k++;
 2288               }
 2289            }
 2290            kb=k-1; /* last index of hidden sub-trajectory */
 2291            if ((*map).verbose) {
 2292               fprintf(stdout,
 2293                  "%s: Adding hidden subtrajectory from ka=%ld to kb=%ld\n",
 2294                  progname,ka,kb);
 2295            }
 2296            add_subtrajectory(outfile,st,ka,kb,map,"hidden");
 2297         }
 2298         k++;
 2299      }
 2300   }
 2301   
 2302   /*
 2303    * Write visible parts of the trajectory to file.
 2304    */
 2305   void add_visible_subtrajectories(FILE *outfile,stoketraject *st,pmap *map) {
 2306      long int k,ka,kb;
 2307      k=ka=kb=1;
 2308      while (k<=(*st).numcoords) {
 2309         if (point_just_became_visible(st,k)) {
 2310            ka=k; /* first index of visible sub-trajectory */
 2311            kb=ka-1; /* initialize kb<ka as check for the following while loop */
 2312            while ((k<=(*st).numcoords)&&(kb<ka)) {
 2313               if (point_just_became_hidden(st,k)) {
 2314                  kb=k-1;
 2315               } else {
 2316                  k++;
 2317               }
 2318            }
 2319            kb=k-1; /* last index of visible sub-trajectory */
 2320            if ((*map).verbose) {
 2321               fprintf(stdout,
 2322                  "%s: Adding visible subtrajectory from ka=%ld to kb=%ld\n",
 2323                  progname,ka,kb);
 2324            }
 2325            /*
 2326             * We allow the ends of visible parts to go one coordinate step
 2327             * behind the sphere so as to give a smooth connection to the
 2328             * hidden parts.
 2329             */
 2330            if (ka>1) ka--;
 2331            if (kb<(*st).numcoords) kb++;
 2332            add_subtrajectory(outfile,st,ka,kb,map,"visible");
 2333         }
 2334         k++;
 2335      }
 2336   }
 2337   
 2338   void add_scanned_trajectory(FILE *outfileptr,stoketraject *st,pmap *map,
 2339         short viewtype) {
 2340      sort_out_visible_and_hidden(outfileptr,st,map);
 2341      if (viewtype==HIDDEN) {
 2342         add_hidden_subtrajectories(outfileptr,st,map);
 2343      } else if (viewtype==VISIBLE) {
 2344         add_visible_subtrajectories(outfileptr,st,map);
 2345      } else { /* if invalid |viewtype| option */
 2346         fprintf(stderr,"%s: Error in add_scanned_trajectory: ",progname);
 2347         fprintf(stderr,"Invalid viewtype value! (%d)\n",viewtype);
 2348         exit(1);
 2349      }
 2350   }
 2351   
 2352   /*
 2353    * The add_scanned_labels() routines the previously scanned text labels
 2354    * of the trajectory |st| to the output MetaPost source file.
 2355    */
 2356   void add_scanned_labels(FILE *outfileptr,stoketraject *st,pmap *map) {
 2357      long int j,k;
 2358      double x,y;
 2359   
 2360      for (k=1;k<=(MAX_NUM_LABELS+2);k++) {
 2361         if ((*st).labellength[k]>0) {
 2362            if ((*st).labelpos[k]==TOPLABEL) {
 2363               fprintf(outfileptr,"   label.top");
 2364            } else if ((*st).labelpos[k]==UPPERLEFTLABEL) {
 2365               fprintf(outfileptr,"   label.ulft");
 2366            } else if ((*st).labelpos[k]==LEFTLABEL) {
 2367               fprintf(outfileptr,"   label.lft");
 2368            } else if ((*st).labelpos[k]==LOWERLEFTLABEL) {
 2369               fprintf(outfileptr,"   label.llft");
 2370            } else if ((*st).labelpos[k]==BOTTOMLABEL) {
 2371               fprintf(outfileptr,"   label.bot");
 2372            } else if ((*st).labelpos[k]==LOWERRIGHTLABEL) {
 2373               fprintf(outfileptr,"   label.lrt");
 2374            } else if ((*st).labelpos[k]==RIGHTLABEL) {
 2375               fprintf(outfileptr,"   label.rt");
 2376            } else if ((*st).labelpos[k]==UPPERRIGHTLABEL) {
 2377               fprintf(outfileptr,"   label.urt");
 2378            } else {
 2379               fprintf(stderr,
 2380                  "%s: add_scanned_labels: Invalid labelpos (%d) detected ",
 2381                  progname,(*st).labelpos[k]);
 2382               fprintf(stderr,"at label No %ld\n",k);
 2383               fprintf(stderr,
 2384                  "%s: add_scanned_labels: Labelstring is \042",progname);
 2385               for (j=1;j<=(*st).labellength[k];j++)
 2386                  fprintf(stderr,"%c",(*st).labeltext[k][j]);
 2387               fprintf(stderr,"\044\n");
 2388               exit(1);
 2389            }
 2390            get_screen_coordinates(&x,&y,(*st).s1[(*st).label[k]],
 2391               (*st).s2[(*st).label[k]],(*st).s3[(*st).label[k]],map);
 2392            fprintf(outfileptr,"(btex ");
 2393            for (j=1;j<=(*st).labellength[k];j++)
 2394               fprintf(outfileptr,"%c",(*st).labeltext[k][j]);
 2395            fprintf(outfileptr," etex,(%f,%f)*radius);\n",x,y);
 2396         }
 2397      }
 2398   }
 2399   
 2400   void scan_for_tickmark(FILE *infileptr,stoketraject *st,long *linenum) {
 2401      if (tickmark(infileptr)) {
 2402         (*st).numtickmarks++;
 2403         (*st).tickmark[(*st).numtickmarks]=(*st).numcoords;
 2404      }
 2405   }
 2406   
 2407   void scan_for_tickmarklabel(FILE *infileptr,stoketraject *st,pmap *map,
 2408         long *linenum) {
 2409      if (tickmarklabel(infileptr)) {
 2410         (*st).numlabels++;
 2411         (*st).label[(*st).numlabels]=(*st).numcoords;
 2412         scan_label(infileptr,st,(*st).numlabels,linenum,map,
 2413            (*st).label[(*st).numlabels]);
 2414      }
 2415   }
 2416   
 2417   void get_tickmark_screen_coordinates(double *xa,double *ya,
 2418         double *xb,double *yb,long int k,stoketraject *st,pmap *map) {
 2419      double xt,yt,snorm,s1n,s2n,s3n,s1a,s2a,s3a,s1b,s2b,s3b,
 2420         s0,s1,s2,s3,p1,p2,p3,q1,q2,q3;
 2421   
 2422      /* Calculate the normalized approximate tangential vector to path
 2423       * at the tickmark point with index k.
 2424       */
 2425      if ((*st).tickmark[k]==1) {
 2426         q1=(*st).s1[(*st).tickmark[k]+1]-(*st).s1[(*st).tickmark[k]];
 2427         q2=(*st).s2[(*st).tickmark[k]+1]-(*st).s2[(*st).tickmark[k]];
 2428         q3=(*st).s3[(*st).tickmark[k]+1]-(*st).s3[(*st).tickmark[k]];
 2429      } else if ((*st).tickmark[k]==(*st).numcoords) {
 2430         q1=(*st).s1[(*st).tickmark[k]]-(*st).s1[(*st).tickmark[k]-1];
 2431         q2=(*st).s2[(*st).tickmark[k]]-(*st).s2[(*st).tickmark[k]-1];
 2432         q3=(*st).s3[(*st).tickmark[k]]-(*st).s3[(*st).tickmark[k]-1];
 2433      } else if ((1<(*st).tickmark[k])&&((*st).tickmark[k]<(*st).numcoords)) {
 2434         q1=(*st).s1[(*st).tickmark[k]+1]-(*st).s1[(*st).tickmark[k]-1];
 2435         q2=(*st).s2[(*st).tickmark[k]+1]-(*st).s2[(*st).tickmark[k]-1];
 2436         q3=(*st).s3[(*st).tickmark[k]+1]-(*st).s3[(*st).tickmark[k]-1];
 2437      } else {
 2438         fprintf(stderr,"%s: In get_tickmark_screen_coordinates:\n",progname);
 2439         fprintf(stderr,"%s: Index k=%ld out of valid range!\n",
 2440            progname,(*st).tickmark[k]);
 2441         exit(1);
 2442      }
 2443      snorm=sqrt(q1*q1+q2*q2+q3*q3);
 2444      q1=q1/snorm;
 2445      q2=q2/snorm;
 2446      q3=q3/snorm;
 2447   
 2448      /* get normalized (unitary) Stokes vector */
 2449      s1=(*st).s1[(*st).tickmark[k]];
 2450      s2=(*st).s2[(*st).tickmark[k]];
 2451      s3=(*st).s3[(*st).tickmark[k]];
 2452      s0=sqrt(s1*s1+s2*s2+s3*s3);
 2453      s1n=s1/s0;
 2454      s2n=s2/s0;
 2455      s3n=s3/s0;
 2456   
 2457      /* Calculate the normalized vector transverse to the tangent of path,
 2458       * {\bf p}={\bf s}\times{\bf q}/|{\bf s}\times{\bf q}|
 2459       */
 2460      p1=s2n*q3-s3n*q2;
 2461      p2=s3n*q1-s1n*q3;
 2462      p3=s1n*q2-s2n*q1;
 2463      snorm=sqrt(p1*p1+p2*p2+p3*p3);
 2464      p1=p1/snorm;
 2465      p2=p2/snorm;
 2466      p3=p3/snorm;
 2467   
 2468      /* Calculate the 1st endpoint of tick mark in Stokes parameter space */
 2469      s1a=s1n+0.028213*p1;
 2470      s2a=s2n+0.028213*p2;
 2471      s3a=s3n+0.028213*p3;
 2472   
 2473   /*   s1a=s1+(*map).ticksize*p1;
 2474      s2a=s2+(*map).ticksize*p2;
 2475      s3a=s3+(*map).ticksize*p3; */
 2476   
 2477      /* Calculate the 2nd endpoint of tick mark in Stokes parameter space */
 2478      s1b=s1n-0.028213*p1;
 2479      s2b=s2n-0.028213*p2;
 2480      s3b=s3n-0.028213*p3;
 2481   
 2482      /* Get the screen coordinates of the ends of the tick mark */
 2483      get_screen_coordinates(&xt,&yt,s0*s1a,s0*s2a,s0*s3a,map);
 2484      (*xa)=xt;
 2485      (*ya)=yt;
 2486      get_screen_coordinates(&xt,&yt,s0*s1b,s0*s2b,s0*s3b,map);
 2487      (*xb)=xt;
 2488      (*yb)=yt;
 2489   }
 2490   
 2491   void add_scanned_tickmarks(FILE *outfileptr,stoketraject *st,pmap *map,
 2492         short viewtype) {
 2493      long int k;
 2494      double xa,ya,xb,yb;
 2495   
 2496      fprintf(outfileptr,"   pickup pencircle scaled %f pt;\n",
 2497         (*map).paththickness/2.0);
 2498      for (k=1;k<=(*st).numtickmarks;k++) {
 2499         get_tickmark_screen_coordinates(&xa,&ya,&xb,&yb,k,st,map);
 2500         fprintf(outfileptr,"   p:=makepath makepen (%f,%f)--(%f,%f);\n",
 2501            xa,ya,xb,yb);
 2502         if (((*st).visible[(*st).tickmark[k]])&&(viewtype==VISIBLE)) {
 2503            fprintf(outfileptr,"   draw p scaled radius;\n");
 2504         } else if ((!((*st).visible[(*st).tickmark[k]]))&&(viewtype==HIDDEN)) {
 2505            fprintf(outfileptr,"   draw p scaled radius");
 2506            fprintf(outfileptr," withcolor %f [black,white];\n",
 2507               (*map).hiddengraytone);
 2508         }
 2509      }
 2510   }
 2511   
 2512   /*-----------------------------------------------------------------------------
 2513   | (REWRITE THIS COMMENT)
 2514   | Draw the trajectories of Stokes parameters on the Poincare sphere (if a
 2515   | filename containing the path was specified at startup of the program).
 2516   | In order to properly take care of the fact that we might be over-writing
 2517   | parts of the 'visible' parts of the path with 'invisible' parts from the
 2518   | hidden parts of the Poincare sphere, we must traverse the path-drawing
 2519   | routine twice, in the second run only writing the parts that are visible.
 2520   | This is only important when the shade of the front, 'visible' paths are
 2521   | different from the paths in the back, 'invisible' parts of the sphere.
 2522   | Notice that, for example, with hidden parts of the path drawn with black,
 2523   | dashed lines, and visible part drawn with solid black lines, it does not
 2524   | matter whether this two-pass drawing takes place or not.
 2525   |
 2526   | This routine also particularly checks for statments of labels to be applied
 2527   | at the begin and end points of the trajectories, that is to say, checking
 2528   | for statements of type 'b <pos> "TeX label string"' or 'e <pos> "TeX label
 2529   | string"' immediately after the 'p' and 'q' statements of each trajectory.
 2530   | These labels could equally well be supplied as regular labels after the
 2531   | first and last triplets of data, but in some cases it is for example in
 2532   | batch mode operation easier to simply concatenate in the Stokes data without
 2533   | having to manually go in and add these label statements. This is the very
 2534   | reason for this additional feature of begin and end label statements in the
 2535   | data syntax as accepted by the program.
 2536   |
 2537   | Depending on whether the boolean input variable |viewtype| is |HIDDEN| or
 2538   | |VISIBLE|, this routine will either flush out the hidden or visible parts
 2539   | of a Stokes trajectory. This construction is done in order for the main
 2540   | program to call the routine in two passes after each other, one for the
 2541   | hidden parts and one for the visible parts (which we do not wish to have
 2542   | later overwritten by hidden parts of other trajectories scanned later on.
 2543   | This switch is not explicitly present in the algorithm of this routine, but
 2544   | rather just used as passed on to the routines |add_scanned_trajectory()| and
 2545   | |add_scanned_tickmarks()|, which take care of adding the hidden ans visible
 2546   | trajectories with tickmarks in a coherent manner.
 2547   |
 2548   | As being the key routine in the flushing of the Stokes trajectories to
 2549   | file, this routine also takes care of writing any present tick marks or
 2550   | labels, by calling the |add_scanned_tickmarks()| and |add_scanned_labels()|
 2551   | routines.
 2552   |
 2553   | Throughout the parsing and writing of the Stokes parameter trajectory to
 2554   | file, the entire data set is kept in the structure |st| (short for Stokes
 2555   | trajectory), of type |stoketraject| as defined in the definitions section
 2556   | of this program.
 2557   -----------------------------------------------------------------------------*/
 2558   void write_scanned_trajectories(FILE *outfileptr,pmap map,short viewtype) {
 2559      FILE *infileptr = NULL;
 2560      stoketraject st; /* data structure for keeping track of trajectories */
 2561      long int linenum,coord;
 2562      int i=0,pathnum,k;
 2563   
 2564      initialize_stoke_trajectory(&st); /* Allocate memory for arrays etc. */
 2565      reset_stokes_trajectory_struct(&st); /* Make sure all data is cleared */
 2566      pathnum=0; /* used to keep track on the number of paths scanned so far */
 2567      infileptr=open_infile(map); /* open file to read Stokes triplets from */
 2568      /*-------------------------------------------------------------------------
 2569      | If the file pointer |infileptr| after opening the file with the
 2570      | |open_infile()| routine points to NULL, this indicates that the filename
 2571      | was empty, and that now Stokes trajectories hence should be read from
 2572      | file. If there was a problem opening the file for reading, due to read
 2573      | permissions, file existence etc, this would automatically have taken
 2574      | care of by the |open_infile()| routine.
 2575      -------------------------------------------------------------------------*/
 2576      if (infileptr!=NULL) {
 2577         fprintf(outfileptr,"  oldahangle:=ahangle;\n");
 2578         fprintf(outfileptr,"  ahangle:=%f;\n",map.arrowheadangle);
 2579         fprintf(outfileptr,"  pickup pencircle scaled %f pt;\n",
 2580            map.paththickness);
 2581         linenum=1; /* counter for keeping track of line numbers in input file */
 2582         coord=0; /* counter for keeping track of trajectory coordinate numbers */
 2583         while (new_trajectory(infileptr)) {
 2584            if (map.verbose) fprintf(stdout,
 2585               "%s: New trajectory detected at line %ld\n",progname,linenum);
 2586            i++; /* increment trajectory counter */
 2587            readaway_comments_and_blanks(infileptr,&linenum);
 2588            if (beginlabel(infileptr)) { /* check for text label at begin point */
 2589               if (map.verbose)
 2590                  fprintf(stdout, "%s: Begin-point label detected at line %ld\n",
 2591                     progname,linenum);
 2592               scan_beginlabel(infileptr,&st,&linenum,&map,1);
 2593               readaway_comments_and_blanks(infileptr,&linenum);
 2594               if (map.verbose) {
 2595                  fprintf(stdout,"%s: Parsed begin label string '",progname);
 2596                  for (k=1;k<=st.labellength[1];k++)
 2597                     fprintf(stdout,"%c",st.labeltext[1][k]);
 2598                  fprintf(stdout,"' [%d characters]\n",st.labellength[1]);
 2599               }
 2600            }
 2601            if (map.verbose) fprintf(stdout,
 2602               "%s: Scanning Stokes trajectory starting at line %ld.\n",
 2603                  progname,linenum);
 2604            st.numcoords=0; /* reset coordinate counter */
 2605            while (!end_of_trajectory(infileptr)) { /* scan current trajectory */
 2606               scan_for_stokes_triplet(infileptr,&st,&linenum);
 2607               coord++;
 2608               readaway_comments_and_blanks(infileptr,&linenum);
 2609               scan_for_tickmark(infileptr,&st,&linenum);
 2610               readaway_comments_and_blanks(infileptr,&linenum);
 2611               scan_for_tickmarklabel(infileptr,&st,&map,&linenum);
 2612               readaway_comments_and_blanks(infileptr,&linenum);
 2613            }
 2614            if (map.verbose) fprintf(stdout,
 2615               "%s: End of Stokes trajectory detected at line %ld.\n",
 2616                  progname,linenum);
 2617            readaway_comments_and_blanks(infileptr,&linenum);
 2618            if (endlabel(infileptr)) { /* check for text label at end point */
 2619               if (map.verbose) fprintf(stdout,
 2620                  "%s: End-point label detected at line %ld\n",progname,linenum);
 2621               scan_endlabel(infileptr,&st,&linenum,&map,coord);
 2622               readaway_comments_and_blanks(infileptr,&linenum);
 2623               if (map.verbose) {
 2624                  fprintf(stdout,"%s: Parsed end label string '",progname);
 2625                  for (k=1;k<=st.labellength[MAX_NUM_LABELS+2];k++)
 2626                     fprintf(stdout,"%c",st.labeltext[MAX_NUM_LABELS+2][k]);
 2627                  fprintf(stdout,"' [%d characters]\n",
 2628                     st.labellength[MAX_NUM_LABELS+2]);
 2629               }
 2630            }
 2631            add_scanned_trajectory(outfileptr,&st,&map,viewtype);
 2632            add_scanned_tickmarks(outfileptr,&st,&map,viewtype);
 2633            add_scanned_labels(outfileptr,&st,&map);
 2634            reset_stokes_trajectory_struct(&st);
 2635         } /* End of "while (new_trajectory(infileptr)) ..." */
 2636         fclose(infileptr);
 2637         fprintf(outfileptr,"  ahangle:=oldahangle;\n");
 2638      } else {
 2639   /*
 2640    *      fprintf(stderr,
 2641    *         "%s: Sorry, trajectory scanning from stdin not implemented yet.",
 2642    *         progname);
 2643    *      exit(1);
 2644    */
 2645      }
 2646   }
 2647   
 2648   /*-----------------------------------------------------------------------------
 2649   | Draw any additional arrows (if specified by user) onto the Poincare sphere.
 2650   | This feature is useful for example whenever there is a distinct transition
 2651   | between two polarization states that need to be pointed out, or whenever
 2652   | just a principal change is to be pointed out.
 2653   -----------------------------------------------------------------------------*/
 2654   void write_additional_arrows(FILE *outfileptr,pmap map) {
 2655      int i;
 2656      double s1,s2,s3,s1a,s2a,s3a,s1b,s2b,s3b,t,dt=0.02,xtmp,ytmp;
 2657      if (map.numarrows > 0) {
 2658         fprintf(outfileptr,
 2659            "%%\n"
 2660            "%% Draw the paths of the arrows specified by the user.\n"
 2661            "%%\n");
 2662         fprintf(outfileptr,"   pickup pencircle scaled 0.5pt;\n");
 2663         for (i=1;i<=map.numarrows;i++) {
 2664            s1a=map.arrows[1][i];
 2665            s2a=map.arrows[2][i];
 2666            s3a=map.arrows[3][i];
 2667            s1b=map.arrows[4][i];
 2668            s2b=map.arrows[5][i];
 2669            s3b=map.arrows[6][i];
 2670            if (map.use_normalized_stokes_params) {
 2671               xtmp=sqrt(s1a*s1a+s2a*s2a+s3a*s3a);
 2672               s1a=s1a/xtmp;
 2673               s2a=s2a/xtmp;
 2674               s3a=s3a/xtmp;
 2675               xtmp=sqrt(s1b*s1b+s2b*s2b+s3b*s3b);
 2676               s1b=s1b/xtmp;
 2677               s2b=s2b/xtmp;
 2678               s3b=s3b/xtmp;
 2679            }
 2680            fprintf(outfileptr,"   p := makepath makepen\n");
 2681            fprintf(outfileptr,"      ");
 2682            for (t=0.0;t<=0.5+dt/1000.0;t+=dt) {
 2683               if (t > dt) fprintf(outfileptr,"                  ");
 2684               if (t > 0.0) fprintf(outfileptr,"..");
 2685               s1=(1.0 - t)*s1a + t*s1b;
 2686               s2=(1.0 - t)*s2a + t*s2b;
 2687               s3=(1.0 - t)*s3a + t*s3b;
 2688               xtmp=sqrt(s1*s1+s2*s2+s3*s3);
 2689               s1 /= xtmp;
 2690               s2 /= xtmp;
 2691               s3 /= xtmp;
 2692               xtmp= s1*sin(map.rot_psi)+s2*cos(map.rot_psi);
 2693               ytmp=-s1*cos(map.rot_psi)*sin(map.rot_phi) + \
 2694                       s2*sin(map.rot_psi)*sin(map.rot_phi) + \
 2695                       s3*cos(map.rot_phi);
 2696               fprintf(outfileptr,"(%f,%f)",xtmp,ytmp);
 2697               if ((t > (0.0+dt/1000.0)) && (t < (0.5-dt/1000.0)))
 2698                  fprintf(outfileptr,"\n");
 2699            }
 2700            fprintf(outfileptr,";\n");
 2701            if ((-0.5 <= map.arrows[7][i]) && (map.arrows[7][i] < 0.5)) {
 2702               fprintf(outfileptr,\
 2703                  "   drawarrow p scaled radius withcolor %f [white,black];\n",\
 2704                                                             map.arrows[8][i]);
 2705            } else if ((0.5 <= map.arrows[7][i]) && (map.arrows[7][i] < 1.5)) {
 2706               fprintf(outfileptr,\
 2707                  "   drawarrow p scaled radius dashed evenly withcolor %f"
 2708                  " [white,black];\n", map.arrows[8][i]);
 2709            }
 2710            fprintf(outfileptr,"   p := makepath makepen\n");
 2711            fprintf(outfileptr,"      ");
 2712            for (t=0.5;t<=1.0+dt/1000.0;t+=dt) {
 2713               if (t > 0.5+dt) fprintf(outfileptr,"                  ");
 2714               if (t > 0.5) fprintf(outfileptr,"..");
 2715               s1=(1.0 - t)*s1a + t*s1b;
 2716               s2=(1.0 - t)*s2a + t*s2b;
 2717               s3=(1.0 - t)*s3a + t*s3b;
 2718               xtmp=sqrt(s1*s1+s2*s2+s3*s3);
 2719               s1 /= xtmp;
 2720               s2 /= xtmp;
 2721               s3 /= xtmp;
 2722               xtmp= s1*sin(map.rot_psi)+s2*cos(map.rot_psi);
 2723               ytmp=-s1*cos(map.rot_psi)*sin(map.rot_phi) + \
 2724                       s2*sin(map.rot_psi)*sin(map.rot_phi) + \
 2725                       s3*cos(map.rot_phi);
 2726               fprintf(outfileptr,"(%f,%f)",xtmp,ytmp);
 2727               if ((t > (0.5+dt/1000.0)) && (t < (1.0-dt/1000.0)))
 2728                  fprintf(outfileptr,"\n");
 2729            }
 2730            fprintf(outfileptr,";\n");
 2731            if ((-0.5 <= map.arrows[7][i]) && (map.arrows[7][i] < 0.5)) {
 2732               fprintf(outfileptr,\
 2733                  "   draw p scaled radius withcolor %f [white,black];\n",\
 2734                                                              map.arrows[8][i]);
 2735            } else if ((0.5 <= map.arrows[7][i]) && (map.arrows[7][i] < 1.5)) {
 2736               fprintf(outfileptr,\
 2737                  "   draw p scaled radius dashed evenly withcolor %f"
 2738                  " [white,black];\n", map.arrows[8][i]);
 2739            }
 2740         }
 2741      }
 2742   }
 2743   
 2744   /*-----------------------------------------------------------------------------
 2745   | Draw the coordinate axes of the $(S_1,S_2,S_3)$-space.
 2746   -----------------------------------------------------------------------------*/
 2747   void write_coordinate_axes(FILE *outfileptr,pmap map) {
 2748   fprintf(outfileptr,
 2749    "%%\n"
 2750    "%% Draw the $S_1$-, $S_2$- and $S_3$-axis of the Poincare sphere.\n"
 2751    "%% First of all, calculate the transformations of the intersections\n"
 2752    "%% for the unity sphere.\n"
 2753    "%%\n");
 2754   fprintf(outfileptr,
 2755    "%% Used variables:\n"
 2756    "%%\n"
 2757    "%%    behind_distance : Specifies the relative distance of the coordi-\n"
 2758    "%%                      axes to be plotted behind origo (in negative di-\n"
 2759    "%%                      rection of respective axis.\n"
 2760    "%%\n");
 2761   fprintf(outfileptr,
 2762    "%%   outside_distance_s1 : The relative distance from origo to the point\n"
 2763    "%%                         of the arrow head of the coordinate axis S1.\n"
 2764    "%%                         If this is set to 1.0, the arrow head will\n"
 2765    "%%                         point directly at the Poincare sphere.\n"
 2766    "%%\n");
 2767   fprintf(outfileptr,
 2768    "%%   outside_distance_s2 : Same as above, except that this one controls\n"
 2769    "%%                         the S2 coordinate axis instead.\n"
 2770    "%%\n");
 2771   fprintf(outfileptr,
 2772    "%%   outside_distance_s3 : Same as above, except that this one controls\n"
 2773    "%%                         the S3 coordinate axis instead.\n"
 2774    "%%\n");
 2775   fprintf(outfileptr,
 2776    "%%    insidecolval :    Specifies the shade of gray to use for the parts\n"
 2777    "%%                      of the coordinate axes that are inside the Poin-\n"
 2778    "%%                      care sphere. Values must be between 0 and 1,\n"
 2779    "%%                      where:  '0.0' <=> 'white';  '1.0' <=> 'black'\n"
 2780    "%%\n");
 2781   fprintf(outfileptr,
 2782    "   behind_distance_s1  := -%f;\n"
 2783    "   behind_distance_s2  := -%f;\n"
 2784    "   behind_distance_s3  := -%f;\n"
 2785    "   outside_distance_s1 :=  %f;\n"
 2786    "   outside_distance_s2 :=  %f;\n"
 2787    "   outside_distance_s3 :=  %f;\n"
 2788    "   insidecolval := .85;    %% '0.0' <=> 'white';  '1.0' <=> 'black'\n"
 2789    "\n",
 2790      map.neg_axis_length_s1, map.neg_axis_length_s2, map.neg_axis_length_s3,
 2791      map.pos_axis_length_s1, map.pos_axis_length_s2, map.pos_axis_length_s3);
 2792      fprintf(outfileptr,
 2793        "   pickup pencircle scaled %f pt;\n",map.coordaxisthickness);
 2794   
 2795   fprintf(outfileptr,
 2796    "%%\n"
 2797    "%% Start with drawing the x-axis...\n"
 2798    "%%\n"
 2799    "   x_bis_start :=  radius*behind_distance_s1*cosd(rot_psi)*cosd(rot_phi);\n"
 2800    "   y_bis_start :=  radius*behind_distance_s1*sind(rot_psi);\n"
 2801    "   z_bis_start := -radius*behind_distance_s1*cosd(rot_psi)*sind(rot_phi);\n"
 2802    "   x_bis_intersect :=  radius*cosd(rot_psi)*cosd(rot_phi);\n"
 2803    "   y_bis_intersect :=  radius*sind(rot_psi);\n"
 2804    "   z_bis_intersect := -radius*cosd(rot_psi)*sind(rot_phi);\n");
 2805   if (map.draw_axes_inside_sphere) fprintf(outfileptr,
 2806    "   p := makepath makepen \n"
 2807    "             (y_bis_start,z_bis_start)--(y_bis_intersect,z_bis_intersect);\n"
 2808    "   draw p dashed evenly withcolor insidecolval [white,black];\n");
 2809   fprintf(outfileptr,
 2810    "   p := makepath makepen (y_bis_intersect,z_bis_intersect)--\n"
 2811    "             (outside_distance_s1*y_bis_intersect,\n"
 2812    "              outside_distance_s1*z_bis_intersect);\n"
 2813    "   drawarrow p;\n");
 2814   fprintf(outfileptr,
 2815    "   label.%s(btex $%s$ etex,\n"
 2816    "             (outside_distance_s1*y_bis_intersect,\n"
 2817    "              outside_distance_s1*z_bis_intersect));\n"
 2818    "\n",map.axislabelposition_s1,
 2819   (map.user_specified_axislabels ? map.axislabel_s1
 2820      : (map.use_normalized_stokes_params ? "S_1/S_0" : "S_1")));
 2821   
 2822   fprintf(outfileptr,
 2823    "%%\n"
 2824    "%% ... then draw the y-axis ...\n"
 2825    "%%\n"
 2826    "   x_bis_start := -radius*behind_distance_s2*sind(rot_psi)*cosd(rot_phi);\n"
 2827    "   y_bis_start :=  radius*behind_distance_s2*cosd(rot_psi);\n"
 2828    "   z_bis_start :=  radius*behind_distance_s2*sind(rot_psi)*sind(rot_phi);\n"
 2829    "   x_bis_intersect := -radius*sind(rot_psi)*cosd(rot_phi);\n"
 2830    "   y_bis_intersect :=  radius*cosd(rot_psi);\n"
 2831    "   z_bis_intersect :=  radius*sind(rot_psi)*sind(rot_phi);\n");
 2832   if (map.draw_axes_inside_sphere) fprintf(outfileptr,
 2833    "   p := makepath makepen \n"
 2834    "             (y_bis_start,z_bis_start)--(y_bis_intersect,z_bis_intersect);\n"
 2835    "   draw p dashed evenly withcolor insidecolval [white,black];\n");
 2836   fprintf(outfileptr,
 2837    "   p := makepath makepen (y_bis_intersect,z_bis_intersect)--\n"
 2838    "             (outside_distance_s2*y_bis_intersect,\n"
 2839    "              outside_distance_s2*z_bis_intersect);\n"
 2840    "   drawarrow p;\n");
 2841   fprintf(outfileptr,
 2842    "   label.%s(btex $%s$ etex,\n"
 2843    "             (outside_distance_s2*y_bis_intersect,\n"
 2844    "              outside_distance_s2*z_bis_intersect));\n"
 2845    "\n",map.axislabelposition_s2,
 2846   (map.user_specified_axislabels ? map.axislabel_s2
 2847      : (map.use_normalized_stokes_params ? "S_2/S_0" : "S_2")));
 2848   
 2849   fprintf(outfileptr,
 2850    "%%\n"
 2851    "%% ... then, finally, draw the z-axis.\n"
 2852    "%%\n"
 2853    "   x_bis_start := radius*behind_distance_s3*sind(rot_phi);\n"
 2854    "   y_bis_start := 0.0;\n"
 2855    "   z_bis_start := radius*behind_distance_s3*cosd(rot_phi);\n"
 2856    "   x_bis_intersect := radius*sind(rot_phi);\n"
 2857    "   y_bis_intersect := 0.0;\n"
 2858    "   z_bis_intersect := radius*cosd(rot_phi);\n");
 2859   if (map.draw_axes_inside_sphere) fprintf(outfileptr,
 2860    "   p := makepath makepen \n"
 2861    "             (y_bis_start,z_bis_start)--(y_bis_intersect,z_bis_intersect);\n"
 2862    "   draw p dashed evenly withcolor insidecolval [white,black];\n");
 2863   fprintf(outfileptr,
 2864    "   p := makepath makepen (y_bis_intersect,z_bis_intersect)--\n"
 2865    "             (outside_distance_s3*y_bis_intersect,\n"
 2866    "              outside_distance_s3*z_bis_intersect);\n"
 2867    "   drawarrow p;\n");
 2868   fprintf(outfileptr,
 2869    "   label.%s(btex $%s$ etex,\n"  /* For z-label to the right of S_3 axis */
 2870    "             (outside_distance_s3*y_bis_intersect,\n"
 2871    "              outside_distance_s3*z_bis_intersect));\n"
 2872    "\n",map.axislabelposition_s3,
 2873   (map.user_specified_axislabels ? map.axislabel_s3
 2874      : (map.use_normalized_stokes_params ? "S_3/S_0" : "S_3")));
 2875   }
 2876   
 2877   /*-----------------------------------------------------------------------------
 2878   | In case the user has specified that an additional coordinate system should
 2879   | be drawn, draw the coordinate axes of the additional $(x,y,z)$-space.
 2880   | The rule is that only coordinate axes that have explicitly declared
 2881   | axis labels should be drawn.  The reason for that I have chosen this
 2882   | way of operation is that the rotated, additional coordinate system might
 2883   | still have one axis in common with the original one, potentially making
 2884   | ugly double-drawn axes with corresponding labels.
 2885   -----------------------------------------------------------------------------*/
 2886   void write_additional_coordinate_axes(FILE *outfileptr,pmap map) {
 2887   if (map.user_specified_additional_coordinate_system)  /* whoaaaou again ... */
 2888   {
 2889   fprintf(outfileptr,
 2890    "%%\n"
 2891    "%% Draw the $S_1$-, $S_2$- and $S_3$-axis of the Poincare sphere.\n"
 2892    "%% First of all, calculate the transformations of the intersections\n"
 2893    "%% for the unity sphere.\n"
 2894    "%%\n");
 2895   fprintf(outfileptr,
 2896    "%% Used variables are similar to the ones described for\n"
 2897    "%% drawing the original coordinate system.\n"
 2898    "%%\n");
 2899   fprintf(outfileptr,
 2900    "   xtra_behind_distance_x  := -%f;\n"
 2901    "   xtra_behind_distance_y  := -%f;\n"
 2902    "   xtra_behind_distance_z  := -%f;\n",
 2903      map.xtra_neg_axis_length_x,
 2904      map.xtra_neg_axis_length_y,
 2905      map.xtra_neg_axis_length_z);
 2906   fprintf(outfileptr,
 2907    "   xtra_outside_distance_x :=  %f;\n"
 2908    "   xtra_outside_distance_y :=  %f;\n"
 2909    "   xtra_outside_distance_z :=  %f;\n",
 2910      map.xtra_pos_axis_length_x,
 2911      map.xtra_pos_axis_length_y,
 2912      map.xtra_pos_axis_length_z);
 2913   fprintf(outfileptr,
 2914    "   insidecolval := .85;    %% '0.0' <=> 'white';  '1.0' <=> 'black'\n\n");
 2915   
 2916   if (map.user_specified_xtra_axislabel_x)
 2917   {
 2918   fprintf(outfileptr,
 2919    "%%\n"
 2920    "%% Start with drawing the x-axis...\n"
 2921    "%%\n");
 2922   fprintf(outfileptr,
 2923    "   x_bis_start :=  radius * xtra_behind_distance_x\n"
 2924    "                          * cosd(rot_psi + delta_rot_psi)\n"
 2925    "                          * cosd(rot_phi + delta_rot_phi);\n"
 2926    "   y_bis_start :=  radius * xtra_behind_distance_x\n"
 2927    "                          * sind(rot_psi + delta_rot_psi);\n"
 2928    "   z_bis_start := -radius * xtra_behind_distance_x\n"
 2929    "                          * cosd(rot_psi + delta_rot_psi)\n"
 2930    "                          * sind(rot_phi + delta_rot_phi);\n");
 2931   fprintf(outfileptr,
 2932    "   x_bis_intersect :=  radius * cosd(rot_psi + delta_rot_psi)\n"
 2933    "                              * cosd(rot_phi + delta_rot_phi);\n"
 2934    "   y_bis_intersect :=  radius * sind(rot_psi + delta_rot_psi);\n"
 2935    "   z_bis_intersect := -radius * cosd(rot_psi + delta_rot_psi)\n"
 2936    "                              * sind(rot_phi + delta_rot_phi);\n");
 2937   if (map.draw_axes_inside_sphere) fprintf(outfileptr,
 2938    "   p := makepath makepen \n"
 2939    "             (y_bis_start,z_bis_start)--(y_bis_intersect,z_bis_intersect);\n"
 2940    "   draw p dashed evenly withcolor insidecolval [white,black];\n");
 2941   fprintf(outfileptr,
 2942    "   p := makepath makepen (y_bis_intersect,z_bis_intersect)--\n"
 2943    "             (xtra_outside_distance_x * y_bis_intersect,\n"
 2944    "              xtra_outside_distance_x * z_bis_intersect);\n"
 2945    "   drawarrow p;\n");
 2946   fprintf(outfileptr,
 2947    "   label.bot(btex $%s$ etex,\n"
 2948    "             (xtra_outside_distance_x * y_bis_intersect,\n"
 2949    "              xtra_outside_distance_x * z_bis_intersect));\n"
 2950    "\n", map.xtra_axislabel_x);
 2951   }
 2952   
 2953   if (map.user_specified_xtra_axislabel_y)
 2954   {
 2955   fprintf(outfileptr,
 2956    "%%\n"
 2957    "%% ... then draw the y-axis ...\n"
 2958    "%%\n");
 2959   fprintf(outfileptr,
 2960    "   x_bis_start := -radius * xtra_behind_distance_y\n"
 2961    "                          * sind(rot_psi + delta_rot_psi)\n"
 2962    "                          * cosd(rot_phi + delta_rot_phi);\n"
 2963    "   y_bis_start :=  radius * xtra_behind_distance_y\n"
 2964    "                          * cosd(rot_psi + delta_rot_psi);\n"
 2965    "   z_bis_start :=  radius * xtra_behind_distance_y\n"
 2966    "                          * sind(rot_psi + delta_rot_psi)\n"
 2967    "                          * sind(rot_phi + delta_rot_phi);\n");
 2968   fprintf(outfileptr,
 2969    "   x_bis_intersect := -radius * sind(rot_psi + delta_rot_psi)\n"
 2970    "                              * cosd(rot_phi + delta_rot_phi);\n"
 2971    "   y_bis_intersect :=  radius * cosd(rot_psi + delta_rot_psi);\n"
 2972    "   z_bis_intersect :=  radius * sind(rot_psi + delta_rot_psi)\n"
 2973    "                              * sind(rot_phi + delta_rot_phi);\n");
 2974   if (map.draw_axes_inside_sphere) fprintf(outfileptr,
 2975    "   p := makepath makepen \n"
 2976    "             (y_bis_start,z_bis_start)--(y_bis_intersect,z_bis_intersect);\n"
 2977    "   draw p dashed evenly withcolor insidecolval [white,black];\n");
 2978   fprintf(outfileptr,
 2979    "   p := makepath makepen (y_bis_intersect,z_bis_intersect)--\n"
 2980    "             (xtra_outside_distance_y * y_bis_intersect,\n"
 2981    "              xtra_outside_distance_y * z_bis_intersect);\n"
 2982    "   drawarrow p;\n");
 2983   fprintf(outfileptr,
 2984    "   label.bot(btex $%s$ etex,\n"
 2985    "             (xtra_outside_distance_y * y_bis_intersect,\n"
 2986    "              xtra_outside_distance_y * z_bis_intersect));\n"
 2987    "\n", map.xtra_axislabel_y);
 2988   }
 2989   
 2990   if (map.user_specified_xtra_axislabel_z)
 2991   {
 2992   fprintf(outfileptr,
 2993    "%%\n"
 2994    "%% ... then, finally, draw the z-axis.\n"
 2995    "%%\n"
 2996    "   x_bis_start := radius * xtra_behind_distance_z\n"
 2997    "                         * sind(rot_phi + delta_rot_phi);\n"
 2998    "   y_bis_start := 0.0;\n"
 2999    "   z_bis_start := radius * xtra_behind_distance_z\n"
 3000    "                         * cosd(rot_phi + delta_rot_phi);\n"
 3001    "   x_bis_intersect := radius * sind(rot_phi + delta_rot_phi);\n"
 3002    "   y_bis_intersect := 0.0;\n"
 3003    "   z_bis_intersect := radius * cosd(rot_phi + delta_rot_phi);\n");
 3004   if (map.draw_axes_inside_sphere) fprintf(outfileptr,
 3005    "   p := makepath makepen \n"
 3006    "             (y_bis_start,z_bis_start)--(y_bis_intersect,z_bis_intersect);\n"
 3007    "   draw p dashed evenly withcolor insidecolval [white,black];\n");
 3008   fprintf(outfileptr,
 3009    "   p := makepath makepen (y_bis_intersect,z_bis_intersect)--\n"
 3010    "             (xtra_outside_distance_z * y_bis_intersect,\n"
 3011    "              xtra_outside_distance_z * z_bis_intersect);\n"
 3012    "   drawarrow p;\n");
 3013   fprintf(outfileptr,
 3014    "   label.top(btex $%s$ etex,\n"
 3015    "             (xtra_outside_distance_z * y_bis_intersect,\n"
 3016    "              xtra_outside_distance_z * z_bis_intersect));\n"
 3017    "\n", map.xtra_axislabel_z);
 3018   }
 3019   }
 3020   }
 3021   
 3022   void write_included_auxiliary_source(FILE *outfileptr,pmap map) {
 3023      if (map.user_specified_auxfile) {
 3024         fprintf(outfileptr,
 3025            "%%\n"
 3026            "%% The following external file is included"
 3027            " (using the --auxsource option): %s  [MetaPost source]\n"
 3028            "%%\n"
 3029            "   input %s\n",map.auxfilename, map.auxfilename);
 3030      }
 3031      fprintf(outfileptr,
 3032         "   endfig;\n"
 3033         "end\n"
 3034      );
 3035   }
 3036   
 3037   /*-----------------------------------------------------------------------
 3038   | Generate Encapsulated PostScript output from the MetaPost-source.
 3039   -----------------------------------------------------------------------*/
 3040   void generate_eps_image(pmap map) {
 3041      char tmpstr[256];
 3042      long int llx,lly,urx,ury;
 3043   
 3044      /*--------------------------------------------------------------------
 3045      | Compile the MetaPost code into EPS with control codes for TeX.
 3046      --------------------------------------------------------------------*/
 3047      sprintf(tmpstr,"mp -job-name %s %s;",map.epsjobname,map.outfilename);
 3048      if (map.verbose)
 3049         fprintf(stdout,"%s: Executing system command: %s\n",progname,tmpstr);
 3050      system(tmpstr);
 3051      /*--------------------------------------------------------------------
 3052      | Use TeX for generating a self-containing DVI output of the figure.
 3053      --------------------------------------------------------------------*/
 3054      sprintf(tmpstr,"tex -job-name %s \'\\input epsf\\nopagenumbers"
 3055         "\\centerline{\\epsfbox{%s.1}}\\bye\';",map.epsjobname,map.epsjobname);
 3056      if (map.verbose)
 3057         fprintf(stdout,"%s: Executing system command: %s\n",progname,tmpstr);
 3058      system(tmpstr);
 3059   
 3060      /*--------------------------------------------------------------------
 3061      | Use DVIPS for generating a self-containing EPS output with a tight
 3062      | bounding box from the DVI.
 3063      --------------------------------------------------------------------*/
 3064      sprintf(tmpstr,"dvips -D1200 -E %s.dvi -o %s.eps",
 3065         map.epsjobname,map.epsjobname);
 3066      if (map.verbose)
 3067         fprintf(stdout,"%s: Executing system command: %s\n",progname,tmpstr);
 3068      system(tmpstr);
 3069   
 3070      /*--------------------------------------------------------------------
 3071      | Finally, extract the bounding box of the generated Encapsulated
 3072      | PostScript file, in order to get an idea of the natural physical
 3073      | size of the figure. This is useful in order to enter correct
 3074      | settings as the resulting image is included in for example TeX
 3075      | documents, in which case the typical inclusion yields
 3076      |    \input epsf
 3077      |    \centerline{\epsfxsize=<reported x-size>\epsfbox{<filename>.eps}}
 3078      --------------------------------------------------------------------*/
 3079      sprintf(tmpstr,"%s.eps",map.epsjobname);
 3080      scan_for_boundingbox(tmpstr,&llx,&lly,&urx,&ury);
 3081      fprintf(stdout,"%s: Bounding box of %s:\n"
 3082         "     width=%-4.2f mm (%ld pts), height=%-4.2f mm (%ld pts)\n",
 3083         progname,tmpstr,(urx-llx)*(25.4/72.27),urx-llx,
 3084                  (ury-lly)*(25.4/72.27),ury-lly);
 3085   }
 3086   
 3087   int main(int argc, char *argv[]) {
 3088      pmap map;              /* The data structure containing input parameters */
 3089      FILE *outfileptr=NULL; /* The destination file for MetaPost code */
 3090   
 3091      map=parse_command_line(argc,argv);
 3092      if (map.verbose) show_banner();
 3093      display_arrow_specs(map);
 3094      outfileptr=open_outfile(map);
 3095      write_header(outfileptr,map,argc,argv);
 3096      write_euler_angle_specs(outfileptr,map);
 3097      write_sphere_shading_specs(outfileptr,map);
 3098      write_shaded_sphere(outfileptr,map); /* Generate the background sphere */
 3099      write_equators(outfileptr,map); /* Generate the equators S_k=0, k=1,2,3 */
 3100      write_scanned_trajectories(outfileptr,map,HIDDEN);
 3101      write_scanned_trajectories(outfileptr,map,VISIBLE);
 3102      write_additional_arrows(outfileptr,map);
 3103      write_coordinate_axes(outfileptr,map);
 3104      write_additional_coordinate_axes(outfileptr,map);
 3105      write_included_auxiliary_source(outfileptr,map);
 3106      fclose(outfileptr);
 3107      if (map.generate_eps_output) generate_eps_image(map);
 3108      return(0);  /* exit clean from error codes if successful execution */
 3109   }
 3110   

Return to previous page

Generated by ::viewsrc::

Last modified Thursday 23 Jan 2014