Search:

Return to previous page

Contents of file 'boxcount/boxcount.w':




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

Return to previous page

Generated by ::viewsrc::

Last modified Saturday 10 Dec 2011