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