\input cwebmac
% File:        boxcount.w [CWEB source code]
% Created:     May 8, 2006 [v.1.0]
% Last change: November 26, 2011 [v.1.6]
% Author:      Fredrik Jonsson
% Description: The CWEB source code for the BOXCOUNT simulator of nonlinear
%              magneto-optical Bragg gratings. For information on the CWEB
%              programming language, see http://www.literateprogramming.com.
% Compilation: Compile this program by using the enclosed Makefile, or use
%              the blocks of the Makefile as listed in the documentation
%              (file boxcount.ps or boxcount.pdf). The C source code (as
%              generated from this CWEB code) conforms to the ANSI standard
%              for the C programming language (which is equivalent to the
%              ISO C90 standard for C).
%
% Copyright (C) 2006-2011, Fredrik Jonsson. Non-commercial copying welcome.
%
\input epsf
\def\version{1.6}
\def\lastrevdate{November 26, 2011}
\font\eightcmr=cmr8
\font\tensc=cmcsc10
\font\eightcmssq=cmssq8
\font\eightcmssqi=cmssqi8
\font\twentycmcsc=cmcsc10 at 20 truept
\def\boxcount{{\eightcmr BOXCOUNT\spacefactor1000}}
\def\poincare{{\eightcmr POINCARE\spacefactor1000}}
\def\ANSI{{\eightcmr ANSI\spacefactor1000}} % The language standard we stick to
\def\SI{{\eightcmr SI\spacefactor1000}}   % Another standard for physical units
\def\GNU{{\eightcmr GNU\spacefactor1000}}       % GNU is Not Unix
\def\GCC{{\eightcmr GCC\spacefactor1000}}       % The GNU C-compiler
\def\CEE{{\eightcmr C\spacefactor1000}}         % The C programming language
\def\ANSICEE{{\eightcmr ANSI~C\spacefactor1000}}% The ANSI-C standard language
\def\CWEB{{\eightcmr CWEB\spacefactor1000}}     % The CWEB programming language
\def\MATLAB{{\eightcmr MATLAB\spacefactor1000}} % The MATLAB ditto
\def\UNIX{{\eightcmr UNIX\spacefactor1000}}
\def\CYGWIN{{\eightcmr CYGWIN\spacefactor1000}}
\def\CTANGLE{{\eightcmr CTANGLE\spacefactor1000}}
\def\CWEAVE{{\eightcmr CWEAVE\spacefactor1000}}
\def\OSX{{OS\,X}}
\def\re{\mathop{\rm Re}\nolimits} % The real part of a complex number
\def\im{\mathop{\rm Im}\nolimits} % The imaginary part of a complex number
\def\dollar{\char'044}             % The `$' character
\def\tothepower{\char'136}         % The `^' character
\def\onehalf{{\textstyle{{1}\over{2}}}}
\def\threefourth{{\textstyle{{3}\over{4}}}}
\def\endalg{\vrule height 1.4ex width .6ex depth .4ex} % Rectangular bullet
\def\ie{i.\thinspace{e.}~\ignorespaces}
\def\eg{e.\thinspace{g.}~\ignorespaces}
\let\,\thinspace
%% \def\subsection[#1]{\goodbreak\noindent{\it #1}\par\nobreak\noindent}
%
% Define a handy macro for listing the steps of an algorithm.
%
\newdimen\aitemindent \aitemindent=26pt
\newdimen\aitemleftskip \aitemleftskip=36pt
\def\aitem[#1]{\smallbreak\noindent\hbox to 10pt{}%
  \hbox to\aitemindent{\bf #1\hfill}%
  \hangindent\aitemleftskip\ignorespaces}
%
% Define a handy macro for bibliographic references.
%
\newdimen\refitemindent \refitemindent=18pt
\def\refitem[#1]{\smallbreak\noindent%
  \hbox to\refitemindent{[#1]\hfill}%
  \hangindent\refitemindent\ignorespaces}
\def\references{\medskip\goodbreak\noindent{\bf References}\nobreak\smallskip}
%
% Define a handy macro for nicely typeset variable descriptions.
%
\newdimen\varitemindent \varitemindent=100pt
\def\varitem[#1]#2{\smallbreak\noindent\hbox to 20pt{}%
  \hbox to\varitemindent{#1\hfill}%
  \hangindent 120pt\ignorespaces#2\par}
%
% Define a handy macro for nicely typeset descriptions of command line options.
%
\newdimen\optitemindent \optitemindent=80pt
\def\optitem[#1]#2{\smallbreak\noindent\hbox to\hsize{\hskip40pt{#1}\hfill}
  \hangindent\optitemindent\ignorespaces{\hskip\optitemindent{#2}}\smallskip}
%
% Define a handy macro for the list of program revisions.
%
\newdimen\citemindent \citemindent=80pt
\newdimen\citemleftskip \citemleftskip=90pt
\def\citem[#1]{\smallbreak\noindent\hbox to 10pt{}%
  \hbox to\citemindent{\bf #1\hfill}%
  \hangindent\citemleftskip\ignorespaces}
\def\citindent{\smallbreak\noindent\hbox to 10pt{}%
  \hbox to\citemindent{\hfil}%
  \hangindent\citemleftskip\ignorespaces}
%
% Define the \beginvrulealign and \endvrulealign macros as described in
% Donald Knuth's The TeXbook, Appendix D: Dirty Tricks. These macros are
% used in typesetting nicely looking tables.
%
\def\beginvrulealign{\setbox0=\vbox\bgroup}
\def\endvrulealign{\egroup % now \box0 holds the entire alignment
  \setbox0=\vbox{\setbox2=\hbox{\vrule height\ht0 depth\dp0 width0pt}
    \unvbox0 \setbox0=\lastbox % now \box0 is the bottom row
    \nointerlineskip \copy0 % put it back
    \global\setbox1=\hbox{} % initialize box that will contain rules
    \setbox4=\hbox{\unhbox0 % now open up the bottom row
      \loop \skip0=\lastskip \unskip % remove tabskip glue
      \advance\skip0 by-.4pt % rules are .4 pt wide
      \divide\skip0 by 2
      \global\setbox1=\hbox{\hskip\skip0\vrule\hskip\skip0
        \unhbox2\unhbox1}%
      \setbox2=\lastbox % remove alignment entry
      \ifhbox2 \setbox2=\hbox{\kern\wd2}\repeat}}%
  \hbox{\rlap{\box0}\box1}} % superimpose the alignment on the rules
%
% Make sure that the METAPOST logo can be loaded even in plain TeX.
%
\ifx\MP\undefined
    \ifx\manfnt\undefined
            \font\manfnt=logo10
    \fi
    \ifx\manfntsl\undefined
            \font\manfntsl=logosl10
    \fi
    \def\MP{{\ifdim\fontdimen1\font>0pt \let\manfnt = \manfntsl \fi
      {\manfnt META}\-{\manfnt POST}}\spacefactor1000 }%
\fi
\ifx\METAPOST\undefined \let\METAPOST=\MP \fi

\datethis


\N{1}{1}Introduction.
\vskip 100pt
\centerline{\twentycmcsc Boxcount}
\vskip 20pt
\centerline{A program for calculating box-counting estimates to the fractal
dimension of curves in the plane.}
\vskip 2pt
\centerline{(Version \version\ of \lastrevdate)}
\vskip 10pt
\centerline{Written by Fredrik Jonsson}
\vskip 10pt
\centerline{\epsfxsize=88mm\epsfbox{figures/figure-05.1}}
\vskip 10pt
\noindent
This \CWEB\footnote{${}^\dagger$}{For information on the \CWEB\ programming
language by Donald E.~Knuth, as well as samples of \CWEB\ programs, see
{\tt http://www-cs-faculty.stanford.edu/\~\ \kern -5pt knuth/cweb.html}.
For general information on literate programming, see
{\tt http://www.literateprogramming.com}.}
computer program calculates box-counting estimates of the fractal dimension
of curves in the two-dimensional plane.

In the box-counting estimate to the fractal dimension of a graph in the domain
$\{x,y:x_{\rm min}\le x\le x_{\rm max},y_{\rm min}\le y\le y_{\rm max}\}$, a
grid of squares, each of horizontal dimension $(x_{\rm max}-x_{\rm min})/2^m$
and vertical dimension $(y_{\rm max}-y_{\rm min})/2^m$, is superimposed onto
the graph for integer numbers $m$.
By counting the total number of such squares $N_m$ needed to cover the entire
graph at a given $m$ (hence the term ``box counting''), an estimate $D_m$ to
the fractal dimension $D$ (or Hausdorff dimension) is obtained as
$D_m=\ln(N_m)/\ln(2^m)$.
This procedure may be repeated may times, with $D_m\to D$ as $m\to\infty$
for real fractal sets. However, for finite-depth fractals (as generated by
a computer), some limit on $m$ is necessary in order to prevent trivial
convergence towards $D_m\to1$.

In addition to mere numerical calculation, the program also generates graphs
of the box distributions, in form of \METAPOST\ code which can be
post-processed by other programs.
\bigskip
\noindent Copyright \copyright\ Fredrik Jonsson, 2006--2011.
All rights reserved.
\vfill

\fi

\N{1}{2}The CWEB programming language.
For the reader who might not be familiar with the concept of the
\CWEB\ programming language, the following citations hopefully will
be useful. For further information, as well as freeware compilers for
compiling \CWEB\ source code, see {\tt http://www.literateprogramming.com}.
\bigskip

{\narrower\narrower\narrower\narrower\eightcmssqi\noindent
I believe that the time is ripe for significantly better documentation of
programs, and that we can best achieve this by considering programs to be
works of literature. Hence, my title: `Literate Programming.'

Let us change our traditional attitude to the construction of programs:
Instead of imagining that our main task is to instruct a computer what to
do, let us concentrate rather on explaining to human beings what we want
a computer to do.

The practitioner of literate programming can be regarded as an essayist,
whose main concern is with exposition and excellence of style. Such an
author, with thesaurus in hand, chooses the names of variables carefully
and explains what each variable means. He or she strives for a program
that is comprehensible because its concepts have been introduced in an
order that is best for human understanding, using a mixture of formal and
informal methods that reinforce each other.
\smallskip
{\eightcmssq --\,Donald Knuth,}
{\eightcmssqi The CWEB System of Structured Documentation}
{\eightcmssq (Addison-Wesley, Massachusetts, 1994)}
}
\bigskip

{\narrower\narrower\narrower\narrower\eightcmssqi\noindent
The philosophy behind CWEB is that an experienced system programmer, who
wants to provide the best possible documentation of his or her software
products, needs two things simultaneously: a language like \TeX\ for
formatting, and a language like C for programming. Neither type of language
can provide the best documentation by itself; but when both are appropriately
combined, we obtain a system that is much more useful than either language
separately.

The structure of a software program may be thought of as a `WEB' that is
made up of many interconnected pieces. To document such a program we want to
explain each individual part of the web and how it relates to its neighbors.
The typographic tools provided by \TeX\ give us an opportunity to explain the
local structure of each part by making that structure visible, and the
programming tools provided by languages like C make it possible for us to
specify the algorithms formally and unambiguously. By combining the two,
we can develop a style of programming that maximizes our ability to perceive
the structure of a complex piece of software, and at the same time the
documented programs can be mechanically translated into a working software
system that matches the documentation.

Besides providing a documentation tool, CWEB enhances the C language by
providing the ability to permute pieces of the program text, so that a
large system can be understood entirely in terms of small sections and
their local interrelationships. The CTANGLE program is so named because
it takes a given web and moves the sections from their web structure into
the order required by C; the advantage of programming in CWEB is that the
algorithms can be expressed in ``untangled'' form, with each section
explained separately. The CWEAVE program is so named because it takes a
given web and intertwines the \TeX\ and C portions contained in each
section, then it knits the whole fabric into a structured document.
\smallskip
{\eightcmssq --\,Donald Knuth, ``Literate Programming'', in}
{\eightcmssqi Literate Programming}
{\eightcmssq (CSLI Lecture Notes, Stanford, 1992)}
}

\fi

\N{1}{3}Revision history of the program.
\medskip

\citem[2006-05-08]{[v.1.0]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
First properly working version of the \boxcount\ program, written in \CWEB\
and (\ANSI-conformant) \CEE. I have so far compiled the code with \GCC\ using
the {\tt --pedantic} option and verified that the box-covering routine
\PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}(\,)} and its more low-level
core engine
\PB{\\{box\_intersection}(\,)} both work properly, by direct inspection of the
compiled
\METAPOST\ graphs generated by this routine.
That the numbers of counted boxes are correct has also been verified and the
only remaining blocks to be added are related to the extraction of the fractal
dimension as such.
This first version of the \boxcount\ program consists of 52671 bytes
(1292 lines) of \CWEB\ code, under \CYGWIN\ generating an executable file
of 22561 bytes and 33 pages of program documentation in 10~pt typeface.

\citem[2006-05-09]{[v.1.1]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
Changed the block for the estimate of the fractal dimension, so that the
estimate now is obtained as the average of $\ln(N_{\rm boxes})/\ln(2^n)$
for a set of $n$ such that $N_{\rm min}\le n\le N_{\rm max}$, rather than
performing a linear curve fit to the data.
In order to sample statistical information on the estimate, such as standard
deviation, average deviation and skewness, I incorporated a slightly modified
version of the routine \PB{\\{moment}(\,)} from {\sl Numerical Recipes in C}.
Also added a preliminary section describing the test application of \boxcount\
to the Koch fractal, being a simple test case which is easily implemented by
means of a recursive generation of the data set for the input trajectory.
As of today, the \boxcount\ program consists of 59788 bytes
(1444 lines) of \CWEB\ code, under \CYGWIN\ generating an executable file
of 24253 bytes and 38 pages of program documentation in 10~pt typeface.

\citem[2006-05-14]{[v.1.2]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
Added a section on the command-line options for supplying the \boxcount\
program with input parameters. Polished the section on the example of
estimation of the box-counting dimension of the Koch fractal, and in
particular changed the example from the Koch curve to the Koch snowflake
instead, just for the sake of visual beauty.
As of today, the \boxcount\ program consists of 72560 bytes (1699 lines) of
\CWEB\ code, under \OSX\ generating an executable file of 24996 bytes and 41
pages of program documentation in 10~pt typeface.

\citem[2006-05-17]{[v.1.3]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
Added documentation on the \PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}(%
\,)} routine
and generally cleaned up in the documentation of the program.
As of today, the \boxcount\ program consists of 82251 bytes
(1844 lines) of \CWEB\ code, under \CYGWIN\ generating an executable file
of 29152 bytes and 40 pages of program documentation in 10~pt typeface.

\citem[2006-06-17]{[v.1.4]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
Added the {\tt --graphsize} option, in order to override the default graph
size. Also changed the inclusion of the input trajectory in the boxmaps, so
that the \boxcount\ program rather than using a \MP\ call of the form
{\tt gdraw "input.dat"} now chops the trajectory into proper pieces and
includes the entire trajectory explicitly in the generated graph.
This way the output \MP\ code naturally increases considerably in size, but
is now at least self-sustaining even is separated from the original data file
containing the input trajectory.

\citem[2006-10-25]{[v.1.5]} {\tt <fj@phys.soton.ac.uk>}\hfill\break
Added two pages of text on the boxcounting estimate of the fractal dimension
of the Koch snowflake fractal in the example section.

\citem[2011-11-26]{[v.1.6]} {\tt <http://jonsson.eu>}\hfill\break
Updated \.{Makefile}:s for the generation of figures. Also corrected a rather
stupid way of removing preceeding paths of file names.

\fi

\N{1}{4}Compiling the source code. The program is written in \CWEB, generating
\ANSICEE\ (ISO C90) conforming source code and documentation as plain
\TeX-source, and is to be compiled using the sequences as outlined in the
{\tt Makefile} listed below.
\bigskip
{\obeylines\obeyspaces\tt
\#
\# Makefile designed for use with ctangle, cweave, gcc, and plain TeX.
\#
\# Copyright (C) 2002-2011, Fredrik Jonsson <http://jonsson.eu>
\#
\# The CTANGLE program converts a CWEB source document into a C program which
\# may be compiled in the usual way. The output file includes \#line specifica-
\# tions so that debugging can be done in terms of the CWEB source file.
\#
\# The CWEAVE program converts the same CWEB file into a TeX file that may be
\# formatted and printed in the usual way. It takes appropriate care of typo-
\# graphic details like page layout and the use of indentation, italics,
\# boldface, etc., and it supplies extensive cross-index information that it
\# gathers automatically.
\#
\# CWEB allows you to prepare a single document containing all the information
\# that is needed both to produce a compilable C program and to produce a well-
\# formatted document describing the program in as much detail as the writer
\# may desire.  The user of CWEB ought to be familiar with TeX as well as C.
\#
PROJECT  = boxcount
CTANGLE  = ctangle
CWEAVE   = cweave
CC       = gcc
CCOPTS   = -O2 -Wall -ansi -std=iso9899:1990 -pedantic
LNOPTS   = -lm
TEX      = tex
DVIPS    = dvips
DVIPSOPT = -ta4 -D1200
METAPOST = mp
PS2PDF   = ps2pdf
~   ~
all: \$(PROJECT) \$(PROJECT).pdf
~        @echo
"==============================================================="
~        @echo " To verify the execution performance of the BOXCOUNT program"
~        @echo " on the Koch snowflake fractal, run 'make verification'."
~        @echo
"==============================================================="
~   ~
\$(PROJECT): \$(PROJECT).o
~        \$(CC) \$(CCOPTS) -o \$(PROJECT) \$(PROJECT).o \$(LNOPTS)
~   ~
\$(PROJECT).o: \$(PROJECT).w
~        \$(CTANGLE) \$(PROJECT)
~        \$(CC) \$(CCOPTS) -c \$(PROJECT).c
~   ~
\$(PROJECT).pdf: \$(PROJECT).ps
~        \$(PS2PDF) \$(PROJECT).ps \$(PROJECT).pdf
~   ~
\$(PROJECT).ps: \$(PROJECT).dvi
~        \$(DVIPS) \$(DVIPSOPT) \$(PROJECT).dvi -o \$(PROJECT).ps
~   ~
\$(PROJECT).dvi: \$(PROJECT).w
~        @make -C figures/
~        @make -C kochxmpl/
~        @make verification
~        \$(CWEAVE) \$(PROJECT)
~        \$(TEX) \$(PROJECT).tex
~   ~
verification:
~        @echo
"==============================================================="
~        @echo " Verifying the performance of the \$(PROJECT) program on the
Koch"
~        @echo " snowflake fractal of iteration order 11."
~        @echo
"==============================================================="
~        make koch -C koch/
~        make kochgraphs -C koch/
~        make fractaldimension -C koch/
~   ~
tidy:
~        make -ik clean -C figures/
~        make -ik clean -C koch/
~        make -ik clean -C kochxmpl/
~        -rm -Rf *~ *.o *.exe *.dat *.tgz *.trj *.aux *.log *.idx *.scn *.dvi
~   ~
clean:
~        make -ik tidy
~        -rm -Rf \$(PROJECT) *.c *.pdf *mp *.toc *.tex *.ps
~   ~
archive:
~        make -ik tidy
~        tar --gzip --directory=../ -cf ../\$(PROJECT).tgz \$(PROJECT)
}
\bigskip
This {\tt Makefile} essentially executes two major calls. First, the \CTANGLE\
program parses the \CWEB\ source document {\tt boxcount.w} to extract a \CEE\
source file {\tt boxcount.c} which may be compiled in the usual way using any
\ANSICEE\ conformant compiler. The output source file includes {\tt \#line}
specifications so that any debugging can be done conveniently in terms of
the original \CWEB\ source file.
Second, the \CWEAVE\ program parses the same \CWEB\ source file to extract a
plain \TeX\ source file {\tt boxcount.tex} which may be compiled in the usual
way.
It takes appropriate care of typographic details like page layout and the use
of indentation, italics, boldface, and so on, and it supplies extensive
cross-index information that it gathers automatically.

After having executed {\tt make} in the same catalogue where the files
{\tt boxcount.w} and {\tt Makefile} are located, one is left with an
executable file {\tt boxcount}, being the ready-to-use compiled program,
and a PostScript file {\tt boxcount.ps}
which contains the full documentation of the program, that is to say the
document you currently are reading.
Notice that on platforms running Windows NT, Windows 2000, Windows ME, or any
other operating system by Microsoft, the executable file will instead
automatically be called {\tt boxcount.exe}.
This convention also applies to programs compiled under the \UNIX-like
environment \CYGWIN.

\fi

\N{1}{5}Running the program. The program is entirely controlled by the command
line options supplied when invoking the program. The syntax for executing the
program is\par
\medskip
\line{\hskip 40pt{\tt boxcount [options]}\hfill}\par
\medskip
\noindent
where {\tt options} include the following, given in their long as well as
their short forms (with prefixes `{\tt --}' and `{\tt -}', respectively):
\medskip
\optitem[{\tt --trajectoryfile}, {\tt -i} $\langle${\it trajectory filename}%
$\rangle$]{Specifies the input trajectory of the graph whose fractal
dimension is to be estimated. The input file should describe the trajectory
as two columns of $x$- and $y$-coordinates of the nodes, between which
straight lines will interpolate the trajectory. Unless the boundary of the
computational window is explicitly stated using the {\tt -w} or
{\tt --computationwindow} options, the minimum and maximum $x$- and
$y$-values of the specified trajectory will be used for the automatic
internal computation of the proper computational domain boundaries.}
\optitem[{\tt --outputfile}, {\tt -o [append{\char'174}overwrite]}
$\langle${\it output filename}$\rangle$]{Specifies the base name of the file
to which the calculated output data will be written.
If the {\tt --outputfile} or {\tt -o} option is followed by ``{\tt append}''
the estimate for the fractal dimension will be appended to the file named
$\langle${\it output filename}$\rangle$.dat, which will be created if it
does not exist. If the following word instead is ``{\tt overwrite}'' the
file will, of course, instead be overwritten.}
\optitem[{\tt -w}, {\tt --computationwindow llx} $\langle x_{\rm LL}\rangle$
{\tt lly} $\langle y_{\rm LL}\rangle$ {\tt urx} $\langle x_{\rm UR}\rangle$
{\tt ury} $\langle y_{\rm UR}\rangle$]{This option explicitly specifies the
domain over which the box-counting algorithm will be applied.
By specifying this option, any automatic calculation of computational window
will be neglected.}
\optitem[{\tt --verbose}, {\tt -v}]{Use this option to toggle verbose mode of
the program execution, controlling the amount of information written to
standard terminal output. (Default is off, that is to say quiet mode.)}
\optitem[{\tt --boxmaps}, {\tt -m}]{If this option is present, the \boxcount\
program will generate \METAPOST\ code for maps of the distribution of boxes,
so-called ``boxmaps''. In doing so, also the input trajectory is included in
the graphs. The convention for the naming of the output map files is that
they are saved as $\langle${\it output filename}$\rangle$.$N$.{\tt dat},
where $\langle${\it output filename}$\rangle$ is the base filename as
specified using the {\tt -o} or {\tt --outputfile} option, $N$ is the
automatically appended current level of resolution refinement in the
box-counting (that is to say, indicating the calculation performed for a
$[2^{N}\times2^{N}]$-grid of boxes), and where {\tt dat} is a file suffix
as automatically appended by the program. This option is useful for tracking
the performance of the program, and for verification of the box counting
algorithm.}
\optitem[{\tt --graphsize} $\langle w\rangle$ $\langle h\rangle$]{If the
{\tt -m} or {\tt --boxmaps} option is present at the command line, then
the {\tt --graphsize} option specifies the physical width $w$ and height
$h$ in millimetres of the generated graphs, overriding the default sizes
$w=80\,{\rm mm}$ and $h=80\,{\rm mm}$.}
\optitem[{\tt --minlevel}, {\tt -Nmin} $\langle N_{\rm min}\rangle$]{Specifies
the minimum level $N_{\rm min}$ of grid refinement that will be used in the
evaluation of the estimate of the fractal dimension.
With this option specified, the coarsest level used in the box-counting will
be a $[2^{N_{\rm min}}\times2^{N_{\rm min}}]$-grid of boxes.}
\optitem[{\tt --maxlevel}, {\tt -Nmax} $\langle N_{\rm max}\rangle$]{This
option is similar to the {\tt --minlevel} or {\tt -Nmin} options, with the
difference that it instead specifies the maximum level $N_{\rm max}$ of grid
refinement that will be used in the evaluation of the estimate of the fractal
dimension.
With this option specified, the finest level used in the box-counting will
be a $2^{N_{\rm max}}\times2^{N_{\rm max}}$-grid of boxes.
If this option is omitted, a default value of $N_{\rm max}=10$ will be
used as default.}


\fi

\N{1}{6}Application example: The Koch fractal.
The Koch curve is one of the earliest described fractal curves, appearing
in the 1904 article {\sl Sur une courbe continue sans tangente, obtenue par
une construction g\'eom\'etrique \'el\'ementaire}, written by the Swedish
mathematician Helge von Koch (1870--1924).
In this application example, we employ the ``snowflake'' variant of the Koch
fractal (merely for the sake of its beauty).
The Koch snowflake fractal is constructed recursively using the following
algorithm.
\medskip
\noindent{\bf Algorithm A} ({\it Construction of the Koch snowflake fractal}).
\aitem[{\bf A1.}]{[Create initiator.] Draw three line segments of equal length
so that they form an equilateral triangle.}
\aitem[{\bf A2.}]{[Line division.] Divide each of the line segments into three
segments of equal length.}
\aitem[{\bf A3.}]{[Replace mid segment by triangle.] Draw an equilateral
triangle that has the middle segment from step
one as its base.}
\aitem[{\bf A4.}]{[Remove base of triangle.] Remove the line segment that is
the base of the triangle from step A3.}
\aitem[{\bf A5.}]{[Recursion step.] For each of the line segments remaining,
repeat steps A2 through A4.}
\quad\endalg
\medskip

The triangle resulting after step A1 is called the {\sl initiator} of the
fractal.
After the first iteration of steps A1--A4, the result should be a shape similar
to the Star of David; this is called the {generator} of the fractal.
The Koch fractal resulting of the infinite iteration of the algorithm of
construction has an infinite length, since each time the steps above are
performed on each line segment of the figure there are four times as many line
segments, the length of each being one-third the length of the segments in the
previous stage.
Hence, the total length of the perimeter of the fractal increases by $4/3$
at each step, and for an initiator of total length $L$ the total length
of the perimeter at the $n$th step of iteration will be $(4/3)^nL$.
The fractal dimension is hence $D=\ln4/\ln3\approx1.26$, being greater
than the dimension of a line ($D=1$) but less than, for example, Peano's
space-filling curve ($D=2$).
The Koch fractal is an example of a curve which is continuous, but not
differentiable anywhere.
The area of the Koch snowflake is $8/5$ that of the initial triangle,
so an infinite perimeter encloses a finite area.
The stepwise construction of the snowflake fractal is illustrated in Fig.~1.

\bigskip
\centerline{\epsfxsize=54mm\epsfbox{koch/kochgraph-1.1}\hskip2mm
\epsfxsize=54mm\epsfbox{koch/kochgraph-2.1}\hskip2mm
\epsfxsize=54mm\epsfbox{koch/kochgraph-3.1}}
\centerline{\epsfxsize=54mm\epsfbox{koch/kochgraph-4.1}\hskip2mm
\epsfxsize=54mm\epsfbox{koch/kochgraph-5.1}\hskip2mm
\epsfxsize=54mm\epsfbox{koch/kochgraph-6.1}}
\medskip\noindent
{\bf Figure 1.} Construction of the Koch fractal of the ``snowflake'' type,
in this case as inscribed in the unitary circle. For this case, the length
of the initiator is $L=3\sqrt{3}$ while its area is
$A=L^2/(6\sqrt{3})=3\sqrt{3}/2$. For each fractal recursion depth $n$, the
trajectory consists of a total of $3\times4^{(n-1)}$ connected linear segments.

\vfill\eject

The data set for the Koch fractal is straightforward to generate by means of
recursion, as for example by using the following compact program (which in
fact was used for the generation of the data sets for the Koch fractals on
the previous page):
\bigskip
{\obeylines\obeyspaces\tt
/*-------------------------------------------------------------------------
{\char'174} File: koch.c [ANSI-C conforming source code]
{\char'174} Created: May 8, 2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
{\char'174} Last modified: May 8, 2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
{\char'174} Compile with: gcc -O2 -g -Wall -pedantic -ansi koch.c -o koch -lm
{\char'174} Description: The KOCH program creates data sets corresponding to
{\char'174}    the Koch fractal, for the purpose of acting as test objects for
{\char'174}    the BOXCOUNT program. The KOCH program is simply executed by
{\char'174}    'koch <N>', where <N> is an integer describing the maximum
{\char'174}    depth of recursion in the generation of the fractal data set.
{\char'174}    If invoked without any arguments, <N>=6 is used as default.
{\char'174}    The generated trajectory is written to standard output.
{\char'174} Copyright (C) 2006 Fredrik Jonsson <fj@phys.soton.ac.uk>
=========================================================================*/
\#include <math.h>
\#include <stdio.h>
extern char *optarg;
\hskip\lineskip
void kochsegment(double xa,double ya,double xb,double yb,
\   int depth,int maxdepth) {\char'173}
\   double xca,yca,xcb,ycb,xcc,ycc;
\   if (depth==maxdepth) {\char'173}
\      fprintf(stdout,"\%2.8f \%2.8f{\char'134}n",xb,yb);
\   {\char'175} else {\char'173}
\      xca=xa+(xb-xa)/3.0;
\      yca=ya+(yb-ya)/3.0;
\      xcb=xb-(xb-xa)/3.0;
\      ycb=yb-(yb-ya)/3.0;
\      xcc=(xa+xb)/2.0-(yb-ya)/(2.0*sqrt(3.0));
\      ycc=(ya+yb)/2.0+(xb-xa)/(2.0*sqrt(3.0));
\      kochsegment(xa,ya,xca,yca,depth+1,maxdepth);
\      kochsegment(xca,yca,xcc,ycc,depth+1,maxdepth);
\      kochsegment(xcc,ycc,xcb,ycb,depth+1,maxdepth);
\      kochsegment(xcb,ycb,xb,yb,depth+1,maxdepth);
\   {\char'175}
{\char'175}
\hskip\lineskip
int main(int argc, char *argv[]) {\char'173}
\   int maxdepth=6;
\   if (argc>1) sscanf(argv[1],"\%d",{\char'046}maxdepth);
\   if (maxdepth>0) {\char'173}
\      fprintf(stdout,"\%2.8f \%2.8f{\char'134}n",0.0,1.0);
\      kochsegment(0.0,1.0,sqrt(3.0)/2.0,-0.5,1,maxdepth);
\      kochsegment(sqrt(3.0)/2.0,-0.5,-sqrt(3.0)/2.0,-0.5,1,maxdepth);
\      kochsegment(-sqrt(3.0)/2.0,-0.5,0.0,1.0,1,maxdepth);
\   {\char'175}
\   return(0);
{\char'175}
}\bigskip
\vfill\eject

The boxcounting dimension of the Koch snowflake fractal can now be investigated
with assistance of the \boxcount\ program. In the analysis as here presented,
this is done using the following Makefile:
\bigskip
{\obeylines\obeyspaces\tt
\#
\# Makefile designed for use with gcc, MetaPost and plain TeX.
\#
\# Copyright (C) 2002-2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
\#
CC       = gcc
CCOPTS   = -O2 -Wall -ansi -std=iso9899:1990 -pedantic
LNOPTS   = -lm
TEX      = tex
DVIPS    = dvips
METAPOST = mpost

\#
\# Define path and executable file for the BOXCOUNT program, used for
\# calculating estimates of the box-counting fractal dimension of a
\# trajectory in the plane.
\#
BOXCOUNTPATH = ../
BOXCOUNT = \$(BOXCOUNTPATH)/boxcount

\#
\# Define path and executable file for the KOCH program, used for generating
\# the trajectory of the Koch snowflake fractal.
\#
KOCHPATH = ../koch/
KOCH = \$(KOCHPATH)/koch

all: koch kochgen kochtab

koch:
\        make -C ../koch/

kochgen:
\        \$(KOCH) 7 > koch.trj
\        \$(BOXCOUNT) --verbose --boxmaps --graphsize 42.0 42.0 {\char'134}
\                --computationwindow llx -1.1 lly -1.1 urx 1.1 ury 1.1 {%
\char'134}
\                --minlevel 3 --maxlevel 8 {\char'134}
\                --trajectoryfile koch.trj --outputfile overwrite koch
\        for k in 03 04 05 06 07 08; do {\char'134}
\        \$(METAPOST) koch.\$\$k.mp ;{\char'134}
\        \$(TEX) -jobname=koch.\$\$k '{\char'134}input epsf{%
\char'134}nopagenumbers{\char'134}
\           {\char'134}centerline{{\char'134}epsfxsize=120mm{\char'134}%
epsfbox{koch.'\$\$k'.1}}{\char'134}bye';{\char'134}
\        \$(DVIPS) -D1200 -E koch.\$\$k -o koch.\$\$k.eps;{\char'134}
\        done
\hskip\lineskip
kochtab:
\        \$(KOCH) 7 > koch.trj
\        \$(BOXCOUNT) --verbose --minlevel 3 --maxlevel 14 {\char'134}
\                --computationwindow llx -1.1 lly -1.1 urx 1.1 ury 1.1 {%
\char'134}
\                --trajectoryfile koch.trj --outputfile overwrite koch
\hskip\lineskip
clean:
\        -rm -Rf koch *~ *.o *.exe *.dat *.mp *.mpx *.trj
\        -rm -Rf *.tex *.aux *.log *.toc *.idx *.scn *.dvi *.1 *.eps
}\bigskip
\vfill\eject

Having executed the {\tt Makefile} as displayed in the previous page, where a
recursion depth of $n=7$ is used for the generation of the Koch fractal, we are
left with a set of images of the consecutively refined grids in the boxcounting
algorithm, and a table containing the estimates of the boxcounting dimension
of the Koch snowflake fractal.
In Fig.~2 the resulting maps of the boxes used in the boxcounting algorithm
for refinement levels $m=3,4,\ldots,8$ are shown, and as the grid refinement
is finer and finer, the boxes finally will be barely visible in the limited
resolution of computer graphics, on screen as well as in the generated
Encapsulated PostScript code.

\bigskip
\centerline{\epsfxsize=54mm\epsfbox{kochxmpl/koch-03.1}\hskip2mm
\epsfxsize=54mm\epsfbox{kochxmpl/koch-04.1}\hskip2mm
\epsfxsize=54mm\epsfbox{kochxmpl/koch-05.1}}
\centerline{\epsfxsize=54mm\epsfbox{kochxmpl/koch-06.1}\hskip2mm
\epsfxsize=54mm\epsfbox{kochxmpl/koch-07.1}\hskip2mm
\epsfxsize=54mm\epsfbox{kochxmpl/koch-08.1}}
\medskip\noindent
{\bf Figure 2.} Consecutive steps of refinement $m=3,4,\ldots,8$ of the grid
of boxes covering the input trajectory, in this case a Koch snowflake fractal
of recursion depth $n=7$ (see Fig.~1). At each step $m$ of refinement, a
virtual grid of $2^m\times2^m$ boxes is superimposed on the trajectory
(fractal) and the number $N_m$ of boxes covering the trajectory are counted.
For the Koch snowflake fractal of recursion depth $n=7$, the trajectory
consists of a total of $3\times4^{(7-1)}=12288$ connected linear segments.
\vfill\eject

The estimated boxcounting dimension $D_m=\ln(N_m)/\ln(2^m)$, with $N_m$ as
previously denoting the number of boxes in a $2^m\times2^m$-grid required to
cover the entire curve, is displayed in Table~A.1.
The values for the estimates could be compared with the analytically obtained
value of $D=\ln4/\ln3\approx1.26$ for the fractal dimension.
However, it should be emphasized that the box counting dimension here just
is an estimate of one definition of the fractal dimension, which in many cases
do not equal to other measures, and that we in the computation of the estimate
always will arrive at a truncated result due to a limited precision and a
limited amount of memory resources and computational time.
As can be seen in the table, the initial estimates at lower resolutions are
pretty crude, but in the final estimate we nevertheless end up with the box
counting estimate of the fractal dimension as 1.29, which is reasonably close
to the analytically obtained value of $D\approx1.26$.
Of course, further refinements such as Richardson extrapolation could be
applied to increase the accuracy, but this is outside of the scope of the
\boxcount\ program, which only serves to provide the raw, basic algorithm
of boxcounting.\footnote{$\dagger$}{For discussions on
different definitions of the fractal dimension, see
the English Wikipedia section on the Minkowski--Bouligand dimension,
{\tt http://en.wikipedia.org/wiki/Minkowski-Bouligand{\char'137}dimension}.}
\bigskip
\noindent{\bf Table A.1} Boxcounting estimates of the fractal dimension of
the Koch snowflake fractal of recursion order $n=7$. In the table, $m$ is the
refinement order as indicated in the graphs in Fig.~2, $N_m$ is the number of
covering boxes counted at refinement level $m$, and $D_m=\ln(N_m)/\ln(2^m)$
is the estimate of the boxcounting dimension.
\smallskip
\noindent{\hrule width 328pt}\vskip1pt
\beginvrulealign
% \tabskip=0pt
\halign{
\hbox to 72pt{\hfil#\hfil}% first column centered
&\hbox to 128pt{\hfil#\hfil}% second column is centered
&\hbox to 128pt{\hfil#\hfil}\cr% third column is centered
\noalign{{\hrule width 328pt}\vskip 2pt}
$m$ & $N_m$ & $D_m=\ln(N_m)/\ln(2^m)$\cr
\noalign{\vskip 1pt{\hrule width 328pt}\vskip 1.4pt}
3 &     44 & 1.8198\cr
4 &     96 & 1.6462\cr
5 &    196 & 1.5229\cr
6 &    504 & 1.4962\cr
7 &   1180 & 1.4578\cr
8 &   2856 & 1.4350\cr
9 &   6844 & 1.4156\cr
10 &  15620 & 1.3931\cr
11 &  32320 & 1.3618\cr
12 &  66200 & 1.3345\cr
13 & 133600 & 1.3098\cr
14 & 268804 & 1.2883\cr
}\endvrulealign
\vskip-9.55pt
{\hskip-.4pt\hrule width 328pt}\vskip 1pt
\noindent{\hskip-.2pt\hrule width 328pt}


\fi

\N{1}{7}The main program. Here follows the general outline of the main program.

\Y\B\X8:Library inclusions\X\6
\X9:Global definitions\X\6
\X10:Global variables\X\6
\X23:Subroutines\X\7
\&{int} \\{main}(\&{int} \\{argc}${},\39{}$\&{char} ${}{*}\\{argv}[\,]){}$\1\1%
\2\2\6
${}\{{}$\1\6
\X11:Declaration of local variables\X\6
\X12:Initialize variables\X\6
\X13:Parse command line for parameters\X\6
\X14:Display starting time of program execution\X\6
\X15:Load input trajectory from file\X\6
\X16:Open file for output of logarithmic estimate\X\6
\X17:Extract boundary of global window of computation from input trajectory\X\6
\X18:Get number of boxes covering the trajectory for all levels of refinement
in resolution\X\6
\X19:Compute the logarithmic estimate of the fractal dimension\X\6
\X20:Save or append the logarithmic estimate to output file\X\6
\X21:Close file for output of logarithmic estimate\X\6
\X22:Display elapsed execution time\X\6
\&{return} (\.{SUCCESS});\6
\4${}\}{}$\2\par
\fi

\M{8}Library dependencies.
The standard \ANSICEE\ libraries included in this program are:\par
\varitem[{\tt math.h}]{For access to common mathematical functions.}
\varitem[{\tt stdio.h}]{For file access and any block involving \PB{%
\\{fprintf}}.}
\varitem[{\tt stdlib.h}]{For access to the \PB{\\{exit}} function.}
\varitem[{\tt string.h}]{For string manipulation, \PB{\\{strcpy}}, \PB{%
\\{strcmp}} etc.}
\varitem[{\tt ctype.h}]{For access to the \PB{\\{isalnum}} function.}

\Y\B\4\X8:Library inclusions\X${}\E{}$\6
\8\#\&{include} \.{<math.h>}\6
\8\#\&{include} \.{<time.h>}\6
\8\#\&{include} \.{<stdio.h>}\6
\8\#\&{include} \.{<stdlib.h>}\6
\8\#\&{include} \.{<string.h>}\6
\8\#\&{include} \.{<ctype.h>}\par
\U7.\fi

\M{9}Global definitions.
There are just a few global definitions present in the \boxcount\ program:\par
\varitem[{\tt VERSION}]{The current program revision number.}
\varitem[{\tt COPYRIGHT}]{The copyright banner.}
\varitem[{\tt SUCCESS}]{The return code for successful program termination.}
\varitem[{\tt FAILURE}]{The return code for unsuccessful program termination.}
\varitem[{\tt NCHMAX}]{The maximum number of characters allowed in strings for
storing file names, including path.}
\varitem[{\tt APPEND}]{Code for the flag \PB{\\{output\_mode}}, used to
determine if
output should append to an existing file.}
\varitem[{\tt OVERWRITE}]{Code for the flag \PB{\\{output\_mode}}, used to
determine
if output should overwrite an existing file.}

\Y\B\4\X9:Global definitions\X${}\E{}$\6
\8\#\&{define} \.{VERSION} \5\.{"1.6"}\6
\8\#\&{define} \.{COPYRIGHT} \5\.{"Copyright\ (C)\ 2006-}\)\.{2011,\ Fredrik\
Jonsso}\)\.{n"}\6
\8\#\&{define} \.{SUCCESS} \5(\T{0})\6
\8\#\&{define} \.{FAILURE} \5(\T{1})\6
\8\#\&{define} \.{NCHMAX} \5(\T{256})\6
\8\#\&{define} \.{APPEND} \5\T{1}\6
\8\#\&{define} \.{OVERWRITE} \5\T{2}\par
\U7.\fi

\M{10}Declaration of global variables. The only global variables allowed in
my programs are \PB{\\{optarg}}, which is a pointer to the the string of
characters
that specified the call from the command line, and \PB{\\{progname}}, which
simply
is a pointer to the string containing the name of the program, as it was
invoked from the command line.

\Y\B\4\X10:Global variables\X${}\E{}$\6
\&{extern} \&{char} ${}{*}\\{optarg};{}$\6
\&{char} ${}{*}\\{progname}{}$;\par
\U7.\fi

\M{11}Declaration of local variables of the \PB{\\{main}} program.
In \CWEB\ one has the option of adding variables along the program, for example
by locally adding temporary variables related to a given sub-block of code.
However, the philosophy in the \boxcount\ program is to keep all variables of
the \PB{\\{main}} section collected together, so as to simplify tasks as, for
example,
tracking down a given variable type definition.
The local variables of the program are as follows:
\medskip
\varitem[\PB{\\{verbose}}]{Boolean flag which, if being nonzero, tells the
program
to display information at terminal output during program execution.}
\varitem[\PB{\\{user\_set\_compwin}}]{Boolean flag which indicates whether the
user has
explicitly stated a window of computation or not.}
\varitem[\PB{\\{output\_mode}}]{Variable which states whether the estimate of
fractal dimension should be appended to an existing file which is
created if it does not exist (for $\PB{\\{output\_mode}}=\PB{\.{APPEND}}$), or
if
the data should overwrite already existing data in the file (for
$\PB{\\{output\_mode}}=\PB{\.{OVERWRITE}}$).}
\varitem[\PB{\\{make\_boxmaps}}]{If nonzero, then graphs showing the
distribution
of covering boxes will be generated.}
\varitem[\PB{${*}\\{num\_boxes}$}]{Pointer to array holding the number of boxes
at
various depths of division.}
\varitem[\PB{\\{initime}}]{The time at which the \boxcount\ program was
initialized.}
\varitem[\PB{\\{now}}]{Dummy variable for extraction of current time from the
system
clock.}
\varitem[\PB{\\{mm}}]{The number of points $M$ in the input trajectory. This is
the
length of the vectors \PB{\\{x\_traj}} and \PB{\\{y\_traj}} as described
below.}
\varitem[\PB{\\{nn}}]{Gives the number of boxes in the $x$- or $y$-directions
as
$2^N$.}
\varitem[\PB{\\{nnmax}}]{The maximum refinement depth $N_{\rm max}$ of $N$.}
\varitem[\PB{\\{nnmin}}]{The minimum refinement depth $N_{\rm min}$ of $N$.}
\varitem[\PB{\\{global\_llx}}, \PB{\\{global\_lly}}]{Lower-left coordinates of
global window
of computation.}
\varitem[\PB{\\{global\_urx}}, \PB{\\{global\_ury}}]{Upper-right coordinates of
global window
of computation.}
\varitem[\PB{${*}\\{x\_traj}$},\PB{${*}\\{y\_traj}$}]{Vectors containing the
$x$- and $y$-values of
the coordinates along the input trajectory.}
\varitem[\PB{${*}\|x$},\PB{${*}\|y$}]{Variables for keeping the refinement and
number of boxes.
(This needs to be changed as they are easily confused with $x$ and $y$
coordinates of the trajectory.)}
\varitem[\PB{${*}\\{trajectory\_file}$}]{Input file pointer, for reading the
trajectory
whose fractal dimension is to be estimated.}
\varitem[\PB{${*}\\{frac\_estimate\_file}$}]{Output file pointer, for saving
the estimated
fractal dimension of the input trajectory.}
\varitem[\PB{${*}\\{boxmap\_file}$}]{Output file pointer, for saving maps of
the spatial
locations of the boxes covering the trajectory.}
\varitem[\PB{\\{boxgraph\_width}}]{The physical width in millimetres of the
generated
\MP\ boxmap graphs.}
\varitem[\PB{\\{boxgraph\_height}}]{The physical height in millimetres of the
generated
\MP\ boxmap graphs.}
\varitem[\PB{\\{trajectory\_filename}}]{String for keeping the name of the file
containing the input trajectory.}
\varitem[\PB{\\{output\_filename}}]{String for keeping the base name of the set
of
output files.}
\varitem[\PB{\\{frac\_estimate\_filename}}]{String keeping the name of the file
to which the estimate of the fractal dimension will be saved.}
\varitem[\PB{\\{boxmap\_filename}}]{String for keeping the name of the file to
which \METAPOST\ code for maps of the spatial distribution of boxes will be
written.}
\varitem[\PB{\\{no\_arg}}]{Variable for extracting the number of input
arguments
present on the command line as the program is executed.}
\varitem[\PB{\|i}]{Dummy variable used in loops when reading the input
trajectory.}
\varitem[\PB{${*}\\{fracdimen\_estimates}$}]{Vector keeping the values
characterizing
estimated fractal dimension for various orders of $N$.}
\varitem[\PB{$\\{ave},\\{adev},\\{sdev}$}]{The average, average deviation, and
standard
deviation of the estimated fractal dimension for various orders of
refinement $N$, as stored in \PB{\\{fracdimen\_estimates}}.}
\varitem[\PB{$\\{var},\\{skew},\\{curt}$}]{The variance, skewness, and kurtosis
of the
estimated fractal dimension for various orders of refinement $N$, as
stored in \PB{\\{fracdimen\_estimates}}.}
\medskip
\noindent
Generally in this program, the maximum number of characters a file name
string can contain is \PB{\.{NCHMAX}}, as defined in the definitions section
of the program.

\Y\B\4\X11:Declaration of local variables\X${}\E{}$\6
\&{short} \\{verbose}${},{}$ \\{user\_set\_compwin}${},{}$ \\{output%
\_mode}${},{}$ \\{make\_boxmaps};\6
\&{long} \&{int} ${}{*}\\{num\_boxes};{}$\6
\&{time\_t} \\{initime};\6
\&{time\_t} \\{now};\6
\&{long} \\{mm}${},{}$ \\{nn}${},{}$ \\{nnmin}${},{}$ \\{nnmax};\6
\&{double} \\{global\_llx}${},{}$ \\{global\_lly}${},{}$ \\{global\_urx}${},{}$
\\{global\_ury};\6
\&{double} ${}{*}\\{x\_traj},{}$ ${}{*}\\{y\_traj},{}$ ${}{*}\|x,{}$ ${}{*}%
\|y;{}$\6
\&{FILE} ${}{*}\\{trajectory\_file},{}$ ${}{*}\\{frac\_estimate\_file},{}$
${}{*}\\{boxmap\_file};{}$\6
\&{char} \\{trajectory\_filename}[\.{NCHMAX}]${},{}$ \\{output\_filename}[%
\.{NCHMAX}]${},{}$ \\{frac\_estimate\_filename}[\.{NCHMAX}];\6
\&{char} \\{boxmap\_filename}[\.{NCHMAX}];\6
\&{double} \\{boxgraph\_width}${},{}$ \\{boxgraph\_height};\6
\&{int} \\{no\_arg};\6
\&{long} \|i;\6
\&{double} ${}{*}\\{fracdimen\_estimates},{}$ \\{ave}${},{}$ \\{adev}${},{}$ %
\\{sdev}${},{}$ \\{var}${},{}$ \\{skew}${},{}$ \\{curt};\par
\U7.\fi

\M{12}Initialization of variables.

\Y\B\4\X12:Initialize variables\X${}\E{}$\6
$\\{verbose}\K\T{0}{}$;\C{ Verbose mode is off by default }\6
${}\\{user\_set\_compwin}\K\T{0}{}$;\C{ Before the command-line is parsed,
nothing is known of these settings }\6
${}\\{output\_mode}\K\.{OVERWRITE}{}$;\C{ default mode is to overwrite existing
files }\6
${}\\{make\_boxmaps}\K\T{0}{}$;\C{ Default is to not create graphs of the box
distributions }\6
${}\\{nnmin}\K\T{0}{}$;\C{ Default value for $N_{\rm min}$ }\6
${}\\{nnmax}\K\T{10}{}$;\C{ Default value for $N_{\rm max}$ }\6
${}\\{strcpy}(\\{output\_filename},\39\.{"out"}){}$;\C{ Default output file
basename. }\6
${}\\{strcpy}(\\{trajectory\_filename},\39\.{""}){}$;\C{ To indicate that no
filename has been set yet. }\6
${}\\{boxgraph\_width}\K\T{80.0}{}$;\C{ Default graph width in millimetres }\6
${}\\{boxgraph\_height}\K\T{56.0}{}$;\C{ Default graph height in millimetres }\6
${}\\{trajectory\_file}\K\NULL;{}$\6
${}\\{frac\_estimate\_file}\K\NULL;{}$\6
${}\\{boxmap\_file}\K\NULL;{}$\6
${}\\{initime}\K\\{time}(\NULL){}$;\C{ Time of initialization of the program. }%
\par
\U7.\fi

\M{13}Parsing command line options. All input parameters are passed to the
program through command line options and arguments to the program.
The syntax of command line options is listed
whenever the program is invoked without any options, or whenever any of the
{\tt --help} or {\tt -h} options are specified at startup.

\Y\B\4\X13:Parse command line for parameters\X${}\E{}$\6
${}\{{}$\1\6
${}\\{progname}\K\\{strip\_away\_path}(\\{argv}[\T{0}]);{}$\6
${}\\{fprintf}(\\{stdout},\39\.{"This\ is\ \%s\ v.\%s.\ \%s}\)\.{\\n"},\39%
\\{progname},\39\.{VERSION},\39\.{COPYRIGHT});{}$\6
${}\\{no\_arg}\K\\{argc};{}$\6
\&{while} ${}(\MM\\{argc}){}$\5
${}\{{}$\1\6
\&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-o"})\V\R%
\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--outputfile"})){}$\5
${}\{{}$\1\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"append"})\V\R%
\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"a"})){}$\5
${}\{{}$\1\6
${}\\{output\_mode}\K\.{APPEND};{}$\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
\.{"overwrite"})\V\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"o"})){}$\5
${}\{{}$\1\6
${}\\{output\_mode}\K\.{OVERWRITE};{}$\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ '-o'\ o}\)\.{r\
'--outputfile'\ opt}\)\.{ion!"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\MM\\{argc};{}$\6
${}\\{strcpy}(\\{output\_filename},\39\\{argv}[\\{no\_arg}-\\{argc}]);{}$\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-w"})\V%
\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--computationwindow}\)%
\.{"})){}$\5
${}\{{}$\1\6
${}\\{user\_set\_compwin}\K\T{1};{}$\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"llx"})){}$\5
${}\{{}$\1\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
\\{global\_llx})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ parsin}\)\.{g\ lower-left\
x-value}\)\.{.\\n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ comput}\)\.{ation\ window\
option\\}\)\.{n"},\39\\{progname});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Expecting\ 'llx'}\)\.{\\n"},\39%
\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"lly"})){}$\5
${}\{{}$\1\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
\\{global\_lly})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ parsin}\)\.{g\ lower-left\
y-value}\)\.{.\\n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ comput}\)\.{ation\ window\
option\\}\)\.{n"},\39\\{progname});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Expecting\ 'lly'}\)\.{\\n"},\39%
\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"urx"})){}$\5
${}\{{}$\1\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
\\{global\_urx})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ parsin}\)\.{g\ lower-left\
x-value}\)\.{.\\n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ comput}\)\.{ation\ window\
option\\}\)\.{n"},\39\\{progname});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Expecting\ 'urx'}\)\.{\\n"},\39%
\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"ury"})){}$\5
${}\{{}$\1\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
\\{global\_ury})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ parsin}\)\.{g\ lower-left\
y-value}\)\.{.\\n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ comput}\)\.{ation\ window\
option\\}\)\.{n"},\39\\{progname});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Expecting\ 'ury'}\)\.{\\n"},\39%
\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-i"})\V%
\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--trajectoryfile"})){}$\5
${}\{{}$\1\6
${}\MM\\{argc};{}$\6
${}\\{strcpy}(\\{trajectory\_filename},\39\\{argv}[\\{no\_arg}-\\{argc}]);{}$\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-v"})\V%
\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--verbose"})){}$\5
${}\{{}$\1\6
${}\\{verbose}\K(\\{verbose}\?\T{0}:\T{1});{}$\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-h"})\V%
\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--help"})){}$\5
${}\{{}$\1\6
\\{showsomehelp}(\,);\6
\\{exit}(\.{SUCCESS});\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"-m"})\V%
\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"--boxmaps"})){}$\5
${}\{{}$\1\6
${}\\{make\_boxmaps}\K\T{1};{}$\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
\.{"--graphsize"})){}$\5
${}\{{}$\1\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
\\{boxgraph\_width})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ width\ }\)\.{of\
'--graphsize'\ opt}\)\.{ion.\\n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%lf"},\39{\AND}%
\\{boxgraph\_height})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ height}\)\.{\ of\
'--graphsize'\ op}\)\.{tion.\\n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
\.{"--minlevel"})\V\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
\.{"-Nmin"})){}$\5
${}\{{}$\1\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%ld"},\39{\AND}%
\\{nnmin})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ '--min}\)\.{level'\ or\
'-Nmin'\ op}\)\.{tion.\\n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
\.{"--maxlevel"})\V\R\\{strcmp}(\\{argv}[\\{no\_arg}-\\{argc}],\39%
\.{"-Nmax"})){}$\5
${}\{{}$\1\6
${}\MM\\{argc};{}$\6
\&{if} ${}(\R\\{sscanf}(\\{argv}[\\{no\_arg}-\\{argc}],\39\.{"\%ld"},\39{\AND}%
\\{nnmax})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ '--max}\)\.{level'\ or\
'-Nmax'\ op}\)\.{tion.\\n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Specified\ optio}\)\.{n\ '\%s'\ invalid!%
\\n"},\39\\{progname},\39\\{argv}[\\{no\_arg}-\\{argc}]);{}$\6
\\{showsomehelp}(\,);\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U7.\fi

\M{14}Display starting time of program execution.

\Y\B\4\X14:Display starting time of program execution\X${}\E{}$\6
${}\{{}$\1\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Program\ executi}\)\.{on\ started\ %
\%s"},\39\\{progname},\39\\{ctime}({\AND}\\{initime}));{}$\6
\4${}\}{}$\2\par
\U7.\fi

\M{15}Load input trajectory from file. This is the section where the trajectory
to be analyzed is loaded into the memory, and where the number~$M$ of input
coordinates present in the input data is determined by a single call to
the subroutine \PB{\\{num\_coordinate\_pairs}(\,)}.

\Y\B\4\X15:Load input trajectory from file\X${}\E{}$\6
${}\{{}$\1\6
\&{if} ${}(\R\\{strcmp}(\\{trajectory\_filename},\39\.{""})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ No\ input\ trajec}\)\.{tory\ specified!%
\\n"},\39\\{progname});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Please\ use\ the\ }\)%
\.{'--trajectoryfile'\ o}\)\.{ption.\\n"},\39\\{progname});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Use\ '--help'\ op}\)\.{tion\ to\ display%
\ a\ he}\)\.{lp\ message.\\n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\&{if} ${}((\\{trajectory\_file}\K\\{fopen}(\\{trajectory\_filename},\39%
\.{"r"}))\E\NULL){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Could\ not\ open\ }\)\.{\%s\ for\
loading\ traje}\)\.{ctory!\\n"},\39\\{progname},\39\\{trajectory%
\_filename});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\\{mm}\K\\{num\_coordinate\_pairs}(\\{trajectory\_file});{}$\6
${}\\{x\_traj}\K\\{dvector}(\T{1},\39\\{mm});{}$\6
${}\\{y\_traj}\K\\{dvector}(\T{1},\39\\{mm});{}$\6
\&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{mm};{}$ ${}\|i\PP){}$\5
${}\{{}$\1\6
${}\\{fscanf}(\\{trajectory\_file},\39\.{"\%lf"},\39{\AND}\\{x\_traj}[\|i]){}$;%
\C{ scan $x$-coordinate }\6
${}\\{fscanf}(\\{trajectory\_file},\39\.{"\%lf"},\39{\AND}\\{y\_traj}[\|i]){}$;%
\C{ scan $y$-coordinate }\6
\4${}\}{}$\2\6
\\{fclose}(\\{trajectory\_file});\6
\&{if} (\\{verbose})\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Loaded\ \%ld\ coor}\)\.{dinate\ pairs\
from\ fi}\)\.{le\ '\%s'.\\n"},\39\\{progname},\39\\{mm},\39\\{trajectory%
\_filename});{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U7.\fi

\M{16}Opening files for output by the program.

\Y\B\4\X16:Open file for output of logarithmic estimate\X${}\E{}$\6
${}\{{}$\1\6
\&{if} ${}(\R\\{strcmp}(\\{output\_filename},\39\.{""})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ No\ output\ base\ }\)\.{name\ specified!%
\\n"},\39\\{progname});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Please\ use\ the\ }\)\.{'--outputfile'\
optio}\)\.{n.\\n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\\{sprintf}(\\{frac\_estimate\_filename},\39\.{"\%s.dat"},\39\\{output%
\_filename});{}$\6
\&{if} ${}(\\{output\_mode}\E\.{APPEND}){}$\5
${}\{{}$\1\6
\&{if} ${}((\\{frac\_estimate\_file}\K\\{fopen}(\\{frac\_estimate\_filename},%
\39\.{"a"}))\E\NULL){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Could\ not\ open\ }\)\.{\%s\ for\
output!\\n"},\39\\{progname},\39\\{frac\_estimate\_filename});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\\{fseek}(\\{frac\_estimate\_file},\39\T{0\$L},\39\.{SEEK\_END}){}$;\C{ set
file pointer to the end of the file }\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\\{output\_mode}\E\.{OVERWRITE}){}$\5
${}\{{}$\1\6
\&{if} ${}((\\{frac\_estimate\_file}\K\\{fopen}(\\{frac\_estimate\_filename},%
\39\.{"w"}))\E\NULL){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Could\ not\ open\ }\)\.{\%s\ for\
loading\ traje}\)\.{ctory!\\n"},\39\\{progname},\39\\{frac\_estimate%
\_filename});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\\{fseek}(\\{frac\_estimate\_file},\39\T{0\$L},\39\.{SEEK\_SET}){}$;\C{ set
file pointer to the beginning of the file }\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error:\ Output\ m}\)\.{ode\ (output%
\_mode)\ un}\)\.{defined!"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U7.\fi

\M{17}Extract global window of computation. This block will only be executed if
there was no explicit computational window stated at startup of the program
(that is to say, with the {\tt -w} or {\tt --computationwindow} option).
In order to determine the minimum area covering the entire input trajectory,
this section scans sequentially through the trajectory and finds the minimum
and maximum of the $x$- and $y$-coordinates.
These values then form the lower-left and upper-right corner coordinates
$(\PB{\\{global\_llx}},\PB{\\{global\_lly}})$ and $(\PB{\\{global\_urx}},\PB{%
\\{global\_ury}})$.

\Y\B\4\X17:Extract boundary of global window of computation from input
trajectory\X${}\E{}$\6
${}\{{}$\1\6
\&{if} ${}(\R\\{user\_set\_compwin}){}$\5
${}\{{}$\1\6
${}\\{global\_llx}\K\\{x\_traj}[\T{1}];{}$\6
${}\\{global\_lly}\K\\{y\_traj}[\T{1}];{}$\6
${}\\{global\_urx}\K\\{global\_llx};{}$\6
${}\\{global\_ury}\K\\{global\_lly};{}$\6
\&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{mm};{}$ ${}\|i\PP){}$\5
${}\{{}$\1\6
\&{if} ${}(\\{x\_traj}[\|i]>\\{global\_urx}){}$\1\5
${}\\{global\_urx}\K\\{x\_traj}[\|i];{}$\2\6
\&{if} ${}(\\{y\_traj}[\|i]>\\{global\_ury}){}$\1\5
${}\\{global\_ury}\K\\{y\_traj}[\|i];{}$\2\6
\&{if} ${}(\\{x\_traj}[\|i]<\\{global\_llx}){}$\1\5
${}\\{global\_llx}\K\\{x\_traj}[\|i];{}$\2\6
\&{if} ${}(\\{y\_traj}[\|i]<\\{global\_lly}){}$\1\5
${}\\{global\_lly}\K\\{y\_traj}[\|i];{}$\2\6
\4${}\}{}$\2\6
\&{if} (\\{verbose})\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Global\ box-coun}\)\.{ting\ window:%
\\n"},\39\\{progname});{}$\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ \ \ (llx,lly)=(\%2}\)\.{.8f,\%2.8f)%
\\n"},\39\\{progname},\39\\{global\_llx},\39\\{global\_lly});{}$\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ \ \ (urx,ury)=(\%2}\)\.{.8f,\%2.8f)%
\\n"},\39\\{progname},\39\\{global\_urx},\39\\{global\_ury});{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U7.\fi

\M{18}Get the number of boxes covering the trajectory for all levels of
refinement.
This is the main loop where the program iterates over all the refinement levels
and meanwhile gather how many boxes are needed to cover the trajectory as
function of the refinement index $N$ for $N_{\rm min}\le N\le N_{\rm max}$.
In the loop, the program starts with a grid of $[2^{N_{\rm min}}\times
2^{N_{\rm min}}]$ and ends with a grid of $[2^{N_{\rm max}}\times
2^{N_{\rm max}}]$ boxes, at each step of refinement increasing the number of
boxes by a factor of four.
At each level of increasing resolution, the program relies on the routines
\PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}(\,)} or \PB{\\{get\_num%
\_covering\_boxes}(\,)} to
perform the actual box-counting.
The results of the box counting are stored in the vector \PB{\\{num\_boxes}},
which
is used later on in the actual extraction of the estimate of fractal dimension.

\Y\B\4\X18:Get number of boxes covering the trajectory for all levels of
refinement in resolution\X${}\E{}$\6
${}\{{}$\1\6
${}\\{num\_boxes}\K\\{livector}(\T{1},\39\\{nnmax});{}$\6
${}\\{nn}\K\T{1};{}$\6
\&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{nnmin}-\T{1};{}$ ${}\|i\PP){}$\1\5
${}\\{nn}\K\T{2}*\\{nn}{}$;\C{ This leaves \PB{\\{nn}} as $2^{(N_{\rm min}-1)}$
}\2\6
\&{for} ${}(\|i\K\\{nnmin};{}$ ${}\|i\Z\\{nnmax};{}$ ${}\|i\PP){}$\5
${}\{{}$\C{ For all levels in refinement of grid resolution }\1\6
${}\\{nn}\K\T{2}*\\{nn};{}$\6
\&{if} (\\{make\_boxmaps})\5
${}\{{}$\C{ do we wish to generate \METAPOST\ graphs of the box distribution? }%
\1\6
${}\\{sprintf}(\\{boxmap\_filename},\39\.{"\%s-\%02ld.mp"},\39\\{output%
\_filename},\39\|i);{}$\6
${}\\{num\_boxes}[\|i]\K\\{get\_num\_covering\_boxes\_with\_boxmaps}(\\{x%
\_traj},\39\\{y\_traj},\39\\{mm},\39\|i,\39\\{global\_llx},\39\\{global\_lly},%
\39\\{global\_urx},\39\\{global\_ury},\39\\{boxmap\_filename},\39\\{boxgraph%
\_width},\39\\{boxgraph\_height},\39\\{trajectory\_filename});{}$\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\C{ if not, just use the regular number-crunching routine }\1\6
${}\\{num\_boxes}[\|i]\K\\{get\_num\_covering\_boxes}(\\{x\_traj},\39\\{y%
\_traj},\39\\{mm},\39\|i,\39\\{global\_llx},\39\\{global\_lly},\39\\{global%
\_urx},\39\\{global\_ury});{}$\6
\4${}\}{}$\2\6
\&{if} (\\{verbose})\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ N=\%ld\ (\%ldx\%ld-}\)\.{grid\ of\
boxes):\ "},\39\\{progname},\39\|i,\39\\{nn},\39\\{nn});{}$\6
${}\\{fprintf}(\\{stdout},\39\.{"Trajectory\ covered\ }\)\.{by\ \%ld\ boxes%
\\n"},\39\\{num\_boxes}[\|i]);{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U7.\fi

\M{19}Compute the logarithmic estimate of the fractal dimension. Having
completed
the task of actually calculating the number of boxes necessary to cover the
trajectory, for varying box dimensions of width
$(\PB{\\{global\_urx}}-\PB{\\{global\_llx}})/2^n$ and height $(\PB{\\{global%
\_ury}}-\PB{\\{global\_lly}})/2^n$,
for $n=1,2,\ldots,N$, the only remaining arithmetics is to actually calculate
the estimate for the fractal dimension.
This is done by performing the fit of the linear function $y=ax+b$ to the data
set obtained with $2\ln(n)$ as abscissa and $\ln(\PB{\\{num\_boxes}}(n))$ as
ordinata.

Also in this block, we analyze whether the option \PB{\\{make\_boxmaps}} is set
as true
(1) or not (0), determining whether a graph of the number of boxes as function
of division depth should be written to file as well.

In the fitting of the linear function $y=ax+b$ to the data set $(x_1,y_1)$,
$(x_2,y_2)$,$\ldots$,$(x_N,y_N)$, the parameters $a$ and $b$ can be explicitly
obtained from the summation formulas
$$
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),
\qquad
b={{1}\over{N}}\Big(\sum^N_{k=1}y_k-a\sum^N_{k=1}x_k\Big),
$$
where
$$
D\equiv N\sum^N_{k=1}x^2_k-\Big(\sum^N_{k=1}x_k\Big)^2.
$$

\Y\B\4\X19:Compute the logarithmic estimate of the fractal dimension\X${}\E{}$\6
${}\{{}$\1\6
${}\|x\K\\{dvector}(\\{nnmin},\39\\{nnmax});{}$\6
${}\|y\K\\{dvector}(\\{nnmin},\39\\{nnmax});{}$\6
${}\\{fracdimen\_estimates}\K\\{dvector}(\\{nnmin},\39\\{nnmax});{}$\6
${}\\{nn}\K\T{1};{}$\6
\&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{nnmax};{}$ ${}\|i\PP){}$\5
${}\{{}$\1\6
${}\\{nn}\K\T{2}*\\{nn};{}$\6
\&{if} ${}(\\{nnmin}\Z\|i){}$\1\5
${}\|x[\|i]\K\\{log}{}$((\&{double}) \\{nn});\2\6
\4${}\}{}$\2\6
\&{for} ${}(\|i\K\\{nnmin};{}$ ${}\|i\Z\\{nnmax};{}$ ${}\|i\PP){}$\1\5
${}\|y[\|i]\K\\{log}(\\{num\_boxes}[\|i]);{}$\2\6
\&{for} ${}(\|i\K\\{nnmin};{}$ ${}\|i\Z\\{nnmax};{}$ ${}\|i\PP){}$\1\5
${}\\{fracdimen\_estimates}[\|i]\K\|y[\|i]/\|x[\|i];{}$\2\6
\&{if} (\\{verbose})\5
${}\{{}$\1\6
\&{for} ${}(\|i\K\\{nnmin};{}$ ${}\|i\Z\\{nnmax};{}$ ${}\|i\PP){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ N=\%ld,\ fractal\ }\)\.{dimension\
estimate=\%}\)\.{f\\n"},\39\\{progname},\39\|i,\39\\{fracdimen\_estimates}[%
\|i]);{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
${}\\{moment}(\\{fracdimen\_estimates},\39\\{nnmin},\39\\{nnmax},\39{\AND}%
\\{ave},\39{\AND}\\{adev},\39{\AND}\\{sdev},\39{\AND}\\{var},\39{\AND}\\{skew},%
\39{\AND}\\{curt});{}$\6
\&{if} (\\{verbose})\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Estimate\ of\ fra}\)\.{ctal\ dimension:\
\%f\\n}\)\.{"},\39\\{progname},\39\\{ave});{}$\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Average\ deviati}\)\.{on:\ \%f\\n"},\39%
\\{progname},\39\\{adev});{}$\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Standard\ deviat}\)\.{ion:\ \%f\\n"},\39%
\\{progname},\39\\{sdev});{}$\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Skewness:\ \%f\\n"},\39\\{progname},\39%
\\{skew});{}$\6
\4${}\}{}$\2\6
${}\\{free\_livector}(\\{num\_boxes},\39\T{1},\39\\{nnmax}){}$;\C{ release the
memory occupied by \PB{\\{num\_boxes}} }\6
${}\\{free\_dvector}(\\{fracdimen\_estimates},\39\\{nnmin},\39\\{nnmax});{}$\6
${}\\{free\_dvector}(\|x,\39\\{nnmin},\39\\{nnmax}){}$;\C{ release the memory
occupied by \PB{\|x} }\6
${}\\{free\_dvector}(\|y,\39\\{nnmin},\39\\{nnmax}){}$;\C{ release the memory
occupied by \PB{\|y} }\6
\4${}\}{}$\2\par
\U7.\fi

\M{20}Save the fractal dimension to file.

\Y\B\4\X20:Save or append the logarithmic estimate to output file\X${}\E{}$\6
${}\{{}$\1\6
\&{if} ${}(\\{output\_mode}\E\.{APPEND}){}$\5
${}\{{}$\1\6
${}\\{fseek}(\\{frac\_estimate\_file},\39\T{0\$L},\39\.{SEEK\_END});{}$\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\\{output\_mode}\E\.{OVERWRITE}){}$\5
${}\{{}$\1\6
${}\\{fseek}(\\{frac\_estimate\_file},\39\T{0\$L},\39\.{SEEK\_SET});{}$\6
\4${}\}{}$\2\6
${}\\{fprintf}(\\{frac\_estimate\_file},\39\.{"\%f\ \%f\ \%f\\n"},\39\\{ave},%
\39\\{sdev},\39\\{skew});{}$\6
\4${}\}{}$\2\par
\U7.\fi

\M{21}Close any open output files.

\Y\B\4\X21:Close file for output of logarithmic estimate\X${}\E{}$\6
${}\{{}$\1\6
\\{fclose}(\\{frac\_estimate\_file});\6
\4${}\}{}$\2\par
\U7.\fi

\M{22}Display elapsed execution time.

\Y\B\4\X22:Display elapsed execution time\X${}\E{}$\6
${}\{{}$\1\6
${}\\{now}\K\\{time}(\NULL);{}$\6
\&{if} (\\{verbose})\1\5
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Total\ execution}\)\.{\ time:\ \%d\ s%
\\n"},\39\\{progname},\39{}$((\&{int}) \\{difftime}${}(\\{now},\39%
\\{initime})));{}$\2\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Program\ executi}\)\.{on\ closed\ \%s"},%
\39\\{progname},\39\\{ctime}({\AND}\\{now}));{}$\6
\4${}\}{}$\2\par
\U7.\fi

\N{1}{23}Subroutines. In this section, all subroutines as used by the main
program
are listed.

\Y\B\4\X23:Subroutines\X${}\E{}$\6
\X24:Routine for computation of average, average deviation and standard
deviation\X\6
\X25:Routine for obtaining the number of coordinate pairs in a file\X\6
\X26:Routines for removing preceding path of filenames\X\6
\X29:Routines for memory allocation of vectors and matrices\X\6
\X38:Routine for displaying help message\X\6
\X39:Routine for determining whether two lines intersect or not\X\6
\X40:Routine for determining whether a line and a box intersect or not\X\6
\X41:Routines for calculation of number of boxes covering the trajectory\X\par
\U7.\fi

\M{24}Routine for computation of average, average deviation and standard
deviation.
This routine is adopted from {\sl Numerical Recipes in C}.

\Y\B\4\X24:Routine for computation of average, average deviation and standard
deviation\X${}\E{}$\6
\&{void} \\{moment}(\&{double} \\{data}[\,]${},\39{}$\&{int} \\{nnmin}${},%
\39{}$\&{int} \\{nnmax}${},\39{}$\&{double} ${}{*}\\{ave},\39{}$\&{double}
${}{*}\\{adev},\39{}$\&{double} ${}{*}\\{sdev},\39{}$\&{double} ${}{*}\\{var},%
\39{}$\&{double} ${}{*}\\{skew},\39{}$\&{double} ${}{*}\\{curt}){}$\1\1\2\2\6
${}\{{}$\1\6
\&{int} \|j;\6
\&{double} \\{ep}${}\K\T{0.0},{}$ \|s${},{}$ \|p;\7
\&{if} ${}(\\{nnmax}-\\{nnmin}\Z\T{1}){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ routin}\)\.{e\ moment()!\ "},%
\39\\{progname});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"(nnmax-nnmin>1\ is\ a}\)\.{\ requirement)%
\\n"});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\|s\K\T{0.0};{}$\6
\&{for} ${}(\|j\K\\{nnmin};{}$ ${}\|j\Z\\{nnmax};{}$ ${}\|j\PP){}$\1\5
${}\|s\MRL{+{\K}}\\{data}[\|j];{}$\2\6
${}{*}\\{ave}\K\|s/(\\{nnmax}-\\{nnmin}+\T{1});{}$\6
${}{*}\\{adev}\K({*}\\{var})\K({*}\\{skew})\K({*}\\{curt})\K\T{0.0};{}$\6
\&{for} ${}(\|j\K\\{nnmin};{}$ ${}\|j\Z\\{nnmax};{}$ ${}\|j\PP){}$\5
${}\{{}$\1\6
${}{*}\\{adev}\MRL{+{\K}}\\{fabs}(\|s\K\\{data}[\|j]-({*}\\{ave}));{}$\6
${}\\{ep}\MRL{+{\K}}\|s;{}$\6
${}\\{ep}\MRL{+{\K}}\|s;{}$\6
${}{*}\\{var}\MRL{+{\K}}(\|p\K\|s*\|s);{}$\6
${}{*}\\{skew}\MRL{+{\K}}(\|p\MRL{*{\K}}\|s);{}$\6
${}{*}\\{curt}\MRL{+{\K}}(\|p\MRL{*{\K}}\|s);{}$\6
\4${}\}{}$\2\6
${}{*}\\{adev}\MRL{{/}{\K}}(\\{nnmax}-\\{nnmin}+\T{1});{}$\6
${}{*}\\{var}\K({*}\\{var}-\\{ep}*\\{ep}/(\\{nnmax}-\\{nnmin}+\T{1}))/(%
\\{nnmax}-\\{nnmin});{}$\6
${}{*}\\{sdev}\K\\{sqrt}({*}\\{var});{}$\6
\&{if} ${}({*}\\{var}){}$\5
${}\{{}$\1\6
${}{*}\\{skew}\MRL{{/}{\K}}((\\{nnmax}-\\{nnmin}+\T{1})*({*}\\{var})*({*}%
\\{sdev}));{}$\6
${}{*}\\{curt}\K({*}\\{curt})/((\\{nnmax}-\\{nnmin}+\T{1})*({*}\\{var})*({*}%
\\{var}))-\T{3.0};{}$\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error\ in\ routin}\)\.{e\ moment()!\ "},%
\39\\{progname});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"No\ skew/kurtosis\ fo}\)\.{r\ zero\ variance%
\\n"});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U23.\fi

\M{25}Routine for obtaining the number of coordinate pairs in a file. This
routine is called prior to loading the trajectory, in order to get the
size needed for allocating the memory for the trajectory vector. As the
number of coordinates~$M$ has been established, two vectors \PB{\\{x\_traj}[%
\T{1..}\|M]}
and \PB{\\{y\_traj}[\T{1..}\|M]}, containing the coordinates $(x_m,y_m)$ of the
trajectory.

\Y\B\4\X25:Routine for obtaining the number of coordinate pairs in a file\X${}%
\E{}$\6
\&{long} \&{int} \\{num\_coordinate\_pairs}(\&{FILE} ${}{*}\\{trajectory%
\_file}){}$\1\1\2\2\6
${}\{{}$\1\6
\&{double} \\{tmp};\6
\&{int} \\{tmpch};\6
\&{long} \&{int} \\{mm}${}\K\T{0};{}$\7
${}\\{fseek}(\\{trajectory\_file},\39\T{0\$L},\39\.{SEEK\_SET}){}$;\C{ rewind
file to beginning }\6
\&{while} ${}((\\{tmpch}\K\\{getc}(\\{trajectory\_file}))\I\.{EOF}){}$\5
${}\{{}$\1\6
${}\\{ungetc}(\\{tmpch},\39\\{trajectory\_file});{}$\6
${}\\{fscanf}(\\{trajectory\_file},\39\.{"\%lf"},\39{\AND}\\{tmp}){}$;\C{ Read
away the $x$ coordinate }\6
${}\\{fscanf}(\\{trajectory\_file},\39\.{"\%lf"},\39{\AND}\\{tmp}){}$;\C{ Read
away the $y$ coordinate }\6
${}\\{mm}\PP;{}$\6
${}\\{tmpch}\K\\{getc}(\\{trajectory\_file}){}$;\C{ Read away blanks and
linefeeds }\6
\&{while} ${}((\\{tmpch}\I\.{EOF})\W(\R\\{isdigit}(\\{tmpch}))){}$\1\5
${}\\{tmpch}\K\\{getc}(\\{trajectory\_file});{}$\2\6
\&{if} ${}(\\{tmpch}\I\.{EOF}){}$\1\5
${}\\{ungetc}(\\{tmpch},\39\\{trajectory\_file});{}$\2\6
\4${}\}{}$\2\6
${}\\{fseek}(\\{trajectory\_file},\39\T{0\$L},\39\.{SEEK\_SET}){}$;\C{ rewind
file to beginning }\6
\&{return} (\\{mm});\6
\4${}\}{}$\2\par
\U23.\fi

\M{26}Routines for removing preceding path of filenames.
In this block all routines related to removing preceding path strings go.
Not really fancy programming, and no contribution to any increase of numerical
efficiency or precision; just for the sake of keeping a tidy terminal output
of the program. The \PB{\\{strip\_away\_path}(\,)} routine is typically called
when
initializing the program name string \PB{\\{progname}} from the command line
string
\PB{\\{argv}[\T{0}]}, and is typically located in the blocks related to parsing
of the
command line options.

\Y\B\4\X26:Routines for removing preceding path of filenames\X${}\E{}$\6
\X27:Routine for checking valid path characters\X\6
\X28:Routine for stripping away path string\X\par
\U23.\fi

\M{27}Checking for a valid path character.
The \PB{\\{pathcharacter}} routine takes one character \PB{\\{ch}} as argument,
and returns
1 (``true'') if the character is valid character of a path string, otherwise 0
(``false'') is returned.

\Y\B\4\X27:Routine for checking valid path characters\X${}\E{}$\6
\&{short} \\{pathcharacter}(\&{int} \\{ch})\1\1\2\2\6
${}\{{}$\1\6
\&{return} ${}(\\{isalnum}(\\{ch})\V(\\{ch}\E\.{'.'})\V(\\{ch}\E\.{'/'})\V(%
\\{ch}\E\.{'\\\\'})\V(\\{ch}\E\.{'\_'})\V(\\{ch}\E\.{'-'})\V(\\{ch}\E%
\.{'+'}));{}$\6
\4${}\}{}$\2\par
\U26.\fi

\M{28}Routine for stripping away path string of a file name.
The \PB{\\{strip\_away\_path}(\,)} routine takes a character string \PB{%
\\{filename}} as
argument, and returns a pointer to the same string but without any
preceding path segments. This routine is, for example, useful for
removing paths from program names as parsed from the command line.

\Y\B\4\X28:Routine for stripping away path string\X${}\E{}$\6
\&{char} ${}{*}{}$\\{strip\_away\_path}(\&{char} \\{filename}[\,])\1\1\2\2\6
${}\{{}$\1\6
\&{int} \|j${},{}$ \|k${}\K\T{0};{}$\7
\&{while} (\\{pathcharacter}(\\{filename}[\|k]))\1\5
${}\|k\PP;{}$\2\6
${}\|j\K(\MM\|k){}$;\C{ this is the uppermost index of the full path+file
string }\6
\&{while} (\\{isalnum}((\&{int})(\\{filename}[\|j])))\1\5
${}\|j\MM;{}$\2\6
${}\|j\PP{}$;\C{ this is the lowermost index of the stripped file name }\6
\&{return} ${}({\AND}\\{filename}[\|j]);{}$\6
\4${}\}{}$\2\par
\U26.\fi

\M{29}Subroutines for memory allocation. Here follows the routines for memory
allocation and deallocation of double precision real and complex valued
vectors, as used for storing the optical field distribution in the grating,
the refractive index distribution of the grating, etc.

\Y\B\4\X29:Routines for memory allocation of vectors and matrices\X${}\E{}$\6
\X30:Routine for allocation of vectors of double precision\X\6
\X33:Routine for deallocation of vectors of integer precision\X\6
\X32:Routine for allocation of vectors of integer precision\X\6
\X31:Routine for deallocation of vectors of double precision\X\6
\X36:Routine for allocation of matrices of short integer precision\X\6
\X37:Routine for deallocation of matrices of short integer precision\X\par
\U23.\fi

\M{30}The \PB{\\{dvector}} routine allocates a real-valued vector of double
precision,
with vector index ranging from \PB{\\{nl}} to \PB{\\{nh}}.

\Y\B\4\X30:Routine for allocation of vectors of double precision\X${}\E{}$\6
\&{double} ${}{*}{}$\\{dvector}(\&{long} \\{nl}${},\39{}$\&{long} \\{nh})\1\1\2%
\2\6
${}\{{}$\1\6
\&{double} ${}{*}\|v;{}$\7
${}\|v\K{}$(\&{double} ${}{*}){}$ \\{malloc}${}((\&{size\_t})((\\{nh}-\\{nl}+%
\T{2})*\&{sizeof}(\&{double})));{}$\6
\&{if} ${}(\R\|v){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"Error:\ Allocation\ f}\)\.{ailure\ in\
dvector()\\}\)\.{n"});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\&{return} \|v${}-\\{nl}+\T{1};{}$\6
\4${}\}{}$\2\par
\U29.\fi

\M{31}The \PB{\\{free\_dvector}} routine release the memory occupied by the
real-valued vector \PB{$\|v[\\{nl}\MRL{{.}{.}}\\{nh}]$}.

\Y\B\4\X31:Routine for deallocation of vectors of double precision\X${}\E{}$\6
\&{void} \\{free\_dvector}(\&{double} ${}{*}\|v,\39{}$\&{long} \\{nl}${},\39{}$%
\&{long} \\{nh})\1\1\2\2\6
${}\{{}$\1\6
\\{free}((\&{char} ${}{*})(\|v+\\{nl}-\T{1}));{}$\6
\4${}\}{}$\2\par
\U29.\fi

\M{32}The \PB{\\{ivector}} routine allocates a real-valued vector of integer
precision,
with vector index ranging from \PB{\\{nl}} to \PB{\\{nh}}.

\Y\B\4\X32:Routine for allocation of vectors of integer precision\X${}\E{}$\6
\&{int} ${}{*}{}$\\{ivector}(\&{long} \\{nl}${},\39{}$\&{long} \\{nh})\1\1\2\2\6
${}\{{}$\1\6
\&{int} ${}{*}\|v;{}$\7
${}\|v\K{}$(\&{int} ${}{*}){}$ \\{malloc}${}((\&{size\_t})((\\{nh}-\\{nl}+%
\T{2})*\&{sizeof}(\&{int})));{}$\6
\&{if} ${}(\R\|v){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"Error:\ Allocation\ f}\)\.{ailure\ in\
ivector()\\}\)\.{n"});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\&{return} \|v${}-\\{nl}+\T{1};{}$\6
\4${}\}{}$\2\par
\A34.
\U29.\fi

\M{33}The \PB{\\{free\_ivector}} routine release the memory occupied by the
real-valued vector \PB{$\|v[\\{nl}\MRL{{.}{.}}\\{nh}]$}.

\Y\B\4\X33:Routine for deallocation of vectors of integer precision\X${}\E{}$\6
\&{void} \\{free\_ivector}(\&{int} ${}{*}\|v,\39{}$\&{long} \\{nl}${},\39{}$%
\&{long} \\{nh})\1\1\2\2\6
${}\{{}$\1\6
\\{free}((\&{char} ${}{*})(\|v+\\{nl}-\T{1}));{}$\6
\4${}\}{}$\2\par
\A35.
\U29.\fi

\M{34}The \PB{\\{livector}} routine allocates a real-valued vector of long
integer
precision, with vector index ranging from \PB{\\{nl}} to \PB{\\{nh}}.

\Y\B\4\X32:Routine for allocation of vectors of integer precision\X${}\mathrel+%
\E{}$\6
\&{long} \&{int} ${}{*}{}$\\{livector}(\&{long} \\{nl}${},\39{}$\&{long} %
\\{nh})\1\1\2\2\6
${}\{{}$\1\6
\&{long} \&{int} ${}{*}\|v;{}$\7
${}\|v\K{}$(\&{long} \&{int} ${}{*}){}$ \\{malloc}${}((\&{size\_t})((\\{nh}-%
\\{nl}+\T{2})*{}$\&{sizeof}(\&{long} \&{int})));\6
\&{if} ${}(\R\|v){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"Error:\ Allocation\ f}\)\.{ailure\ in\
livector()}\)\.{\\n"});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\&{return} \|v${}-\\{nl}+\T{1};{}$\6
\4${}\}{}$\2\par
\fi

\M{35}The \PB{\\{free\_livector}} routine release the memory occupied by the
real-valued vector \PB{$\|v[\\{nl}\MRL{{.}{.}}\\{nh}]$}.

\Y\B\4\X33:Routine for deallocation of vectors of integer precision\X${}%
\mathrel+\E{}$\6
\&{void} \\{free\_livector}(\&{long} \&{int} ${}{*}\|v,\39{}$\&{long} %
\\{nl}${},\39{}$\&{long} \\{nh})\1\1\2\2\6
${}\{{}$\1\6
\\{free}((\&{char} ${}{*})(\|v+\\{nl}-\T{1}));{}$\6
\4${}\}{}$\2\par
\fi

\M{36}The \PB{\\{simatrix}} routine allocates an array of short integer
precision,
with array row index ranging from \PB{\\{nrl}} to \PB{\\{nrh}} and column index
ranging
from \PB{\\{ncl}} to \PB{\\{nch}}.

\Y\B\4\X36:Routine for allocation of matrices of short integer precision\X${}%
\E{}$\6
\&{short} \&{int} ${}{*}{*}{}$\\{simatrix}(\&{long} \\{nrl}${},\39{}$\&{long} %
\\{nrh}${},\39{}$\&{long} \\{ncl}${},\39{}$\&{long} \\{nch})\1\1\2\2\6
${}\{{}$\1\6
\&{long} \|i${},{}$ \\{nrow}${}\K\\{nrh}-\\{nrl}+\T{1},{}$ \\{ncol}${}\K%
\\{nch}-\\{ncl}+\T{1};{}$\6
\&{short} \&{int} ${}{*}{*}\|m;{}$\7
${}\|m\K{}$(\&{short} \&{int} ${}{*}{*}){}$ \\{malloc}${}((\&{size\_t})((%
\\{nrow}+\T{1})*{}$\&{sizeof}(\&{short} \&{int} ${}{*})));{}$\6
\&{if} ${}(\R\|m){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Allocation\ fail}\)\.{ure\ 1\ in\
simatrix()\\}\)\.{n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\|m\MRL{+{\K}}\T{1};{}$\6
${}\|m\MRL{-{\K}}\\{nrl};{}$\6
${}\|m[\\{nrl}]\K{}$(\&{short} \&{int} ${}{*}){}$ \\{malloc}${}((\&{size\_t})((%
\\{nrow}*\\{ncol}+\T{1})*{}$\&{sizeof}(\&{short} \&{int})));\6
\&{if} ${}(\R\|m[\\{nrl}]){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Allocation\ fail}\)\.{ure\ 2\ in\
simatrix()\\}\)\.{n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\|m[\\{nrl}]\MRL{+{\K}}\T{1};{}$\6
${}\|m[\\{nrl}]\MRL{-{\K}}\\{ncl};{}$\6
\&{for} ${}(\|i\K\\{nrl}+\T{1};{}$ ${}\|i\Z\\{nrh};{}$ ${}\|i\PP){}$\1\5
${}\|m[\|i]\K\|m[\|i-\T{1}]+\\{ncol};{}$\2\6
\&{return} \|m;\6
\4${}\}{}$\2\par
\U29.\fi

\M{37}The \PB{\\{free\_simatrix}} routine releases the memory occupied by the
short integer matrix \PB{$\|v[\\{nl}\MRL{{.}{.}}\\{nh}]$}, as allocated by \PB{%
\\{simatrix}(\,)}.

\Y\B\4\X37:Routine for deallocation of matrices of short integer precision\X${}%
\E{}$\6
\&{void} \\{free\_simatrix}(\&{short} \&{int} ${}{*}{*}\|m,\39{}$\&{long} %
\\{nrl}${},\39{}$\&{long} \\{nrh}${},\39{}$\&{long} \\{ncl}${},\39{}$\&{long} %
\\{nch})\1\1\2\2\6
${}\{{}$\1\6
\\{free}((\&{char} ${}{*})(\|m[\\{nrl}]+\\{ncl}-\T{1}));{}$\6
\\{free}((\&{char} ${}{*})(\|m+\\{nrl}-\T{1}));{}$\6
\4${}\}{}$\2\par
\U29.\fi

\M{38}Routine for displaying help message to standard terminal output.

\Y\B\4\X38:Routine for displaying help message\X${}\E{}$\6
\&{void} \\{showsomehelp}(\&{void})\1\1\2\2\6
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"Usage:\ \%s\ [options]}\)\.{\\n"},\39%
\\{progname});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"Options:\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ -h,\ --help\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ Display\ this\ he}\)\.{lp\ message\
and\ exit\ }\)\.{clean.\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ -v,\ --verbose\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ Toggle\ verbose\ }\)\.{mode\ on/off.%
\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ -o,\ --outputfile\ <}\)\.{str>\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ Specifies\ the\ b}\)\.{ase\ name\ of\
the\ outp}\)\.{ut\ files\ where\ the\ p}\)\.{rogram\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ is\ to\ save\ the\ }\)\.{calculated\
data.\ If\ }\)\.{the\ --outputfile\ or\ }\)\.{-o\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ option\ is\ follo}\)\.{wed\ by\
'append'\ the\ }\)\.{estimate\ for\ the\ fra}\)\.{ctal\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ dimension\ will\ }\)\.{be\ appended\
to\ the\ f}\)\.{ile\ named\ <str>.dat,}\)\.{\ which\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ will\ be\ created}\)\.{\ if\ it\ does\
not\ exis}\)\.{t.\ If\ the\ following\ }\)\.{word\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ instead\ is\ `ove}\)\.{rwrite'\ the\
file\ wil}\)\.{l\ instead\ be\ overwri}\)\.{tten.\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ -i,\ --trajectoryfi}\)\.{le\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ Specifies\ the\ i}\)\.{nput\
trajectory\ of\ t}\)\.{he\ graph\ whose\ fract}\)\.{al\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ dimension\ is\ to}\)\.{\ be\
estimated.\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ -w,\ --computationw}\)\.{indow\ llx\ <num>\
lly\ }\)\.{<num>\ urx\ <num>\ ury\ }\)\.{<num>\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ This\ option\ exp}\)\.{licitly\
specifies\ th}\)\.{e\ domain\ over\ which\ }\)\.{the\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ box-counting\ al}\)\.{gorithm\ will\
be\ appl}\)\.{ied,\ in\ terms\ of\ the}\)\.{\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ lower-left\ and\ }\)\.{upper-right\
corners\ }\)\.{(llx,lly)\ and\ (urx,u}\)\.{ry),\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ respectively.\ B}\)\.{y\ specifying\
this\ op}\)\.{tion,\ any\ automatic\\}\)\.{n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ calculation\ of\ }\)\.{computational\
window}\)\.{\ will\ be\ neglected.\\}\)\.{n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ -m,\ --boxmaps\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ If\ this\ option\ }\)\.{is\ present,\
the\ prog}\)\.{ram\ will\ generate\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ MetaPost\ code\ f}\)\.{or\ maps\ of\
the\ distr}\)\.{ibution\ of\ boxes.\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ In\ doing\ so,\ al}\)\.{so\ the\ input%
\ traject}\)\.{ory\ is\ included\ in\\n}\)\.{"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ the\ graphs.\ The}\)\.{\ convention\
for\ the\ }\)\.{naming\ of\ the\ output}\)\.{\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ map\ files\ is\ th}\)\.{at\ they\ are\
saved\ as}\)\.{\ <str>.<N>.dat,\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ where\ <str>\ is\ }\)\.{the\ base\
filename\ as}\)\.{\ specified\ using\ the}\)\.{\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ -o\ or\ --outputf}\)\.{ile\ option,\
<N>\ is\ t}\)\.{he\ automatically\ app}\)\.{ended\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ current\ level\ o}\)\.{f\ resolution\
refinem}\)\.{ent\ in\ the\ box-count}\)\.{ing,\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ and\ where\ '.dat}\)\.{'\ is\ the\
file\ suffix}\)\.{\ as\ automatically\ ap}\)\.{pended\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ by\ the\ program.}\)\.{\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ --graphsize\ <width}\)\.{\ in\ mm>\ <height\
in\ m}\)\.{m>\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ If\ the\ -m\ or\ --}\)\.{boxmaps\
option\ is\ pr}\)\.{esent\ at\ the\ command}\)\.{\ line,\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ then\ the\ --grap}\)\.{hsize\ option\
will\ ov}\)\.{erride\ the\ default\ g}\)\.{raph\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ size\ of\ the\ gen}\)\.{erated\
boxmaps.\ (Def}\)\.{ault\ graph\ size\ is\ 8}\)\.{0\ mm\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ width\ and\ 56\ mm}\)\.{\ height.)%
\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ -Nmin,\ --minlevel\ }\)\.{<num>\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ Specifies\ the\ m}\)\.{inimum\ level\
Nmin\ of}\)\.{\ grid\ refinement\ tha}\)\.{t\ \\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ will\ be\ used\ in}\)\.{\ the\
evaluation\ of\ t}\)\.{he\ estimate\ of\ the\ f}\)\.{ractal\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ dimension.\ With}\)\.{\ this\ option\
specifi}\)\.{ed,\ the\ coarsest\ lev}\)\.{el\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ used\ in\ the\ box}\)\.{-counting\
will\ be\ a\ }\)\.{[(2\^Nmin)x(2\^Nmin)]-}\)\.{grid\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ of\ boxes.\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ -Nmax,\ --maxlevel\ }\)\.{<num>\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ This\ option\ is\ }\)\.{similar\ to\
the\ --min}\)\.{level\ or\ -Nmin\ optio}\)\.{ns,\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ with\ the\ differ}\)\.{ence\ that\ it\
instead}\)\.{\ specifies\ the\ maxim}\)\.{um\\n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ level\ Nmax\ of\ g}\)\.{rid\
refinement\ that\ }\)\.{will\ be\ used\ in\ the\\}\)\.{n"});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\ \ \ \ evaluation\ of\ t}\)\.{he\ estimate\
of\ the\ f}\)\.{ractal\ dimension.\\n"});{}$\6
\4${}\}{}$\2\par
\U23.\fi

\M{39}The \PB{$\\{lines\_intersect}(\\{p1x},\\{p1y},\\{q1x},\\{q1y},\\{p2x},%
\\{p2y},\\{q2x},\\{q2y})$} routine takes the
start and end points of two line segments, the first between $(\PB{\\{p1x}},%
\PB{\\{p1y}})$
and $(\PB{\\{q1x}},\PB{\\{q1y}})$ and the second between $(\PB{\\{p2x}},\PB{%
\\{p2y}})$ and $(\PB{\\{q2x}},\PB{\\{q2y}})$,
and returns $1$ (`true') if they are found to intersect, and $0$ (`false')
otherwise.

For a brief sumnmary of the algorithm behind this routine, consider two line
segments in the plane, the first one between the points ${\bf p}_1$ and
${\bf q}_1$ and the second one between ${\bf p}_2$ and ${\bf q}_1$.
In general, these segments can be written in parametric forms as
${\bf r}_1={\bf p}_1+t_1({\bf q}_1-{\bf p}_1)$ and
${\bf r}_2={\bf p}_2+t_2({\bf q}_2-{\bf p}_2)$
for $0\le t_1\le1$ and  $0\le t_2\le1$.
These line segments intersect each other if they for these intervals for the
parameters $t_1$ and $t_2$ share a common point, or equivalently if the
solutions to
$$
{\bf p}_1+t_1({\bf q}_1-{\bf p}_1)={\bf p}_2+t_2({\bf q}_2-{\bf p}_2)
\qquad\Leftrightarrow\qquad
({\bf q}_1-{\bf p}_1)t_1+({\bf p}_2-{\bf q}_2)t_2={\bf p}_2-{\bf p}_1
$$
both simultaneously satisfy $0\le t_1\le1$ and $0\le t_2\le1$.
This vectorial equation can equivalently be written in component form as
$$
\eqalign{
(q_{1x}-p_{1x})t_1+(p_{2x}-q_{2x})t_2&=p_{2x}-p_{1x},\cr
(q_{1y}-p_{1y})t_1+(p_{2y}-q_{2y})t_2&=p_{2y}-p_{1y},\cr
}
$$
which after some straightforward algebra gives the explicit solutions for the
parameters $t_1$ and $t_2$ as
$$
t_1={{ed-bf}\over{ad-bc}},\qquad t_2={{af-ec}\over{ad-bc}},
$$
where
$$
\eqalign{
a\equiv(q_{1x}-p_{1x}),\qquad
b\equiv(p_{2x}-q_{2x}),\qquad
c\equiv(q_{1y}-p_{1y}),\cr
d\equiv(p_{2y}-q_{2y}),\qquad
e\equiv(p_{2x}-p_{1x}),\qquad
f\equiv(p_{2x}-p_{1x}).\cr
}
$$
Hence, the two line segments intersect if and only if
$$
0\le{{ed-bf}\over{ad-bc}}\le1,\qquad
{\rm and}\qquad 0\le{{af-ec}\over{ad-bc}}\le1.
$$
By observing that their denominators are equal, the evaluation of these quotes
involves in total 6 floating-point multiplications, 2 divisions and 3
subtractions.
Notice that whenever $ad-bd=0$, the two lines are parallell and will never
intersect, regardless of the values of $t_1$ and $t_2$.

\Y\B\4\X39:Routine for determining whether two lines intersect or not\X${}\E{}$%
\6
\&{short} \&{int} \\{lines\_intersect}(\&{double} \\{p1x}${},\39{}$\&{double} %
\\{p1y}${},\39{}$\&{double} \\{q1x}${},\39{}$\&{double} \\{q1y}${},\39{}$%
\&{double} \\{p2x}${},\39{}$\&{double} \\{p2y}${},\39{}$\&{double} \\{q2x}${},%
\39{}$\&{double} \\{q2y})\1\1\2\2\6
${}\{{}$\1\6
\&{double} \|a${},{}$ \|b${},{}$ \|c${},{}$ \|d${},{}$ \|e${},{}$ \|f${},{}$ %
\\{det}${},{}$ \\{tmp1}${},{}$ \\{tmp2};\6
\&{short} \&{int} \\{intersect};\7
${}\|a\K\\{q1x}-\\{p1x};{}$\6
${}\|b\K\\{p2x}-\\{q2x};{}$\6
${}\|c\K\\{q1y}-\\{p1y};{}$\6
${}\|d\K\\{p2y}-\\{q2y};{}$\6
${}\|e\K\\{p2x}-\\{p1x};{}$\6
${}\|f\K\\{p2y}-\\{p1y};{}$\6
${}\\{det}\K\|a*\|d-\|b*\|c;{}$\6
${}\\{tmp1}\K\|e*\|d-\|b*\|f;{}$\6
${}\\{tmp2}\K\|a*\|f-\|e*\|c;{}$\6
${}\\{intersect}\K\T{0};{}$\6
\&{if} ${}(\\{det}>\T{0}){}$\5
${}\{{}$\1\6
\&{if} ${}(((\T{0.0}\Z\\{tmp1})\W(\\{tmp1}\Z\\{det}))\W((\T{0.0}\Z\\{tmp2})\W(%
\\{tmp2}\Z\\{det}))){}$\1\5
${}\\{intersect}\K\T{1};{}$\2\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\\{det}<\T{0}){}$\5
${}\{{}$\1\6
\&{if} ${}(((\\{det}\Z\\{tmp1})\W(\\{tmp1}\Z\T{0.0}))\W((\\{det}\Z\\{tmp2})\W(%
\\{tmp2}\Z\T{0.0}))){}$\1\5
${}\\{intersect}\K\T{1};{}$\2\6
\4${}\}{}$\2\6
\&{return} (\\{intersect});\6
\4${}\}{}$\2\par
\U23.\fi

\M{40}Routine for determining whether a line and a box intersect or not.
The \PB{\\{box\_intersection}(\,)} routine simply divides the input box, being
characterized by its lower-left and upper-right corners $(\PB{\\{llx}},\PB{%
\\{lly}})$
and $(\PB{\\{urx}},\PB{\\{ury}})$, into the four line segments corresponding to
its four
edges, followed by calling the routine \PB{\\{lines\_intersect}(\,)} to
determine if
any of these edges intersect the line segment. If an intersection is found,
the \PB{\\{box\_intersection}(\,)} routine returns 1 (true), otherwise 0
(false).
\medskip\noindent
Input variables:
\varitem[\PB{\\{px}}, \PB{\\{py}}]{The coordinates of the first end point
${\bf p}=(p_x,p_y)$ of the line segment.}
\varitem[\PB{\\{qx}}, \PB{\\{qy}}]{The coordinates of the second end point
${\bf q}=(q_x,q_y)$ of the line segment.}
\varitem[\PB{\\{llx}}, \PB{\\{lly}}]{Coordinates of the lower-left corner of
the box.}
\varitem[\PB{\\{urx}}, \PB{\\{ury}}]{Coordinates of the upper-right corner of
the box.}
\medskip\noindent
Output variables:
\varitem[]{On exit, the routine returns 1 if an intersection is found,
otherwise 0 is returned, in either case the value are returned as
integers of \PB{\&{short}} precision.}
\medskip

\Y\B\4\X40:Routine for determining whether a line and a box intersect or not%
\X${}\E{}$\6
\&{short} \&{int} \\{box\_intersection}(\&{double} \\{px}${},\39{}$\&{double} %
\\{py}${},\39{}$\&{double} \\{qx}${},\39{}$\&{double} \\{qy}${},\39{}$%
\&{double} \\{llx}${},\39{}$\&{double} \\{lly}${},\39{}$\&{double} \\{urx}${},%
\39{}$\&{double} \\{ury})\1\1\2\2\6
${}\{{}$\1\6
\&{if} ${}(\\{lines\_intersect}(\\{px},\39\\{py},\39\\{qx},\39\\{qy},\39%
\\{llx},\39\\{lly},\39\\{urx},\39\\{lly})){}$\1\5
\&{return} (\T{1});\C{ intersection with bottom edge }\2\6
\&{if} ${}(\\{lines\_intersect}(\\{px},\39\\{py},\39\\{qx},\39\\{qy},\39%
\\{urx},\39\\{lly},\39\\{urx},\39\\{ury})){}$\1\5
\&{return} (\T{1});\C{ intersection with right edge }\2\6
\&{if} ${}(\\{lines\_intersect}(\\{px},\39\\{py},\39\\{qx},\39\\{qy},\39%
\\{urx},\39\\{ury},\39\\{llx},\39\\{ury})){}$\1\5
\&{return} (\T{1});\C{ intersection with top edge }\2\6
\&{if} ${}(\\{lines\_intersect}(\\{px},\39\\{py},\39\\{qx},\39\\{qy},\39%
\\{llx},\39\\{ury},\39\\{llx},\39\\{lly})){}$\1\5
\&{return} (\T{1});\C{ intersection with left edge }\2\6
\&{return} (\T{0});\C{ this happens only if no edge is intersecting the line
segment }\6
\4${}\}{}$\2\par
\U23.\fi

\M{41}Routines for calculation of number of boxes covering the trajectory.
There
are two almost identical routines for the calculation of the number of boxes
covering the input trajectory at a given level of subdivision of the box sizes.
The first routine, \PB{\\{get\_num\_covering\_boxes}(\,)} simply performs this
task
without caring of documenting the box distributions as graphs, or ``box maps'',
while the second one, \PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}(\,)}
also includes
the generation of these maps, with output in terms of \METAPOST\ code.

\Y\B\4\X41:Routines for calculation of number of boxes covering the trajectory%
\X${}\E{}$\6
\X42:Routine for calculation of number of boxes covering the trajectory\X\6
\X52:Simplified routine for calculation of number of boxes covering the
trajectory\X\par
\U23.\fi

\M{42}Routine for calculation of number of boxes covering the trajectory, also
including the generation of documenting graphs of the box distributions.
The \PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}(\,)} routine takes a
trajectory as
described by a discrete set of coordinates as input, and for a given grid
refinement $N$ returns the number of boxes needed to entirely cover the
trajectory. The grid refinement factor $N$ indicates that the domain of
computation is divided into a $[2^N\times2^N]$-grid of boxes.

The computational domain in which the box counting is to be performed is
explicitly stated by the coordinates of its lower-left and upper-right
corners, $(global_llx,global_lly)$ and $(global_urx,global_ury)$, respectively.
Parts of the trajectory which are outside of this domain are not included in
the box-counting.
If the entire input trajectory is to be used in the box counting, simply use
$(\PB{\\{global\_llx}},\PB{\\{global\_lly}})=({\rm min}[\PB{\\{x\_traj}}],{\rm
min}[\PB{\\{y\_traj}}])$ and
$(\PB{\\{global\_urx}},\PB{\\{global\_ury}})=({\rm max}[\PB{\\{x\_traj}}],{\rm
max}[\PB{\\{y\_traj}}])$ for the
specification of the computational domain.
\medskip\noindent
Input variables:
\varitem[\PB{\\{mm}}]{The total number of coordinates forming the input
trajectory,
or equivalently the length of the vectors \PB{${*}\\{x\_traj}$} and \PB{${*}%
\\{y\_traj}$}.}
\varitem[\PB{${*}\\{x\_traj}$}, \PB{${*}\\{y\_traj}$}]{Vectors of length \PB{%
\\{mm}} containing the $x$- and
$y$-coordinates of the input trajectory.}
\varitem[\PB{\\{resolution}}]{The grid refinement factor $N$.}
\varitem[\PB{\\{global\_llx}}, \PB{\\{global\_lly}}]{Coordinates of the
lower-left corner
of the computational domain in which the box-counting is to be performed.}
\varitem[\PB{\\{global\_urx}}, \PB{\\{global\_ury}}]{Coordinates of the
upper-right corner
of the computational domain in which the box-counting is to be performed.}
\varitem[\PB{\\{trajectory\_filename}}]{String containing the name of the file
containing the input trajectory.}
\varitem[\PB{\\{boxgraph\_filename}}]{String which, if non-empty, states the
file name
to which the map of the spatial distribution of the covering boxes is to be
written, as \METAPOST\ code.}
\medskip\noindent
Output variables:
\varitem[]{On exit, the routine returns the number of covering boxes
as an integer of \PB{\&{long} \&{unsigned}} precision.}
\medskip\noindent
Internal variables:
\varitem[\PB{\\{px}},\PB{\\{py}},\PB{\\{qx}},\PB{\\{qy}}]{Keeps track of the
$x$- and $y$-coordinates of
the start and end points of line segments, between ${\bf p}=(p_x,p_y)$ and
${\bf q}=(q_x,q_y)$.}
\medskip

\Y\B\4\X42:Routine for calculation of number of boxes covering the trajectory%
\X${}\E{}$\6
\&{long} \&{unsigned} \&{int} \\{get\_num\_covering\_boxes\_with\_boxmaps}(%
\&{double} ${}{*}\\{x\_traj},\39{}$\&{double} ${}{*}\\{y\_traj},\39{}$\&{long} %
\&{int} \\{mm}${},\39{}$\&{int} \\{resolution}${},\39{}$\&{double} \\{global%
\_llx}${},\39{}$\&{double} \\{global\_lly}${},\39{}$\&{double} \\{global%
\_urx}${},\39{}$\&{double} \\{global\_ury}${},\39{}$\&{char} \\{boxgraph%
\_filename}[\,]${},\39{}$\&{double} \\{boxgraph\_width}${},\39{}$\&{double} %
\\{boxgraph\_height}${},\39{}$\&{char} \\{trajectory\_filename}[\,])\1\1\2\2\6
${}\{{}$\1\6
\&{short} \&{int} ${}{*}{*}\\{box};{}$\6
\&{long} \&{unsigned} \&{int} \|i${},{}$ \|j${},{}$ \|m${},{}$ \\{nn}${},{}$ %
\\{istart}${},{}$ \\{jstart}${},{}$ \\{istop}${},{}$ \\{jstop}${},{}$ %
\\{iincr}${},{}$ \\{jincr}${},{}$ \\{num\_boxes};\6
\&{double} ${}{*}\\{x\_box},{}$ ${}{*}\\{y\_box}{}$;\C{ Keeps track of the
lower-left coordinates of                             the boxes. }\6
\&{double} \\{px}${},{}$ \\{py}${},{}$ \\{qx}${},{}$ \\{qy};\6
\&{FILE} ${}{*}\\{boxgraph\_file};{}$\7
\X43:Set up the grid of $2^N\times2^N$ boxes covering the entire global window
of computation\X\6
\X44:Find indices $(i_a,j_a)$ of the box covering the first coordinate of the
trajectory\X\6
\&{for} ${}(\|m\K\T{1};{}$ ${}\|m\Z\\{mm}-\T{1};{}$ ${}\|m\PP){}$\5
${}\{{}$\1\6
\X45:Find indices $(i_b,j_b)$ of the box covering the end point of $m$th
trajectory segment\X\6
\X46:Scan the $m$th trajectory segment for intersecting boxes\X\6
\4${}\}{}$\2\6
\X47:Open file for output of box distribution graph\X\6
\X48:Write the input trajectory to the box distribution graph\X\6
\X49:Count the total number of boxes \PB{\\{num\_boxes}} intersecting the
trajectory\X\6
\X50:Write closing blocks the box distribution graph\X\6
\X51:Close any open file for output of box distribution graph\X\6
\&{return} (\\{num\_boxes});\6
\4${}\}{}$\2\par
\U41.\fi

\M{43}Set up the grid of $2^N\times2^N$ boxes covering the entire global window
of computation.
In this block, the resolution of the grid of boxes is defined. Notice that in
many cases, only a certain subset of boxes will actally intersect the input
trajectory. However, we do not \'a priori know this number of boxes, and in
order to safeguard against the possible danger of running out of allocated
memory, with the need for a dynamic memory allocation depending on the state
of computation, a matrix of short integer precision is allocated covering the
entire computational window.
In order to keep track of the coordinates of the boxes, two vectors
$\PB{\\{x\_box}}[1\ldots 2^N]$ and $\PB{\\{y\_box}}[1\ldots 2^N]$ are allocated
to contain
the $x$- and $y$-coordinates of the lower-left corners of the boxes, with the
last elements $\PB{\\{x\_box}}[2^N]$ and $\PB{\\{y\_box}}[2^N]$ containing the
upper-right
corner coordinates of the upper-right-most box of the global window.

\Y\B\4\X43:Set up the grid of $2^N\times2^N$ boxes covering the entire global
window of computation\X${}\E{}$\6
${}\{{}$\1\6
${}\\{nn}\K\T{1};{}$\6
\&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{resolution};{}$ ${}\|i\PP){}$\1\5
${}\\{nn}\K\T{2}*\\{nn};{}$\2\6
${}\\{box}\K\\{simatrix}(\T{1},\39\\{nn},\39\T{1},\39\\{nn});{}$\6
\&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{nn};{}$ ${}\|i\PP){}$\1\6
\&{for} ${}(\|j\K\T{1};{}$ ${}\|j\Z\\{nn};{}$ ${}\|j\PP){}$\1\5
${}\\{box}[\|i][\|j]\K\T{0};{}$\2\2\6
${}\\{x\_box}\K\\{dvector}(\T{1},\39\\{nn}+\T{1});{}$\6
${}\\{y\_box}\K\\{dvector}(\T{1},\39\\{nn}+\T{1});{}$\6
\&{for} ${}(\|m\K\T{1};{}$ ${}\|m\Z\\{nn}+\T{1};{}$ ${}\|m\PP){}$\5
${}\{{}$\1\6
${}\\{x\_box}[\|m]\K\\{global\_llx}+((\&{double})(\|m-\T{1}))*(\\{global\_urx}-%
\\{global\_llx})/((\&{double})(\\{nn}));{}$\6
${}\\{y\_box}[\|m]\K\\{global\_lly}+((\&{double})(\|m-\T{1}))*(\\{global\_ury}-%
\\{global\_lly})/((\&{double})(\\{nn}));{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U42.\fi

\M{44}Find indices $(i_a,j_a)$ of the box covering the first coordinate of the
trajectory.

\Y\B\4\X44:Find indices $(i_a,j_a)$ of the box covering the first coordinate of
the trajectory\X${}\E{}$\6
${}\{{}$\1\6
${}\\{istart}\K\T{0};{}$\6
${}\\{jstart}\K\T{0};{}$\6
\&{for} ${}(\|m\K\T{1};{}$ ${}\|m\Z\\{nn};{}$ ${}\|m\PP){}$\5
${}\{{}$\1\6
\&{if} ${}((\\{x\_box}[\|m]\Z\\{x\_traj}[\T{1}])\W(\\{x\_traj}[\T{1}]\Z\\{x%
\_box}[\|m+\T{1}])){}$\1\5
${}\\{istart}\K\|m;{}$\2\6
\&{if} ${}((\\{y\_box}[\|m]\Z\\{y\_traj}[\T{1}])\W(\\{y\_traj}[\T{1}]\Z\\{y%
\_box}[\|m+\T{1}])){}$\1\5
${}\\{jstart}\K\|m;{}$\2\6
\4${}\}{}$\2\6
\&{if} ${}((\\{istart}\E\T{0})\V(\\{jstart}\E\T{0})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Error!\ Cannot\ f}\)\.{ind\ box\ indices%
\ of\ 1}\)\.{st\ coordinate!\\n"},\39\\{progname});{}$\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Please\ check\ da}\)\.{ta\ for\ input\
trajeto}\)\.{ry.\\n"},\39\\{progname});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U42.\fi

\M{45}Find indices $(i_b,j_b)$ of the box covering the end point of the $m$th
trajectory segment.

\Y\B\4\X45:Find indices $(i_b,j_b)$ of the box covering the end point of $m$th
trajectory segment\X${}\E{}$\6
${}\{{}$\1\6
${}\\{px}\K\\{x\_traj}[\|m];{}$\6
${}\\{py}\K\\{y\_traj}[\|m];{}$\6
${}\\{qx}\K\\{x\_traj}[\|m+\T{1}];{}$\6
${}\\{qy}\K\\{y\_traj}[\|m+\T{1}];{}$\6
${}\\{istop}\K\\{istart};{}$\6
${}\\{jstop}\K\\{jstart};{}$\6
\&{if} ${}(\\{px}<\\{qx}){}$\5
${}\{{}$\1\6
${}\\{iincr}\K\T{1};{}$\6
\&{while} ${}(\\{x\_box}[\\{istop}+\T{1}]<\\{qx}){}$\1\5
${}\\{istop}\PP;{}$\2\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
${}\\{iincr}\K{-}\T{1};{}$\6
\&{while} ${}(\\{x\_box}[\\{istop}]>\\{qx}){}$\1\5
${}\\{istop}\MM;{}$\2\6
\4${}\}{}$\2\6
\&{if} ${}(\\{py}<\\{qy}){}$\5
${}\{{}$\1\6
${}\\{jincr}\K\T{1};{}$\6
\&{while} ${}(\\{y\_box}[\\{jstop}+\T{1}]<\\{qy}){}$\1\5
${}\\{jstop}\PP;{}$\2\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
${}\\{jincr}\K{-}\T{1};{}$\6
\&{while} ${}(\\{y\_box}[\\{jstop}]>\\{qy}){}$\1\5
${}\\{jstop}\MM;{}$\2\6
\4${}\}{}$\2\6
\&{if} ${}(\T{0}\E\T{1}){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ Endpoint\ box\ in}\)\.{dices:\ (i=%
\%ld,j=\%ld)}\)\.{\\n"},\39\\{progname},\39\\{istop},\39\\{jstop});{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U42.\fi

\M{46}Scan the $m$th trajectory segment for intersecting boxes.
As the indices of the boxes covering the start and end points of the $m$th
trajectory segment have been previously determined, one may now use this
information in order to reduce the search for intersecting boxes to the
domain as covered by box indices $i_{\rm start}\le i\le i_{\rm stop}$ and
$j_{\rm start}\le j\le j_{\rm stop}$.
This way, the computational time needed is greatly reduced as compared to if
the entire window would be scanned for every segment.

\Y\B\4\X46:Scan the $m$th trajectory segment for intersecting boxes\X${}\E{}$\6
${}\{{}$\1\6
\&{for} ${}(\|i\K\\{istart};{}$ ${}\|i\I(\\{istop}+\\{iincr});{}$ ${}\|i\MRL{+{%
\K}}\\{iincr}){}$\5
${}\{{}$\1\6
\&{for} ${}(\|j\K\\{jstart};{}$ ${}\|j\I(\\{jstop}+\\{jincr});{}$ ${}\|j\MRL{+{%
\K}}\\{jincr}){}$\5
${}\{{}$\1\6
\&{if} ${}(\\{box\_intersection}(\\{px},\39\\{py},\39\\{qx},\39\\{qy},\39\\{x%
\_box}[\|i],\39\\{y\_box}[\|j],\39\\{x\_box}[\|i+\T{1}],\39\\{y\_box}[\|j+%
\T{1}])){}$\5
${}\{{}$\1\6
${}\\{box}[\|i][\|j]\K\T{1};{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
${}\\{istart}\K\\{istop};{}$\6
${}\\{jstart}\K\\{jstop};{}$\6
\4${}\}{}$\2\par
\U42.\fi

\M{47}Open file for output of box distribution graph. In this block the
preamble
of the output \METAPOST\ code is written to file. This preamble contains a
macro {\tt box(i,j)} which simply is a parameter-specific \METAPOST\ subroutine
which is used in order to reduce the final size of the code to be compiled
into a graph over the distribution of covering boxes.

\Y\B\4\X47:Open file for output of box distribution graph\X${}\E{}$\6
${}\{{}$\1\6
\&{if} ${}(\R\\{strcmp}(\\{boxgraph\_filename},\39\.{""})){}$\5
${}\{{}$\C{ is \PB{\\{boxgraph\_filename}} an empty string? }\1\6
${}\\{boxgraph\_file}\K\NULL;{}$\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
\&{if} ${}((\\{boxgraph\_file}\K\\{fopen}(\\{boxgraph\_filename},\39\.{"w"}))\E%
\NULL){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"\%s:\ Could\ not\ open\ }\)\.{\%s\ for\ box\
graphs!\\n}\)\.{"},\39\\{progname},\39\\{boxgraph\_filename});{}$\6
\\{exit}(\.{FAILURE});\6
\4${}\}{}$\2\6
${}\\{fseek}(\\{boxgraph\_file},\39\T{0\$L},\39\.{SEEK\_SET});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"input\ graph;\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"def\ box(expr\ i,j)=\\}\)\.{n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"begingroup\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"llx:=\%2.8f+(i-1)*\%2}\)\.{.8f;\\n"},%
\39\\{global\_llx},\39(\\{global\_urx}-\\{global\_llx})/(\\{nn}));{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"lly:=\%2.8f+(j-1)*\%2}\)\.{.8f;\\n"},%
\39\\{global\_lly},\39(\\{global\_ury}-\\{global\_lly})/(\\{nn}));{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"urx:=\%2.8f+(i)*\%2.8}\)\.{f;\\n"},\39%
\\{global\_llx},\39(\\{global\_urx}-\\{global\_llx})/(\\{nn}));{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"ury:=\%2.8f+(j)*\%2.8}\)\.{f;\\n"},\39%
\\{global\_lly},\39(\\{global\_ury}-\\{global\_lly})/(\\{nn}));{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (llx,lly)--(u}\)\.{rx,lly);%
\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (urx,lly)--(u}\)\.{rx,ury);%
\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (urx,ury)--(l}\)\.{lx,ury);%
\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (llx,ury)--(l}\)\.{lx,lly);%
\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"endgroup\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"enddef;\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"beginfig(1);\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"w:=\%2.4fmm;\ h:=\%2.4}\)\.{fmm;\\n"},%
\39\\{boxgraph\_width},\39\\{boxgraph\_height});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"draw\ begingraph(w,h}\)\.{);\\n"});{}$%
\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"pickup\ pencircle\ sc}\)\.{aled\ .3pt;%
\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"setrange(\%2.6f,\%2.6}\)\.{f,\%2.6f,%
\%2.6f);\\n"},\39\\{global\_llx},\39\\{global\_lly},\39\\{global\_urx},\39%
\\{global\_ury});{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U42.\fi

\M{48}Write the input trajectory to the box distribution graph. Here there are
two
possible choices for how the input trajectory is to be included in the box
graph.

First, we may use \MP\ to automatically scan and draw the trajectory
for us, simply by using a statment like {\tt gdraw input.dat}, assuming that
the file {\tt input.dat} contains properly formatted columnwise data.
However, this choice would have two major drawbacks, namely that the generated
code would be dependent on that the original input file always is in place,
hence not allowing the \MP\ code to be exported as a standalone code as such,
and also that this would put a limitation on the number of nodes allowed in the
input trajectory, as the {\tt graph.mp} macro package of \MP\ only accepts
roughly 4000 points before it cuts the mapping.

The other choice is to include the input trajectory directly into the generated
code, preferrably by chopping the trajectory into pieces small enough to easily
be mapped by \MP\ without reaching a critical limit of exhaustion.
This choice of course significantly increases the file size of the generated
code, but this is a price we will have to accept in order to get stand-alone
output.
In the \boxcount\ program, the second alternative was naturtally chosen, simply
because the author is a fan of self-sustaining code which can be exported for
later use.

In this block, the status of the pointer \PB{\\{boxgraph\_file}} is used to
determine
whether to write the trajectory to file or not. If \PB{\\{boxgraph\_file}}
equals to
\PB{$\NULL$} (NULL), then the \boxcount\ program will not attempt to write the
input
trajectory to file.

\Y\B\4\X48:Write the input trajectory to the box distribution graph\X${}\E{}$\6
${}\{{}$\1\6
\&{if} ${}(\\{boxgraph\_file}\I\NULL){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"pickup\ pencircle\ sc}\)\.{aled\ .5pt;%
\\n"});{}$\6
${}\|i\K\T{0};{}$\6
\&{for} ${}(\|m\K\T{1};{}$ ${}\|m\Z\\{mm};{}$ ${}\|m\PP){}$\5
${}\{{}$\1\6
\&{if} ${}(\|i\E\T{0}){}$\5
${}\{{}$\1\6
\&{if} ${}(\|m\E\T{1}){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (\%2.4f,\%2.4f)}\)\.{"},\39\\{x%
\_traj}[\|m],\39\\{y\_traj}[\|m]);{}$\6
\4${}\}{}$\2\6
\&{else} \&{if} ${}(\T{2}<\\{mm}){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"gdraw\ (\%2.4f,\%2.4f)}\)\.{--(\%2.4f,%
\%2.4f)"},\39\\{x\_traj}[\|m-\T{1}],\39\\{y\_traj}[\|m-\T{1}],\39\\{x\_traj}[%
\|m],\39\\{y\_traj}[\|m]);{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\&{else}\5
${}\{{}$\1\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"--(\%2.4f,\%2.4f)"},\39\\{x\_traj}[%
\|m],\39\\{y\_traj}[\|m]);{}$\6
\4${}\}{}$\2\6
${}\|i\PP;{}$\6
\&{if} ${}((\|i\E\T{5})\V(\|m\E\\{mm})){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{";\\n"});{}$\6
${}\|i\K\T{0};{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U42.\fi

\M{49}Count the total number of boxes \PB{\\{num\_boxes}} intersecting the
trajectory.
Having traversed the entire trajectory at the current depth of resolution,
the only remaining task is to sum up the total number of boxes needed to
cover the entire trajectory.
This is simply done by extracting the number of set elements in the \PB{%
\\{box}}
matrix.
Having extracted the total number of boxes, this block also takes care of
releasing the memory occupied by the \PB{\\{box}} matrix, as this memory might
be
needed for iterations to come in which an even finer mesh of boxes is used.

\Y\B\4\X49:Count the total number of boxes \PB{\\{num\_boxes}} intersecting the
trajectory\X${}\E{}$\6
${}\{{}$\1\6
${}\\{num\_boxes}\K\T{0};{}$\6
\&{for} ${}(\|i\K\T{1};{}$ ${}\|i\Z\\{nn};{}$ ${}\|i\PP){}$\5
${}\{{}$\1\6
\&{for} ${}(\|j\K\T{1};{}$ ${}\|j\Z\\{nn};{}$ ${}\|j\PP){}$\5
${}\{{}$\1\6
\&{if} ${}(\\{box}[\|i][\|j]\E\T{1}){}$\5
${}\{{}$\1\6
${}\\{num\_boxes}\PP;{}$\6
\&{if} ${}(\\{boxgraph\_file}\I\NULL){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"box(\%ld,\%ld);\\n"},\39\|i,\39%
\|j);{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
\4${}\}{}$\2\6
${}\\{free\_simatrix}(\\{box},\39\T{1},\39\\{nn},\39\T{1},\39\\{nn});{}$\6
\4${}\}{}$\2\par
\U42.\fi

\M{50}Write closing blocks the box distribution graph. Here follows the
specification of tick marking (set as automatic for the sake of simplicity)
and axis labels, just before the closing statements of the \METAPOST\ code
for the graphs of box distributions.

\Y\B\4\X50:Write closing blocks the box distribution graph\X${}\E{}$\6
${}\{{}$\1\6
\&{if} ${}(\\{boxgraph\_file}\I\NULL){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"autogrid(itick\ bot,}\)\.{itick\ lft);%
\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"glabel.bot(btex\ \$x\$}\)\.{\
etex,OUT);\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"glabel.lft(btex\ \$y\$}\)\.{\
etex,OUT);\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"endgraph;\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"endfig;\\n"});{}$\6
${}\\{fprintf}(\\{boxgraph\_file},\39\.{"end\\n"});{}$\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U42.\fi

\M{51}Close any open file for output of box distribution graph.

\Y\B\4\X51:Close any open file for output of box distribution graph\X${}\E{}$\6
${}\{{}$\1\6
\&{if} ${}(\\{boxgraph\_file}\I\NULL){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stdout},\39\.{"\%s:\ MetaPost\ output}\)\.{\ box\
distribution\ gr}\)\.{aph\ written\ to\ \%s.\\n}\)\.{"},\39\\{progname},\39%
\\{boxgraph\_filename});{}$\6
\\{fclose}(\\{boxgraph\_file});\6
\4${}\}{}$\2\6
\4${}\}{}$\2\par
\U42.\fi

\M{52}Routine for calculation of number of boxes covering the trajectory. This
routine provides a simplified interface to the general boxcountng routine
\PB{\\{get\_num\_covering\_boxes\_with\_boxmaps}}, with the only difference
that no
graph of the box distribution over the trajectory is generated.
The \PB{\\{get\_num\_covering\_boxes}(\,)} routine takes a trajectory as
described by a discrete set of coordinates as input, and for a given grid
refinement $N$ returns the number of boxes needed to entirely cover the
trajectory. The grid refinement factor $N$ indicates that the domain of
computation is divided into a $[2^N\times2^N]$-grid of boxes.

The computational domain in which the box counting is to be performed is
explicitly stated by the coordinates of its lower-left and upper-right
corners, $(global_llx,global_lly)$ and $(global_urx,global_ury)$, respectively.
Parts of the trajectory which are outside of this domain are not included in
the box-counting.
If the entire input trajectory is to be used in the box counting, simply use
$(\PB{\\{global\_llx}},\PB{\\{global\_lly}})=({\rm min}[\PB{\\{x\_traj}}],{\rm
min}[\PB{\\{y\_traj}}])$ and
$(\PB{\\{global\_urx}},\PB{\\{global\_ury}})=({\rm max}[\PB{\\{x\_traj}}],{\rm
max}[\PB{\\{y\_traj}}])$ for the
specification of the computational domain.
\medskip\noindent
Input variables:
\varitem[\PB{\\{mm}}]{The total number of coordinates forming the input
trajectory,
or equivalently the length of the vectors \PB{${*}\\{x\_traj}$} and \PB{${*}%
\\{y\_traj}$}.}
\varitem[\PB{${*}\\{x\_traj}$}, \PB{${*}\\{y\_traj}$}]{Vectors of length \PB{%
\\{mm}} containing the $x$- and
$y$-coordinates of the input trajectory.}
\varitem[\PB{\\{resolution}}]{The grid refinement factor $N$.}
\varitem[\PB{\\{global\_llx}}, \PB{\\{global\_lly}}]{Coordinates of the
lower-left corner
of the computational domain in which the box-counting is to be performed.}
\varitem[\PB{\\{global\_urx}}, \PB{\\{global\_ury}}]{Coordinates of the
upper-right corner
of the computational domain in which the box-counting is to be performed.}
\medskip\noindent
Output variables:
\varitem[]{On exit, the routine returns the number of covering boxes
as an integer of \PB{\&{long} \&{unsigned}} precision.}
\medskip

\Y\B\4\X52:Simplified routine for calculation of number of boxes covering the
trajectory\X${}\E{}$\6
\&{long} \&{unsigned} \&{int} \\{get\_num\_covering\_boxes}(\&{double} ${}{*}%
\\{x\_traj},\39{}$\&{double} ${}{*}\\{y\_traj},\39{}$\&{long} \&{int} %
\\{mm}${},\39{}$\&{int} \|i${},\39{}$\&{double} \\{global\_llx}${},\39{}$%
\&{double} \\{global\_lly}${},\39{}$\&{double} \\{global\_urx}${},\39{}$%
\&{double} \\{global\_ury})\1\1\2\2\6
${}\{{}$\1\6
\&{return} ${}(\\{get\_num\_covering\_boxes\_with\_boxmaps}(\\{x\_traj},\39\\{y%
\_traj},\39\\{mm},\39\|i,\39\\{global\_llx},\39\\{global\_lly},\39\\{global%
\_urx},\39\\{global\_ury},\39\.{""},\39\T{0.0},\39\T{0.0},\39\.{""}));{}$\6
\4${}\}{}$\2\par
\U41.\fi

\N{1}{53}Index.
\fi

\inx
\fin
\con
