Search:

Return to previous page

Contents of file 'boxcount/boxcount.tex':




Deprecated: Function split() is deprecated in /storage/content/45/2011745/jonsson.eu/public_html/php/htmlicise.php on line 46
    1   \input cwebmac
    2   % File:        boxcount.w [CWEB source code]
    3   % Created:     May 8, 2006 [v.1.0]
    4   % Last change: November 26, 2011 [v.1.6]
    5   % Author:      Fredrik Jonsson
    6   % Description: The CWEB source code for the BOXCOUNT simulator of nonlinear
    7   %              magneto-optical Bragg gratings. For information on the CWEB
    8   %              programming language, see http://www.literateprogramming.com.
    9   % Compilation: Compile this program by using the enclosed Makefile, or use
   10   %              the blocks of the Makefile as listed in the documentation
   11   %              (file boxcount.ps or boxcount.pdf). The C source code (as
   12   %              generated from this CWEB code) conforms to the ANSI standard
   13   %              for the C programming language (which is equivalent to the
   14   %              ISO C90 standard for C).
   15   %
   16   % Copyright (C) 2006-2011, Fredrik Jonsson. Non-commercial copying welcome.
   17   %
   18   \input epsf
   19   \def\version{1.6}
   20   \def\lastrevdate{November 26, 2011}
   21   \font\eightcmr=cmr8
   22   \font\tensc=cmcsc10
   23   \font\eightcmssq=cmssq8
   24   \font\eightcmssqi=cmssqi8
   25   \font\twentycmcsc=cmcsc10 at 20 truept
   26   \def\boxcount{{\eightcmr BOXCOUNT\spacefactor1000}}
   27   \def\poincare{{\eightcmr POINCARE\spacefactor1000}}
   28   \def\ANSI{{\eightcmr ANSI\spacefactor1000}} % The language standard we stick to
   29   \def\SI{{\eightcmr SI\spacefactor1000}}   % Another standard for physical units
   30   \def\GNU{{\eightcmr GNU\spacefactor1000}}       % GNU is Not Unix
   31   \def\GCC{{\eightcmr GCC\spacefactor1000}}       % The GNU C-compiler
   32   \def\CEE{{\eightcmr C\spacefactor1000}}         % The C programming language
   33   \def\ANSICEE{{\eightcmr ANSI~C\spacefactor1000}}% The ANSI-C standard language
   34   \def\CWEB{{\eightcmr CWEB\spacefactor1000}}     % The CWEB programming language
   35   \def\MATLAB{{\eightcmr MATLAB\spacefactor1000}} % The MATLAB ditto
   36   \def\UNIX{{\eightcmr UNIX\spacefactor1000}}
   37   \def\CYGWIN{{\eightcmr CYGWIN\spacefactor1000}}
   38   \def\CTANGLE{{\eightcmr CTANGLE\spacefactor1000}}
   39   \def\CWEAVE{{\eightcmr CWEAVE\spacefactor1000}}
   40   \def\OSX{{OS\,X}}
   41   \def\re{\mathop{\rm Re}\nolimits} % The real part of a complex number
   42   \def\im{\mathop{\rm Im}\nolimits} % The imaginary part of a complex number
   43   \def\dollar{\char'044}             % The `$' character
   44   \def\tothepower{\char'136}         % The `^' character
   45   \def\onehalf{{\textstyle{{1}\over{2}}}}
   46   \def\threefourth{{\textstyle{{3}\over{4}}}}
   47   \def\endalg{\vrule height 1.4ex width .6ex depth .4ex} % Rectangular bullet
   48   \def\ie{i.\thinspace{e.}~\ignorespaces}
   49   \def\eg{e.\thinspace{g.}~\ignorespaces}
   50   \let\,\thinspace
   51   %% \def\subsection[#1]{\goodbreak\noindent{\it #1}\par\nobreak\noindent}
   52   %
   53   % Define a handy macro for listing the steps of an algorithm.
   54   %
   55   \newdimen\aitemindent \aitemindent=26pt
   56   \newdimen\aitemleftskip \aitemleftskip=36pt
   57   \def\aitem[#1]{\smallbreak\noindent\hbox to 10pt{}%
   58     \hbox to\aitemindent{\bf #1\hfill}%
   59     \hangindent\aitemleftskip\ignorespaces}
   60   %
   61   % Define a handy macro for bibliographic references.
   62   %
   63   \newdimen\refitemindent \refitemindent=18pt
   64   \def\refitem[#1]{\smallbreak\noindent%
   65     \hbox to\refitemindent{[#1]\hfill}%
   66     \hangindent\refitemindent\ignorespaces}
   67   \def\references{\medskip\goodbreak\noindent{\bf References}\nobreak\smallskip}
   68   %
   69   % Define a handy macro for nicely typeset variable descriptions.
   70   %
   71   \newdimen\varitemindent \varitemindent=100pt
   72   \def\varitem[#1]#2{\smallbreak\noindent\hbox to 20pt{}%
   73     \hbox to\varitemindent{#1\hfill}%
   74     \hangindent 120pt\ignorespaces#2\par}
   75   %
   76   % Define a handy macro for nicely typeset descriptions of command line options.
   77   %
   78   \newdimen\optitemindent \optitemindent=80pt
   79   \def\optitem[#1]#2{\smallbreak\noindent\hbox to\hsize{\hskip40pt{#1}\hfill}
   80     \hangindent\optitemindent\ignorespaces{\hskip\optitemindent{#2}}\smallskip}
   81   %
   82   % Define a handy macro for the list of program revisions.
   83   %
   84   \newdimen\citemindent \citemindent=80pt
   85   \newdimen\citemleftskip \citemleftskip=90pt
   86   \def\citem[#1]{\smallbreak\noindent\hbox to 10pt{}%
   87     \hbox to\citemindent{\bf #1\hfill}%
   88     \hangindent\citemleftskip\ignorespaces}
   89   \def\citindent{\smallbreak\noindent\hbox to 10pt{}%
   90     \hbox to\citemindent{\hfil}%
   91     \hangindent\citemleftskip\ignorespaces}
   92   %
   93   % Define the \beginvrulealign and \endvrulealign macros as described in
   94   % Donald Knuth's The TeXbook, Appendix D: Dirty Tricks. These macros are
   95   % used in typesetting nicely looking tables.
   96   %
   97   \def\beginvrulealign{\setbox0=\vbox\bgroup}
   98   \def\endvrulealign{\egroup % now \box0 holds the entire alignment
   99     \setbox0=\vbox{\setbox2=\hbox{\vrule height\ht0 depth\dp0 width0pt}
  100       \unvbox0 \setbox0=\lastbox % now \box0 is the bottom row
  101       \nointerlineskip \copy0 % put it back
  102       \global\setbox1=\hbox{} % initialize box that will contain rules
  103       \setbox4=\hbox{\unhbox0 % now open up the bottom row
  104         \loop \skip0=\lastskip \unskip % remove tabskip glue
  105         \advance\skip0 by-.4pt % rules are .4 pt wide
  106         \divide\skip0 by 2
  107         \global\setbox1=\hbox{\hskip\skip0\vrule\hskip\skip0
  108           \unhbox2\unhbox1}%
  109         \setbox2=\lastbox % remove alignment entry
  110         \ifhbox2 \setbox2=\hbox{\kern\wd2}\repeat}}%
  111     \hbox{\rlap{\box0}\box1}} % superimpose the alignment on the rules
  112   %
  113   % Make sure that the METAPOST logo can be loaded even in plain TeX.
  114   %
  115   \ifx\MP\undefined
  116       \ifx\manfnt\undefined
  117               \font\manfnt=logo10
  118       \fi
  119       \ifx\manfntsl\undefined
  120               \font\manfntsl=logosl10
  121       \fi
  122       \def\MP{{\ifdim\fontdimen1\font>0pt \let\manfnt = \manfntsl \fi
  123         {\manfnt META}\-{\manfnt POST}}\spacefactor1000 }%
  124   \fi
  125   \ifx\METAPOST\undefined \let\METAPOST=\MP \fi
  126   
  127   \datethis
  128   
  129   
  130   \N{1}{1}Introduction.
  131   \vskip 100pt
  132   \centerline{\twentycmcsc Boxcount}
  133   \vskip 20pt
  134   \centerline{A program for calculating box-counting estimates to the fractal
  135   dimension of curves in the plane.}
  136   \vskip 2pt
  137   \centerline{(Version \version\ of \lastrevdate)}
  138   \vskip 10pt
  139   \centerline{Written by Fredrik Jonsson}
  140   \vskip 10pt
  141   \centerline{\epsfxsize=88mm\epsfbox{figures/figure-05.1}}
  142   \vskip 10pt
  143   \noindent
  144   This \CWEB\footnote{${}^\dagger$}{For information on the \CWEB\ programming
  145   language by Donald E.~Knuth, as well as samples of \CWEB\ programs, see
  146   {\tt http://www-cs-faculty.stanford.edu/\~\ \kern -5pt knuth/cweb.html}.
  147   For general information on literate programming, see
  148   {\tt http://www.literateprogramming.com}.}
  149   computer program calculates box-counting estimates of the fractal dimension
  150   of curves in the two-dimensional plane.
  151   
  152   In the box-counting estimate to the fractal dimension of a graph in the domain
  153   $\{x,y:x_{\rm min}\le x\le x_{\rm max},y_{\rm min}\le y\le y_{\rm max}\}$, a
  154   grid of squares, each of horizontal dimension $(x_{\rm max}-x_{\rm min})/2^m$
  155   and vertical dimension $(y_{\rm max}-y_{\rm min})/2^m$, is superimposed onto
  156   the graph for integer numbers $m$.
  157   By counting the total number of such squares $N_m$ needed to cover the entire
  158   graph at a given $m$ (hence the term ``box counting''), an estimate $D_m$ to
  159   the fractal dimension $D$ (or Hausdorff dimension) is obtained as
  160   $D_m=\ln(N_m)/\ln(2^m)$.
  161   This procedure may be repeated may times, with $D_m\to D$ as $m\to\infty$
  162   for real fractal sets. However, for finite-depth fractals (as generated by
  163   a computer), some limit on $m$ is necessary in order to prevent trivial
  164   convergence towards $D_m\to1$.
  165   
  166   In addition to mere numerical calculation, the program also generates graphs
  167   of the box distributions, in form of \METAPOST\ code which can be
  168   post-processed by other programs.
  169   \bigskip
  170   \noindent Copyright \copyright\ Fredrik Jonsson, 2006--2011.
  171   All rights reserved.
  172   \vfill
  173   
  174   \fi
  175   
  176   \N{1}{2}The CWEB programming language.
  177   For the reader who might not be familiar with the concept of the
  178   \CWEB\ programming language, the following citations hopefully will
  179   be useful. For further information, as well as freeware compilers for
  180   compiling \CWEB\ source code, see {\tt http://www.literateprogramming.com}.
  181   \bigskip
  182   
  183   {\narrower\narrower\narrower\narrower\eightcmssqi\noindent
  184   I believe that the time is ripe for significantly better documentation of
  185   programs, and that we can best achieve this by considering programs to be
  186   works of literature. Hence, my title: `Literate Programming.'
  187   
  188   Let us change our traditional attitude to the construction of programs:
  189   Instead of imagining that our main task is to instruct a computer what to
  190   do, let us concentrate rather on explaining to human beings what we want
  191   a computer to do.
  192   
  193   The practitioner of literate programming can be regarded as an essayist,
  194   whose main concern is with exposition and excellence of style. Such an
  195   author, with thesaurus in hand, chooses the names of variables carefully
  196   and explains what each variable means. He or she strives for a program
  197   that is comprehensible because its concepts have been introduced in an
  198   order that is best for human understanding, using a mixture of formal and
  199   informal methods that reinforce each other.
  200   \smallskip
  201   {\eightcmssq --\,Donald Knuth,}
  202   {\eightcmssqi The CWEB System of Structured Documentation}
  203   {\eightcmssq (Addison-Wesley, Massachusetts, 1994)}
  204   }
  205   \bigskip
  206   
  207   {\narrower\narrower\narrower\narrower\eightcmssqi\noindent
  208   The philosophy behind CWEB is that an experienced system programmer, who
  209   wants to provide the best possible documentation of his or her software
  210   products, needs two things simultaneously: a language like \TeX\ for
  211   formatting, and a language like C for programming. Neither type of language
  212   can provide the best documentation by itself; but when both are appropriately
  213   combined, we obtain a system that is much more useful than either language
  214   separately.
  215   
  216   The structure of a software program may be thought of as a `WEB' that is
  217   made up of many interconnected pieces. To document such a program we want to
  218   explain each individual part of the web and how it relates to its neighbors.
  219   The typographic tools provided by \TeX\ give us an opportunity to explain the
  220   local structure of each part by making that structure visible, and the
  221   programming tools provided by languages like C make it possible for us to
  222   specify the algorithms formally and unambiguously. By combining the two,
  223   we can develop a style of programming that maximizes our ability to perceive
  224   the structure of a complex piece of software, and at the same time the
  225   documented programs can be mechanically translated into a working software
  226   system that matches the documentation.
  227   
  228   Besides providing a documentation tool, CWEB enhances the C language by
  229   providing the ability to permute pieces of the program text, so that a
  230   large system can be understood entirely in terms of small sections and
  231   their local interrelationships. The CTANGLE program is so named because
  232   it takes a given web and moves the sections from their web structure into
  233   the order required by C; the advantage of programming in CWEB is that the
  234   algorithms can be expressed in ``untangled'' form, with each section
  235   explained separately. The CWEAVE program is so named because it takes a
  236   given web and intertwines the \TeX\ and C portions contained in each
  237   section, then it knits the whole fabric into a structured document.
  238   \smallskip
  239   {\eightcmssq --\,Donald Knuth, ``Literate Programming'', in}
  240   {\eightcmssqi Literate Programming}
  241   {\eightcmssq (CSLI Lecture Notes, Stanford, 1992)}
  242   }
  243   
  244   \fi
  245   
  246   \N{1}{3}Revision history of the program.
  247   \medskip
  248   
  249   \citem[2006-05-08]{[v.1.0]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
  250   First properly working version of the \boxcount\ program, written in \CWEB\
  251   and (\ANSI-conformant) \CEE. I have so far compiled the code with \GCC\ using
  252   the {\tt --pedantic} option and verified that the box-covering routine
  253   \PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}(\,)} and its more low-level
  254   core engine
  255   \PB{\\{box\_intersection}(\,)} both work properly, by direct inspection of the
  256   compiled
  257   \METAPOST\ graphs generated by this routine.
  258   That the numbers of counted boxes are correct has also been verified and the
  259   only remaining blocks to be added are related to the extraction of the fractal
  260   dimension as such.
  261   This first version of the \boxcount\ program consists of 52671 bytes
  262   (1292 lines) of \CWEB\ code, under \CYGWIN\ generating an executable file
  263   of 22561 bytes and 33 pages of program documentation in 10~pt typeface.
  264   
  265   \citem[2006-05-09]{[v.1.1]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
  266   Changed the block for the estimate of the fractal dimension, so that the
  267   estimate now is obtained as the average of $\ln(N_{\rm boxes})/\ln(2^n)$
  268   for a set of $n$ such that $N_{\rm min}\le n\le N_{\rm max}$, rather than
  269   performing a linear curve fit to the data.
  270   In order to sample statistical information on the estimate, such as standard
  271   deviation, average deviation and skewness, I incorporated a slightly modified
  272   version of the routine \PB{\\{moment}(\,)} from {\sl Numerical Recipes in C}.
  273   Also added a preliminary section describing the test application of \boxcount\
  274   to the Koch fractal, being a simple test case which is easily implemented by
  275   means of a recursive generation of the data set for the input trajectory.
  276   As of today, the \boxcount\ program consists of 59788 bytes
  277   (1444 lines) of \CWEB\ code, under \CYGWIN\ generating an executable file
  278   of 24253 bytes and 38 pages of program documentation in 10~pt typeface.
  279   
  280   \citem[2006-05-14]{[v.1.2]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
  281   Added a section on the command-line options for supplying the \boxcount\
  282   program with input parameters. Polished the section on the example of
  283   estimation of the box-counting dimension of the Koch fractal, and in
  284   particular changed the example from the Koch curve to the Koch snowflake
  285   instead, just for the sake of visual beauty.
  286   As of today, the \boxcount\ program consists of 72560 bytes (1699 lines) of
  287   \CWEB\ code, under \OSX\ generating an executable file of 24996 bytes and 41
  288   pages of program documentation in 10~pt typeface.
  289   
  290   \citem[2006-05-17]{[v.1.3]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
  291   Added documentation on the \PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}(%
  292   \,)} routine
  293   and generally cleaned up in the documentation of the program.
  294   As of today, the \boxcount\ program consists of 82251 bytes
  295   (1844 lines) of \CWEB\ code, under \CYGWIN\ generating an executable file
  296   of 29152 bytes and 40 pages of program documentation in 10~pt typeface.
  297   
  298   \citem[2006-06-17]{[v.1.4]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
  299   Added the {\tt --graphsize} option, in order to override the default graph
  300   size. Also changed the inclusion of the input trajectory in the boxmaps, so
  301   that the \boxcount\ program rather than using a \MP\ call of the form
  302   {\tt gdraw "input.dat"} now chops the trajectory into proper pieces and
  303   includes the entire trajectory explicitly in the generated graph.
  304   This way the output \MP\ code naturally increases considerably in size, but
  305   is now at least self-sustaining even is separated from the original data file
  306   containing the input trajectory.
  307   
  308   \citem[2006-10-25]{[v.1.5]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
  309   Added two pages of text on the boxcounting estimate of the fractal dimension
  310   of the Koch snowflake fractal in the example section.
  311   
  312   \citem[2011-11-26]{[v.1.6]} {\tt <http://jonsson.eu>}\hfill\break
  313   Updated \.{Makefile}:s for the generation of figures. Also corrected a rather
  314   stupid way of removing preceeding paths of file names.
  315   
  316   \fi
  317   
  318   \N{1}{4}Compiling the source code. The program is written in \CWEB, generating
  319   \ANSICEE\ (ISO C90) conforming source code and documentation as plain
  320   \TeX-source, and is to be compiled using the sequences as outlined in the
  321   {\tt Makefile} listed below.
  322   \bigskip
  323   {\obeylines\obeyspaces\tt
  324   \#
  325   \# Makefile designed for use with ctangle, cweave, gcc, and plain TeX.
  326   \#
  327   \# Copyright (C) 2002-2011, Fredrik Jonsson <http://jonsson.eu>
  328   \#
  329   \# The CTANGLE program converts a CWEB source document into a C program which
  330   \# may be compiled in the usual way. The output file includes \#line specifica-
  331   \# tions so that debugging can be done in terms of the CWEB source file.
  332   \#
  333   \# The CWEAVE program converts the same CWEB file into a TeX file that may be
  334   \# formatted and printed in the usual way. It takes appropriate care of typo-
  335   \# graphic details like page layout and the use of indentation, italics,
  336   \# boldface, etc., and it supplies extensive cross-index information that it
  337   \# gathers automatically.
  338   \#
  339   \# CWEB allows you to prepare a single document containing all the information
  340   \# that is needed both to produce a compilable C program and to produce a well-
  341   \# formatted document describing the program in as much detail as the writer
  342   \# may desire.  The user of CWEB ought to be familiar with TeX as well as C.
  343   \#
  344   PROJECT  = boxcount
  345   CTANGLE  = ctangle
  346   CWEAVE   = cweave
  347   CC       = gcc
  348   CCOPTS   = -O2 -Wall -ansi -std=iso9899:1990 -pedantic
  349   LNOPTS   = -lm
  350   TEX      = tex
  351   DVIPS    = dvips
  352   DVIPSOPT = -ta4 -D1200
  353   METAPOST = mp
  354   PS2PDF   = ps2pdf
  355   ~   ~
  356   all: \$(PROJECT) \$(PROJECT).pdf
  357   ~        @echo
  358   "==============================================================="
  359   ~        @echo " To verify the execution performance of the BOXCOUNT program"
  360   ~        @echo " on the Koch snowflake fractal, run 'make verification'."
  361   ~        @echo
  362   "==============================================================="
  363   ~   ~
  364   \$(PROJECT): \$(PROJECT).o
  365   ~        \$(CC) \$(CCOPTS) -o \$(PROJECT) \$(PROJECT).o \$(LNOPTS)
  366   ~   ~
  367   \$(PROJECT).o: \$(PROJECT).w
  368   ~        \$(CTANGLE) \$(PROJECT)
  369   ~        \$(CC) \$(CCOPTS) -c \$(PROJECT).c
  370   ~   ~
  371   \$(PROJECT).pdf: \$(PROJECT).ps
  372   ~        \$(PS2PDF) \$(PROJECT).ps \$(PROJECT).pdf
  373   ~   ~
  374   \$(PROJECT).ps: \$(PROJECT).dvi
  375   ~        \$(DVIPS) \$(DVIPSOPT) \$(PROJECT).dvi -o \$(PROJECT).ps
  376   ~   ~
  377   \$(PROJECT).dvi: \$(PROJECT).w
  378   ~        @make -C figures/
  379   ~        @make -C kochxmpl/
  380   ~        @make verification
  381   ~        \$(CWEAVE) \$(PROJECT)
  382   ~        \$(TEX) \$(PROJECT).tex
  383   ~   ~
  384   verification:
  385   ~        @echo
  386   "==============================================================="
  387   ~        @echo " Verifying the performance of the \$(PROJECT) program on the
  388   Koch"
  389   ~        @echo " snowflake fractal of iteration order 11."
  390   ~        @echo
  391   "==============================================================="
  392   ~        make koch -C koch/
  393   ~        make kochgraphs -C koch/
  394   ~        make fractaldimension -C koch/
  395   ~   ~
  396   tidy:
  397   ~        make -ik clean -C figures/
  398   ~        make -ik clean -C koch/
  399   ~        make -ik clean -C kochxmpl/
  400   ~        -rm -Rf *~ *.o *.exe *.dat *.tgz *.trj *.aux *.log *.idx *.scn *.dvi
  401   ~   ~
  402   clean:
  403   ~        make -ik tidy
  404   ~        -rm -Rf \$(PROJECT) *.c *.pdf *mp *.toc *.tex *.ps
  405   ~   ~
  406   archive:
  407   ~        make -ik tidy
  408   ~        tar --gzip --directory=../ -cf ../\$(PROJECT).tgz \$(PROJECT)
  409   }
  410   \bigskip
  411   This {\tt Makefile} essentially executes two major calls. First, the \CTANGLE\
  412   program parses the \CWEB\ source document {\tt boxcount.w} to extract a \CEE\
  413   source file {\tt boxcount.c} which may be compiled in the usual way using any
  414   \ANSICEE\ conformant compiler. The output source file includes {\tt \#line}
  415   specifications so that any debugging can be done conveniently in terms of
  416   the original \CWEB\ source file.
  417   Second, the \CWEAVE\ program parses the same \CWEB\ source file to extract a
  418   plain \TeX\ source file {\tt boxcount.tex} which may be compiled in the usual
  419   way.
  420   It takes appropriate care of typographic details like page layout and the use
  421   of indentation, italics, boldface, and so on, and it supplies extensive
  422   cross-index information that it gathers automatically.
  423   
  424   After having executed {\tt make} in the same catalogue where the files
  425   {\tt boxcount.w} and {\tt Makefile} are located, one is left with an
  426   executable file {\tt boxcount}, being the ready-to-use compiled program,
  427   and a PostScript file {\tt boxcount.ps}
  428   which contains the full documentation of the program, that is to say the
  429   document you currently are reading.
  430   Notice that on platforms running Windows NT, Windows 2000, Windows ME, or any
  431   other operating system by Microsoft, the executable file will instead
  432   automatically be called {\tt boxcount.exe}.
  433   This convention also applies to programs compiled under the \UNIX-like
  434   environment \CYGWIN.
  435   
  436   \fi
  437   
  438   \N{1}{5}Running the program. The program is entirely controlled by the command
  439   line options supplied when invoking the program. The syntax for executing the
  440   program is\par
  441   \medskip
  442   \line{\hskip 40pt{\tt boxcount [options]}\hfill}\par
  443   \medskip
  444   \noindent
  445   where {\tt options} include the following, given in their long as well as
  446   their short forms (with prefixes `{\tt --}' and `{\tt -}', respectively):
  447   \medskip
  448   \optitem[{\tt --trajectoryfile}, {\tt -i} $\langle${\it trajectory filename}%
  449   $\rangle$]{Specifies the input trajectory of the graph whose fractal
  450   dimension is to be estimated. The input file should describe the trajectory
  451   as two columns of $x$- and $y$-coordinates of the nodes, between which
  452   straight lines will interpolate the trajectory. Unless the boundary of the
  453   computational window is explicitly stated using the {\tt -w} or
  454   {\tt --computationwindow} options, the minimum and maximum $x$- and
  455   $y$-values of the specified trajectory will be used for the automatic
  456   internal computation of the proper computational domain boundaries.}
  457   \optitem[{\tt --outputfile}, {\tt -o [append{\char'174}overwrite]}
  458   $\langle${\it output filename}$\rangle$]{Specifies the base name of the file
  459   to which the calculated output data will be written.
  460   If the {\tt --outputfile} or {\tt -o} option is followed by ``{\tt append}''
  461   the estimate for the fractal dimension will be appended to the file named
  462   $\langle${\it output filename}$\rangle$.dat, which will be created if it
  463   does not exist. If the following word instead is ``{\tt overwrite}'' the
  464   file will, of course, instead be overwritten.}
  465   \optitem[{\tt -w}, {\tt --computationwindow llx} $\langle x_{\rm LL}\rangle$
  466   {\tt lly} $\langle y_{\rm LL}\rangle$ {\tt urx} $\langle x_{\rm UR}\rangle$
  467   {\tt ury} $\langle y_{\rm UR}\rangle$]{This option explicitly specifies the
  468   domain over which the box-counting algorithm will be applied.
  469   By specifying this option, any automatic calculation of computational window
  470   will be neglected.}
  471   \optitem[{\tt --verbose}, {\tt -v}]{Use this option to toggle verbose mode of
  472   the program execution, controlling the amount of information written to
  473   standard terminal output. (Default is off, that is to say quiet mode.)}
  474   \optitem[{\tt --boxmaps}, {\tt -m}]{If this option is present, the \boxcount\
  475   program will generate \METAPOST\ code for maps of the distribution of boxes,
  476   so-called ``boxmaps''. In doing so, also the input trajectory is included in
  477   the graphs. The convention for the naming of the output map files is that
  478   they are saved as $\langle${\it output filename}$\rangle$.$N$.{\tt dat},
  479   where $\langle${\it output filename}$\rangle$ is the base filename as
  480   specified using the {\tt -o} or {\tt --outputfile} option, $N$ is the
  481   automatically appended current level of resolution refinement in the
  482   box-counting (that is to say, indicating the calculation performed for a
  483   $[2^{N}\times2^{N}]$-grid of boxes), and where {\tt dat} is a file suffix
  484   as automatically appended by the program. This option is useful for tracking
  485   the performance of the program, and for verification of the box counting
  486   algorithm.}
  487   \optitem[{\tt --graphsize} $\langle w\rangle$ $\langle h\rangle$]{If the
  488   {\tt -m} or {\tt --boxmaps} option is present at the command line, then
  489   the {\tt --graphsize} option specifies the physical width $w$ and height
  490   $h$ in millimetres of the generated graphs, overriding the default sizes
  491   $w=80\,{\rm mm}$ and $h=80\,{\rm mm}$.}
  492   \optitem[{\tt --minlevel}, {\tt -Nmin} $\langle N_{\rm min}\rangle$]{Specifies
  493   the minimum level $N_{\rm min}$ of grid refinement that will be used in the
  494   evaluation of the estimate of the fractal dimension.
  495   With this option specified, the coarsest level used in the box-counting will
  496   be a $[2^{N_{\rm min}}\times2^{N_{\rm min}}]$-grid of boxes.}
  497   \optitem[{\tt --maxlevel}, {\tt -Nmax} $\langle N_{\rm max}\rangle$]{This
  498   option is similar to the {\tt --minlevel} or {\tt -Nmin} options, with the
  499   difference that it instead specifies the maximum level $N_{\rm max}$ of grid
  500   refinement that will be used in the evaluation of the estimate of the fractal
  501   dimension.
  502   With this option specified, the finest level used in the box-counting will
  503   be a $2^{N_{\rm max}}\times2^{N_{\rm max}}$-grid of boxes.
  504   If this option is omitted, a default value of $N_{\rm max}=10$ will be
  505   used as default.}
  506   
  507   
  508   \fi
  509   
  510   \N{1}{6}Application example: The Koch fractal.
  511   The Koch curve is one of the earliest described fractal curves, appearing
  512   in the 1904 article {\sl Sur une courbe continue sans tangente, obtenue par
  513   une construction g\'eom\'etrique \'el\'ementaire}, written by the Swedish
  514   mathematician Helge von Koch (1870--1924).
  515   In this application example, we employ the ``snowflake'' variant of the Koch
  516   fractal (merely for the sake of its beauty).
  517   The Koch snowflake fractal is constructed recursively using the following
  518   algorithm.
  519   \medskip
  520   \noindent{\bf Algorithm A} ({\it Construction of the Koch snowflake fractal}).
  521   \aitem[{\bf A1.}]{[Create initiator.] Draw three line segments of equal length
  522   so that they form an equilateral triangle.}
  523   \aitem[{\bf A2.}]{[Line division.] Divide each of the line segments into three
  524   segments of equal length.}
  525   \aitem[{\bf A3.}]{[Replace mid segment by triangle.] Draw an equilateral
  526   triangle that has the middle segment from step
  527   one as its base.}
  528   \aitem[{\bf A4.}]{[Remove base of triangle.] Remove the line segment that is
  529   the base of the triangle from step A3.}
  530   \aitem[{\bf A5.}]{[Recursion step.] For each of the line segments remaining,
  531   repeat steps A2 through A4.}
  532   \quad\endalg
  533   \medskip
  534   
  535   The triangle resulting after step A1 is called the {\sl initiator} of the
  536   fractal.
  537   After the first iteration of steps A1--A4, the result should be a shape similar
  538   to the Star of David; this is called the {generator} of the fractal.
  539   The Koch fractal resulting of the infinite iteration of the algorithm of
  540   construction has an infinite length, since each time the steps above are
  541   performed on each line segment of the figure there are four times as many line
  542   segments, the length of each being one-third the length of the segments in the
  543   previous stage.
  544   Hence, the total length of the perimeter of the fractal increases by $4/3$
  545   at each step, and for an initiator of total length $L$ the total length
  546   of the perimeter at the $n$th step of iteration will be $(4/3)^nL$.
  547   The fractal dimension is hence $D=\ln4/\ln3\approx1.26$, being greater
  548   than the dimension of a line ($D=1$) but less than, for example, Peano's
  549   space-filling curve ($D=2$).
  550   The Koch fractal is an example of a curve which is continuous, but not
  551   differentiable anywhere.
  552   The area of the Koch snowflake is $8/5$ that of the initial triangle,
  553   so an infinite perimeter encloses a finite area.
  554   The stepwise construction of the snowflake fractal is illustrated in Fig.~1.
  555   
  556   \bigskip
  557   \centerline{\epsfxsize=54mm\epsfbox{koch/kochgraph-1.1}\hskip2mm
  558   \epsfxsize=54mm\epsfbox{koch/kochgraph-2.1}\hskip2mm
  559   \epsfxsize=54mm\epsfbox{koch/kochgraph-3.1}}
  560   \centerline{\epsfxsize=54mm\epsfbox{koch/kochgraph-4.1}\hskip2mm
  561   \epsfxsize=54mm\epsfbox{koch/kochgraph-5.1}\hskip2mm
  562   \epsfxsize=54mm\epsfbox{koch/kochgraph-6.1}}
  563   \medskip\noindent
  564   {\bf Figure 1.} Construction of the Koch fractal of the ``snowflake'' type,
  565   in this case as inscribed in the unitary circle. For this case, the length
  566   of the initiator is $L=3\sqrt{3}$ while its area is
  567   $A=L^2/(6\sqrt{3})=3\sqrt{3}/2$. For each fractal recursion depth $n$, the
  568   trajectory consists of a total of $3\times4^{(n-1)}$ connected linear segments.
  569   
  570   \vfill\eject
  571   
  572   The data set for the Koch fractal is straightforward to generate by means of
  573   recursion, as for example by using the following compact program (which in
  574   fact was used for the generation of the data sets for the Koch fractals on
  575   the previous page):
  576   \bigskip
  577   {\obeylines\obeyspaces\tt
  578   /*-------------------------------------------------------------------------
  579   {\char'174} File: koch.c [ANSI-C conforming source code]
  580   {\char'174} Created: May 8, 2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
  581   {\char'174} Last modified: May 8, 2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
  582   {\char'174} Compile with: gcc -O2 -g -Wall -pedantic -ansi koch.c -o koch -lm
  583   {\char'174} Description: The KOCH program creates data sets corresponding to
  584   {\char'174}    the Koch fractal, for the purpose of acting as test objects for
  585   {\char'174}    the BOXCOUNT program. The KOCH program is simply executed by
  586   {\char'174}    'koch <N>', where <N> is an integer describing the maximum
  587   {\char'174}    depth of recursion in the generation of the fractal data set.
  588   {\char'174}    If invoked without any arguments, <N>=6 is used as default.
  589   {\char'174}    The generated trajectory is written to standard output.
  590   {\char'174} Copyright (C) 2006 Fredrik Jonsson <fj@phys.soton.ac.uk>
  591   =========================================================================*/
  592   \#include <math.h>
  593   \#include <stdio.h>
  594   extern char *optarg;
  595   \hskip\lineskip
  596   void kochsegment(double xa,double ya,double xb,double yb,
  597   \   int depth,int maxdepth) {\char'173}
  598   \   double xca,yca,xcb,ycb,xcc,ycc;
  599   \   if (depth==maxdepth) {\char'173}
  600   \      fprintf(stdout,"\%2.8f \%2.8f{\char'134}n",xb,yb);
  601   \   {\char'175} else {\char'173}
  602   \      xca=xa+(xb-xa)/3.0;
  603   \      yca=ya+(yb-ya)/3.0;
  604   \      xcb=xb-(xb-xa)/3.0;
  605   \      ycb=yb-(yb-ya)/3.0;
  606   \      xcc=(xa+xb)/2.0-(yb-ya)/(2.0*sqrt(3.0));
  607   \      ycc=(ya+yb)/2.0+(xb-xa)/(2.0*sqrt(3.0));
  608   \      kochsegment(xa,ya,xca,yca,depth+1,maxdepth);
  609   \      kochsegment(xca,yca,xcc,ycc,depth+1,maxdepth);
  610   \      kochsegment(xcc,ycc,xcb,ycb,depth+1,maxdepth);
  611   \      kochsegment(xcb,ycb,xb,yb,depth+1,maxdepth);
  612   \   {\char'175}
  613   {\char'175}
  614   \hskip\lineskip
  615   int main(int argc, char *argv[]) {\char'173}
  616   \   int maxdepth=6;
  617   \   if (argc>1) sscanf(argv[1],"\%d",{\char'046}maxdepth);
  618   \   if (maxdepth>0) {\char'173}
  619   \      fprintf(stdout,"\%2.8f \%2.8f{\char'134}n",0.0,1.0);
  620   \      kochsegment(0.0,1.0,sqrt(3.0)/2.0,-0.5,1,maxdepth);
  621   \      kochsegment(sqrt(3.0)/2.0,-0.5,-sqrt(3.0)/2.0,-0.5,1,maxdepth);
  622   \      kochsegment(-sqrt(3.0)/2.0,-0.5,0.0,1.0,1,maxdepth);
  623   \   {\char'175}
  624   \   return(0);
  625   {\char'175}
  626   }\bigskip
  627   \vfill\eject
  628   
  629   The boxcounting dimension of the Koch snowflake fractal can now be investigated
  630   with assistance of the \boxcount\ program. In the analysis as here presented,
  631   this is done using the following Makefile:
  632   \bigskip
  633   {\obeylines\obeyspaces\tt
  634   \#
  635   \# Makefile designed for use with gcc, MetaPost and plain TeX.
  636   \#
  637   \# Copyright (C) 2002-2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
  638   \#
  639   CC       = gcc
  640   CCOPTS   = -O2 -Wall -ansi -std=iso9899:1990 -pedantic
  641   LNOPTS   = -lm
  642   TEX      = tex
  643   DVIPS    = dvips
  644   METAPOST = mpost
  645   
  646   \#
  647   \# Define path and executable file for the BOXCOUNT program, used for
  648   \# calculating estimates of the box-counting fractal dimension of a
  649   \# trajectory in the plane.
  650   \#
  651   BOXCOUNTPATH = ../
  652   BOXCOUNT = \$(BOXCOUNTPATH)/boxcount
  653   
  654   \#
  655   \# Define path and executable file for the KOCH program, used for generating
  656   \# the trajectory of the Koch snowflake fractal.
  657   \#
  658   KOCHPATH = ../koch/
  659   KOCH = \$(KOCHPATH)/koch
  660   
  661   all: koch kochgen kochtab
  662   
  663   koch:
  664   \        make -C ../koch/
  665   
  666   kochgen:
  667   \        \$(KOCH) 7 > koch.trj
  668   \        \$(BOXCOUNT) --verbose --boxmaps --graphsize 42.0 42.0 {\char'134}
  669   \                --computationwindow llx -1.1 lly -1.1 urx 1.1 ury 1.1 {%
  670   \char'134}
  671   \                --minlevel 3 --maxlevel 8 {\char'134}
  672   \                --trajectoryfile koch.trj --outputfile overwrite koch
  673   \        for k in 03 04 05 06 07 08; do {\char'134}
  674   \        \$(METAPOST) koch.\$\$k.mp ;{\char'134}
  675   \        \$(TEX) -jobname=koch.\$\$k '{\char'134}input epsf{%
  676   \char'134}nopagenumbers{\char'134}
  677   \           {\char'134}centerline{{\char'134}epsfxsize=120mm{\char'134}%
  678   epsfbox{koch.'\$\$k'.1}}{\char'134}bye';{\char'134}
  679   \        \$(DVIPS) -D1200 -E koch.\$\$k -o koch.\$\$k.eps;{\char'134}
  680   \        done
  681   \hskip\lineskip
  682   kochtab:
  683   \        \$(KOCH) 7 > koch.trj
  684   \        \$(BOXCOUNT) --verbose --minlevel 3 --maxlevel 14 {\char'134}
  685   \                --computationwindow llx -1.1 lly -1.1 urx 1.1 ury 1.1 {%
  686   \char'134}
  687   \                --trajectoryfile koch.trj --outputfile overwrite koch
  688   \hskip\lineskip
  689   clean:
  690   \        -rm -Rf koch *~ *.o *.exe *.dat *.mp *.mpx *.trj
  691   \        -rm -Rf *.tex *.aux *.log *.toc *.idx *.scn *.dvi *.1 *.eps
  692   }\bigskip
  693   \vfill\eject
  694   
  695   Having executed the {\tt Makefile} as displayed in the previous page, where a
  696   recursion depth of $n=7$ is used for the generation of the Koch fractal, we are
  697   left with a set of images of the consecutively refined grids in the boxcounting
  698   algorithm, and a table containing the estimates of the boxcounting dimension
  699   of the Koch snowflake fractal.
  700   In Fig.~2 the resulting maps of the boxes used in the boxcounting algorithm
  701   for refinement levels $m=3,4,\ldots,8$ are shown, and as the grid refinement
  702   is finer and finer, the boxes finally will be barely visible in the limited
  703   resolution of computer graphics, on screen as well as in the generated
  704   Encapsulated PostScript code.
  705   
  706   \bigskip
  707   \centerline{\epsfxsize=54mm\epsfbox{kochxmpl/koch-03.1}\hskip2mm
  708   \epsfxsize=54mm\epsfbox{kochxmpl/koch-04.1}\hskip2mm
  709   \epsfxsize=54mm\epsfbox{kochxmpl/koch-05.1}}
  710   \centerline{\epsfxsize=54mm\epsfbox{kochxmpl/koch-06.1}\hskip2mm
  711   \epsfxsize=54mm\epsfbox{kochxmpl/koch-07.1}\hskip2mm
  712   \epsfxsize=54mm\epsfbox{kochxmpl/koch-08.1}}
  713   \medskip\noindent
  714   {\bf Figure 2.} Consecutive steps of refinement $m=3,4,\ldots,8$ of the grid
  715   of boxes covering the input trajectory, in this case a Koch snowflake fractal
  716   of recursion depth $n=7$ (see Fig.~1). At each step $m$ of refinement, a
  717   virtual grid of $2^m\times2^m$ boxes is superimposed on the trajectory
  718   (fractal) and the number $N_m$ of boxes covering the trajectory are counted.
  719   For the Koch snowflake fractal of recursion depth $n=7$, the trajectory
  720   consists of a total of $3\times4^{(7-1)}=12288$ connected linear segments.
  721   \vfill\eject
  722   
  723   The estimated boxcounting dimension $D_m=\ln(N_m)/\ln(2^m)$, with $N_m$ as
  724   previously denoting the number of boxes in a $2^m\times2^m$-grid required to
  725   cover the entire curve, is displayed in Table~A.1.
  726   The values for the estimates could be compared with the analytically obtained
  727   value of $D=\ln4/\ln3\approx1.26$ for the fractal dimension.
  728   However, it should be emphasized that the box counting dimension here just
  729   is an estimate of one definition of the fractal dimension, which in many cases
  730   do not equal to other measures, and that we in the computation of the estimate
  731   always will arrive at a truncated result due to a limited precision and a
  732   limited amount of memory resources and computational time.
  733   As can be seen in the table, the initial estimates at lower resolutions are
  734   pretty crude, but in the final estimate we nevertheless end up with the box
  735   counting estimate of the fractal dimension as 1.29, which is reasonably close
  736   to the analytically obtained value of $D\approx1.26$.
  737   Of course, further refinements such as Richardson extrapolation could be
  738   applied to increase the accuracy, but this is outside of the scope of the
  739   \boxcount\ program, which only serves to provide the raw, basic algorithm
  740   of boxcounting.\footnote{$\dagger$}{For discussions on
  741   different definitions of the fractal dimension, see
  742   the English Wikipedia section on the Minkowski--Bouligand dimension,
  743   {\tt http://en.wikipedia.org/wiki/Minkowski-Bouligand{\char'137}dimension}.}
  744   \bigskip
  745   \noindent{\bf Table A.1} Boxcounting estimates of the fractal dimension of
  746   the Koch snowflake fractal of recursion order $n=7$. In the table, $m$ is the
  747   refinement order as indicated in the graphs in Fig.~2, $N_m$ is the number of
  748   covering boxes counted at refinement level $m$, and $D_m=\ln(N_m)/\ln(2^m)$
  749   is the estimate of the boxcounting dimension.
  750   \smallskip
  751   \noindent{\hrule width 328pt}\vskip1pt
  752   \beginvrulealign
  753   % \tabskip=0pt
  754   \halign{
  755   \hbox to 72pt{\hfil#\hfil}% first column centered
  756   &\hbox to 128pt{\hfil#\hfil}% second column is centered
  757   &\hbox to 128pt{\hfil#\hfil}\cr% third column is centered
  758   \noalign{{\hrule width 328pt}\vskip 2pt}
  759   $m$ & $N_m$ & $D_m=\ln(N_m)/\ln(2^m)$\cr
  760   \noalign{\vskip 1pt{\hrule width 328pt}\vskip 1.4pt}
  761   3 &     44 & 1.8198\cr
  762   4 &     96 & 1.6462\cr
  763   5 &    196 & 1.5229\cr
  764   6 &    504 & 1.4962\cr
  765   7 &   1180 & 1.4578\cr
  766   8 &   2856 & 1.4350\cr
  767   9 &   6844 & 1.4156\cr
  768   10 &  15620 & 1.3931\cr
  769   11 &  32320 & 1.3618\cr
  770   12 &  66200 & 1.3345\cr
  771   13 & 133600 & 1.3098\cr
  772   14 & 268804 & 1.2883\cr
  773   }\endvrulealign
  774   \vskip-9.55pt
  775   {\hskip-.4pt\hrule width 328pt}\vskip 1pt
  776   \noindent{\hskip-.2pt\hrule width 328pt}
  777   
  778   
  779   \fi
  780   
  781   \N{1}{7}The main program. Here follows the general outline of the main program.
  782   
  783   \Y\B\X8:Library inclusions\X\6
  784   \X9:Global definitions\X\6
  785   \X10:Global variables\X\6
  786   \X23:Subroutines\X\7
  787   \&{int} \\{main}(\&{int} \\{argc}${},\39{}$\&{char} ${}{*}\\{argv}[\,]){}$\1\1%
  788   \2\2\6
  789   ${}\{{}$\1\6
  790   \X11:Declaration of local variables\X\6
  791   \X12:Initialize variables\X\6
  792   \X13:Parse command line for parameters\X\6
  793   \X14:Display starting time of program execution\X\6
  794   \X15:Load input trajectory from file\X\6
  795   \X16:Open file for output of logarithmic estimate\X\6
  796   \X17:Extract boundary of global window of computation from input trajectory\X\6
  797   \X18:Get number of boxes covering the trajectory for all levels of refinement
  798   in resolution\X\6
  799   \X19:Compute the logarithmic estimate of the fractal dimension\X\6
  800   \X20:Save or append the logarithmic estimate to output file\X\6
  801   \X21:Close file for output of logarithmic estimate\X\6
  802   \X22:Display elapsed execution time\X\6
  803   \&{return} (\.{SUCCESS});\6
  804   \4${}\}{}$\2\par
  805   \fi
  806   
  807   \M{8}Library dependencies.
  808   The standard \ANSICEE\ libraries included in this program are:\par
  809   \varitem[{\tt math.h}]{For access to common mathematical functions.}
  810   \varitem[{\tt stdio.h}]{For file access and any block involving \PB{%
  811   \\{fprintf}}.}
  812   \varitem[{\tt stdlib.h}]{For access to the \PB{\\{exit}} function.}
  813   \varitem[{\tt string.h}]{For string manipulation, \PB{\\{strcpy}}, \PB{%
  814   \\{strcmp}} etc.}
  815   \varitem[{\tt ctype.h}]{For access to the \PB{\\{isalnum}} function.}
  816   
  817   \Y\B\4\X8:Library inclusions\X${}\E{}$\6
  818   \8\#\&{include} \.{<math.h>}\6
  819   \8\#\&{include} \.{<time.h>}\6
  820   \8\#\&{include} \.{<stdio.h>}\6
  821   \8\#\&{include} \.{<stdlib.h>}\6
  822   \8\#\&{include} \.{<string.h>}\6
  823   \8\#\&{include} \.{<ctype.h>}\par
  824   \U7.\fi
  825   
  826   \M{9}Global definitions.
  827   There are just a few global definitions present in the \boxcount\ program:\par
  828   \varitem[{\tt VERSION}]{The current program revision number.}
  829   \varitem[{\tt COPYRIGHT}]{The copyright banner.}
  830   \varitem[{\tt SUCCESS}]{The return code for successful program termination.}
  831   \varitem[{\tt FAILURE}]{The return code for unsuccessful program termination.}
  832   \varitem[{\tt NCHMAX}]{The maximum number of characters allowed in strings for
  833   storing file names, including path.}
  834   \varitem[{\tt APPEND}]{Code for the flag \PB{\\{output\_mode}}, used to
  835   determine if
  836   output should append to an existing file.}
  837   \varitem[{\tt OVERWRITE}]{Code for the flag \PB{\\{output\_mode}}, used to
  838   determine
  839   if output should overwrite an existing file.}
  840   
  841   \Y\B\4\X9:Global definitions\X${}\E{}$\6
  842   \8\#\&{define} \.{VERSION} \5\.{"1.6"}\6
  843   \8\#\&{define} \.{COPYRIGHT} \5\.{"Copyright\ (C)\ 2006-}\)\.{2011,\ Fredrik\
  844   Jonsso}\)\.{n"}\6
  845   \8\#\&{define} \.{SUCCESS} \5(\T{0})\6
  846   \8\#\&{define} \.{FAILURE} \5(\T{1})\6
  847   \8\#\&{define} \.{NCHMAX} \5(\T{256})\6
  848   \8\#\&{define} \.{APPEND} \5\T{1}\6
  849   \8\#\&{define} \.{OVERWRITE} \5\T{2}\par
  850   \U7.\fi
  851   
  852   \M{10}Declaration of global variables. The only global variables allowed in
  853   my programs are \PB{\\{optarg}}, which is a pointer to the the string of
  854   characters
  855   that specified the call from the command line, and \PB{\\{progname}}, which
  856   simply
  857   is a pointer to the string containing the name of the program, as it was
  858   invoked from the command line.
  859   
  860   \Y\B\4\X10:Global variables\X${}\E{}$\6
  861   \&{extern} \&{char} ${}{*}\\{optarg};{}$\6
  862   \&{char} ${}{*}\\{progname}{}$;\par
  863   \U7.\fi
  864   
  865   \M{11}Declaration of local variables of the \PB{\\{main}} program.
  866   In \CWEB\ one has the option of adding variables along the program, for example
  867   by locally adding temporary variables related to a given sub-block of code.
  868   However, the philosophy in the \boxcount\ program is to keep all variables of
  869   the \PB{\\{main}} section collected together, so as to simplify tasks as, for
  870   example,
  871   tracking down a given variable type definition.
  872   The local variables of the program are as follows:
  873   \medskip
  874   \varitem[\PB{\\{verbose}}]{Boolean flag which, if being nonzero, tells the
  875   program
  876   to display information at terminal output during program execution.}
  877   \varitem[\PB{\\{user\_set\_compwin}}]{Boolean flag which indicates whether the
  878   user has
  879   explicitly stated a window of computation or not.}
  880   \varitem[\PB{\\{output\_mode}}]{Variable which states whether the estimate of
  881   fractal dimension should be appended to an existing file which is
  882   created if it does not exist (for $\PB{\\{output\_mode}}=\PB{\.{APPEND}}$), or
  883   if
  884   the data should overwrite already existing data in the file (for
  885   $\PB{\\{output\_mode}}=\PB{\.{OVERWRITE}}$).}
  886   \varitem[\PB{\\{make\_boxmaps}}]{If nonzero, then graphs showing the
  887   distribution
  888   of covering boxes will be generated.}
  889   \varitem[\PB{${*}\\{num\_boxes}$}]{Pointer to array holding the number of boxes
  890   at
  891   various depths of division.}
  892   \varitem[\PB{\\{initime}}]{The time at which the \boxcount\ program was
  893   initialized.}
  894   \varitem[\PB{\\{now}}]{Dummy variable for extraction of current time from the
  895   system
  896   clock.}
  897   \varitem[\PB{\\{mm}}]{The number of points $M$ in the input trajectory. This is
  898   the
  899   length of the vectors \PB{\\{x\_traj}} and \PB{\\{y\_traj}} as described
  900   below.}
  901   \varitem[\PB{\\{nn}}]{Gives the number of boxes in the $x$- or $y$-directions
  902   as
  903   $2^N$.}
  904   \varitem[\PB{\\{nnmax}}]{The maximum refinement depth $N_{\rm max}$ of $N$.}
  905   \varitem[\PB{\\{nnmin}}]{The minimum refinement depth $N_{\rm min}$ of $N$.}
  906   \varitem[\PB{\\{global\_llx}}, \PB{\\{global\_lly}}]{Lower-left coordinates of
  907   global window
  908   of computation.}
  909   \varitem[\PB{\\{global\_urx}}, \PB{\\{global\_ury}}]{Upper-right coordinates of
  910   global window
  911   of computation.}
  912   \varitem[\PB{${*}\\{x\_traj}$},\PB{${*}\\{y\_traj}$}]{Vectors containing the
  913   $x$- and $y$-values of
  914   the coordinates along the input trajectory.}
  915   \varitem[\PB{${*}\|x$},\PB{${*}\|y$}]{Variables for keeping the refinement and
  916   number of boxes.
  917   (This needs to be changed as they are easily confused with $x$ and $y$
  918   coordinates of the trajectory.)}
  919   \varitem[\PB{${*}\\{trajectory\_file}$}]{Input file pointer, for reading the
  920   trajectory
  921   whose fractal dimension is to be estimated.}
  922   \varitem[\PB{${*}\\{frac\_estimate\_file}$}]{Output file pointer, for saving
  923   the estimated
  924   fractal dimension of the input trajectory.}
  925   \varitem[\PB{${*}\\{boxmap\_file}$}]{Output file pointer, for saving maps of
  926   the spatial
  927   locations of the boxes covering the trajectory.}
  928   \varitem[\PB{\\{boxgraph\_width}}]{The physical width in millimetres of the
  929   generated
  930   \MP\ boxmap graphs.}
  931   \varitem[\PB{\\{boxgraph\_height}}]{The physical height in millimetres of the
  932   generated
  933   \MP\ boxmap graphs.}
  934   \varitem[\PB{\\{trajectory\_filename}}]{String for keeping the name of the file
  935   containing the input trajectory.}
  936   \varitem[\PB{\\{output\_filename}}]{String for keeping the base name of the set
  937   of
  938   output files.}
  939   \varitem[\PB{\\{frac\_estimate\_filename}}]{String keeping the name of the file
  940   to which the estimate of the fractal dimension will be saved.}
  941   \varitem[\PB{\\{boxmap\_filename}}]{String for keeping the name of the file to
  942   which \METAPOST\ code for maps of the spatial distribution of boxes will be
  943   written.}
  944   \varitem[\PB{\\{no\_arg}}]{Variable for extracting the number of input
  945   arguments
  946   present on the command line as the program is executed.}
  947   \varitem[\PB{\|i}]{Dummy variable used in loops when reading the input
  948   trajectory.}
  949   \varitem[\PB{${*}\\{fracdimen\_estimates}$}]{Vector keeping the values
  950   characterizing
  951   estimated fractal dimension for various orders of $N$.}
  952   \varitem[\PB{$\\{ave},\\{adev},\\{sdev}$}]{The average, average deviation, and
  953   standard
  954   deviation of the estimated fractal dimension for various orders of
  955   refinement $N$, as stored in \PB{\\{fracdimen\_estimates}}.}
  956   \varitem[\PB{$\\{var},\\{skew},\\{curt}$}]{The variance, skewness, and kurtosis
  957   of the
  958   estimated fractal dimension for various orders of refinement $N$, as
  959   stored in \PB{\\{fracdimen\_estimates}}.}
  960   \medskip
  961   \noindent
  962   Generally in this program, the maximum number of characters a file name
  963   string can contain is \PB{\.{NCHMAX}}, as defined in the definitions section
  964   of the program.
  965   
  966   \Y\B\4\X11:Declaration of local variables\X${}\E{}$\6
  967   \&{short} \\{verbose}${},{}$ \\{user\_set\_compwin}${},{}$ \\{output%
  968   \_mode}${},{}$ \\{make\_boxmaps};\6
  969   \&{long} \&{int} ${}{*}\\{num\_boxes};{}$\6
  970   \&{time\_t} \\{initime};\6
  971   \&{time\_t} \\{now};\6
  972   \&{long} \\{mm}${},{}$ \\{nn}${},{}$ \\{nnmin}${},{}$ \\{nnmax};\6
  973   \&{double} \\{global\_llx}${},{}$ \\{global\_lly}${},{}$ \\{global\_urx}${},{}$
  974   \\{global\_ury};\6
  975   \&{double} ${}{*}\\{x\_traj},{}$ ${}{*}\\{y\_traj},{}$ ${}{*}\|x,{}$ ${}{*}%
  976   \|y;{}$\6
  977   \&{FILE} ${}{*}\\{trajectory\_file},{}$ ${}{*}\\{frac\_estimate\_file},{}$
  978   ${}{*}\\{boxmap\_file};{}$\6
  979   \&{char} \\{trajectory\_filename}[\.{NCHMAX}]${},{}$ \\{output\_filename}[%
  980   \.{NCHMAX}]${},{}$ \\{frac\_estimate\_filename}[\.{NCHMAX}];\6
  981   \&{char} \\{boxmap\_filename}[\.{NCHMAX}];\6
  982   \&{double} \\{boxgraph\_width}${},{}$ \\{boxgraph\_height};\6
  983   \&{int} \\{no\_arg};\6
  984   \&{long} \|i;\6
  985   \&{double} ${}{*}\\{fracdimen\_estimates},{}$ \\{ave}${},{}$ \\{adev}${},{}$ %
  986   \\{sdev}${},{}$ \\{var}${},{}$ \\{skew}${},{}$ \\{curt};\par
  987   \U7.\fi
  988   
  989   \M{12}Initialization of variables.
  990   
  991   \Y\B\4\X12:Initialize variables\X${}\E{}$\6
  992   $\\{verbose}\K\T{0}{}$;\C{ Verbose mode is off by default }\6
  993   ${}\\{user\_set\_compwin}\K\T{0}{}$;\C{ Before the command-line is parsed,
  994   nothing is known of these settings }\6
  995   ${}\\{output\_mode}\K\.{OVERWRITE}{}$;\C{ default mode is to overwrite existing
  996   files }\6
  997   ${}\\{make\_boxmaps}\K\T{0}{}$;\C{ Default is to not create graphs of the box
  998   distributions }\6
  999   ${}\\{nnmin}\K\T{0}{}$;\C{ Default value for $N_{\rm min}$ }\6
 1000   ${}\\{nnmax}\K\T{10}{}$;\C{ Default value for $N_{\rm max}$ }\6
 1001   ${}\\{strcpy}(\\{output\_filename},\39\.{"out"}){}$;\C{ Default output file
 1002   basename. }\6
 1003   ${}\\{strcpy}(\\{trajectory\_filename},\39\.{""}){}$;\C{ To indicate that no
 1004   filename has been set yet. }\6
 1005   ${}\\{boxgraph\_width}\K\T{80.0}{}$;\C{ Default graph width in millimetres }\6
 1006   ${}\\{boxgraph\_height}\K\T{56.0}{}$;\C{ Default graph height in millimetres }\6
 1007   ${}\\{trajectory\_file}\K\NULL;{}$\6
 1008   ${}\\{frac\_estimate\_file}\K\NULL;{}$\6
 1009   ${}\\{boxmap\_file}\K\NULL;{}$\6
 1010   ${}\\{initime}\K\\{time}(\NULL){}$;\C{ Time of initialization of the program. }%
 1011   \par
 1012   \U7.\fi
 1013   
 1014   \M{13}Parsing command line options. All input parameters are passed to the
 1015   program through command line options and arguments to the program.
 1016   The syntax of command line options is listed
 1017   whenever the program is invoked without any options, or whenever any of the
 1018   {\tt --help} or {\tt -h} options are specified at startup.
 1019   
 1020   \Y\B\4\X13:Parse command line for parameters\X${}\E{}$\6
 1021   ${}\{{}$\1\6
 1022   ${}\\{progname}\K\\{strip\_away\_path}(\\{argv}[\T{0}]);{}$\6
 1023   ${}\\{fprintf}(\\{stdout},\39\.{"This\ is\ \%s\ v.\%s.\ \%s}\)\.{\\n"},\39%
 1024   \\{progname},\39\.{VERSION},\39\.{COPYRIGHT});{}$\6
 1025   ${}\\{no\_arg}\K\\{argc};{}$\6
 1026   \&{while} ${}(\MM\\{argc}){}$\5
 1027   ${}\{{}$\1\6
 1028   \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-o"})\V\R%
 1029   \\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--outputfile"})){}$\5
 1030   ${}\{{}$\1\6
 1031   ${}\MM\\{argc};{}$\6
 1032   \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"append"})\V\R%
 1033   \\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"a"})){}$\5
 1034   ${}\{{}$\1\6
 1035   ${}\\{output\_mode}\K\.{APPEND};{}$\6
 1036   \4${}\}{}$\2\6
 1037   \&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
 1038   \.{"overwrite"})\V\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"o"})){}$\5
 1039   ${}\{{}$\1\6
 1040   ${}\\{output\_mode}\K\.{OVERWRITE};{}$\6
 1041   \4${}\}{}$\2\6
 1042   \&{else}\5
 1043   ${}\{{}$\1\6
 1044   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ '-o'\ o}\)\.{r\
 1045   '--outputfile'\ opt}\)\.{ion!"},\39\\{progname});{}$\6
 1046   \\{exit}(\.{FAILURE});\6
 1047   \4${}\}{}$\2\6
 1048   ${}\MM\\{argc};{}$\6
 1049   ${}\\{strcpy}(\\{output\_filename},\39\\{argv}[\\{no\_arg}-\\{argc}]);{}$\6
 1050   \4${}\}{}$\2\6
 1051   \&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-w"})\V%
 1052   \R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--computationwindow}\)%
 1053   \.{"})){}$\5
 1054   ${}\{{}$\1\6
 1055   ${}\\{user\_set\_compwin}\K\T{1};{}$\6
 1056   ${}\MM\\{argc};{}$\6
 1057   \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"llx"})){}$\5
 1058   ${}\{{}$\1\6
 1059   ${}\MM\\{argc};{}$\6
 1060   \&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
 1061   \\{global\_llx})){}$\5
 1062   ${}\{{}$\1\6
 1063   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ parsin}\)\.{g\ lower-left\
 1064   x-value}\)\.{.\\n"},\39\\{progname});{}$\6
 1065   \\{exit}(\.{FAILURE});\6
 1066   \4${}\}{}$\2\6
 1067   \4${}\}{}$\2\6
 1068   \&{else}\5
 1069   ${}\{{}$\1\6
 1070   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ comput}\)\.{ation\ window\
 1071   option\\}\)\.{n"},\39\\{progname});{}$\6
 1072   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Expecting\ 'llx'}\)\.{\\n"},\39%
 1073   \\{progname});{}$\6
 1074   \\{exit}(\.{FAILURE});\6
 1075   \4${}\}{}$\2\6
 1076   ${}\MM\\{argc};{}$\6
 1077   \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"lly"})){}$\5
 1078   ${}\{{}$\1\6
 1079   ${}\MM\\{argc};{}$\6
 1080   \&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
 1081   \\{global\_lly})){}$\5
 1082   ${}\{{}$\1\6
 1083   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ parsin}\)\.{g\ lower-left\
 1084   y-value}\)\.{.\\n"},\39\\{progname});{}$\6
 1085   \\{exit}(\.{FAILURE});\6
 1086   \4${}\}{}$\2\6
 1087   \4${}\}{}$\2\6
 1088   \&{else}\5
 1089   ${}\{{}$\1\6
 1090   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ comput}\)\.{ation\ window\
 1091   option\\}\)\.{n"},\39\\{progname});{}$\6
 1092   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Expecting\ 'lly'}\)\.{\\n"},\39%
 1093   \\{progname});{}$\6
 1094   \\{exit}(\.{FAILURE});\6
 1095   \4${}\}{}$\2\6
 1096   ${}\MM\\{argc};{}$\6
 1097   \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"urx"})){}$\5
 1098   ${}\{{}$\1\6
 1099   ${}\MM\\{argc};{}$\6
 1100   \&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
 1101   \\{global\_urx})){}$\5
 1102   ${}\{{}$\1\6
 1103   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ parsin}\)\.{g\ lower-left\
 1104   x-value}\)\.{.\\n"},\39\\{progname});{}$\6
 1105   \\{exit}(\.{FAILURE});\6
 1106   \4${}\}{}$\2\6
 1107   \4${}\}{}$\2\6
 1108   \&{else}\5
 1109   ${}\{{}$\1\6
 1110   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ comput}\)\.{ation\ window\
 1111   option\\}\)\.{n"},\39\\{progname});{}$\6
 1112   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Expecting\ 'urx'}\)\.{\\n"},\39%
 1113   \\{progname});{}$\6
 1114   \\{exit}(\.{FAILURE});\6
 1115   \4${}\}{}$\2\6
 1116   ${}\MM\\{argc};{}$\6
 1117   \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"ury"})){}$\5
 1118   ${}\{{}$\1\6
 1119   ${}\MM\\{argc};{}$\6
 1120   \&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
 1121   \\{global\_ury})){}$\5
 1122   ${}\{{}$\1\6
 1123   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ parsin}\)\.{g\ lower-left\
 1124   y-value}\)\.{.\\n"},\39\\{progname});{}$\6
 1125   \\{exit}(\.{FAILURE});\6
 1126   \4${}\}{}$\2\6
 1127   \4${}\}{}$\2\6
 1128   \&{else}\5
 1129   ${}\{{}$\1\6
 1130   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ comput}\)\.{ation\ window\
 1131   option\\}\)\.{n"},\39\\{progname});{}$\6
 1132   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Expecting\ 'ury'}\)\.{\\n"},\39%
 1133   \\{progname});{}$\6
 1134   \\{exit}(\.{FAILURE});\6
 1135   \4${}\}{}$\2\6
 1136   \4${}\}{}$\2\6
 1137   \&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-i"})\V%
 1138   \R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--trajectoryfile"})){}$\5
 1139   ${}\{{}$\1\6
 1140   ${}\MM\\{argc};{}$\6
 1141   ${}\\{strcpy}(\\{trajectory\_filename},\39\\{argv}[\\{no\_arg}-\\{argc}]);{}$\6
 1142   \4${}\}{}$\2\6
 1143   \&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-v"})\V%
 1144   \R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--verbose"})){}$\5
 1145   ${}\{{}$\1\6
 1146   ${}\\{verbose}\K(\\{verbose}\?\T{0}:\T{1});{}$\6
 1147   \4${}\}{}$\2\6
 1148   \&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-h"})\V%
 1149   \R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--help"})){}$\5
 1150   ${}\{{}$\1\6
 1151   \\{showsomehelp}(\,);\6
 1152   \\{exit}(\.{SUCCESS});\6
 1153   \4${}\}{}$\2\6
 1154   \&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-m"})\V%
 1155   \R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--boxmaps"})){}$\5
 1156   ${}\{{}$\1\6
 1157   ${}\\{make\_boxmaps}\K\T{1};{}$\6
 1158   \4${}\}{}$\2\6
 1159   \&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
 1160   \.{"--graphsize"})){}$\5
 1161   ${}\{{}$\1\6
 1162   ${}\MM\\{argc};{}$\6
 1163   \&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
 1164   \\{boxgraph\_width})){}$\5
 1165   ${}\{{}$\1\6
 1166   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ width\ }\)\.{of\
 1167   '--graphsize'\ opt}\)\.{ion.\\n"},\39\\{progname});{}$\6
 1168   \\{exit}(\.{FAILURE});\6
 1169   \4${}\}{}$\2\6
 1170   ${}\MM\\{argc};{}$\6
 1171   \&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
 1172   \\{boxgraph\_height})){}$\5
 1173   ${}\{{}$\1\6
 1174   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ height}\)\.{\ of\
 1175   '--graphsize'\ op}\)\.{tion.\\n"},\39\\{progname});{}$\6
 1176   \\{exit}(\.{FAILURE});\6
 1177   \4${}\}{}$\2\6
 1178   \4${}\}{}$\2\6
 1179   \&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
 1180   \.{"--minlevel"})\V\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
 1181   \.{"-Nmin"})){}$\5
 1182   ${}\{{}$\1\6
 1183   ${}\MM\\{argc};{}$\6
 1184   \&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%ld"},\39{\AND}%
 1185   \\{nnmin})){}$\5
 1186   ${}\{{}$\1\6
 1187   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ '--min}\)\.{level'\ or\
 1188   '-Nmin'\ op}\)\.{tion.\\n"},\39\\{progname});{}$\6
 1189   \\{exit}(\.{FAILURE});\6
 1190   \4${}\}{}$\2\6
 1191   \4${}\}{}$\2\6
 1192   \&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
 1193   \.{"--maxlevel"})\V\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
 1194   \.{"-Nmax"})){}$\5
 1195   ${}\{{}$\1\6
 1196   ${}\MM\\{argc};{}$\6
 1197   \&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%ld"},\39{\AND}%
 1198   \\{nnmax})){}$\5
 1199   ${}\{{}$\1\6
 1200   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ '--max}\)\.{level'\ or\
 1201   '-Nmax'\ op}\)\.{tion.\\n"},\39\\{progname});{}$\6
 1202   \\{exit}(\.{FAILURE});\6
 1203   \4${}\}{}$\2\6
 1204   \4${}\}{}$\2\6
 1205   \&{else}\5
 1206   ${}\{{}$\1\6
 1207   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Specified\ optio}\)\.{n\ '\%s'\ invalid!%
 1208   \\n"},\39\\{progname},\39\\{argv}[\\{no\_arg}-\\{argc}]);{}$\6
 1209   \\{showsomehelp}(\,);\6
 1210   \\{exit}(\.{FAILURE});\6
 1211   \4${}\}{}$\2\6
 1212   \4${}\}{}$\2\6
 1213   \4${}\}{}$\2\par
 1214   \U7.\fi
 1215   
 1216   \M{14}Display starting time of program execution.
 1217   
 1218   \Y\B\4\X14:Display starting time of program execution\X${}\E{}$\6
 1219   ${}\{{}$\1\6
 1220   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Program\ executi}\)\.{on\ started\ %
 1221   \%s"},\39\\{progname},\39\\{ctime}({\AND}\\{initime}));{}$\6
 1222   \4${}\}{}$\2\par
 1223   \U7.\fi
 1224   
 1225   \M{15}Load input trajectory from file. This is the section where the trajectory
 1226   to be analyzed is loaded into the memory, and where the number~$M$ of input
 1227   coordinates present in the input data is determined by a single call to
 1228   the subroutine \PB{\\{num\_coordinate\_pairs}(\,)}.
 1229   
 1230   \Y\B\4\X15:Load input trajectory from file\X${}\E{}$\6
 1231   ${}\{{}$\1\6
 1232   \&{if} ${}(\R\\{strcmp}(\\{trajectory\_filename},\39\.{""})){}$\5
 1233   ${}\{{}$\1\6
 1234   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ No\ input\ trajec}\)\.{tory\ specified!%
 1235   \\n"},\39\\{progname});{}$\6
 1236   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Please\ use\ the\ }\)%
 1237   \.{'--trajectoryfile'\ o}\)\.{ption.\\n"},\39\\{progname});{}$\6
 1238   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Use\ '--help'\ op}\)\.{tion\ to\ display%
 1239   \ a\ he}\)\.{lp\ message.\\n"},\39\\{progname});{}$\6
 1240   \\{exit}(\.{FAILURE});\6
 1241   \4${}\}{}$\2\6
 1242   \&{if} ${}((\\{trajectory\_file}\K\\{fopen}(\\{trajectory\_filename},\39%
 1243   \.{"r"}))\E\NULL){}$\5
 1244   ${}\{{}$\1\6
 1245   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Could\ not\ open\ }\)\.{\%s\ for\
 1246   loading\ traje}\)\.{ctory!\\n"},\39\\{progname},\39\\{trajectory%
 1247   \_filename});{}$\6
 1248   \\{exit}(\.{FAILURE});\6
 1249   \4${}\}{}$\2\6
 1250   ${}\\{mm}\K\\{num\_coordinate\_pairs}(\\{trajectory\_file});{}$\6
 1251   ${}\\{x\_traj}\K\\{dvector}(\T{1},\39\\{mm});{}$\6
 1252   ${}\\{y\_traj}\K\\{dvector}(\T{1},\39\\{mm});{}$\6
 1253   \&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{mm};{}$ ${}\|i\PP){}$\5
 1254   ${}\{{}$\1\6
 1255   ${}\\{fscanf}(\\{trajectory\_file},\39\.{"\%lf"},\39{\AND}\\{x\_traj}[\|i]){}$;%
 1256   \C{ scan $x$-coordinate }\6
 1257   ${}\\{fscanf}(\\{trajectory\_file},\39\.{"\%lf"},\39{\AND}\\{y\_traj}[\|i]){}$;%
 1258   \C{ scan $y$-coordinate }\6
 1259   \4${}\}{}$\2\6
 1260   \\{fclose}(\\{trajectory\_file});\6
 1261   \&{if} (\\{verbose})\5
 1262   ${}\{{}$\1\6
 1263   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Loaded\ \%ld\ coor}\)\.{dinate\ pairs\
 1264   from\ fi}\)\.{le\ '\%s'.\\n"},\39\\{progname},\39\\{mm},\39\\{trajectory%
 1265   \_filename});{}$\6
 1266   \4${}\}{}$\2\6
 1267   \4${}\}{}$\2\par
 1268   \U7.\fi
 1269   
 1270   \M{16}Opening files for output by the program.
 1271   
 1272   \Y\B\4\X16:Open file for output of logarithmic estimate\X${}\E{}$\6
 1273   ${}\{{}$\1\6
 1274   \&{if} ${}(\R\\{strcmp}(\\{output\_filename},\39\.{""})){}$\5
 1275   ${}\{{}$\1\6
 1276   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ No\ output\ base\ }\)\.{name\ specified!%
 1277   \\n"},\39\\{progname});{}$\6
 1278   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Please\ use\ the\ }\)\.{'--outputfile'\
 1279   optio}\)\.{n.\\n"},\39\\{progname});{}$\6
 1280   \\{exit}(\.{FAILURE});\6
 1281   \4${}\}{}$\2\6
 1282   ${}\\{sprintf}(\\{frac\_estimate\_filename},\39\.{"\%s.dat"},\39\\{output%
 1283   \_filename});{}$\6
 1284   \&{if} ${}(\\{output\_mode}\E\.{APPEND}){}$\5
 1285   ${}\{{}$\1\6
 1286   \&{if} ${}((\\{frac\_estimate\_file}\K\\{fopen}(\\{frac\_estimate\_filename},%
 1287   \39\.{"a"}))\E\NULL){}$\5
 1288   ${}\{{}$\1\6
 1289   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Could\ not\ open\ }\)\.{\%s\ for\
 1290   output!\\n"},\39\\{progname},\39\\{frac\_estimate\_filename});{}$\6
 1291   \\{exit}(\.{FAILURE});\6
 1292   \4${}\}{}$\2\6
 1293   ${}\\{fseek}(\\{frac\_estimate\_file},\39\T{0\$L},\39\.{SEEK\_END}){}$;\C{ set
 1294   file pointer to the end of the file }\6
 1295   \4${}\}{}$\2\6
 1296   \&{else} \&{if} ${}(\\{output\_mode}\E\.{OVERWRITE}){}$\5
 1297   ${}\{{}$\1\6
 1298   \&{if} ${}((\\{frac\_estimate\_file}\K\\{fopen}(\\{frac\_estimate\_filename},%
 1299   \39\.{"w"}))\E\NULL){}$\5
 1300   ${}\{{}$\1\6
 1301   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Could\ not\ open\ }\)\.{\%s\ for\
 1302   loading\ traje}\)\.{ctory!\\n"},\39\\{progname},\39\\{frac\_estimate%
 1303   \_filename});{}$\6
 1304   \\{exit}(\.{FAILURE});\6
 1305   \4${}\}{}$\2\6
 1306   ${}\\{fseek}(\\{frac\_estimate\_file},\39\T{0\$L},\39\.{SEEK\_SET}){}$;\C{ set
 1307   file pointer to the beginning of the file }\6
 1308   \4${}\}{}$\2\6
 1309   \&{else}\5
 1310   ${}\{{}$\1\6
 1311   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error:\ Output\ m}\)\.{ode\ (output%
 1312   \_mode)\ un}\)\.{defined!"},\39\\{progname});{}$\6
 1313   \\{exit}(\.{FAILURE});\6
 1314   \4${}\}{}$\2\6
 1315   \4${}\}{}$\2\par
 1316   \U7.\fi
 1317   
 1318   \M{17}Extract global window of computation. This block will only be executed if
 1319   there was no explicit computational window stated at startup of the program
 1320   (that is to say, with the {\tt -w} or {\tt --computationwindow} option).
 1321   In order to determine the minimum area covering the entire input trajectory,
 1322   this section scans sequentially through the trajectory and finds the minimum
 1323   and maximum of the $x$- and $y$-coordinates.
 1324   These values then form the lower-left and upper-right corner coordinates
 1325   $(\PB{\\{global\_llx}},\PB{\\{global\_lly}})$ and $(\PB{\\{global\_urx}},\PB{%
 1326   \\{global\_ury}})$.
 1327   
 1328   \Y\B\4\X17:Extract boundary of global window of computation from input
 1329   trajectory\X${}\E{}$\6
 1330   ${}\{{}$\1\6
 1331   \&{if} ${}(\R\\{user\_set\_compwin}){}$\5
 1332   ${}\{{}$\1\6
 1333   ${}\\{global\_llx}\K\\{x\_traj}[\T{1}];{}$\6
 1334   ${}\\{global\_lly}\K\\{y\_traj}[\T{1}];{}$\6
 1335   ${}\\{global\_urx}\K\\{global\_llx};{}$\6
 1336   ${}\\{global\_ury}\K\\{global\_lly};{}$\6
 1337   \&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{mm};{}$ ${}\|i\PP){}$\5
 1338   ${}\{{}$\1\6
 1339   \&{if} ${}(\\{x\_traj}[\|i]>\\{global\_urx}){}$\1\5
 1340   ${}\\{global\_urx}\K\\{x\_traj}[\|i];{}$\2\6
 1341   \&{if} ${}(\\{y\_traj}[\|i]>\\{global\_ury}){}$\1\5
 1342   ${}\\{global\_ury}\K\\{y\_traj}[\|i];{}$\2\6
 1343   \&{if} ${}(\\{x\_traj}[\|i]<\\{global\_llx}){}$\1\5
 1344   ${}\\{global\_llx}\K\\{x\_traj}[\|i];{}$\2\6
 1345   \&{if} ${}(\\{y\_traj}[\|i]<\\{global\_lly}){}$\1\5
 1346   ${}\\{global\_lly}\K\\{y\_traj}[\|i];{}$\2\6
 1347   \4${}\}{}$\2\6
 1348   \&{if} (\\{verbose})\5
 1349   ${}\{{}$\1\6
 1350   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Global\ box-coun}\)\.{ting\ window:%
 1351   \\n"},\39\\{progname});{}$\6
 1352   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ \ \ (llx,lly)=(\%2}\)\.{.8f,\%2.8f)%
 1353   \\n"},\39\\{progname},\39\\{global\_llx},\39\\{global\_lly});{}$\6
 1354   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ \ \ (urx,ury)=(\%2}\)\.{.8f,\%2.8f)%
 1355   \\n"},\39\\{progname},\39\\{global\_urx},\39\\{global\_ury});{}$\6
 1356   \4${}\}{}$\2\6
 1357   \4${}\}{}$\2\6
 1358   \4${}\}{}$\2\par
 1359   \U7.\fi
 1360   
 1361   \M{18}Get the number of boxes covering the trajectory for all levels of
 1362   refinement.
 1363   This is the main loop where the program iterates over all the refinement levels
 1364   and meanwhile gather how many boxes are needed to cover the trajectory as
 1365   function of the refinement index $N$ for $N_{\rm min}\le N\le N_{\rm max}$.
 1366   In the loop, the program starts with a grid of $[2^{N_{\rm min}}\times
 1367   2^{N_{\rm min}}]$ and ends with a grid of $[2^{N_{\rm max}}\times
 1368   2^{N_{\rm max}}]$ boxes, at each step of refinement increasing the number of
 1369   boxes by a factor of four.
 1370   At each level of increasing resolution, the program relies on the routines
 1371   \PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}(\,)} or \PB{\\{get\_num%
 1372   \_covering\_boxes}(\,)} to
 1373   perform the actual box-counting.
 1374   The results of the box counting are stored in the vector \PB{\\{num\_boxes}},
 1375   which
 1376   is used later on in the actual extraction of the estimate of fractal dimension.
 1377   
 1378   \Y\B\4\X18:Get number of boxes covering the trajectory for all levels of
 1379   refinement in resolution\X${}\E{}$\6
 1380   ${}\{{}$\1\6
 1381   ${}\\{num\_boxes}\K\\{livector}(\T{1},\39\\{nnmax});{}$\6
 1382   ${}\\{nn}\K\T{1};{}$\6
 1383   \&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{nnmin}-\T{1};{}$ ${}\|i\PP){}$\1\5
 1384   ${}\\{nn}\K\T{2}*\\{nn}{}$;\C{ This leaves \PB{\\{nn}} as $2^{(N_{\rm min}-1)}$
 1385   }\2\6
 1386   \&{for} ${}(\|i\K\\{nnmin};{}$ ${}\|i\Z\\{nnmax};{}$ ${}\|i\PP){}$\5
 1387   ${}\{{}$\C{ For all levels in refinement of grid resolution }\1\6
 1388   ${}\\{nn}\K\T{2}*\\{nn};{}$\6
 1389   \&{if} (\\{make\_boxmaps})\5
 1390   ${}\{{}$\C{ do we wish to generate \METAPOST\ graphs of the box distribution? }%
 1391   \1\6
 1392   ${}\\{sprintf}(\\{boxmap\_filename},\39\.{"\%s-\%02ld.mp"},\39\\{output%
 1393   \_filename},\39\|i);{}$\6
 1394   ${}\\{num\_boxes}[\|i]\K\\{get\_num\_covering\_boxes\_with\_boxmaps}(\\{x%
 1395   \_traj},\39\\{y\_traj},\39\\{mm},\39\|i,\39\\{global\_llx},\39\\{global\_lly},%
 1396   \39\\{global\_urx},\39\\{global\_ury},\39\\{boxmap\_filename},\39\\{boxgraph%
 1397   \_width},\39\\{boxgraph\_height},\39\\{trajectory\_filename});{}$\6
 1398   \4${}\}{}$\2\6
 1399   \&{else}\5
 1400   ${}\{{}$\C{ if not, just use the regular number-crunching routine }\1\6
 1401   ${}\\{num\_boxes}[\|i]\K\\{get\_num\_covering\_boxes}(\\{x\_traj},\39\\{y%
 1402   \_traj},\39\\{mm},\39\|i,\39\\{global\_llx},\39\\{global\_lly},\39\\{global%
 1403   \_urx},\39\\{global\_ury});{}$\6
 1404   \4${}\}{}$\2\6
 1405   \&{if} (\\{verbose})\5
 1406   ${}\{{}$\1\6
 1407   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ N=\%ld\ (\%ldx\%ld-}\)\.{grid\ of\
 1408   boxes):\ "},\39\\{progname},\39\|i,\39\\{nn},\39\\{nn});{}$\6
 1409   ${}\\{fprintf}(\\{stdout},\39\.{"Trajectory\ covered\ }\)\.{by\ \%ld\ boxes%
 1410   \\n"},\39\\{num\_boxes}[\|i]);{}$\6
 1411   \4${}\}{}$\2\6
 1412   \4${}\}{}$\2\6
 1413   \4${}\}{}$\2\par
 1414   \U7.\fi
 1415   
 1416   \M{19}Compute the logarithmic estimate of the fractal dimension. Having
 1417   completed
 1418   the task of actually calculating the number of boxes necessary to cover the
 1419   trajectory, for varying box dimensions of width
 1420   $(\PB{\\{global\_urx}}-\PB{\\{global\_llx}})/2^n$ and height $(\PB{\\{global%
 1421   \_ury}}-\PB{\\{global\_lly}})/2^n$,
 1422   for $n=1,2,\ldots,N$, the only remaining arithmetics is to actually calculate
 1423   the estimate for the fractal dimension.
 1424   This is done by performing the fit of the linear function $y=ax+b$ to the data
 1425   set obtained with $2\ln(n)$ as abscissa and $\ln(\PB{\\{num\_boxes}}(n))$ as
 1426   ordinata.
 1427   
 1428   Also in this block, we analyze whether the option \PB{\\{make\_boxmaps}} is set
 1429   as true
 1430   (1) or not (0), determining whether a graph of the number of boxes as function
 1431   of division depth should be written to file as well.
 1432   
 1433   In the fitting of the linear function $y=ax+b$ to the data set $(x_1,y_1)$,
 1434   $(x_2,y_2)$,$\ldots$,$(x_N,y_N)$, the parameters $a$ and $b$ can be explicitly
 1435   obtained from the summation formulas
 1436   $$
 1437   a={{1}\over{D}}\Big(N\sum^N_{k=1}x_ky_k-\sum^N_{k=1}x_k\sum^N_{k=1}y_k\Big),
 1438   \qquad
 1439   b={{1}\over{N}}\Big(\sum^N_{k=1}y_k-a\sum^N_{k=1}x_k\Big),
 1440   $$
 1441   where
 1442   $$
 1443   D\equiv N\sum^N_{k=1}x^2_k-\Big(\sum^N_{k=1}x_k\Big)^2.
 1444   $$
 1445   
 1446   \Y\B\4\X19:Compute the logarithmic estimate of the fractal dimension\X${}\E{}$\6
 1447   ${}\{{}$\1\6
 1448   ${}\|x\K\\{dvector}(\\{nnmin},\39\\{nnmax});{}$\6
 1449   ${}\|y\K\\{dvector}(\\{nnmin},\39\\{nnmax});{}$\6
 1450   ${}\\{fracdimen\_estimates}\K\\{dvector}(\\{nnmin},\39\\{nnmax});{}$\6
 1451   ${}\\{nn}\K\T{1};{}$\6
 1452   \&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{nnmax};{}$ ${}\|i\PP){}$\5
 1453   ${}\{{}$\1\6
 1454   ${}\\{nn}\K\T{2}*\\{nn};{}$\6
 1455   \&{if} ${}(\\{nnmin}\Z\|i){}$\1\5
 1456   ${}\|x[\|i]\K\\{log}{}$((\&{double}) \\{nn});\2\6
 1457   \4${}\}{}$\2\6
 1458   \&{for} ${}(\|i\K\\{nnmin};{}$ ${}\|i\Z\\{nnmax};{}$ ${}\|i\PP){}$\1\5
 1459   ${}\|y[\|i]\K\\{log}(\\{num\_boxes}[\|i]);{}$\2\6
 1460   \&{for} ${}(\|i\K\\{nnmin};{}$ ${}\|i\Z\\{nnmax};{}$ ${}\|i\PP){}$\1\5
 1461   ${}\\{fracdimen\_estimates}[\|i]\K\|y[\|i]/\|x[\|i];{}$\2\6
 1462   \&{if} (\\{verbose})\5
 1463   ${}\{{}$\1\6
 1464   \&{for} ${}(\|i\K\\{nnmin};{}$ ${}\|i\Z\\{nnmax};{}$ ${}\|i\PP){}$\5
 1465   ${}\{{}$\1\6
 1466   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ N=\%ld,\ fractal\ }\)\.{dimension\
 1467   estimate=\%}\)\.{f\\n"},\39\\{progname},\39\|i,\39\\{fracdimen\_estimates}[%
 1468   \|i]);{}$\6
 1469   \4${}\}{}$\2\6
 1470   \4${}\}{}$\2\6
 1471   ${}\\{moment}(\\{fracdimen\_estimates},\39\\{nnmin},\39\\{nnmax},\39{\AND}%
 1472   \\{ave},\39{\AND}\\{adev},\39{\AND}\\{sdev},\39{\AND}\\{var},\39{\AND}\\{skew},%
 1473   \39{\AND}\\{curt});{}$\6
 1474   \&{if} (\\{verbose})\5
 1475   ${}\{{}$\1\6
 1476   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Estimate\ of\ fra}\)\.{ctal\ dimension:\
 1477   \%f\\n}\)\.{"},\39\\{progname},\39\\{ave});{}$\6
 1478   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Average\ deviati}\)\.{on:\ \%f\\n"},\39%
 1479   \\{progname},\39\\{adev});{}$\6
 1480   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Standard\ deviat}\)\.{ion:\ \%f\\n"},\39%
 1481   \\{progname},\39\\{sdev});{}$\6
 1482   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Skewness:\ \%f\\n"},\39\\{progname},\39%
 1483   \\{skew});{}$\6
 1484   \4${}\}{}$\2\6
 1485   ${}\\{free\_livector}(\\{num\_boxes},\39\T{1},\39\\{nnmax}){}$;\C{ release the
 1486   memory occupied by \PB{\\{num\_boxes}} }\6
 1487   ${}\\{free\_dvector}(\\{fracdimen\_estimates},\39\\{nnmin},\39\\{nnmax});{}$\6
 1488   ${}\\{free\_dvector}(\|x,\39\\{nnmin},\39\\{nnmax}){}$;\C{ release the memory
 1489   occupied by \PB{\|x} }\6
 1490   ${}\\{free\_dvector}(\|y,\39\\{nnmin},\39\\{nnmax}){}$;\C{ release the memory
 1491   occupied by \PB{\|y} }\6
 1492   \4${}\}{}$\2\par
 1493   \U7.\fi
 1494   
 1495   \M{20}Save the fractal dimension to file.
 1496   
 1497   \Y\B\4\X20:Save or append the logarithmic estimate to output file\X${}\E{}$\6
 1498   ${}\{{}$\1\6
 1499   \&{if} ${}(\\{output\_mode}\E\.{APPEND}){}$\5
 1500   ${}\{{}$\1\6
 1501   ${}\\{fseek}(\\{frac\_estimate\_file},\39\T{0\$L},\39\.{SEEK\_END});{}$\6
 1502   \4${}\}{}$\2\6
 1503   \&{else} \&{if} ${}(\\{output\_mode}\E\.{OVERWRITE}){}$\5
 1504   ${}\{{}$\1\6
 1505   ${}\\{fseek}(\\{frac\_estimate\_file},\39\T{0\$L},\39\.{SEEK\_SET});{}$\6
 1506   \4${}\}{}$\2\6
 1507   ${}\\{fprintf}(\\{frac\_estimate\_file},\39\.{"\%f\ \%f\ \%f\\n"},\39\\{ave},%
 1508   \39\\{sdev},\39\\{skew});{}$\6
 1509   \4${}\}{}$\2\par
 1510   \U7.\fi
 1511   
 1512   \M{21}Close any open output files.
 1513   
 1514   \Y\B\4\X21:Close file for output of logarithmic estimate\X${}\E{}$\6
 1515   ${}\{{}$\1\6
 1516   \\{fclose}(\\{frac\_estimate\_file});\6
 1517   \4${}\}{}$\2\par
 1518   \U7.\fi
 1519   
 1520   \M{22}Display elapsed execution time.
 1521   
 1522   \Y\B\4\X22:Display elapsed execution time\X${}\E{}$\6
 1523   ${}\{{}$\1\6
 1524   ${}\\{now}\K\\{time}(\NULL);{}$\6
 1525   \&{if} (\\{verbose})\1\5
 1526   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Total\ execution}\)\.{\ time:\ \%d\ s%
 1527   \\n"},\39\\{progname},\39{}$((\&{int}) \\{difftime}${}(\\{now},\39%
 1528   \\{initime})));{}$\2\6
 1529   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Program\ executi}\)\.{on\ closed\ \%s"},%
 1530   \39\\{progname},\39\\{ctime}({\AND}\\{now}));{}$\6
 1531   \4${}\}{}$\2\par
 1532   \U7.\fi
 1533   
 1534   \N{1}{23}Subroutines. In this section, all subroutines as used by the main
 1535   program
 1536   are listed.
 1537   
 1538   \Y\B\4\X23:Subroutines\X${}\E{}$\6
 1539   \X24:Routine for computation of average, average deviation and standard
 1540   deviation\X\6
 1541   \X25:Routine for obtaining the number of coordinate pairs in a file\X\6
 1542   \X26:Routines for removing preceding path of filenames\X\6
 1543   \X29:Routines for memory allocation of vectors and matrices\X\6
 1544   \X38:Routine for displaying help message\X\6
 1545   \X39:Routine for determining whether two lines intersect or not\X\6
 1546   \X40:Routine for determining whether a line and a box intersect or not\X\6
 1547   \X41:Routines for calculation of number of boxes covering the trajectory\X\par
 1548   \U7.\fi
 1549   
 1550   \M{24}Routine for computation of average, average deviation and standard
 1551   deviation.
 1552   This routine is adopted from {\sl Numerical Recipes in C}.
 1553   
 1554   \Y\B\4\X24:Routine for computation of average, average deviation and standard
 1555   deviation\X${}\E{}$\6
 1556   \&{void} \\{moment}(\&{double} \\{data}[\,]${},\39{}$\&{int} \\{nnmin}${},%
 1557   \39{}$\&{int} \\{nnmax}${},\39{}$\&{double} ${}{*}\\{ave},\39{}$\&{double}
 1558   ${}{*}\\{adev},\39{}$\&{double} ${}{*}\\{sdev},\39{}$\&{double} ${}{*}\\{var},%
 1559   \39{}$\&{double} ${}{*}\\{skew},\39{}$\&{double} ${}{*}\\{curt}){}$\1\1\2\2\6
 1560   ${}\{{}$\1\6
 1561   \&{int} \|j;\6
 1562   \&{double} \\{ep}${}\K\T{0.0},{}$ \|s${},{}$ \|p;\7
 1563   \&{if} ${}(\\{nnmax}-\\{nnmin}\Z\T{1}){}$\5
 1564   ${}\{{}$\1\6
 1565   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ routin}\)\.{e\ moment()!\ "},%
 1566   \39\\{progname});{}$\6
 1567   ${}\\{fprintf}(\\{stderr},\39\.{"(nnmax-nnmin>1\ is\ a}\)\.{\ requirement)%
 1568   \\n"});{}$\6
 1569   \\{exit}(\.{FAILURE});\6
 1570   \4${}\}{}$\2\6
 1571   ${}\|s\K\T{0.0};{}$\6
 1572   \&{for} ${}(\|j\K\\{nnmin};{}$ ${}\|j\Z\\{nnmax};{}$ ${}\|j\PP){}$\1\5
 1573   ${}\|s\MRL{+{\K}}\\{data}[\|j];{}$\2\6
 1574   ${}{*}\\{ave}\K\|s/(\\{nnmax}-\\{nnmin}+\T{1});{}$\6
 1575   ${}{*}\\{adev}\K({*}\\{var})\K({*}\\{skew})\K({*}\\{curt})\K\T{0.0};{}$\6
 1576   \&{for} ${}(\|j\K\\{nnmin};{}$ ${}\|j\Z\\{nnmax};{}$ ${}\|j\PP){}$\5
 1577   ${}\{{}$\1\6
 1578   ${}{*}\\{adev}\MRL{+{\K}}\\{fabs}(\|s\K\\{data}[\|j]-({*}\\{ave}));{}$\6
 1579   ${}\\{ep}\MRL{+{\K}}\|s;{}$\6
 1580   ${}\\{ep}\MRL{+{\K}}\|s;{}$\6
 1581   ${}{*}\\{var}\MRL{+{\K}}(\|p\K\|s*\|s);{}$\6
 1582   ${}{*}\\{skew}\MRL{+{\K}}(\|p\MRL{*{\K}}\|s);{}$\6
 1583   ${}{*}\\{curt}\MRL{+{\K}}(\|p\MRL{*{\K}}\|s);{}$\6
 1584   \4${}\}{}$\2\6
 1585   ${}{*}\\{adev}\MRL{{/}{\K}}(\\{nnmax}-\\{nnmin}+\T{1});{}$\6
 1586   ${}{*}\\{var}\K({*}\\{var}-\\{ep}*\\{ep}/(\\{nnmax}-\\{nnmin}+\T{1}))/(%
 1587   \\{nnmax}-\\{nnmin});{}$\6
 1588   ${}{*}\\{sdev}\K\\{sqrt}({*}\\{var});{}$\6
 1589   \&{if} ${}({*}\\{var}){}$\5
 1590   ${}\{{}$\1\6
 1591   ${}{*}\\{skew}\MRL{{/}{\K}}((\\{nnmax}-\\{nnmin}+\T{1})*({*}\\{var})*({*}%
 1592   \\{sdev}));{}$\6
 1593   ${}{*}\\{curt}\K({*}\\{curt})/((\\{nnmax}-\\{nnmin}+\T{1})*({*}\\{var})*({*}%
 1594   \\{var}))-\T{3.0};{}$\6
 1595   \4${}\}{}$\2\6
 1596   \&{else}\5
 1597   ${}\{{}$\1\6
 1598   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ routin}\)\.{e\ moment()!\ "},%
 1599   \39\\{progname});{}$\6
 1600   ${}\\{fprintf}(\\{stderr},\39\.{"No\ skew/kurtosis\ fo}\)\.{r\ zero\ variance%
 1601   \\n"});{}$\6
 1602   \\{exit}(\.{FAILURE});\6
 1603   \4${}\}{}$\2\6
 1604   \4${}\}{}$\2\par
 1605   \U23.\fi
 1606   
 1607   \M{25}Routine for obtaining the number of coordinate pairs in a file. This
 1608   routine is called prior to loading the trajectory, in order to get the
 1609   size needed for allocating the memory for the trajectory vector. As the
 1610   number of coordinates~$M$ has been established, two vectors \PB{\\{x\_traj}[%
 1611   \T{1..}\|M]}
 1612   and \PB{\\{y\_traj}[\T{1..}\|M]}, containing the coordinates $(x_m,y_m)$ of the
 1613   trajectory.
 1614   
 1615   \Y\B\4\X25:Routine for obtaining the number of coordinate pairs in a file\X${}%
 1616   \E{}$\6
 1617   \&{long} \&{int} \\{num\_coordinate\_pairs}(\&{FILE} ${}{*}\\{trajectory%
 1618   \_file}){}$\1\1\2\2\6
 1619   ${}\{{}$\1\6
 1620   \&{double} \\{tmp};\6
 1621   \&{int} \\{tmpch};\6
 1622   \&{long} \&{int} \\{mm}${}\K\T{0};{}$\7
 1623   ${}\\{fseek}(\\{trajectory\_file},\39\T{0\$L},\39\.{SEEK\_SET}){}$;\C{ rewind
 1624   file to beginning }\6
 1625   \&{while} ${}((\\{tmpch}\K\\{getc}(\\{trajectory\_file}))\I\.{EOF}){}$\5
 1626   ${}\{{}$\1\6
 1627   ${}\\{ungetc}(\\{tmpch},\39\\{trajectory\_file});{}$\6
 1628   ${}\\{fscanf}(\\{trajectory\_file},\39\.{"\%lf"},\39{\AND}\\{tmp}){}$;\C{ Read
 1629   away the $x$ coordinate }\6
 1630   ${}\\{fscanf}(\\{trajectory\_file},\39\.{"\%lf"},\39{\AND}\\{tmp}){}$;\C{ Read
 1631   away the $y$ coordinate }\6
 1632   ${}\\{mm}\PP;{}$\6
 1633   ${}\\{tmpch}\K\\{getc}(\\{trajectory\_file}){}$;\C{ Read away blanks and
 1634   linefeeds }\6
 1635   \&{while} ${}((\\{tmpch}\I\.{EOF})\W(\R\\{isdigit}(\\{tmpch}))){}$\1\5
 1636   ${}\\{tmpch}\K\\{getc}(\\{trajectory\_file});{}$\2\6
 1637   \&{if} ${}(\\{tmpch}\I\.{EOF}){}$\1\5
 1638   ${}\\{ungetc}(\\{tmpch},\39\\{trajectory\_file});{}$\2\6
 1639   \4${}\}{}$\2\6
 1640   ${}\\{fseek}(\\{trajectory\_file},\39\T{0\$L},\39\.{SEEK\_SET}){}$;\C{ rewind
 1641   file to beginning }\6
 1642   \&{return} (\\{mm});\6
 1643   \4${}\}{}$\2\par
 1644   \U23.\fi
 1645   
 1646   \M{26}Routines for removing preceding path of filenames.
 1647   In this block all routines related to removing preceding path strings go.
 1648   Not really fancy programming, and no contribution to any increase of numerical
 1649   efficiency or precision; just for the sake of keeping a tidy terminal output
 1650   of the program. The \PB{\\{strip\_away\_path}(\,)} routine is typically called
 1651   when
 1652   initializing the program name string \PB{\\{progname}} from the command line
 1653   string
 1654   \PB{\\{argv}[\T{0}]}, and is typically located in the blocks related to parsing
 1655   of the
 1656   command line options.
 1657   
 1658   \Y\B\4\X26:Routines for removing preceding path of filenames\X${}\E{}$\6
 1659   \X27:Routine for checking valid path characters\X\6
 1660   \X28:Routine for stripping away path string\X\par
 1661   \U23.\fi
 1662   
 1663   \M{27}Checking for a valid path character.
 1664   The \PB{\\{pathcharacter}} routine takes one character \PB{\\{ch}} as argument,
 1665   and returns
 1666   1 (``true'') if the character is valid character of a path string, otherwise 0
 1667   (``false'') is returned.
 1668   
 1669   \Y\B\4\X27:Routine for checking valid path characters\X${}\E{}$\6
 1670   \&{short} \\{pathcharacter}(\&{int} \\{ch})\1\1\2\2\6
 1671   ${}\{{}$\1\6
 1672   \&{return} ${}(\\{isalnum}(\\{ch})\V(\\{ch}\E\.{'.'})\V(\\{ch}\E\.{'/'})\V(%
 1673   \\{ch}\E\.{'\\\\'})\V(\\{ch}\E\.{'\_'})\V(\\{ch}\E\.{'-'})\V(\\{ch}\E%
 1674   \.{'+'}));{}$\6
 1675   \4${}\}{}$\2\par
 1676   \U26.\fi
 1677   
 1678   \M{28}Routine for stripping away path string of a file name.
 1679   The \PB{\\{strip\_away\_path}(\,)} routine takes a character string \PB{%
 1680   \\{filename}} as
 1681   argument, and returns a pointer to the same string but without any
 1682   preceding path segments. This routine is, for example, useful for
 1683   removing paths from program names as parsed from the command line.
 1684   
 1685   \Y\B\4\X28:Routine for stripping away path string\X${}\E{}$\6
 1686   \&{char} ${}{*}{}$\\{strip\_away\_path}(\&{char} \\{filename}[\,])\1\1\2\2\6
 1687   ${}\{{}$\1\6
 1688   \&{int} \|j${},{}$ \|k${}\K\T{0};{}$\7
 1689   \&{while} (\\{pathcharacter}(\\{filename}[\|k]))\1\5
 1690   ${}\|k\PP;{}$\2\6
 1691   ${}\|j\K(\MM\|k){}$;\C{ this is the uppermost index of the full path+file
 1692   string }\6
 1693   \&{while} (\\{isalnum}((\&{int})(\\{filename}[\|j])))\1\5
 1694   ${}\|j\MM;{}$\2\6
 1695   ${}\|j\PP{}$;\C{ this is the lowermost index of the stripped file name }\6
 1696   \&{return} ${}({\AND}\\{filename}[\|j]);{}$\6
 1697   \4${}\}{}$\2\par
 1698   \U26.\fi
 1699   
 1700   \M{29}Subroutines for memory allocation. Here follows the routines for memory
 1701   allocation and deallocation of double precision real and complex valued
 1702   vectors, as used for storing the optical field distribution in the grating,
 1703   the refractive index distribution of the grating, etc.
 1704   
 1705   \Y\B\4\X29:Routines for memory allocation of vectors and matrices\X${}\E{}$\6
 1706   \X30:Routine for allocation of vectors of double precision\X\6
 1707   \X33:Routine for deallocation of vectors of integer precision\X\6
 1708   \X32:Routine for allocation of vectors of integer precision\X\6
 1709   \X31:Routine for deallocation of vectors of double precision\X\6
 1710   \X36:Routine for allocation of matrices of short integer precision\X\6
 1711   \X37:Routine for deallocation of matrices of short integer precision\X\par
 1712   \U23.\fi
 1713   
 1714   \M{30}The \PB{\\{dvector}} routine allocates a real-valued vector of double
 1715   precision,
 1716   with vector index ranging from \PB{\\{nl}} to \PB{\\{nh}}.
 1717   
 1718   \Y\B\4\X30:Routine for allocation of vectors of double precision\X${}\E{}$\6
 1719   \&{double} ${}{*}{}$\\{dvector}(\&{long} \\{nl}${},\39{}$\&{long} \\{nh})\1\1\2%
 1720   \2\6
 1721   ${}\{{}$\1\6
 1722   \&{double} ${}{*}\|v;{}$\7
 1723   ${}\|v\K{}$(\&{double} ${}{*}){}$ \\{malloc}${}((\&{size\_t})((\\{nh}-\\{nl}+%
 1724   \T{2})*\&{sizeof}(\&{double})));{}$\6
 1725   \&{if} ${}(\R\|v){}$\5
 1726   ${}\{{}$\1\6
 1727   ${}\\{fprintf}(\\{stderr},\39\.{"Error:\ Allocation\ f}\)\.{ailure\ in\
 1728   dvector()\\}\)\.{n"});{}$\6
 1729   \\{exit}(\.{FAILURE});\6
 1730   \4${}\}{}$\2\6
 1731   \&{return} \|v${}-\\{nl}+\T{1};{}$\6
 1732   \4${}\}{}$\2\par
 1733   \U29.\fi
 1734   
 1735   \M{31}The \PB{\\{free\_dvector}} routine release the memory occupied by the
 1736   real-valued vector \PB{$\|v[\\{nl}\MRL{{.}{.}}\\{nh}]$}.
 1737   
 1738   \Y\B\4\X31:Routine for deallocation of vectors of double precision\X${}\E{}$\6
 1739   \&{void} \\{free\_dvector}(\&{double} ${}{*}\|v,\39{}$\&{long} \\{nl}${},\39{}$%
 1740   \&{long} \\{nh})\1\1\2\2\6
 1741   ${}\{{}$\1\6
 1742   \\{free}((\&{char} ${}{*})(\|v+\\{nl}-\T{1}));{}$\6
 1743   \4${}\}{}$\2\par
 1744   \U29.\fi
 1745   
 1746   \M{32}The \PB{\\{ivector}} routine allocates a real-valued vector of integer
 1747   precision,
 1748   with vector index ranging from \PB{\\{nl}} to \PB{\\{nh}}.
 1749   
 1750   \Y\B\4\X32:Routine for allocation of vectors of integer precision\X${}\E{}$\6
 1751   \&{int} ${}{*}{}$\\{ivector}(\&{long} \\{nl}${},\39{}$\&{long} \\{nh})\1\1\2\2\6
 1752   ${}\{{}$\1\6
 1753   \&{int} ${}{*}\|v;{}$\7
 1754   ${}\|v\K{}$(\&{int} ${}{*}){}$ \\{malloc}${}((\&{size\_t})((\\{nh}-\\{nl}+%
 1755   \T{2})*\&{sizeof}(\&{int})));{}$\6
 1756   \&{if} ${}(\R\|v){}$\5
 1757   ${}\{{}$\1\6
 1758   ${}\\{fprintf}(\\{stderr},\39\.{"Error:\ Allocation\ f}\)\.{ailure\ in\
 1759   ivector()\\}\)\.{n"});{}$\6
 1760   \\{exit}(\.{FAILURE});\6
 1761   \4${}\}{}$\2\6
 1762   \&{return} \|v${}-\\{nl}+\T{1};{}$\6
 1763   \4${}\}{}$\2\par
 1764   \A34.
 1765   \U29.\fi
 1766   
 1767   \M{33}The \PB{\\{free\_ivector}} routine release the memory occupied by the
 1768   real-valued vector \PB{$\|v[\\{nl}\MRL{{.}{.}}\\{nh}]$}.
 1769   
 1770   \Y\B\4\X33:Routine for deallocation of vectors of integer precision\X${}\E{}$\6
 1771   \&{void} \\{free\_ivector}(\&{int} ${}{*}\|v,\39{}$\&{long} \\{nl}${},\39{}$%
 1772   \&{long} \\{nh})\1\1\2\2\6
 1773   ${}\{{}$\1\6
 1774   \\{free}((\&{char} ${}{*})(\|v+\\{nl}-\T{1}));{}$\6
 1775   \4${}\}{}$\2\par
 1776   \A35.
 1777   \U29.\fi
 1778   
 1779   \M{34}The \PB{\\{livector}} routine allocates a real-valued vector of long
 1780   integer
 1781   precision, with vector index ranging from \PB{\\{nl}} to \PB{\\{nh}}.
 1782   
 1783   \Y\B\4\X32:Routine for allocation of vectors of integer precision\X${}\mathrel+%
 1784   \E{}$\6
 1785   \&{long} \&{int} ${}{*}{}$\\{livector}(\&{long} \\{nl}${},\39{}$\&{long} %
 1786   \\{nh})\1\1\2\2\6
 1787   ${}\{{}$\1\6
 1788   \&{long} \&{int} ${}{*}\|v;{}$\7
 1789   ${}\|v\K{}$(\&{long} \&{int} ${}{*}){}$ \\{malloc}${}((\&{size\_t})((\\{nh}-%
 1790   \\{nl}+\T{2})*{}$\&{sizeof}(\&{long} \&{int})));\6
 1791   \&{if} ${}(\R\|v){}$\5
 1792   ${}\{{}$\1\6
 1793   ${}\\{fprintf}(\\{stderr},\39\.{"Error:\ Allocation\ f}\)\.{ailure\ in\
 1794   livector()}\)\.{\\n"});{}$\6
 1795   \\{exit}(\.{FAILURE});\6
 1796   \4${}\}{}$\2\6
 1797   \&{return} \|v${}-\\{nl}+\T{1};{}$\6
 1798   \4${}\}{}$\2\par
 1799   \fi
 1800   
 1801   \M{35}The \PB{\\{free\_livector}} routine release the memory occupied by the
 1802   real-valued vector \PB{$\|v[\\{nl}\MRL{{.}{.}}\\{nh}]$}.
 1803   
 1804   \Y\B\4\X33:Routine for deallocation of vectors of integer precision\X${}%
 1805   \mathrel+\E{}$\6
 1806   \&{void} \\{free\_livector}(\&{long} \&{int} ${}{*}\|v,\39{}$\&{long} %
 1807   \\{nl}${},\39{}$\&{long} \\{nh})\1\1\2\2\6
 1808   ${}\{{}$\1\6
 1809   \\{free}((\&{char} ${}{*})(\|v+\\{nl}-\T{1}));{}$\6
 1810   \4${}\}{}$\2\par
 1811   \fi
 1812   
 1813   \M{36}The \PB{\\{simatrix}} routine allocates an array of short integer
 1814   precision,
 1815   with array row index ranging from \PB{\\{nrl}} to \PB{\\{nrh}} and column index
 1816   ranging
 1817   from \PB{\\{ncl}} to \PB{\\{nch}}.
 1818   
 1819   \Y\B\4\X36:Routine for allocation of matrices of short integer precision\X${}%
 1820   \E{}$\6
 1821   \&{short} \&{int} ${}{*}{*}{}$\\{simatrix}(\&{long} \\{nrl}${},\39{}$\&{long} %
 1822   \\{nrh}${},\39{}$\&{long} \\{ncl}${},\39{}$\&{long} \\{nch})\1\1\2\2\6
 1823   ${}\{{}$\1\6
 1824   \&{long} \|i${},{}$ \\{nrow}${}\K\\{nrh}-\\{nrl}+\T{1},{}$ \\{ncol}${}\K%
 1825   \\{nch}-\\{ncl}+\T{1};{}$\6
 1826   \&{short} \&{int} ${}{*}{*}\|m;{}$\7
 1827   ${}\|m\K{}$(\&{short} \&{int} ${}{*}{*}){}$ \\{malloc}${}((\&{size\_t})((%
 1828   \\{nrow}+\T{1})*{}$\&{sizeof}(\&{short} \&{int} ${}{*})));{}$\6
 1829   \&{if} ${}(\R\|m){}$\5
 1830   ${}\{{}$\1\6
 1831   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Allocation\ fail}\)\.{ure\ 1\ in\
 1832   simatrix()\\}\)\.{n"},\39\\{progname});{}$\6
 1833   \\{exit}(\.{FAILURE});\6
 1834   \4${}\}{}$\2\6
 1835   ${}\|m\MRL{+{\K}}\T{1};{}$\6
 1836   ${}\|m\MRL{-{\K}}\\{nrl};{}$\6
 1837   ${}\|m[\\{nrl}]\K{}$(\&{short} \&{int} ${}{*}){}$ \\{malloc}${}((\&{size\_t})((%
 1838   \\{nrow}*\\{ncol}+\T{1})*{}$\&{sizeof}(\&{short} \&{int})));\6
 1839   \&{if} ${}(\R\|m[\\{nrl}]){}$\5
 1840   ${}\{{}$\1\6
 1841   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Allocation\ fail}\)\.{ure\ 2\ in\
 1842   simatrix()\\}\)\.{n"},\39\\{progname});{}$\6
 1843   \\{exit}(\.{FAILURE});\6
 1844   \4${}\}{}$\2\6
 1845   ${}\|m[\\{nrl}]\MRL{+{\K}}\T{1};{}$\6
 1846   ${}\|m[\\{nrl}]\MRL{-{\K}}\\{ncl};{}$\6
 1847   \&{for} ${}(\|i\K\\{nrl}+\T{1};{}$ ${}\|i\Z\\{nrh};{}$ ${}\|i\PP){}$\1\5
 1848   ${}\|m[\|i]\K\|m[\|i-\T{1}]+\\{ncol};{}$\2\6
 1849   \&{return} \|m;\6
 1850   \4${}\}{}$\2\par
 1851   \U29.\fi
 1852   
 1853   \M{37}The \PB{\\{free\_simatrix}} routine releases the memory occupied by the
 1854   short integer matrix \PB{$\|v[\\{nl}\MRL{{.}{.}}\\{nh}]$}, as allocated by \PB{%
 1855   \\{simatrix}(\,)}.
 1856   
 1857   \Y\B\4\X37:Routine for deallocation of matrices of short integer precision\X${}%
 1858   \E{}$\6
 1859   \&{void} \\{free\_simatrix}(\&{short} \&{int} ${}{*}{*}\|m,\39{}$\&{long} %
 1860   \\{nrl}${},\39{}$\&{long} \\{nrh}${},\39{}$\&{long} \\{ncl}${},\39{}$\&{long} %
 1861   \\{nch})\1\1\2\2\6
 1862   ${}\{{}$\1\6
 1863   \\{free}((\&{char} ${}{*})(\|m[\\{nrl}]+\\{ncl}-\T{1}));{}$\6
 1864   \\{free}((\&{char} ${}{*})(\|m+\\{nrl}-\T{1}));{}$\6
 1865   \4${}\}{}$\2\par
 1866   \U29.\fi
 1867   
 1868   \M{38}Routine for displaying help message to standard terminal output.
 1869   
 1870   \Y\B\4\X38:Routine for displaying help message\X${}\E{}$\6
 1871   \&{void} \\{showsomehelp}(\&{void})\1\1\2\2\6
 1872   ${}\{{}$\1\6
 1873   ${}\\{fprintf}(\\{stderr},\39\.{"Usage:\ \%s\ [options]}\)\.{\\n"},\39%
 1874   \\{progname});{}$\6
 1875   ${}\\{fprintf}(\\{stderr},\39\.{"Options:\\n"});{}$\6
 1876   ${}\\{fprintf}(\\{stderr},\39\.{"\ -h,\ --help\\n"});{}$\6
 1877   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ Display\ this\ he}\)\.{lp\ message\
 1878   and\ exit\ }\)\.{clean.\\n"});{}$\6
 1879   ${}\\{fprintf}(\\{stderr},\39\.{"\ -v,\ --verbose\\n"});{}$\6
 1880   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ Toggle\ verbose\ }\)\.{mode\ on/off.%
 1881   \\n"});{}$\6
 1882   ${}\\{fprintf}(\\{stderr},\39\.{"\ -o,\ --outputfile\ <}\)\.{str>\\n"});{}$\6
 1883   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ Specifies\ the\ b}\)\.{ase\ name\ of\
 1884   the\ outp}\)\.{ut\ files\ where\ the\ p}\)\.{rogram\\n"});{}$\6
 1885   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ is\ to\ save\ the\ }\)\.{calculated\
 1886   data.\ If\ }\)\.{the\ --outputfile\ or\ }\)\.{-o\\n"});{}$\6
 1887   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ option\ is\ follo}\)\.{wed\ by\
 1888   'append'\ the\ }\)\.{estimate\ for\ the\ fra}\)\.{ctal\\n"});{}$\6
 1889   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ dimension\ will\ }\)\.{be\ appended\
 1890   to\ the\ f}\)\.{ile\ named\ <str>.dat,}\)\.{\ which\\n"});{}$\6
 1891   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ will\ be\ created}\)\.{\ if\ it\ does\
 1892   not\ exis}\)\.{t.\ If\ the\ following\ }\)\.{word\\n"});{}$\6
 1893   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ instead\ is\ `ove}\)\.{rwrite'\ the\
 1894   file\ wil}\)\.{l\ instead\ be\ overwri}\)\.{tten.\\n"});{}$\6
 1895   ${}\\{fprintf}(\\{stderr},\39\.{"\ -i,\ --trajectoryfi}\)\.{le\\n"});{}$\6
 1896   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ Specifies\ the\ i}\)\.{nput\
 1897   trajectory\ of\ t}\)\.{he\ graph\ whose\ fract}\)\.{al\\n"});{}$\6
 1898   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ dimension\ is\ to}\)\.{\ be\
 1899   estimated.\\n"});{}$\6
 1900   ${}\\{fprintf}(\\{stderr},\39\.{"\ -w,\ --computationw}\)\.{indow\ llx\ <num>\
 1901   lly\ }\)\.{<num>\ urx\ <num>\ ury\ }\)\.{<num>\\n"});{}$\6
 1902   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ This\ option\ exp}\)\.{licitly\
 1903   specifies\ th}\)\.{e\ domain\ over\ which\ }\)\.{the\\n"});{}$\6
 1904   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ box-counting\ al}\)\.{gorithm\ will\
 1905   be\ appl}\)\.{ied,\ in\ terms\ of\ the}\)\.{\\n"});{}$\6
 1906   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ lower-left\ and\ }\)\.{upper-right\
 1907   corners\ }\)\.{(llx,lly)\ and\ (urx,u}\)\.{ry),\\n"});{}$\6
 1908   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ respectively.\ B}\)\.{y\ specifying\
 1909   this\ op}\)\.{tion,\ any\ automatic\\}\)\.{n"});{}$\6
 1910   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ calculation\ of\ }\)\.{computational\
 1911   window}\)\.{\ will\ be\ neglected.\\}\)\.{n"});{}$\6
 1912   ${}\\{fprintf}(\\{stderr},\39\.{"\ -m,\ --boxmaps\\n"});{}$\6
 1913   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ If\ this\ option\ }\)\.{is\ present,\
 1914   the\ prog}\)\.{ram\ will\ generate\\n"});{}$\6
 1915   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ MetaPost\ code\ f}\)\.{or\ maps\ of\
 1916   the\ distr}\)\.{ibution\ of\ boxes.\\n"});{}$\6
 1917   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ In\ doing\ so,\ al}\)\.{so\ the\ input%
 1918   \ traject}\)\.{ory\ is\ included\ in\\n}\)\.{"});{}$\6
 1919   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ the\ graphs.\ The}\)\.{\ convention\
 1920   for\ the\ }\)\.{naming\ of\ the\ output}\)\.{\\n"});{}$\6
 1921   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ map\ files\ is\ th}\)\.{at\ they\ are\
 1922   saved\ as}\)\.{\ <str>.<N>.dat,\\n"});{}$\6
 1923   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ where\ <str>\ is\ }\)\.{the\ base\
 1924   filename\ as}\)\.{\ specified\ using\ the}\)\.{\\n"});{}$\6
 1925   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ -o\ or\ --outputf}\)\.{ile\ option,\
 1926   <N>\ is\ t}\)\.{he\ automatically\ app}\)\.{ended\\n"});{}$\6
 1927   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ current\ level\ o}\)\.{f\ resolution\
 1928   refinem}\)\.{ent\ in\ the\ box-count}\)\.{ing,\\n"});{}$\6
 1929   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ and\ where\ '.dat}\)\.{'\ is\ the\
 1930   file\ suffix}\)\.{\ as\ automatically\ ap}\)\.{pended\\n"});{}$\6
 1931   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ by\ the\ program.}\)\.{\\n"});{}$\6
 1932   ${}\\{fprintf}(\\{stderr},\39\.{"\ --graphsize\ <width}\)\.{\ in\ mm>\ <height\
 1933   in\ m}\)\.{m>\\n"});{}$\6
 1934   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ If\ the\ -m\ or\ --}\)\.{boxmaps\
 1935   option\ is\ pr}\)\.{esent\ at\ the\ command}\)\.{\ line,\\n"});{}$\6
 1936   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ then\ the\ --grap}\)\.{hsize\ option\
 1937   will\ ov}\)\.{erride\ the\ default\ g}\)\.{raph\\n"});{}$\6
 1938   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ size\ of\ the\ gen}\)\.{erated\
 1939   boxmaps.\ (Def}\)\.{ault\ graph\ size\ is\ 8}\)\.{0\ mm\\n"});{}$\6
 1940   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ width\ and\ 56\ mm}\)\.{\ height.)%
 1941   \\n"});{}$\6
 1942   ${}\\{fprintf}(\\{stderr},\39\.{"\ -Nmin,\ --minlevel\ }\)\.{<num>\\n"});{}$\6
 1943   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ Specifies\ the\ m}\)\.{inimum\ level\
 1944   Nmin\ of}\)\.{\ grid\ refinement\ tha}\)\.{t\ \\n"});{}$\6
 1945   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ will\ be\ used\ in}\)\.{\ the\
 1946   evaluation\ of\ t}\)\.{he\ estimate\ of\ the\ f}\)\.{ractal\\n"});{}$\6
 1947   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ dimension.\ With}\)\.{\ this\ option\
 1948   specifi}\)\.{ed,\ the\ coarsest\ lev}\)\.{el\\n"});{}$\6
 1949   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ used\ in\ the\ box}\)\.{-counting\
 1950   will\ be\ a\ }\)\.{[(2\^Nmin)x(2\^Nmin)]-}\)\.{grid\\n"});{}$\6
 1951   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ of\ boxes.\\n"});{}$\6
 1952   ${}\\{fprintf}(\\{stderr},\39\.{"\ -Nmax,\ --maxlevel\ }\)\.{<num>\\n"});{}$\6
 1953   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ This\ option\ is\ }\)\.{similar\ to\
 1954   the\ --min}\)\.{level\ or\ -Nmin\ optio}\)\.{ns,\\n"});{}$\6
 1955   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ with\ the\ differ}\)\.{ence\ that\ it\
 1956   instead}\)\.{\ specifies\ the\ maxim}\)\.{um\\n"});{}$\6
 1957   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ level\ Nmax\ of\ g}\)\.{rid\
 1958   refinement\ that\ }\)\.{will\ be\ used\ in\ the\\}\)\.{n"});{}$\6
 1959   ${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ evaluation\ of\ t}\)\.{he\ estimate\
 1960   of\ the\ f}\)\.{ractal\ dimension.\\n"});{}$\6
 1961   \4${}\}{}$\2\par
 1962   \U23.\fi
 1963   
 1964   \M{39}The \PB{$\\{lines\_intersect}(\\{p1x},\\{p1y},\\{q1x},\\{q1y},\\{p2x},%
 1965   \\{p2y},\\{q2x},\\{q2y})$} routine takes the
 1966   start and end points of two line segments, the first between $(\PB{\\{p1x}},%
 1967   \PB{\\{p1y}})$
 1968   and $(\PB{\\{q1x}},\PB{\\{q1y}})$ and the second between $(\PB{\\{p2x}},\PB{%
 1969   \\{p2y}})$ and $(\PB{\\{q2x}},\PB{\\{q2y}})$,
 1970   and returns $1$ (`true') if they are found to intersect, and $0$ (`false')
 1971   otherwise.
 1972   
 1973   For a brief sumnmary of the algorithm behind this routine, consider two line
 1974   segments in the plane, the first one between the points ${\bf p}_1$ and
 1975   ${\bf q}_1$ and the second one between ${\bf p}_2$ and ${\bf q}_1$.
 1976   In general, these segments can be written in parametric forms as
 1977   ${\bf r}_1={\bf p}_1+t_1({\bf q}_1-{\bf p}_1)$ and
 1978   ${\bf r}_2={\bf p}_2+t_2({\bf q}_2-{\bf p}_2)$
 1979   for $0\le t_1\le1$ and  $0\le t_2\le1$.
 1980   These line segments intersect each other if they for these intervals for the
 1981   parameters $t_1$ and $t_2$ share a common point, or equivalently if the
 1982   solutions to
 1983   $$
 1984   {\bf p}_1+t_1({\bf q}_1-{\bf p}_1)={\bf p}_2+t_2({\bf q}_2-{\bf p}_2)
 1985   \qquad\Leftrightarrow\qquad
 1986   ({\bf q}_1-{\bf p}_1)t_1+({\bf p}_2-{\bf q}_2)t_2={\bf p}_2-{\bf p}_1
 1987   $$
 1988   both simultaneously satisfy $0\le t_1\le1$ and $0\le t_2\le1$.
 1989   This vectorial equation can equivalently be written in component form as
 1990   $$
 1991   \eqalign{
 1992   (q_{1x}-p_{1x})t_1+(p_{2x}-q_{2x})t_2&=p_{2x}-p_{1x},\cr
 1993   (q_{1y}-p_{1y})t_1+(p_{2y}-q_{2y})t_2&=p_{2y}-p_{1y},\cr
 1994   }
 1995   $$
 1996   which after some straightforward algebra gives the explicit solutions for the
 1997   parameters $t_1$ and $t_2$ as
 1998   $$
 1999   t_1={{ed-bf}\over{ad-bc}},\qquad t_2={{af-ec}\over{ad-bc}},
 2000   $$
 2001   where
 2002   $$
 2003   \eqalign{
 2004   a\equiv(q_{1x}-p_{1x}),\qquad
 2005   b\equiv(p_{2x}-q_{2x}),\qquad
 2006   c\equiv(q_{1y}-p_{1y}),\cr
 2007   d\equiv(p_{2y}-q_{2y}),\qquad
 2008   e\equiv(p_{2x}-p_{1x}),\qquad
 2009   f\equiv(p_{2x}-p_{1x}).\cr
 2010   }
 2011   $$
 2012   Hence, the two line segments intersect if and only if
 2013   $$
 2014   0\le{{ed-bf}\over{ad-bc}}\le1,\qquad
 2015   {\rm and}\qquad 0\le{{af-ec}\over{ad-bc}}\le1.
 2016   $$
 2017   By observing that their denominators are equal, the evaluation of these quotes
 2018   involves in total 6 floating-point multiplications, 2 divisions and 3
 2019   subtractions.
 2020   Notice that whenever $ad-bd=0$, the two lines are parallell and will never
 2021   intersect, regardless of the values of $t_1$ and $t_2$.
 2022   
 2023   \Y\B\4\X39:Routine for determining whether two lines intersect or not\X${}\E{}$%
 2024   \6
 2025   \&{short} \&{int} \\{lines\_intersect}(\&{double} \\{p1x}${},\39{}$\&{double} %
 2026   \\{p1y}${},\39{}$\&{double} \\{q1x}${},\39{}$\&{double} \\{q1y}${},\39{}$%
 2027   \&{double} \\{p2x}${},\39{}$\&{double} \\{p2y}${},\39{}$\&{double} \\{q2x}${},%
 2028   \39{}$\&{double} \\{q2y})\1\1\2\2\6
 2029   ${}\{{}$\1\6
 2030   \&{double} \|a${},{}$ \|b${},{}$ \|c${},{}$ \|d${},{}$ \|e${},{}$ \|f${},{}$ %
 2031   \\{det}${},{}$ \\{tmp1}${},{}$ \\{tmp2};\6
 2032   \&{short} \&{int} \\{intersect};\7
 2033   ${}\|a\K\\{q1x}-\\{p1x};{}$\6
 2034   ${}\|b\K\\{p2x}-\\{q2x};{}$\6
 2035   ${}\|c\K\\{q1y}-\\{p1y};{}$\6
 2036   ${}\|d\K\\{p2y}-\\{q2y};{}$\6
 2037   ${}\|e\K\\{p2x}-\\{p1x};{}$\6
 2038   ${}\|f\K\\{p2y}-\\{p1y};{}$\6
 2039   ${}\\{det}\K\|a*\|d-\|b*\|c;{}$\6
 2040   ${}\\{tmp1}\K\|e*\|d-\|b*\|f;{}$\6
 2041   ${}\\{tmp2}\K\|a*\|f-\|e*\|c;{}$\6
 2042   ${}\\{intersect}\K\T{0};{}$\6
 2043   \&{if} ${}(\\{det}>\T{0}){}$\5
 2044   ${}\{{}$\1\6
 2045   \&{if} ${}(((\T{0.0}\Z\\{tmp1})\W(\\{tmp1}\Z\\{det}))\W((\T{0.0}\Z\\{tmp2})\W(%
 2046   \\{tmp2}\Z\\{det}))){}$\1\5
 2047   ${}\\{intersect}\K\T{1};{}$\2\6
 2048   \4${}\}{}$\2\6
 2049   \&{else} \&{if} ${}(\\{det}<\T{0}){}$\5
 2050   ${}\{{}$\1\6
 2051   \&{if} ${}(((\\{det}\Z\\{tmp1})\W(\\{tmp1}\Z\T{0.0}))\W((\\{det}\Z\\{tmp2})\W(%
 2052   \\{tmp2}\Z\T{0.0}))){}$\1\5
 2053   ${}\\{intersect}\K\T{1};{}$\2\6
 2054   \4${}\}{}$\2\6
 2055   \&{return} (\\{intersect});\6
 2056   \4${}\}{}$\2\par
 2057   \U23.\fi
 2058   
 2059   \M{40}Routine for determining whether a line and a box intersect or not.
 2060   The \PB{\\{box\_intersection}(\,)} routine simply divides the input box, being
 2061   characterized by its lower-left and upper-right corners $(\PB{\\{llx}},\PB{%
 2062   \\{lly}})$
 2063   and $(\PB{\\{urx}},\PB{\\{ury}})$, into the four line segments corresponding to
 2064   its four
 2065   edges, followed by calling the routine \PB{\\{lines\_intersect}(\,)} to
 2066   determine if
 2067   any of these edges intersect the line segment. If an intersection is found,
 2068   the \PB{\\{box\_intersection}(\,)} routine returns 1 (true), otherwise 0
 2069   (false).
 2070   \medskip\noindent
 2071   Input variables:
 2072   \varitem[\PB{\\{px}}, \PB{\\{py}}]{The coordinates of the first end point
 2073   ${\bf p}=(p_x,p_y)$ of the line segment.}
 2074   \varitem[\PB{\\{qx}}, \PB{\\{qy}}]{The coordinates of the second end point
 2075   ${\bf q}=(q_x,q_y)$ of the line segment.}
 2076   \varitem[\PB{\\{llx}}, \PB{\\{lly}}]{Coordinates of the lower-left corner of
 2077   the box.}
 2078   \varitem[\PB{\\{urx}}, \PB{\\{ury}}]{Coordinates of the upper-right corner of
 2079   the box.}
 2080   \medskip\noindent
 2081   Output variables:
 2082   \varitem[]{On exit, the routine returns 1 if an intersection is found,
 2083   otherwise 0 is returned, in either case the value are returned as
 2084   integers of \PB{\&{short}} precision.}
 2085   \medskip
 2086   
 2087   \Y\B\4\X40:Routine for determining whether a line and a box intersect or not%
 2088   \X${}\E{}$\6
 2089   \&{short} \&{int} \\{box\_intersection}(\&{double} \\{px}${},\39{}$\&{double} %
 2090   \\{py}${},\39{}$\&{double} \\{qx}${},\39{}$\&{double} \\{qy}${},\39{}$%
 2091   \&{double} \\{llx}${},\39{}$\&{double} \\{lly}${},\39{}$\&{double} \\{urx}${},%
 2092   \39{}$\&{double} \\{ury})\1\1\2\2\6
 2093   ${}\{{}$\1\6
 2094   \&{if} ${}(\\{lines\_intersect}(\\{px},\39\\{py},\39\\{qx},\39\\{qy},\39%
 2095   \\{llx},\39\\{lly},\39\\{urx},\39\\{lly})){}$\1\5
 2096   \&{return} (\T{1});\C{ intersection with bottom edge }\2\6
 2097   \&{if} ${}(\\{lines\_intersect}(\\{px},\39\\{py},\39\\{qx},\39\\{qy},\39%
 2098   \\{urx},\39\\{lly},\39\\{urx},\39\\{ury})){}$\1\5
 2099   \&{return} (\T{1});\C{ intersection with right edge }\2\6
 2100   \&{if} ${}(\\{lines\_intersect}(\\{px},\39\\{py},\39\\{qx},\39\\{qy},\39%
 2101   \\{urx},\39\\{ury},\39\\{llx},\39\\{ury})){}$\1\5
 2102   \&{return} (\T{1});\C{ intersection with top edge }\2\6
 2103   \&{if} ${}(\\{lines\_intersect}(\\{px},\39\\{py},\39\\{qx},\39\\{qy},\39%
 2104   \\{llx},\39\\{ury},\39\\{llx},\39\\{lly})){}$\1\5
 2105   \&{return} (\T{1});\C{ intersection with left edge }\2\6
 2106   \&{return} (\T{0});\C{ this happens only if no edge is intersecting the line
 2107   segment }\6
 2108   \4${}\}{}$\2\par
 2109   \U23.\fi
 2110   
 2111   \M{41}Routines for calculation of number of boxes covering the trajectory.
 2112   There
 2113   are two almost identical routines for the calculation of the number of boxes
 2114   covering the input trajectory at a given level of subdivision of the box sizes.
 2115   The first routine, \PB{\\{get\_num\_covering\_boxes}(\,)} simply performs this
 2116   task
 2117   without caring of documenting the box distributions as graphs, or ``box maps'',
 2118   while the second one, \PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}(\,)}
 2119   also includes
 2120   the generation of these maps, with output in terms of \METAPOST\ code.
 2121   
 2122   \Y\B\4\X41:Routines for calculation of number of boxes covering the trajectory%
 2123   \X${}\E{}$\6
 2124   \X42:Routine for calculation of number of boxes covering the trajectory\X\6
 2125   \X52:Simplified routine for calculation of number of boxes covering the
 2126   trajectory\X\par
 2127   \U23.\fi
 2128   
 2129   \M{42}Routine for calculation of number of boxes covering the trajectory, also
 2130   including the generation of documenting graphs of the box distributions.
 2131   The \PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}(\,)} routine takes a
 2132   trajectory as
 2133   described by a discrete set of coordinates as input, and for a given grid
 2134   refinement $N$ returns the number of boxes needed to entirely cover the
 2135   trajectory. The grid refinement factor $N$ indicates that the domain of
 2136   computation is divided into a $[2^N\times2^N]$-grid of boxes.
 2137   
 2138   The computational domain in which the box counting is to be performed is
 2139   explicitly stated by the coordinates of its lower-left and upper-right
 2140   corners, $(global_llx,global_lly)$ and $(global_urx,global_ury)$, respectively.
 2141   Parts of the trajectory which are outside of this domain are not included in
 2142   the box-counting.
 2143   If the entire input trajectory is to be used in the box counting, simply use
 2144   $(\PB{\\{global\_llx}},\PB{\\{global\_lly}})=({\rm min}[\PB{\\{x\_traj}}],{\rm
 2145   min}[\PB{\\{y\_traj}}])$ and
 2146   $(\PB{\\{global\_urx}},\PB{\\{global\_ury}})=({\rm max}[\PB{\\{x\_traj}}],{\rm
 2147   max}[\PB{\\{y\_traj}}])$ for the
 2148   specification of the computational domain.
 2149   \medskip\noindent
 2150   Input variables:
 2151   \varitem[\PB{\\{mm}}]{The total number of coordinates forming the input
 2152   trajectory,
 2153   or equivalently the length of the vectors \PB{${*}\\{x\_traj}$} and \PB{${*}%
 2154   \\{y\_traj}$}.}
 2155   \varitem[\PB{${*}\\{x\_traj}$}, \PB{${*}\\{y\_traj}$}]{Vectors of length \PB{%
 2156   \\{mm}} containing the $x$- and
 2157   $y$-coordinates of the input trajectory.}
 2158   \varitem[\PB{\\{resolution}}]{The grid refinement factor $N$.}
 2159   \varitem[\PB{\\{global\_llx}}, \PB{\\{global\_lly}}]{Coordinates of the
 2160   lower-left corner
 2161   of the computational domain in which the box-counting is to be performed.}
 2162   \varitem[\PB{\\{global\_urx}}, \PB{\\{global\_ury}}]{Coordinates of the
 2163   upper-right corner
 2164   of the computational domain in which the box-counting is to be performed.}
 2165   \varitem[\PB{\\{trajectory\_filename}}]{String containing the name of the file
 2166   containing the input trajectory.}
 2167   \varitem[\PB{\\{boxgraph\_filename}}]{String which, if non-empty, states the
 2168   file name
 2169   to which the map of the spatial distribution of the covering boxes is to be
 2170   written, as \METAPOST\ code.}
 2171   \medskip\noindent
 2172   Output variables:
 2173   \varitem[]{On exit, the routine returns the number of covering boxes
 2174   as an integer of \PB{\&{long} \&{unsigned}} precision.}
 2175   \medskip\noindent
 2176   Internal variables:
 2177   \varitem[\PB{\\{px}},\PB{\\{py}},\PB{\\{qx}},\PB{\\{qy}}]{Keeps track of the
 2178   $x$- and $y$-coordinates of
 2179   the start and end points of line segments, between ${\bf p}=(p_x,p_y)$ and
 2180   ${\bf q}=(q_x,q_y)$.}
 2181   \medskip
 2182   
 2183   \Y\B\4\X42:Routine for calculation of number of boxes covering the trajectory%
 2184   \X${}\E{}$\6
 2185   \&{long} \&{unsigned} \&{int} \\{get\_num\_covering\_boxes\_with\_boxmaps}(%
 2186   \&{double} ${}{*}\\{x\_traj},\39{}$\&{double} ${}{*}\\{y\_traj},\39{}$\&{long} %
 2187   \&{int} \\{mm}${},\39{}$\&{int} \\{resolution}${},\39{}$\&{double} \\{global%
 2188   \_llx}${},\39{}$\&{double} \\{global\_lly}${},\39{}$\&{double} \\{global%
 2189   \_urx}${},\39{}$\&{double} \\{global\_ury}${},\39{}$\&{char} \\{boxgraph%
 2190   \_filename}[\,]${},\39{}$\&{double} \\{boxgraph\_width}${},\39{}$\&{double} %
 2191   \\{boxgraph\_height}${},\39{}$\&{char} \\{trajectory\_filename}[\,])\1\1\2\2\6
 2192   ${}\{{}$\1\6
 2193   \&{short} \&{int} ${}{*}{*}\\{box};{}$\6
 2194   \&{long} \&{unsigned} \&{int} \|i${},{}$ \|j${},{}$ \|m${},{}$ \\{nn}${},{}$ %
 2195   \\{istart}${},{}$ \\{jstart}${},{}$ \\{istop}${},{}$ \\{jstop}${},{}$ %
 2196   \\{iincr}${},{}$ \\{jincr}${},{}$ \\{num\_boxes};\6
 2197   \&{double} ${}{*}\\{x\_box},{}$ ${}{*}\\{y\_box}{}$;\C{ Keeps track of the
 2198   lower-left coordinates of                             the boxes. }\6
 2199   \&{double} \\{px}${},{}$ \\{py}${},{}$ \\{qx}${},{}$ \\{qy};\6
 2200   \&{FILE} ${}{*}\\{boxgraph\_file};{}$\7
 2201   \X43:Set up the grid of $2^N\times2^N$ boxes covering the entire global window
 2202   of computation\X\6
 2203   \X44:Find indices $(i_a,j_a)$ of the box covering the first coordinate of the
 2204   trajectory\X\6
 2205   \&{for} ${}(\|m\K\T{1};{}$ ${}\|m\Z\\{mm}-\T{1};{}$ ${}\|m\PP){}$\5
 2206   ${}\{{}$\1\6
 2207   \X45:Find indices $(i_b,j_b)$ of the box covering the end point of $m$th
 2208   trajectory segment\X\6
 2209   \X46:Scan the $m$th trajectory segment for intersecting boxes\X\6
 2210   \4${}\}{}$\2\6
 2211   \X47:Open file for output of box distribution graph\X\6
 2212   \X48:Write the input trajectory to the box distribution graph\X\6
 2213   \X49:Count the total number of boxes \PB{\\{num\_boxes}} intersecting the
 2214   trajectory\X\6
 2215   \X50:Write closing blocks the box distribution graph\X\6
 2216   \X51:Close any open file for output of box distribution graph\X\6
 2217   \&{return} (\\{num\_boxes});\6
 2218   \4${}\}{}$\2\par
 2219   \U41.\fi
 2220   
 2221   \M{43}Set up the grid of $2^N\times2^N$ boxes covering the entire global window
 2222   of computation.
 2223   In this block, the resolution of the grid of boxes is defined. Notice that in
 2224   many cases, only a certain subset of boxes will actally intersect the input
 2225   trajectory. However, we do not \'a priori know this number of boxes, and in
 2226   order to safeguard against the possible danger of running out of allocated
 2227   memory, with the need for a dynamic memory allocation depending on the state
 2228   of computation, a matrix of short integer precision is allocated covering the
 2229   entire computational window.
 2230   In order to keep track of the coordinates of the boxes, two vectors
 2231   $\PB{\\{x\_box}}[1\ldots 2^N]$ and $\PB{\\{y\_box}}[1\ldots 2^N]$ are allocated
 2232   to contain
 2233   the $x$- and $y$-coordinates of the lower-left corners of the boxes, with the
 2234   last elements $\PB{\\{x\_box}}[2^N]$ and $\PB{\\{y\_box}}[2^N]$ containing the
 2235   upper-right
 2236   corner coordinates of the upper-right-most box of the global window.
 2237   
 2238   \Y\B\4\X43:Set up the grid of $2^N\times2^N$ boxes covering the entire global
 2239   window of computation\X${}\E{}$\6
 2240   ${}\{{}$\1\6
 2241   ${}\\{nn}\K\T{1};{}$\6
 2242   \&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{resolution};{}$ ${}\|i\PP){}$\1\5
 2243   ${}\\{nn}\K\T{2}*\\{nn};{}$\2\6
 2244   ${}\\{box}\K\\{simatrix}(\T{1},\39\\{nn},\39\T{1},\39\\{nn});{}$\6
 2245   \&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{nn};{}$ ${}\|i\PP){}$\1\6
 2246   \&{for} ${}(\|j\K\T{1};{}$ ${}\|j\Z\\{nn};{}$ ${}\|j\PP){}$\1\5
 2247   ${}\\{box}[\|i][\|j]\K\T{0};{}$\2\2\6
 2248   ${}\\{x\_box}\K\\{dvector}(\T{1},\39\\{nn}+\T{1});{}$\6
 2249   ${}\\{y\_box}\K\\{dvector}(\T{1},\39\\{nn}+\T{1});{}$\6
 2250   \&{for} ${}(\|m\K\T{1};{}$ ${}\|m\Z\\{nn}+\T{1};{}$ ${}\|m\PP){}$\5
 2251   ${}\{{}$\1\6
 2252   ${}\\{x\_box}[\|m]\K\\{global\_llx}+((\&{double})(\|m-\T{1}))*(\\{global\_urx}-%
 2253   \\{global\_llx})/((\&{double})(\\{nn}));{}$\6
 2254   ${}\\{y\_box}[\|m]\K\\{global\_lly}+((\&{double})(\|m-\T{1}))*(\\{global\_ury}-%
 2255   \\{global\_lly})/((\&{double})(\\{nn}));{}$\6
 2256   \4${}\}{}$\2\6
 2257   \4${}\}{}$\2\par
 2258   \U42.\fi
 2259   
 2260   \M{44}Find indices $(i_a,j_a)$ of the box covering the first coordinate of the
 2261   trajectory.
 2262   
 2263   \Y\B\4\X44:Find indices $(i_a,j_a)$ of the box covering the first coordinate of
 2264   the trajectory\X${}\E{}$\6
 2265   ${}\{{}$\1\6
 2266   ${}\\{istart}\K\T{0};{}$\6
 2267   ${}\\{jstart}\K\T{0};{}$\6
 2268   \&{for} ${}(\|m\K\T{1};{}$ ${}\|m\Z\\{nn};{}$ ${}\|m\PP){}$\5
 2269   ${}\{{}$\1\6
 2270   \&{if} ${}((\\{x\_box}[\|m]\Z\\{x\_traj}[\T{1}])\W(\\{x\_traj}[\T{1}]\Z\\{x%
 2271   \_box}[\|m+\T{1}])){}$\1\5
 2272   ${}\\{istart}\K\|m;{}$\2\6
 2273   \&{if} ${}((\\{y\_box}[\|m]\Z\\{y\_traj}[\T{1}])\W(\\{y\_traj}[\T{1}]\Z\\{y%
 2274   \_box}[\|m+\T{1}])){}$\1\5
 2275   ${}\\{jstart}\K\|m;{}$\2\6
 2276   \4${}\}{}$\2\6
 2277   \&{if} ${}((\\{istart}\E\T{0})\V(\\{jstart}\E\T{0})){}$\5
 2278   ${}\{{}$\1\6
 2279   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error!\ Cannot\ f}\)\.{ind\ box\ indices%
 2280   \ of\ 1}\)\.{st\ coordinate!\\n"},\39\\{progname});{}$\6
 2281   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Please\ check\ da}\)\.{ta\ for\ input\
 2282   trajeto}\)\.{ry.\\n"},\39\\{progname});{}$\6
 2283   \\{exit}(\.{FAILURE});\6
 2284   \4${}\}{}$\2\6
 2285   \4${}\}{}$\2\par
 2286   \U42.\fi
 2287   
 2288   \M{45}Find indices $(i_b,j_b)$ of the box covering the end point of the $m$th
 2289   trajectory segment.
 2290   
 2291   \Y\B\4\X45:Find indices $(i_b,j_b)$ of the box covering the end point of $m$th
 2292   trajectory segment\X${}\E{}$\6
 2293   ${}\{{}$\1\6
 2294   ${}\\{px}\K\\{x\_traj}[\|m];{}$\6
 2295   ${}\\{py}\K\\{y\_traj}[\|m];{}$\6
 2296   ${}\\{qx}\K\\{x\_traj}[\|m+\T{1}];{}$\6
 2297   ${}\\{qy}\K\\{y\_traj}[\|m+\T{1}];{}$\6
 2298   ${}\\{istop}\K\\{istart};{}$\6
 2299   ${}\\{jstop}\K\\{jstart};{}$\6
 2300   \&{if} ${}(\\{px}<\\{qx}){}$\5
 2301   ${}\{{}$\1\6
 2302   ${}\\{iincr}\K\T{1};{}$\6
 2303   \&{while} ${}(\\{x\_box}[\\{istop}+\T{1}]<\\{qx}){}$\1\5
 2304   ${}\\{istop}\PP;{}$\2\6
 2305   \4${}\}{}$\2\6
 2306   \&{else}\5
 2307   ${}\{{}$\1\6
 2308   ${}\\{iincr}\K{-}\T{1};{}$\6
 2309   \&{while} ${}(\\{x\_box}[\\{istop}]>\\{qx}){}$\1\5
 2310   ${}\\{istop}\MM;{}$\2\6
 2311   \4${}\}{}$\2\6
 2312   \&{if} ${}(\\{py}<\\{qy}){}$\5
 2313   ${}\{{}$\1\6
 2314   ${}\\{jincr}\K\T{1};{}$\6
 2315   \&{while} ${}(\\{y\_box}[\\{jstop}+\T{1}]<\\{qy}){}$\1\5
 2316   ${}\\{jstop}\PP;{}$\2\6
 2317   \4${}\}{}$\2\6
 2318   \&{else}\5
 2319   ${}\{{}$\1\6
 2320   ${}\\{jincr}\K{-}\T{1};{}$\6
 2321   \&{while} ${}(\\{y\_box}[\\{jstop}]>\\{qy}){}$\1\5
 2322   ${}\\{jstop}\MM;{}$\2\6
 2323   \4${}\}{}$\2\6
 2324   \&{if} ${}(\T{0}\E\T{1}){}$\5
 2325   ${}\{{}$\1\6
 2326   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Endpoint\ box\ in}\)\.{dices:\ (i=%
 2327   \%ld,j=\%ld)}\)\.{\\n"},\39\\{progname},\39\\{istop},\39\\{jstop});{}$\6
 2328   \4${}\}{}$\2\6
 2329   \4${}\}{}$\2\par
 2330   \U42.\fi
 2331   
 2332   \M{46}Scan the $m$th trajectory segment for intersecting boxes.
 2333   As the indices of the boxes covering the start and end points of the $m$th
 2334   trajectory segment have been previously determined, one may now use this
 2335   information in order to reduce the search for intersecting boxes to the
 2336   domain as covered by box indices $i_{\rm start}\le i\le i_{\rm stop}$ and
 2337   $j_{\rm start}\le j\le j_{\rm stop}$.
 2338   This way, the computational time needed is greatly reduced as compared to if
 2339   the entire window would be scanned for every segment.
 2340   
 2341   \Y\B\4\X46:Scan the $m$th trajectory segment for intersecting boxes\X${}\E{}$\6
 2342   ${}\{{}$\1\6
 2343   \&{for} ${}(\|i\K\\{istart};{}$ ${}\|i\I(\\{istop}+\\{iincr});{}$ ${}\|i\MRL{+{%
 2344   \K}}\\{iincr}){}$\5
 2345   ${}\{{}$\1\6
 2346   \&{for} ${}(\|j\K\\{jstart};{}$ ${}\|j\I(\\{jstop}+\\{jincr});{}$ ${}\|j\MRL{+{%
 2347   \K}}\\{jincr}){}$\5
 2348   ${}\{{}$\1\6
 2349   \&{if} ${}(\\{box\_intersection}(\\{px},\39\\{py},\39\\{qx},\39\\{qy},\39\\{x%
 2350   \_box}[\|i],\39\\{y\_box}[\|j],\39\\{x\_box}[\|i+\T{1}],\39\\{y\_box}[\|j+%
 2351   \T{1}])){}$\5
 2352   ${}\{{}$\1\6
 2353   ${}\\{box}[\|i][\|j]\K\T{1};{}$\6
 2354   \4${}\}{}$\2\6
 2355   \4${}\}{}$\2\6
 2356   \4${}\}{}$\2\6
 2357   ${}\\{istart}\K\\{istop};{}$\6
 2358   ${}\\{jstart}\K\\{jstop};{}$\6
 2359   \4${}\}{}$\2\par
 2360   \U42.\fi
 2361   
 2362   \M{47}Open file for output of box distribution graph. In this block the
 2363   preamble
 2364   of the output \METAPOST\ code is written to file. This preamble contains a
 2365   macro {\tt box(i,j)} which simply is a parameter-specific \METAPOST\ subroutine
 2366   which is used in order to reduce the final size of the code to be compiled
 2367   into a graph over the distribution of covering boxes.
 2368   
 2369   \Y\B\4\X47:Open file for output of box distribution graph\X${}\E{}$\6
 2370   ${}\{{}$\1\6
 2371   \&{if} ${}(\R\\{strcmp}(\\{boxgraph\_filename},\39\.{""})){}$\5
 2372   ${}\{{}$\C{ is \PB{\\{boxgraph\_filename}} an empty string? }\1\6
 2373   ${}\\{boxgraph\_file}\K\NULL;{}$\6
 2374   \4${}\}{}$\2\6
 2375   \&{else}\5
 2376   ${}\{{}$\1\6
 2377   \&{if} ${}((\\{boxgraph\_file}\K\\{fopen}(\\{boxgraph\_filename},\39\.{"w"}))\E%
 2378   \NULL){}$\5
 2379   ${}\{{}$\1\6
 2380   ${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Could\ not\ open\ }\)\.{\%s\ for\ box\
 2381   graphs!\\n}\)\.{"},\39\\{progname},\39\\{boxgraph\_filename});{}$\6
 2382   \\{exit}(\.{FAILURE});\6
 2383   \4${}\}{}$\2\6
 2384   ${}\\{fseek}(\\{boxgraph\_file},\39\T{0\$L},\39\.{SEEK\_SET});{}$\6
 2385   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"input\ graph;\\n"});{}$\6
 2386   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"def\ box(expr\ i,j)=\\}\)\.{n"});{}$\6
 2387   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"begingroup\\n"});{}$\6
 2388   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"llx:=\%2.8f+(i-1)*\%2}\)\.{.8f;\\n"},%
 2389   \39\\{global\_llx},\39(\\{global\_urx}-\\{global\_llx})/(\\{nn}));{}$\6
 2390   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"lly:=\%2.8f+(j-1)*\%2}\)\.{.8f;\\n"},%
 2391   \39\\{global\_lly},\39(\\{global\_ury}-\\{global\_lly})/(\\{nn}));{}$\6
 2392   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"urx:=\%2.8f+(i)*\%2.8}\)\.{f;\\n"},\39%
 2393   \\{global\_llx},\39(\\{global\_urx}-\\{global\_llx})/(\\{nn}));{}$\6
 2394   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"ury:=\%2.8f+(j)*\%2.8}\)\.{f;\\n"},\39%
 2395   \\{global\_lly},\39(\\{global\_ury}-\\{global\_lly})/(\\{nn}));{}$\6
 2396   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (llx,lly)--(u}\)\.{rx,lly);%
 2397   \\n"});{}$\6
 2398   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (urx,lly)--(u}\)\.{rx,ury);%
 2399   \\n"});{}$\6
 2400   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (urx,ury)--(l}\)\.{lx,ury);%
 2401   \\n"});{}$\6
 2402   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (llx,ury)--(l}\)\.{lx,lly);%
 2403   \\n"});{}$\6
 2404   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"endgroup\\n"});{}$\6
 2405   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"enddef;\\n"});{}$\6
 2406   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"beginfig(1);\\n"});{}$\6
 2407   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"w:=\%2.4fmm;\ h:=\%2.4}\)\.{fmm;\\n"},%
 2408   \39\\{boxgraph\_width},\39\\{boxgraph\_height});{}$\6
 2409   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"draw\ begingraph(w,h}\)\.{);\\n"});{}$%
 2410   \6
 2411   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"pickup\ pencircle\ sc}\)\.{aled\ .3pt;%
 2412   \\n"});{}$\6
 2413   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"setrange(\%2.6f,\%2.6}\)\.{f,\%2.6f,%
 2414   \%2.6f);\\n"},\39\\{global\_llx},\39\\{global\_lly},\39\\{global\_urx},\39%
 2415   \\{global\_ury});{}$\6
 2416   \4${}\}{}$\2\6
 2417   \4${}\}{}$\2\par
 2418   \U42.\fi
 2419   
 2420   \M{48}Write the input trajectory to the box distribution graph. Here there are
 2421   two
 2422   possible choices for how the input trajectory is to be included in the box
 2423   graph.
 2424   
 2425   First, we may use \MP\ to automatically scan and draw the trajectory
 2426   for us, simply by using a statment like {\tt gdraw input.dat}, assuming that
 2427   the file {\tt input.dat} contains properly formatted columnwise data.
 2428   However, this choice would have two major drawbacks, namely that the generated
 2429   code would be dependent on that the original input file always is in place,
 2430   hence not allowing the \MP\ code to be exported as a standalone code as such,
 2431   and also that this would put a limitation on the number of nodes allowed in the
 2432   input trajectory, as the {\tt graph.mp} macro package of \MP\ only accepts
 2433   roughly 4000 points before it cuts the mapping.
 2434   
 2435   The other choice is to include the input trajectory directly into the generated
 2436   code, preferrably by chopping the trajectory into pieces small enough to easily
 2437   be mapped by \MP\ without reaching a critical limit of exhaustion.
 2438   This choice of course significantly increases the file size of the generated
 2439   code, but this is a price we will have to accept in order to get stand-alone
 2440   output.
 2441   In the \boxcount\ program, the second alternative was naturtally chosen, simply
 2442   because the author is a fan of self-sustaining code which can be exported for
 2443   later use.
 2444   
 2445   In this block, the status of the pointer \PB{\\{boxgraph\_file}} is used to
 2446   determine
 2447   whether to write the trajectory to file or not. If \PB{\\{boxgraph\_file}}
 2448   equals to
 2449   \PB{$\NULL$} (NULL), then the \boxcount\ program will not attempt to write the
 2450   input
 2451   trajectory to file.
 2452   
 2453   \Y\B\4\X48:Write the input trajectory to the box distribution graph\X${}\E{}$\6
 2454   ${}\{{}$\1\6
 2455   \&{if} ${}(\\{boxgraph\_file}\I\NULL){}$\5
 2456   ${}\{{}$\1\6
 2457   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"pickup\ pencircle\ sc}\)\.{aled\ .5pt;%
 2458   \\n"});{}$\6
 2459   ${}\|i\K\T{0};{}$\6
 2460   \&{for} ${}(\|m\K\T{1};{}$ ${}\|m\Z\\{mm};{}$ ${}\|m\PP){}$\5
 2461   ${}\{{}$\1\6
 2462   \&{if} ${}(\|i\E\T{0}){}$\5
 2463   ${}\{{}$\1\6
 2464   \&{if} ${}(\|m\E\T{1}){}$\5
 2465   ${}\{{}$\1\6
 2466   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (\%2.4f,\%2.4f)}\)\.{"},\39\\{x%
 2467   \_traj}[\|m],\39\\{y\_traj}[\|m]);{}$\6
 2468   \4${}\}{}$\2\6
 2469   \&{else} \&{if} ${}(\T{2}<\\{mm}){}$\5
 2470   ${}\{{}$\1\6
 2471   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (\%2.4f,\%2.4f)}\)\.{--(\%2.4f,%
 2472   \%2.4f)"},\39\\{x\_traj}[\|m-\T{1}],\39\\{y\_traj}[\|m-\T{1}],\39\\{x\_traj}[%
 2473   \|m],\39\\{y\_traj}[\|m]);{}$\6
 2474   \4${}\}{}$\2\6
 2475   \4${}\}{}$\2\6
 2476   \&{else}\5
 2477   ${}\{{}$\1\6
 2478   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"--(\%2.4f,\%2.4f)"},\39\\{x\_traj}[%
 2479   \|m],\39\\{y\_traj}[\|m]);{}$\6
 2480   \4${}\}{}$\2\6
 2481   ${}\|i\PP;{}$\6
 2482   \&{if} ${}((\|i\E\T{5})\V(\|m\E\\{mm})){}$\5
 2483   ${}\{{}$\1\6
 2484   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{";\\n"});{}$\6
 2485   ${}\|i\K\T{0};{}$\6
 2486   \4${}\}{}$\2\6
 2487   \4${}\}{}$\2\6
 2488   \4${}\}{}$\2\6
 2489   \4${}\}{}$\2\par
 2490   \U42.\fi
 2491   
 2492   \M{49}Count the total number of boxes \PB{\\{num\_boxes}} intersecting the
 2493   trajectory.
 2494   Having traversed the entire trajectory at the current depth of resolution,
 2495   the only remaining task is to sum up the total number of boxes needed to
 2496   cover the entire trajectory.
 2497   This is simply done by extracting the number of set elements in the \PB{%
 2498   \\{box}}
 2499   matrix.
 2500   Having extracted the total number of boxes, this block also takes care of
 2501   releasing the memory occupied by the \PB{\\{box}} matrix, as this memory might
 2502   be
 2503   needed for iterations to come in which an even finer mesh of boxes is used.
 2504   
 2505   \Y\B\4\X49:Count the total number of boxes \PB{\\{num\_boxes}} intersecting the
 2506   trajectory\X${}\E{}$\6
 2507   ${}\{{}$\1\6
 2508   ${}\\{num\_boxes}\K\T{0};{}$\6
 2509   \&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{nn};{}$ ${}\|i\PP){}$\5
 2510   ${}\{{}$\1\6
 2511   \&{for} ${}(\|j\K\T{1};{}$ ${}\|j\Z\\{nn};{}$ ${}\|j\PP){}$\5
 2512   ${}\{{}$\1\6
 2513   \&{if} ${}(\\{box}[\|i][\|j]\E\T{1}){}$\5
 2514   ${}\{{}$\1\6
 2515   ${}\\{num\_boxes}\PP;{}$\6
 2516   \&{if} ${}(\\{boxgraph\_file}\I\NULL){}$\5
 2517   ${}\{{}$\1\6
 2518   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"box(\%ld,\%ld);\\n"},\39\|i,\39%
 2519   \|j);{}$\6
 2520   \4${}\}{}$\2\6
 2521   \4${}\}{}$\2\6
 2522   \4${}\}{}$\2\6
 2523   \4${}\}{}$\2\6
 2524   ${}\\{free\_simatrix}(\\{box},\39\T{1},\39\\{nn},\39\T{1},\39\\{nn});{}$\6
 2525   \4${}\}{}$\2\par
 2526   \U42.\fi
 2527   
 2528   \M{50}Write closing blocks the box distribution graph. Here follows the
 2529   specification of tick marking (set as automatic for the sake of simplicity)
 2530   and axis labels, just before the closing statements of the \METAPOST\ code
 2531   for the graphs of box distributions.
 2532   
 2533   \Y\B\4\X50:Write closing blocks the box distribution graph\X${}\E{}$\6
 2534   ${}\{{}$\1\6
 2535   \&{if} ${}(\\{boxgraph\_file}\I\NULL){}$\5
 2536   ${}\{{}$\1\6
 2537   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"autogrid(itick\ bot,}\)\.{itick\ lft);%
 2538   \\n"});{}$\6
 2539   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"glabel.bot(btex\ \$x\$}\)\.{\
 2540   etex,OUT);\\n"});{}$\6
 2541   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"glabel.lft(btex\ \$y\$}\)\.{\
 2542   etex,OUT);\\n"});{}$\6
 2543   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"endgraph;\\n"});{}$\6
 2544   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"endfig;\\n"});{}$\6
 2545   ${}\\{fprintf}(\\{boxgraph\_file},\39\.{"end\\n"});{}$\6
 2546   \4${}\}{}$\2\6
 2547   \4${}\}{}$\2\par
 2548   \U42.\fi
 2549   
 2550   \M{51}Close any open file for output of box distribution graph.
 2551   
 2552   \Y\B\4\X51:Close any open file for output of box distribution graph\X${}\E{}$\6
 2553   ${}\{{}$\1\6
 2554   \&{if} ${}(\\{boxgraph\_file}\I\NULL){}$\5
 2555   ${}\{{}$\1\6
 2556   ${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ MetaPost\ output}\)\.{\ box\
 2557   distribution\ gr}\)\.{aph\ written\ to\ \%s.\\n}\)\.{"},\39\\{progname},\39%
 2558   \\{boxgraph\_filename});{}$\6
 2559   \\{fclose}(\\{boxgraph\_file});\6
 2560   \4${}\}{}$\2\6
 2561   \4${}\}{}$\2\par
 2562   \U42.\fi
 2563   
 2564   \M{52}Routine for calculation of number of boxes covering the trajectory. This
 2565   routine provides a simplified interface to the general boxcountng routine
 2566   \PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}}, with the only difference
 2567   that no
 2568   graph of the box distribution over the trajectory is generated.
 2569   The \PB{\\{get\_num\_covering\_boxes}(\,)} routine takes a trajectory as
 2570   described by a discrete set of coordinates as input, and for a given grid
 2571   refinement $N$ returns the number of boxes needed to entirely cover the
 2572   trajectory. The grid refinement factor $N$ indicates that the domain of
 2573   computation is divided into a $[2^N\times2^N]$-grid of boxes.
 2574   
 2575   The computational domain in which the box counting is to be performed is
 2576   explicitly stated by the coordinates of its lower-left and upper-right
 2577   corners, $(global_llx,global_lly)$ and $(global_urx,global_ury)$, respectively.
 2578   Parts of the trajectory which are outside of this domain are not included in
 2579   the box-counting.
 2580   If the entire input trajectory is to be used in the box counting, simply use
 2581   $(\PB{\\{global\_llx}},\PB{\\{global\_lly}})=({\rm min}[\PB{\\{x\_traj}}],{\rm
 2582   min}[\PB{\\{y\_traj}}])$ and
 2583   $(\PB{\\{global\_urx}},\PB{\\{global\_ury}})=({\rm max}[\PB{\\{x\_traj}}],{\rm
 2584   max}[\PB{\\{y\_traj}}])$ for the
 2585   specification of the computational domain.
 2586   \medskip\noindent
 2587   Input variables:
 2588   \varitem[\PB{\\{mm}}]{The total number of coordinates forming the input
 2589   trajectory,
 2590   or equivalently the length of the vectors \PB{${*}\\{x\_traj}$} and \PB{${*}%
 2591   \\{y\_traj}$}.}
 2592   \varitem[\PB{${*}\\{x\_traj}$}, \PB{${*}\\{y\_traj}$}]{Vectors of length \PB{%
 2593   \\{mm}} containing the $x$- and
 2594   $y$-coordinates of the input trajectory.}
 2595   \varitem[\PB{\\{resolution}}]{The grid refinement factor $N$.}
 2596   \varitem[\PB{\\{global\_llx}}, \PB{\\{global\_lly}}]{Coordinates of the
 2597   lower-left corner
 2598   of the computational domain in which the box-counting is to be performed.}
 2599   \varitem[\PB{\\{global\_urx}}, \PB{\\{global\_ury}}]{Coordinates of the
 2600   upper-right corner
 2601   of the computational domain in which the box-counting is to be performed.}
 2602   \medskip\noindent
 2603   Output variables:
 2604   \varitem[]{On exit, the routine returns the number of covering boxes
 2605   as an integer of \PB{\&{long} \&{unsigned}} precision.}
 2606   \medskip
 2607   
 2608   \Y\B\4\X52:Simplified routine for calculation of number of boxes covering the
 2609   trajectory\X${}\E{}$\6
 2610   \&{long} \&{unsigned} \&{int} \\{get\_num\_covering\_boxes}(\&{double} ${}{*}%
 2611   \\{x\_traj},\39{}$\&{double} ${}{*}\\{y\_traj},\39{}$\&{long} \&{int} %
 2612   \\{mm}${},\39{}$\&{int} \|i${},\39{}$\&{double} \\{global\_llx}${},\39{}$%
 2613   \&{double} \\{global\_lly}${},\39{}$\&{double} \\{global\_urx}${},\39{}$%
 2614   \&{double} \\{global\_ury})\1\1\2\2\6
 2615   ${}\{{}$\1\6
 2616   \&{return} ${}(\\{get\_num\_covering\_boxes\_with\_boxmaps}(\\{x\_traj},\39\\{y%
 2617   \_traj},\39\\{mm},\39\|i,\39\\{global\_llx},\39\\{global\_lly},\39\\{global%
 2618   \_urx},\39\\{global\_ury},\39\.{""},\39\T{0.0},\39\T{0.0},\39\.{""}));{}$\6
 2619   \4${}\}{}$\2\par
 2620   \U41.\fi
 2621   
 2622   \N{1}{53}Index.
 2623   \fi
 2624   
 2625   \inx
 2626   \fin
 2627   \con
 2628   

Return to previous page

Generated by ::viewsrc::

Last modified Saturday 10 Dec 2011