% 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 @*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 @*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)} } @*Revision history of the program. \medskip \citem[2006-05-08]{[v.1.0]} {\tt }\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 |get_num_covering_boxes_with_boxmaps()| and its more low-level core engine |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 }\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 |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 }\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 }\hfill\break Added documentation on the |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 }\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 }\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 }\hfill\break Updated \.{Makefile}:s for the generation of figures. Also corrected a rather stupid way of removing preceeding paths of file names. @*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 \# \# 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. @*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.} @*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 {\char'174} Last modified: May 8, 2006, Fredrik Jonsson {\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 ', where 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, =6 is used as default. {\char'174} The generated trajectory is written to standard output. {\char'174} Copyright (C) 2006 Fredrik Jonsson =========================================================================*/ \#include \#include 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 \# 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} @*The main program. Here follows the general outline of the main program. @c @@; @@; @@; @@; int main(int argc, char *argv[]) { @@; @@; @@; @@; @@; @@; @@; @@; @@; @@; @@; @@; return(SUCCESS); } @ 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 |fprintf|.} \varitem[{\tt stdlib.h}]{For access to the |exit| function.} \varitem[{\tt string.h}]{For string manipulation, |strcpy|, |strcmp| etc.} \varitem[{\tt ctype.h}]{For access to the |isalnum| function.} @= #include #include #include #include #include #include @ 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 |output_mode|, used to determine if output should append to an existing file.} \varitem[{\tt OVERWRITE}]{Code for the flag |output_mode|, used to determine if output should overwrite an existing file.} @= #define VERSION "1.6" #define COPYRIGHT "Copyright (C) 2006-2011, Fredrik Jonsson" #define SUCCESS (0) #define FAILURE (1) #define NCHMAX (256) #define APPEND 1 #define OVERWRITE 2 @ Declaration of global variables. The only global variables allowed in my programs are |optarg|, which is a pointer to the the string of characters that specified the call from the command line, and |progname|, which simply is a pointer to the string containing the name of the program, as it was invoked from the command line. @= extern char *optarg; char *progname; @ Declaration of local variables of the |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 |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[|verbose|]{Boolean flag which, if being nonzero, tells the program to display information at terminal output during program execution.} \varitem[|user_set_compwin|]{Boolean flag which indicates whether the user has explicitly stated a window of computation or not.} \varitem[|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 $|output_mode|=|APPEND|$), or if the data should overwrite already existing data in the file (for $|output_mode|=|OVERWRITE|$).} \varitem[|make_boxmaps|]{If nonzero, then graphs showing the distribution of covering boxes will be generated.} \varitem[|*num_boxes|]{Pointer to array holding the number of boxes at various depths of division.} \varitem[|initime|]{The time at which the \boxcount\ program was initialized.} \varitem[|now|]{Dummy variable for extraction of current time from the system clock.} \varitem[|mm|]{The number of points $M$ in the input trajectory. This is the length of the vectors |x_traj| and |y_traj| as described below.} \varitem[|nn|]{Gives the number of boxes in the $x$- or $y$-directions as $2^N$.} \varitem[|nnmax|]{The maximum refinement depth $N_{\rm max}$ of $N$.} \varitem[|nnmin|]{The minimum refinement depth $N_{\rm min}$ of $N$.} \varitem[|global_llx|, |global_lly|]{Lower-left coordinates of global window of computation.} \varitem[|global_urx|, |global_ury|]{Upper-right coordinates of global window of computation.} \varitem[|*x_traj|,|*y_traj|]{Vectors containing the $x$- and $y$-values of the coordinates along the input trajectory.} \varitem[|*x|,|*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[|*trajectory_file|]{Input file pointer, for reading the trajectory whose fractal dimension is to be estimated.} \varitem[|*frac_estimate_file|]{Output file pointer, for saving the estimated fractal dimension of the input trajectory.} \varitem[|*boxmap_file|]{Output file pointer, for saving maps of the spatial locations of the boxes covering the trajectory.} \varitem[|boxgraph_width|]{The physical width in millimetres of the generated \MP\ boxmap graphs.} \varitem[|boxgraph_height|]{The physical height in millimetres of the generated \MP\ boxmap graphs.} \varitem[|trajectory_filename|]{String for keeping the name of the file containing the input trajectory.} \varitem[|output_filename|]{String for keeping the base name of the set of output files.} \varitem[|frac_estimate_filename|]{String keeping the name of the file to which the estimate of the fractal dimension will be saved.} \varitem[|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[|no_arg|]{Variable for extracting the number of input arguments present on the command line as the program is executed.} \varitem[|i|]{Dummy variable used in loops when reading the input trajectory.} \varitem[|*fracdimen_estimates|]{Vector keeping the values characterizing estimated fractal dimension for various orders of $N$.} \varitem[|ave,adev,sdev|]{The average, average deviation, and standard deviation of the estimated fractal dimension for various orders of refinement $N$, as stored in |fracdimen_estimates|.} \varitem[|var,skew,curt|]{The variance, skewness, and kurtosis of the estimated fractal dimension for various orders of refinement $N$, as stored in |fracdimen_estimates|.} \medskip \noindent Generally in this program, the maximum number of characters a file name string can contain is |NCHMAX|, as defined in the definitions section of the program. @= short verbose,user_set_compwin,output_mode,make_boxmaps; long int *num_boxes; time_t initime; time_t now; long mm,nn,nnmin,nnmax; double global_llx,global_lly,global_urx,global_ury; double *x_traj,*y_traj,*x,*y; FILE *trajectory_file,*frac_estimate_file,*boxmap_file; char trajectory_filename[NCHMAX],output_filename[NCHMAX],frac_estimate_filename[NCHMAX]; char boxmap_filename[NCHMAX]; double boxgraph_width,boxgraph_height; int no_arg; long i; double *fracdimen_estimates,ave,adev,sdev,var,skew,curt; @ Initialization of variables. @= verbose=0; /* Verbose mode is off by default */ user_set_compwin=0; /* Before the command-line is parsed, nothing is known of these settings */ output_mode=OVERWRITE; /* default mode is to overwrite existing files */ make_boxmaps=0; /* Default is to not create graphs of the box distributions */ nnmin=0; /* Default value for $N_{\rm min}$ */ nnmax=10; /* Default value for $N_{\rm max}$ */ strcpy(output_filename,"out"); /* Default output file basename. */ strcpy(trajectory_filename,""); /* To indicate that no filename has been set yet. */ boxgraph_width=80.0; /* Default graph width in millimetres */ boxgraph_height=56.0; /* Default graph height in millimetres */ trajectory_file=NULL; frac_estimate_file=NULL; boxmap_file=NULL; initime=time(NULL); /* Time of initialization of the program. */ @ 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. @= { progname=strip_away_path(argv[0]); fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT); no_arg=argc; while(--argc) { if(!strcmp(argv[no_arg-argc],"-o") || !strcmp(argv[no_arg-argc],"--outputfile")) { --argc; if(!strcmp(argv[no_arg-argc],"append") || !strcmp(argv[no_arg-argc],"a")) { output_mode=APPEND; } else if(!strcmp(argv[no_arg-argc],"overwrite") || !strcmp(argv[no_arg-argc],"o")) { output_mode=OVERWRITE; } else { fprintf(stderr,"%s: Error in '-o' or '--outputfile' option!", progname); exit(FAILURE); } --argc; strcpy(output_filename,argv[no_arg-argc]); } else if(!strcmp(argv[no_arg-argc],"-w") || !strcmp(argv[no_arg-argc],"--computationwindow")) { user_set_compwin=1; --argc; if(!strcmp(argv[no_arg-argc],"llx")) { --argc; if(!sscanf(argv[no_arg-argc],"%lf",&global_llx)) { fprintf(stderr,"%s: Error in parsing lower-left x-value.\n", progname); exit(FAILURE); } } else { fprintf(stderr,"%s: Error in computation window option\n", progname); fprintf(stderr,"%s: Expecting 'llx'\n",progname); exit(FAILURE); } --argc; if(!strcmp(argv[no_arg-argc],"lly")) { --argc; if(!sscanf(argv[no_arg-argc],"%lf",&global_lly)) { fprintf(stderr,"%s: Error in parsing lower-left y-value.\n", progname); exit(FAILURE); } } else { fprintf(stderr,"%s: Error in computation window option\n", progname); fprintf(stderr,"%s: Expecting 'lly'\n",progname); exit(FAILURE); } --argc; if(!strcmp(argv[no_arg-argc],"urx")) { --argc; if(!sscanf(argv[no_arg-argc],"%lf",&global_urx)) { fprintf(stderr,"%s: Error in parsing lower-left x-value.\n", progname); exit(FAILURE); } } else { fprintf(stderr,"%s: Error in computation window option\n", progname); fprintf(stderr,"%s: Expecting 'urx'\n",progname); exit(FAILURE); } --argc; if(!strcmp(argv[no_arg-argc],"ury")) { --argc; if(!sscanf(argv[no_arg-argc],"%lf",&global_ury)) { fprintf(stderr,"%s: Error in parsing lower-left y-value.\n", progname); exit(FAILURE); } } else { fprintf(stderr,"%s: Error in computation window option\n", progname); fprintf(stderr,"%s: Expecting 'ury'\n",progname); exit(FAILURE); } } else if(!strcmp(argv[no_arg-argc],"-i") || !strcmp(argv[no_arg-argc],"--trajectoryfile")) { --argc; strcpy(trajectory_filename,argv[no_arg-argc]); } else if(!strcmp(argv[no_arg-argc],"-v") || !strcmp(argv[no_arg-argc],"--verbose")) { verbose=(verbose?0:1); } else if(!strcmp(argv[no_arg-argc],"-h") || !strcmp(argv[no_arg-argc],"--help")) { showsomehelp(); exit(SUCCESS); } else if(!strcmp(argv[no_arg-argc],"-m") || !strcmp(argv[no_arg-argc],"--boxmaps")) { make_boxmaps=1; } else if(!strcmp(argv[no_arg-argc],"--graphsize")) { --argc; if(!sscanf(argv[no_arg-argc],"%lf",&boxgraph_width)) { fprintf(stderr,"%s: Error in width of '--graphsize' option.\n", progname); exit(FAILURE); } --argc; if(!sscanf(argv[no_arg-argc],"%lf",&boxgraph_height)) { fprintf(stderr,"%s: Error in height of '--graphsize' option.\n", progname); exit(FAILURE); } } else if(!strcmp(argv[no_arg-argc],"--minlevel") || !strcmp(argv[no_arg-argc],"-Nmin")) { --argc; if(!sscanf(argv[no_arg-argc],"%ld",&nnmin)) { fprintf(stderr,"%s: Error in '--minlevel' or '-Nmin' option.\n", progname); exit(FAILURE); } } else if(!strcmp(argv[no_arg-argc],"--maxlevel") || !strcmp(argv[no_arg-argc],"-Nmax")) { --argc; if(!sscanf(argv[no_arg-argc],"%ld",&nnmax)) { fprintf(stderr,"%s: Error in '--maxlevel' or '-Nmax' option.\n", progname); exit(FAILURE); } } else { fprintf(stderr,"%s: Specified option '%s' invalid!\n", progname,argv[no_arg-argc]); showsomehelp(); exit(FAILURE); } } } @ Display starting time of program execution. @= { fprintf(stdout,"%s: Program execution started %s",progname,ctime(&initime)); } @ 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 |num_coordinate_pairs()|. @={ if (!strcmp(trajectory_filename,"")) { fprintf(stderr,"%s: No input trajectory specified!\n",progname); fprintf(stderr,"%s: Please use the '--trajectoryfile' option.\n", progname); fprintf(stderr,"%s: Use '--help' option to display a help message.\n", progname); exit(FAILURE); } if ((trajectory_file=fopen(trajectory_filename,"r"))==NULL) { fprintf(stderr,"%s: Could not open %s for loading trajectory!\n", progname,trajectory_filename); exit(FAILURE); } mm=num_coordinate_pairs(trajectory_file); x_traj=dvector(1,mm); y_traj=dvector(1,mm); for (i=1;i<=mm;i++) { fscanf(trajectory_file,"%lf",&x_traj[i]); /* scan $x$-coordinate */ fscanf(trajectory_file,"%lf",&y_traj[i]); /* scan $y$-coordinate */ } fclose(trajectory_file); if (verbose) { fprintf(stdout,"%s: Loaded %ld coordinate pairs from file '%s'.\n", progname,mm,trajectory_filename); } } @ Opening files for output by the program. @={ if (!strcmp(output_filename,"")) { fprintf(stderr,"%s: No output base name specified!\n",progname); fprintf(stderr,"%s: Please use the '--outputfile' option.\n", progname); exit(FAILURE); } sprintf(frac_estimate_filename,"%s.dat",output_filename); if (output_mode==APPEND) { if ((frac_estimate_file=fopen(frac_estimate_filename,"a"))==NULL) { fprintf(stderr,"%s: Could not open %s for output!\n", progname,frac_estimate_filename); exit(FAILURE); } fseek(frac_estimate_file,0L,SEEK_END); /* set file pointer to the end of the file */ } else if (output_mode==OVERWRITE) { if ((frac_estimate_file=fopen(frac_estimate_filename,"w"))==NULL) { fprintf(stderr,"%s: Could not open %s for loading trajectory!\n", progname,frac_estimate_filename); exit(FAILURE); } fseek(frac_estimate_file,0L,SEEK_SET); /* set file pointer to the beginning of the file */ } else { fprintf(stderr,"%s: Error: Output mode (output_mode) undefined!",progname); exit(FAILURE); } } @ 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 $(|global_llx|,|global_lly|)$ and $(|global_urx|,|global_ury|)$. @= { if (!user_set_compwin) { global_llx=x_traj[1]; global_lly=y_traj[1]; global_urx=global_llx; global_ury=global_lly; for (i=1;i<=mm;i++) { if (x_traj[i]>global_urx) global_urx=x_traj[i]; if (y_traj[i]>global_ury) global_ury=y_traj[i]; if (x_traj[i]={ num_boxes=livector(1,nnmax); nn=1; for (i=1;i<=nnmin-1;i++) nn=2*nn; /* This leaves |nn| as $2^{(N_{\rm min}-1)}$ */ for (i=nnmin;i<=nnmax;i++) { /* For all levels in refinement of grid resolution */ nn=2*nn; if (make_boxmaps) { /* do we wish to generate \METAPOST\ graphs of the box distribution? */ sprintf(boxmap_filename,"%s-%02ld.mp",output_filename,i); num_boxes[i]=get_num_covering_boxes_with_boxmaps(x_traj,y_traj,mm,i, global_llx,global_lly,global_urx,global_ury,boxmap_filename, boxgraph_width,boxgraph_height,trajectory_filename); } else { /* if not, just use the regular number-crunching routine */ num_boxes[i]=get_num_covering_boxes(x_traj,y_traj,mm,i, global_llx,global_lly,global_urx,global_ury); } if (verbose) { fprintf(stdout,"%s: N=%ld (%ldx%ld-grid of boxes): ",progname,i,nn,nn); fprintf(stdout,"Trajectory covered by %ld boxes\n",num_boxes[i]); } } } @ 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 $(|global_urx|-|global_llx|)/2^n$ and height $(|global_ury|-|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(|num_boxes|(n))$ as ordinata. Also in this block, we analyze whether the option |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. $$ @={ x=dvector(nnmin,nnmax); y=dvector(nnmin,nnmax); fracdimen_estimates=dvector(nnmin,nnmax); nn=1; for (i=1;i<=nnmax;i++) { nn=2*nn; if (nnmin<=i) x[i]=log((double)nn); } for (i=nnmin;i<=nnmax;i++) y[i]=log(num_boxes[i]); for (i=nnmin;i<=nnmax;i++) fracdimen_estimates[i]=y[i]/x[i]; if (verbose) { for (i=nnmin;i<=nnmax;i++) { fprintf(stdout,"%s: N=%ld, fractal dimension estimate=%f\n",progname, i,fracdimen_estimates[i]); } } moment(fracdimen_estimates,nnmin,nnmax,&ave,&adev,&sdev,&var,&skew,&curt); if (verbose) { fprintf(stdout,"%s: Estimate of fractal dimension: %f\n",progname,ave); fprintf(stdout,"%s: Average deviation: %f\n",progname,adev); fprintf(stdout,"%s: Standard deviation: %f\n",progname,sdev); fprintf(stdout,"%s: Skewness: %f\n",progname,skew); } free_livector(num_boxes,1,nnmax); /* release the memory occupied by |num_boxes| */ free_dvector(fracdimen_estimates,nnmin,nnmax); free_dvector(x,nnmin,nnmax); /* release the memory occupied by |x| */ free_dvector(y,nnmin,nnmax); /* release the memory occupied by |y| */ } @ Save the fractal dimension to file. @={ if (output_mode==APPEND) { fseek(frac_estimate_file,0L,SEEK_END); } else if (output_mode==OVERWRITE) { fseek(frac_estimate_file,0L,SEEK_SET); } fprintf(frac_estimate_file,"%f %f %f\n",ave,sdev,skew); } @ Close any open output files. @={ fclose(frac_estimate_file); } @ Display elapsed execution time. @={ now=time(NULL); if (verbose) fprintf(stdout,"%s: Total execution time: %d s\n",progname, ((int)difftime(now,initime))); fprintf(stdout,"%s: Program execution closed %s",progname,ctime(&now)); } @*Subroutines. In this section, all subroutines as used by the main program are listed. @= @@; @@; @@; @@; @@; @@; @@; @@; @ Routine for computation of average, average deviation and standard deviation. This routine is adopted from {\sl Numerical Recipes in C}. @= void moment(double data[],int nnmin,int nnmax,double *ave,double *adev, double *sdev,double *var,double *skew,double *curt) { int j; double ep=0.0,s,p; if (nnmax-nnmin <= 1) { fprintf(stderr,"%s: Error in routine moment()! ",progname); fprintf(stderr,"(nnmax-nnmin>1 is a requirement)\n"); exit(FAILURE); } s=0.0; for (j=nnmin;j<=nnmax;j++) s += data[j]; *ave=s/(nnmax-nnmin+1); *adev=(*var)=(*skew)=(*curt)=0.0; for (j=nnmin;j<=nnmax;j++) { *adev += fabs(s=data[j]-(*ave)); ep += s; ep += s; *var += (p=s*s); *skew += (p *= s); *curt += (p *= s); } *adev /= (nnmax-nnmin+1); *var=(*var-ep*ep/(nnmax-nnmin+1))/(nnmax-nnmin); *sdev=sqrt(*var); if (*var) { *skew /= ((nnmax-nnmin+1)*(*var)*(*sdev)); *curt=(*curt)/((nnmax-nnmin+1)*(*var)*(*var))-3.0; } else { fprintf(stderr,"%s: Error in routine moment()! ",progname); fprintf(stderr,"No skew/kurtosis for zero variance\n"); exit(FAILURE); } } @ 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 |x_traj[1..M]| and |y_traj[1..M]|, containing the coordinates $(x_m,y_m)$ of the trajectory. @= long int num_coordinate_pairs(FILE *trajectory_file) { double tmp; int tmpch; long int mm=0; fseek(trajectory_file,0L,SEEK_SET); /* rewind file to beginning */ while((tmpch=getc(trajectory_file))!=EOF) { ungetc(tmpch,trajectory_file); fscanf(trajectory_file,"%lf",&tmp); /* Read away the $x$ coordinate */ fscanf(trajectory_file,"%lf",&tmp); /* Read away the $y$ coordinate */ mm++; tmpch=getc(trajectory_file); /* Read away blanks and linefeeds */ while ((tmpch!=EOF)&&(!isdigit(tmpch))) tmpch=getc(trajectory_file); if (tmpch!=EOF) ungetc(tmpch,trajectory_file); } fseek(trajectory_file,0L,SEEK_SET); /* rewind file to beginning */ return(mm); } @ 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 |strip_away_path()| routine is typically called when initializing the program name string |progname| from the command line string |argv[0]|, and is typically located in the blocks related to parsing of the command line options. @= @@; @@; @ Checking for a valid path character. The |pathcharacter| routine takes one character |ch| as argument, and returns 1 (``true'') if the character is valid character of a path string, otherwise 0 (``false'') is returned. @= short pathcharacter(int ch) { return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')|| (ch=='-')||(ch=='+')); } @ Routine for stripping away path string of a file name. The |strip_away_path()| routine takes a character string |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. @= char *strip_away_path(char filename[]) { int j,k=0; while (pathcharacter(filename[k])) k++; j=(--k); /* this is the uppermost index of the full path+file string */ while (isalnum((int)(filename[j]))) j--; j++; /* this is the lowermost index of the stripped file name */ return(&filename[j]); } @ 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. @= @@; @@; @@; @@; @@; @@; @ The |dvector| routine allocates a real-valued vector of double precision, with vector index ranging from |nl| to |nh|. @= double *dvector(long nl, long nh) { double *v; v=(double *)malloc((size_t) ((nh-nl+2)*sizeof(double))); if (!v) { fprintf(stderr,"Error: Allocation failure in dvector()\n"); exit(FAILURE); } return v-nl+1; } @ The |free_dvector| routine release the memory occupied by the real-valued vector |v[nl..nh]|. @= void free_dvector(double *v, long nl, long nh) { free((char*) (v+nl-1)); } @ The |ivector| routine allocates a real-valued vector of integer precision, with vector index ranging from |nl| to |nh|. @= int *ivector(long nl, long nh) { int *v; v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int))); if (!v) { fprintf(stderr,"Error: Allocation failure in ivector()\n"); exit(FAILURE); } return v-nl+1; } @ The |free_ivector| routine release the memory occupied by the real-valued vector |v[nl..nh]|. @= void free_ivector(int *v, long nl, long nh) { free((char*) (v+nl-1)); } @ The |livector| routine allocates a real-valued vector of long integer precision, with vector index ranging from |nl| to |nh|. @= long int *livector(long nl, long nh) { long int *v; v=(long int *)malloc((size_t) ((nh-nl+2)*sizeof(long int))); if (!v) { fprintf(stderr,"Error: Allocation failure in livector()\n"); exit(FAILURE); } return v-nl+1; } @ The |free_livector| routine release the memory occupied by the real-valued vector |v[nl..nh]|. @= void free_livector(long int *v, long nl, long nh) { free((char*) (v+nl-1)); } @ The |simatrix| routine allocates an array of short integer precision, with array row index ranging from |nrl| to |nrh| and column index ranging from |ncl| to |nch|. @= short int **simatrix(long nrl, long nrh, long ncl, long nch) { long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; short int **m; m=(short int **) malloc((size_t)((nrow+1)*sizeof(short int*))); if (!m) { fprintf(stderr,"%s: Allocation failure 1 in simatrix()\n",progname); exit(FAILURE); } m += 1; m -= nrl; m[nrl]=(short int *) malloc((size_t)((nrow*ncol+1)*sizeof(short int))); if (!m[nrl]) { fprintf(stderr,"%s: Allocation failure 2 in simatrix()\n",progname); exit(FAILURE); } m[nrl] += 1; m[nrl] -= ncl; for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; return m; } @ The |free_simatrix| routine releases the memory occupied by the short integer matrix |v[nl..nh]|, as allocated by |simatrix()|. @= void free_simatrix(short int **m,long nrl,long nrh,long ncl,long nch) { free((char*) (m[nrl]+ncl-1)); free((char*) (m+nrl-1)); } @ Routine for displaying help message to standard terminal output. @= void showsomehelp(void) { fprintf(stderr,"Usage: %s [options]\n",progname); fprintf(stderr,"Options:\n"); fprintf(stderr," -h, --help\n"); fprintf(stderr," Display this help message and exit clean.\n"); fprintf(stderr," -v, --verbose\n"); fprintf(stderr," Toggle verbose mode on/off.\n"); fprintf(stderr," -o, --outputfile \n"); fprintf(stderr," Specifies the base name of the output files where the program\n"); fprintf(stderr," is to save the calculated data. If the --outputfile or -o\n"); fprintf(stderr," option is followed by 'append' the estimate for the fractal\n"); fprintf(stderr," dimension will be appended to the file named .dat, which\n"); fprintf(stderr," will be created if it does not exist. If the following word\n"); fprintf(stderr," instead is `overwrite' the file will instead be overwritten.\n"); fprintf(stderr," -i, --trajectoryfile\n"); fprintf(stderr," Specifies the input trajectory of the graph whose fractal\n"); fprintf(stderr," dimension is to be estimated.\n"); fprintf(stderr," -w, --computationwindow llx lly urx ury \n"); fprintf(stderr," This option explicitly specifies the domain over which the\n"); fprintf(stderr," box-counting algorithm will be applied, in terms of the\n"); fprintf(stderr," lower-left and upper-right corners (llx,lly) and (urx,ury),\n"); fprintf(stderr," respectively. By specifying this option, any automatic\n"); fprintf(stderr," calculation of computational window will be neglected.\n"); fprintf(stderr," -m, --boxmaps\n"); fprintf(stderr," If this option is present, the program will generate\n"); fprintf(stderr," MetaPost code for maps of the distribution of boxes.\n"); fprintf(stderr," In doing so, also the input trajectory is included in\n"); fprintf(stderr," the graphs. The convention for the naming of the output\n"); fprintf(stderr," map files is that they are saved as ..dat,\n"); fprintf(stderr," where is the base filename as specified using the\n"); fprintf(stderr," -o or --outputfile option, is the automatically appended\n"); fprintf(stderr," current level of resolution refinement in the box-counting,\n"); fprintf(stderr," and where '.dat' is the file suffix as automatically appended\n"); fprintf(stderr," by the program.\n"); fprintf(stderr," --graphsize \n"); fprintf(stderr," If the -m or --boxmaps option is present at the command line,\n"); fprintf(stderr," then the --graphsize option will override the default graph\n"); fprintf(stderr," size of the generated boxmaps. (Default graph size is 80 mm\n"); fprintf(stderr," width and 56 mm height.)\n"); fprintf(stderr," -Nmin, --minlevel \n"); fprintf(stderr," Specifies the minimum level Nmin of grid refinement that \n"); fprintf(stderr," will be used in the evaluation of the estimate of the fractal\n"); fprintf(stderr," dimension. With this option specified, the coarsest level\n"); fprintf(stderr," used in the box-counting will be a [(2^Nmin)x(2^Nmin)]-grid\n"); fprintf(stderr," of boxes.\n"); fprintf(stderr," -Nmax, --maxlevel \n"); fprintf(stderr," This option is similar to the --minlevel or -Nmin options,\n"); fprintf(stderr," with the difference that it instead specifies the maximum\n"); fprintf(stderr," level Nmax of grid refinement that will be used in the\n"); fprintf(stderr," evaluation of the estimate of the fractal dimension.\n"); } @ The |lines_intersect(p1x,p1y,q1x,q1y,p2x,p2y,q2x,q2y)| routine takes the start and end points of two line segments, the first between $(|p1x|,|p1y|)$ and $(|q1x|,|q1y|)$ and the second between $(|p2x|,|p2y|)$ and $(|q2x|,|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$. @= short int lines_intersect(double p1x,double p1y,double q1x,double q1y, double p2x,double p2y,double q2x,double q2y) { double a,b,c,d,e,f,det,tmp1,tmp2; short int intersect; a=q1x-p1x; b=p2x-q2x; c=q1y-p1y; d=p2y-q2y; e=p2x-p1x; f=p2y-p1y; det=a*d-b*c; tmp1=e*d-b*f; tmp2=a*f-e*c; intersect=0; if (det>0) { if (((0.0<=tmp1)&&(tmp1<=det))&&((0.0<=tmp2)&&(tmp2<=det))) intersect=1; } else if (det<0) { if (((det<=tmp1)&&(tmp1<=0.0))&&((det<=tmp2)&&(tmp2<=0.0))) intersect=1; } return(intersect); } @ Routine for determining whether a line and a box intersect or not. The |box_intersection()| routine simply divides the input box, being characterized by its lower-left and upper-right corners $(|llx|,|lly|)$ and $(|urx|,|ury|)$, into the four line segments corresponding to its four edges, followed by calling the routine |lines_intersect()| to determine if any of these edges intersect the line segment. If an intersection is found, the |box_intersection()| routine returns 1 (true), otherwise 0 (false). \medskip\noindent Input variables: \varitem[|px|, |py|]{The coordinates of the first end point ${\bf p}=(p_x,p_y)$ of the line segment.} \varitem[|qx|, |qy|]{The coordinates of the second end point ${\bf q}=(q_x,q_y)$ of the line segment.} \varitem[|llx|, |lly|]{Coordinates of the lower-left corner of the box.} \varitem[|urx|, |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 |short| precision.} \medskip @= short int box_intersection(double px,double py,double qx,double qy, double llx,double lly,double urx,double ury) { if (lines_intersect(px,py,qx,qy,llx,lly,urx,lly)) return(1); /* intersection with bottom edge */ if (lines_intersect(px,py,qx,qy,urx,lly,urx,ury)) return(1); /* intersection with right edge */ if (lines_intersect(px,py,qx,qy,urx,ury,llx,ury)) return(1); /* intersection with top edge */ if (lines_intersect(px,py,qx,qy,llx,ury,llx,lly)) return(1); /* intersection with left edge */ return(0); /* this happens only if no edge is intersecting the line segment */ } @ 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, |get_num_covering_boxes()| simply performs this task without caring of documenting the box distributions as graphs, or ``box maps'', while the second one, |get_num_covering_boxes_with_boxmaps()| also includes the generation of these maps, with output in terms of \METAPOST\ code. @= @@; @@; @ Routine for calculation of number of boxes covering the trajectory, also including the generation of documenting graphs of the box distributions. The |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 $(|global_llx|,|global_lly|)=({\rm min}[|x_traj|],{\rm min}[|y_traj|])$ and $(|global_urx|,|global_ury|)=({\rm max}[|x_traj|],{\rm max}[|y_traj|])$ for the specification of the computational domain. \medskip\noindent Input variables: \varitem[|mm|]{The total number of coordinates forming the input trajectory, or equivalently the length of the vectors |*x_traj| and |*y_traj|.} \varitem[|*x_traj|, |*y_traj|]{Vectors of length |mm| containing the $x$- and $y$-coordinates of the input trajectory.} \varitem[|resolution|]{The grid refinement factor $N$.} \varitem[|global_llx|, |global_lly|]{Coordinates of the lower-left corner of the computational domain in which the box-counting is to be performed.} \varitem[|global_urx|, |global_ury|]{Coordinates of the upper-right corner of the computational domain in which the box-counting is to be performed.} \varitem[|trajectory_filename|]{String containing the name of the file containing the input trajectory.} \varitem[|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 |long unsigned| precision.} \medskip\noindent Internal variables: \varitem[|px|,|py|,|qx|,|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 @= long unsigned int get_num_covering_boxes_with_boxmaps( double *x_traj,double *y_traj, long int mm,int resolution,double global_llx,double global_lly, double global_urx,double global_ury,char boxgraph_filename[], double boxgraph_width,double boxgraph_height, char trajectory_filename[]) { short int **box; long unsigned int i,j,m,nn,istart,jstart,istop,jstop,iincr,jincr,num_boxes; double *x_box,*y_box; /* Keeps track of the lower-left coordinates of the boxes. */ double px,py,qx,qy; FILE *boxgraph_file; @@; @@; for (m=1;m<=mm-1;m++) { @@; @@; } @@; @@; @@; @@; @@; return(num_boxes); } @ 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 $|x_box|[1\ldots 2^N]$ and $|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 $|x_box|[2^N]$ and $|y_box|[2^N]$ containing the upper-right corner coordinates of the upper-right-most box of the global window. @= { nn=1; for (i=1;i<=resolution;i++) nn=2*nn; box=simatrix(1,nn,1,nn); for (i=1;i<=nn;i++) for (j=1;j<=nn;j++) box[i][j]=0; x_box=dvector(1,nn+1); y_box=dvector(1,nn+1); for (m=1;m<=nn+1;m++) { x_box[m]=global_llx+((double)(m-1))*(global_urx-global_llx)/((double)(nn)); y_box[m]=global_lly+((double)(m-1))*(global_ury-global_lly)/((double)(nn)); } } @ Find indices $(i_a,j_a)$ of the box covering the first coordinate of the trajectory. @= { istart=0; jstart=0; for (m=1;m<=nn;m++) { if ((x_box[m]<=x_traj[1])&&(x_traj[1]<=x_box[m+1])) istart=m; if ((y_box[m]<=y_traj[1])&&(y_traj[1]<=y_box[m+1])) jstart=m; } if ((istart==0)||(jstart==0)) { fprintf(stderr,"%s: Error! Cannot find box indices of 1st coordinate!\n", progname); fprintf(stderr,"%s: Please check data for input trajetory.\n",progname); exit(FAILURE); } } @ Find indices $(i_b,j_b)$ of the box covering the end point of the $m$th trajectory segment. @={ px=x_traj[m]; py=y_traj[m]; qx=x_traj[m+1]; qy=y_traj[m+1]; istop=istart; jstop=jstart; if (pxqx) istop--; } if (pyqy) jstop--; } if (0==1) { fprintf(stdout,"%s: Endpoint box indices: (i=%ld,j=%ld)\n", progname,istop,jstop); } } @ 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. @= { for (i=istart;i!=(istop+iincr);i+=iincr) { for (j=jstart;j!=(jstop+jincr);j+=jincr) { if (box_intersection(px,py,qx,qy, x_box[i],y_box[j],x_box[i+1],y_box[j+1])) { box[i][j]=1; } } } istart=istop; jstart=jstop; } @ 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. @={ if (!strcmp(boxgraph_filename,"")) { /* is |boxgraph_filename| an empty string? */ boxgraph_file=NULL; } else { if ((boxgraph_file=fopen(boxgraph_filename,"w"))==NULL) { fprintf(stderr,"%s: Could not open %s for box graphs!\n", progname,boxgraph_filename); exit(FAILURE); } fseek(boxgraph_file,0L,SEEK_SET); fprintf(boxgraph_file,"input graph;\n"); fprintf(boxgraph_file,"def box(expr i,j)=\n"); fprintf(boxgraph_file,"begingroup\n"); fprintf(boxgraph_file,"llx:=%2.8f+(i-1)*%2.8f;\n", global_llx,(global_urx-global_llx)/(nn)); fprintf(boxgraph_file,"lly:=%2.8f+(j-1)*%2.8f;\n", global_lly,(global_ury-global_lly)/(nn)); fprintf(boxgraph_file,"urx:=%2.8f+(i)*%2.8f;\n", global_llx,(global_urx-global_llx)/(nn)); fprintf(boxgraph_file,"ury:=%2.8f+(j)*%2.8f;\n", global_lly,(global_ury-global_lly)/(nn)); fprintf(boxgraph_file,"gdraw (llx,lly)--(urx,lly);\n"); fprintf(boxgraph_file,"gdraw (urx,lly)--(urx,ury);\n"); fprintf(boxgraph_file,"gdraw (urx,ury)--(llx,ury);\n"); fprintf(boxgraph_file,"gdraw (llx,ury)--(llx,lly);\n"); fprintf(boxgraph_file,"endgroup\n"); fprintf(boxgraph_file,"enddef;\n"); fprintf(boxgraph_file,"beginfig(1);\n"); fprintf(boxgraph_file,"w:=%2.4fmm; h:=%2.4fmm;\n", boxgraph_width,boxgraph_height); fprintf(boxgraph_file,"draw begingraph(w,h);\n"); fprintf(boxgraph_file,"pickup pencircle scaled .3pt;\n"); fprintf(boxgraph_file,"setrange(%2.6f,%2.6f,%2.6f,%2.6f);\n", global_llx,global_lly,global_urx,global_ury); } } @ 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 |boxgraph_file| is used to determine whether to write the trajectory to file or not. If |boxgraph_file| equals to |NULL| (NULL), then the \boxcount\ program will not attempt to write the input trajectory to file. @={ if (boxgraph_file!=NULL) { fprintf(boxgraph_file,"pickup pencircle scaled .5pt;\n"); i=0; for (m=1;m<=mm;m++) { if (i==0) { if (m==1) { fprintf(boxgraph_file,"gdraw (%2.4f,%2.4f)", x_traj[m],y_traj[m]); } else if (2= { num_boxes=0; for (i=1;i<=nn;i++) { for (j=1;j<=nn;j++) { if(box[i][j]==1) { num_boxes++; if (boxgraph_file!=NULL) { fprintf(boxgraph_file,"box(%ld,%ld);\n",i,j); } } } } free_simatrix(box,1,nn,1,nn); } @ 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. @={ if (boxgraph_file!=NULL) { fprintf(boxgraph_file,"autogrid(itick bot,itick lft);\n"); fprintf(boxgraph_file,"glabel.bot(btex $x$ etex,OUT);\n"); fprintf(boxgraph_file,"glabel.lft(btex $y$ etex,OUT);\n"); fprintf(boxgraph_file,"endgraph;\n"); fprintf(boxgraph_file,"endfig;\n"); fprintf(boxgraph_file,"end\n"); } } @ Close any open file for output of box distribution graph. @={ if (boxgraph_file!=NULL) { fprintf(stdout, "%s: MetaPost output box distribution graph written to %s.\n", progname,boxgraph_filename); fclose(boxgraph_file); } } @ Routine for calculation of number of boxes covering the trajectory. This routine provides a simplified interface to the general boxcountng routine |get_num_covering_boxes_with_boxmaps|, with the only difference that no graph of the box distribution over the trajectory is generated. The |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 $(|global_llx|,|global_lly|)=({\rm min}[|x_traj|],{\rm min}[|y_traj|])$ and $(|global_urx|,|global_ury|)=({\rm max}[|x_traj|],{\rm max}[|y_traj|])$ for the specification of the computational domain. \medskip\noindent Input variables: \varitem[|mm|]{The total number of coordinates forming the input trajectory, or equivalently the length of the vectors |*x_traj| and |*y_traj|.} \varitem[|*x_traj|, |*y_traj|]{Vectors of length |mm| containing the $x$- and $y$-coordinates of the input trajectory.} \varitem[|resolution|]{The grid refinement factor $N$.} \varitem[|global_llx|, |global_lly|]{Coordinates of the lower-left corner of the computational domain in which the box-counting is to be performed.} \varitem[|global_urx|, |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 |long unsigned| precision.} \medskip @= long unsigned int get_num_covering_boxes(double *x_traj,double *y_traj,long int mm, int i,double global_llx,double global_lly, double global_urx,double global_ury) { return(get_num_covering_boxes_with_boxmaps(x_traj,y_traj,mm,i, global_llx,global_lly,global_urx,global_ury,"",0.0,0.0,"")); } @*Index.