This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International (CC BY-SA 4.0) License. To view a copy of this license, visit https://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

Preface

With the dramatic increases in home computer performance over the past few years it is not surprising that chaos theory and fractal graphics programming have become extremely popular. People can now create, in their own homes, the programming experiments and images which only 10 years ago were leading edge scientific research. The Amiga is an ideal system for exploring this exciting field, as it can perform the necessary calculations very quickly and the output can be displayed using high resolution graphics.

This book details the methods used to produce most important fractals and also provides a wide variety of ideas for further experimentation, using many of the Amiga’s most powerful features including colour graphics and even sound. The concepts behind these images are discussed using only the most elementary maths and the minimum of technical jargon. Wider aspects of chaos, such as its history and notable applications in the real world, have also been included.

Fully annotated Amiga BASIC (compatible with Hisoft BASIC) and GFA BASIC example programs are given throughout to illustrate the concepts and to provide a base for experimentation. Such examples can even be appreciated by novice programmers who will find this an enjoyable way to learn more about Amiga programming techniques. Appendices provide useful BASIC routines and give hints for converting the programs to other Amiga languages.

The layout of this book is such that it may be used either as a practical introduction to chaos or as a reference text for the experienced reader, where chapters may be skipped or read out of sequence if necessary.

I would like to thank my parents, Ange and Mike, for their support and assistance during the writing of this book. Also thanks to Gemma for distribution of the support disk, and to Russell for his help with the initial experiments.

1. Introduction

Until a few centuries ago it was believed that nature was a totally unpredictable force, and that forecasting or altering natural events was impossible. Religions were born to explain away this unpredictability, and indeed the thought of a great being who arbitrarily decided what was going to happen in the universe fits in well, as it is a known fact that human actions are often unforeseeable.

During the 16th and 17th centuries scientists such as Galileo, Brahe and Kepler combined observational and academic skills in order to determine that the planets of the solar system had cyclic motions, and also devised methods of precisely calculating the future positions of planets. Newton later determined more general laws governing motion, allowing the results of forces and collisions on smaller bodies to be predicted. Naturally the theories of these scientists came into conflict with the religious movements, that still believed the world was flat, and that it lay at the centre of the solar system.

The Church put a lot of effort into silencing these scientists, but gradually their beliefs became widely accepted and society now appreciates that the results of many natural events can be predicted. By the 18th century scientists were becoming complacent, believing that there were very few things on Earth that they could not understand. However, at the beginning of the 20th century scientists were questioning the validity of some of the old theories and equations, and indeed many modifications were made. Even so, scientists were still not content with the accuracy of some of these predictions.

During World War Two the importance of accurate weather forecasts became acutely apparent, and in the following decades scientists were encouraged to become more involved in meteorology, the science of weather and climate. It was popularly hoped that with the new technology available, and with large research funds, accurate long-term weather forecasting would be possible. Despite well financed and well co-ordinated schemes in Europe and the United States this hope was never realised.

At the time people were surprised at the scientists' apparent failure, and could not understand why mundane weather forecasting was so elusive when something as distant and obscure as planetary motion had been predicted 300 years previously. A simple hypothesis for the lack of success is that the weather is of a much more complicated nature than those systems for which rules have been found, and it is influenced by a wide variety of outside events. This is a fair assumption for weather, but unpredictable processes do not have to be complicated. For example the flow of water in a pipe, the path of energetic smoke particles in a container and even the action of a simple pendulum bob are hard to forecast in some circumstances.

An answer closer to the truth is that unlike predictable processes, such as the ebb and flow of the tide, the day-to-day weather is not cyclic, it does not repeat to a constant rhythm. Cyclic processes are easier to forecast mainly because a simple model can be built around data collected from previous cycles. However, examples contradict this again, as useful rules have been found for forecasting the results of all manner of events, from collisions to chemical reactions. Naturally, if it was possible to accurately predict the effects of things like the weather, the paths of oil slicks and fluctuations of the world money markets many problems could be solved, and much effort saved. Currently, however, general trends are all that can be found in these important areas.

Scientists working on meteorology theory, such as Edward Lorenz, began research into why weather prediction was so difficult, using simple computer-based mathematical models. Other researchers looked at different unpredictable elements of nature. For example Robert May concentrated on population dynamics, and Benoit Mandelbrot discovered unpredictability in pure mathematics. As people like these published papers of their work connections began to appear between the different kinds of irregularity being investigated. Common characteristics could be seen to be exhibited in erratic processes taking place in all manner of different subject areas, from economics to chemistry. Processes displaying these certain characteristics became known as chaotic processes due to their pseudo-random nature, and the study of their similarities became the new science: chaos.

Chaos

For hundreds of years scientists had disregarded the intricacies of the non-cyclic, erratic side of nature, believing it to be too difficult to predict and managing instead with generalised behaviour patterns. With the birth of chaos theory it seemed as if a solution to these fundamental problems may at last be in sight. Surely if the cause of unpredictability could be found it would be possible to cure it. Some chaos crusaders even rated it as the third great physics discovery of the 20th century, after quantum mechanics and the theory of relativity.

The backbone of chaos theory is the notion that complex behaviour does not have to come from complex rules, thus inferring that a system as complicated as the weather may be summed up in one simple equation. The hope is that chaos will be able to give us an equation for everything. This does, however, seem far from likely at the present time as the science is still in its infancy. What a lot of people fail to understand is why it has taken so long for scientists to acknowledge the importance of chaos. Surely if it is so fundamental it could have been appreciated centuries ago.

The catalyst for the research has undoubtedly been the computer. The millions of calculations that must be performed in even the simplest of chaotic processes make manual calculation practically impossible. In the beginning it was only well-funded defence establishments and universities that could afford the equipment necessary to carry out such experiments. But since the early 80s the price and size of computers have fallen dramatically, while power and availability have increased, bringing chaos within easy reach of many other people. Never before have so many people been able to become actively involved in such a state-of-the-art science, and this has been a major factor in its increasing popularity.

Many books have been written on the subject at all levels and mathematical periodicals, not normally touched with a barge-pole if possible, strangely went missing from reference libraries across the world. This incredible level of interest has grown ever since, with some chaos books becoming instant best-sellers, and less than five years after the discovery of the Mandelbrot set, popular computing magazines were publishing polished programs to display it on home computers.

Fractals

Fractals are abstract shapes existing in two or more dimensions, whose incredible complexity prevents them from being treated like normal geometrical shapes. Fractals have surprisingly little ordered structure, despite the fact that they originate from very simple rules. Such images have been irreversibly linked to chaos theory because they are such a good example of how immensely complicated output can be created from a very simple process.

The typical example of a computer generated fractal is the Mandelbrot set, whose infinitely complex structure is uncharacteristic of the simple equation from which it originates. Like most things in chaos theory these computer images are not abstract from the real world, frost on a window and the structure of a tree are natural fractals. Fractals don’t always have to be static though. For example, a November 5th sparkler produces a constantly changing fractal structure of sparks.

This is a very generalised definition of a fractal, and a more complete one will be built up as the book progresses. Fractals are most successfully described with the aid of practical experiments and examples, such as those provided for the Amiga in later chapters.

The Amiga

In terms of power, price and software availability the Amiga is undoubtedly one of the best value computers around, now rivalling the less powerful Atari ST. The Amiga is particularly suitable for fractal work because it can perform the necessary calculations very quickly and the output can be displayed in up to 4096 colours, or up to a resolution of 640x512 pixels. The inclusion of the mouse, printer port and disk drive mean that chaos programs can be made very user-friendly, with easy storage for the final output.

The only problem with the Amiga is the language supplied with it, Amiga BASIC, lacks some features (e.g. recursion) which are essential for fractal work. However, because Amiga BASIC is readily available, easy to read and not unworkably slow it has been used for most of the examples in this book. In fact a listing can be assumed to be written in Amiga BASIC unless otherwise stated.

The decision regarding the best language on which to base other parts of the book was a more difficult one. Machine code was ruled out immediately because of its unstructured, often difficult to understand nature. The C language is widely used by chaos researchers due to its speed, mathematical capabilities and structure. However, versions of C for the Amiga generally cost about five times more than comparable BASICs, meaning that few Amiga owners have had a chance to become familiar with them. The only other well supported high-level language for the Amiga is BASIC, so the choice was narrowed down to the four main Amiga versions of the language, shown in Table 1.1.

Language Vendor Approx. price

AMOS

Mandarin

49.99 *

Blitz BASIC

Siren Software

69.99

GFA BASIC 3.5

GFA Data Media

50.00

Power BASIC

Hisoft

49.95 (price includes compiler)

* includes various utilities such as music and sprite editors

Due to its pricing and popularity GFA BASIC appeared to be the most suitable choice. Hisoft BASIC would also have been a good choice, but even though it is compiled it is still slower than interpreted GFA BASIC. AMOS and Blitz BASIC were not considered as they are geared more towards games than mathematical applications. Appendix C gives advice on converting the GFA BASIC and Amiga BASIC example programs to other dialogues of BASIC and to C.

In the GFA BASIC programs only version 3.0 syntax has been used, so the programs should work on both version 3.0 and 3.5 of the language without alteration. If you have the GFA compiler it is recommended that you use it once you are sure that the program has been typed in correctly (i.e. when it works with the interpreter). Note that the GFA BASIC compiler only gives a very slight increase in speed over the interpreted programs.

How to Use this Book

Chaos theory can appear daunting, but the use of simple programming procedures and notational conventions in this book should make it very accessible. The following few pages provide details of the conventions used and give other hints to help you get the most out of both the text and the listings.

Monitor Compatibility

The design of the Amiga allows any compatible monitor to display any of the Amiga’s main screen modes. However, the non-interlace high resolution mode is the best one suitable for all monitors, giving a resolution of 640x200 pixels in Amiga BASIC with up to 32 colours. This specification is perfectly acceptable for fractal work, so this mode is used throughout the book. Details on converting the programs to run in different screen modes are provided in Appendix D.

Format of Equations

Many of the fractals in this book are based on mathematical equations. Such equations are shown using programming syntax wherever possible, instead of normal mathematical notation, to provide a solid link to the example programs and allow non-mathematicians to get to grips with the concepts more readily. A summary of the most commonly used conventions is shown in Table 1.2. If you are an experienced programmer you will probably already be familiar with most of these, but the table should still prove useful when cross-referencing to other chaos literature.

Convention Description Convention

(in this book)

(normal maths)

+

+

-

minus/negative

-

*

multiply

x

/

divide

÷

1/2

a fraction less than 1

frac{1}{2}

23/4

a fraction larger than 1

frac{23}{4}

=

equal to

=

equivalent to

α

proportional to

α

less than or equal to

>=

greater than or equal to

<>

not equal to

2^2

a number raised to a power

2^2

SQR(16)

square root of a number

sqrt(16)

(1+2)

bracketed expression

(1+2)

((1+1)*2)

multi-bracketed expression

[(1+1)x2]

These operators can generally be said to apply to all the major Amiga programming languages. However, the function used to raise a number to a power varies considerably, for example in C the `pow()` function is usually used. It is important to note that the 'equivalent to' and 'proportional to' symbols are not implemented in most languages as their definitions are rather too flexible.

Although variable names will be kept consistent between text and program listings their appearance may vary. This is because some variables will be shown in the text with subscripts in order to differentiate between variables with similar characteristics, for example pnew and pold. Because BASIC doesn’t allow subscripts the variables will be shown in the listings as `pnew` and `pold`.

Program Listings

Great effort has gone into keeping the program listings as short as possible. This has been done to increase clarity, to reduce the possibility of typing errors, and to allow output similar to that shown in the figures to be created with as little effort as possible. The few points discussed below are intended to make entering the programs as painless as possible.

NB: A program can be assumed to be written in Amiga BASIC unless otherwise stated.

Line Length

Lines which are longer than the physical width of the page are wrapped around and right justified. When typing in such lines you should ignore the wrap-around, and simply carry on typing the line to the end (this may require typing past the right hand edge of the editing window if the line is over a certain length).

Remarks

Remarks are used, along with self-explanatory procedure and variable names, to help explain the actions of particular lines or groups of lines. These are always preceded by either REM, an apostrophe (') or, in GFA BASIC, an exclamation mark (!). The ! type of remarks are the only type which GFA BASIC allows at the end of a program line. These will be used in a fashion similar to:

``  COLOR 1	!Change colour``

In Amiga BASIC an apostrophe would be used instead, for example:

``  COLOR 1	'Change colour``

The remarks are only included to aid understanding and may be omitted if this is desirable.

Screen Dumps

Except where specifically stated, all the screen dumps given in this book were produced from a high-resolution interlaced monochrome screen image. This means that these dumps are of a 640x400 pixel resolution. However, some of the figures were produced on a 24 pin printer using programs that I have written specifically for the task. These programs generally produce output with a very high resolution of 1416x1416 pixels. Further information on the 24 pin technique is included in Appendix D, and details of more general screen dump methods are given in Appendix A.

Entering and Debugging the Examples

All the programs contained in this book have been thoroughly tested, and should hopefully contain no errors. The listings were all saved from Amiga BASIC and GFA BASIC as ASCII files after being test run, and then directly imported into the text file. This means that little room was left for human error, and also means that indentation and variable capitalisation are shown exactly as they will appear when you enter them. However, errors may be induced while the programs are being typed in, and although the interpreters will catch some problems, such as syntax errors and incomplete loops, before the program is executed, other bugs can creep through. If you enter a program and it fails to produce the expected output, the first place to look is at lines containing numerical calculations, as a single wrong digit or variable can fundamentally alter the program’s output.

If you find typing too tedious and error prone I recommend that you obtain the support disk mentioned a few pages earlier, as this contains all the examples from this book in a form ready to run immediately from either Amiga or GFA BASIC.

In general it is safe to say that a random process produces random results, and an ordered process produces ordered results. For example, randomly sprinkling something over a flat surface produces an even spread of particles with no structured pattern, but a set of explicit commands or equations (a deterministic process) gives a predictable, ordered result as in a computer program. There are, however, some processes which do not abide by this rule of thumb. The Sierpiński triangle is the result of one such process.

The Sierpiński Triangle

Because it was discovered by mathematicians in the early 1900s, the Sierpiński triangle is one of the oldest fractals covered in this book. The reason for the triangle’s relatively early discovery is that it is based on very simple rules and can be drawn by hand, with little or no calculation. This meant that its creators were not hindered by their lack of computational machinery, the like of which is now indispensable in chaos research. Although the Sierpiński triangle can be drawn by hand it does take a considerable amount of time, so to demonstrate it here we will be using the Amiga. However, if you would like to follow the original method you are quite welcome!

Figure 2.1: Outline of the Sierpiński triangle

The starting point for the triangle is a blank plane, on which a triangle is outlined by three single points, numbered 0, 1 and 2 for reference. Such a situation is shown in Figure 2.1. The Sierpiński triangle can now be created on this plane by applying the following rules:

1. Pick one of the vertices (corners) at random and go directly to it

2. Choose another vertex at random, move halfway towards it and plot a point

3. Repeat from step 2.

This type of process is known as an iterative process, because a single set of rules is applied repeatedly in order to produce the result. Each application of the rules is called an iteration, in our example the plotting of one pixel represents one iteration. Figures 2.2a to 2.2d show the evolution of the triangle after different numbers of iterations.

Figure 2.2: The Sierpiński triangle appears

It is quite natural to expect this random process to yield a random result but, as Figure 2.2 shows, the process eventually gives rise to a rather complex, ordered structure - the Sierpiński triangle. This unusual behaviour is best demonstrated with the aid of a simple program which executes the rules given above. Such a program is shown in Listing 2.1.

Basically it can be broken down into three sections. The first dimensions two arrays to store the horizontal (x) and vertical (y) positions of each of the three vertices of the triangle, and places suitable values in them. The next section corresponds to stage one of the above process, as it initialises the drawing position by randomly setting it to be at one of three vertices. The drawing position can be thought of as the position of the computers 'pen' to give consistency with the manual drawing method. The drawing position is held in the variables <F102>px<F> and <F102>py<F>. The final section of the program is the WHILE…​WEND loop which actually does the drawing of the points, representing stages two and three of the process.

The only relatively difficult area in this program is the part which moves the graphics cursor halfway towards a vertex. This is done by finding the mid-point of the imaginary line between the current position and the relevant vertex, and moving to it. The method of finding the mid-point of the line has been borrowed from co-ordinate geometry, which states that if there are two points with co-ordinates (x,y) and (x1,y1) respectively, the mid point of the line between them is at (x+(x1-x)/2,y+(y2-y)/2). More details on graphs and the (x,y) notation can be found in Appendix B.

``  Unresolved directive in chapter2.adoc - include::../../src/amiga/2.1[]``
Listing 2.1: Amiga BASIC program to produce the Sierpiński triangle

When typing in Listing 2.1 you should ensure that you follow the guide-lines given in Chapter 1. However, as this is the first program in the book, explicit instructions will be given here to help you get the triangle up and running.

First you should load Amiga BASIC as normal and ensure that there are no BASIC programs in memory. Listing 1.1 should then be typed in exactly as shown (remarks may be omitted if desired). Assuming you have corrected any syntax errors you may have made the program can now be executed by clicking with the mouse on the start option of the run menu, or by simultaneously pressing the right Amiga key and <R>. If you have typed in the program correctly you will see the triangle gradually appear, filling most of the screen.

A Fractal

The Sierpiński triangle exhibits two important features which distinguish it as a fractal. The first is the fact that an immensely complex and structured pattern is created by just a few very simple rules. The Mandelbrot set is probably the most famous example of this, as one equation gives rise to an infinitely rich structure.

The second point is that if any part of triangular structure is sufficiently magnified the same general shape can be seen (see Figure 2.3). In this case the shape is an equilateral triangle. This property is called self-similarity, which also makes its most notable appearance in the Mandelbrot set. Note that the image cannot be optically magnified, as the resolution of the program’s output is not good enough, and even if it was possible to produce more detailed output the paper can only be magnified to a certain level. This means that the image must be magnified mathematically by completely recalculating the part of the image that needs to be enlarged.

Figure 2.3: A magnified section of the Sierpiński triangle

Unlike fractals, normal geometrical shapes are not self-similar and lose their identity when magnified enough, for example a circle becomes a straight line. An everyday example of a circle losing its identity is the phenomenon of the Earth appearing to be flat to a person standing on it, because they can only see a tiny fraction of its surface.

Figure 2.4: Magnification of the Sierpiński triangle by enlargement

Self-similarity can easily be seen in the Sierpiński triangle by altering the program in Listing 2.1. The easiest way to do this is to enlarge the triangle so that only a small part of it lies within the screen area, as shown in Figure 2.4. Using this method the program will still effectively move to each point, but will only plot points which lie on the screen. In practice the triangle is enlarged by altering the positions of the three vertices, which involves changing the section of the program which sets the vertex positions to read as follows:

```   x(0)=565
y(0)=-530
x(1)=25
y(1)=190
y(2)=1105
y(2)=190```

When enlarging fractals in this way it is imperative that the ratio of width to height (the aspect ratio) is kept constant. If this ratio is not maintained the fractal will become distorted, making self-similarity hard to identify. The Sierpiński triangle used in this chapter is 50% wider than it is tall and I have ensured that the enlarged triangle also has this property, thereby preserving the aspect ratio.

Something to note while experimenting with magnifications of the Sierpiński triangle, is its lack of substance. This may sound strange but in fact all the areas which appear to consist of solid black lines are actually full of self-similar triangular shaped holes. This continues to be true for consecutive magnifications, which suggests that the lines between the triangles are infinitely thin (hence non-existent), and implies the Sierpiński triangle to be nothing more than a group of holes! Many other objects with similar properties were devised at the same time, most notably the cube shaped Menger sponge, which, because it didn’t really exist, could hold a volume of water equal to 100% of its volume.

Although turn of the century theorists were intrigued by the strange objects that they had discovered, they found it hard to continue their investigations due to the vast number of calculations required. Unlike most classic mathematics the outcome of these types of processes could only be determined by brute force calculation, and were very sensitive to small errors. Because of these problems objects such as the Sierpiński triangle were regarded as nothing more than mathematical oddities for many years, until the advent of chaos, when natural processes were shown to exhibit similar behaviour.

Population Dynamics and the Feigenbaum Diagram

Physicists and mathematicians are lucky. The problems they have to solve are generally based on inanimate objects which can usually be studied, or at least modelled, in a controlled laboratory experiment or on paper. Biologists, on the other hand, often have to study in a much larger laboratory of the Earth, where their experiments are constantly at risk from natural sabotage. The complex interactions between living organisms and their surroundings make it impossible to enact biological studies in the lab without upsetting the results, which means that patterns and predictions are very difficult to establish.

However, to make their lives somewhat easier biologists have developed a branch of their science called ecology, in which individual creatures and events are generalised to produce simple mathematical equations which can predict the behaviour of large groups of organisms over a period of time. One such equation, used to predict simple population changes, is the May equation, developed in 1845 by Verhulst but much studied by the 'mathematical biologist', Robert May. It is this equation which forms the basis of the Feigenbaum diagram, named after Mitchell Feigenbaum, the first person to draw the whole diagram and analyse it in detail.

The Feigenbaum diagram represents the opposite situation to that of the Sierpiński triangle. That is to say the complex, and in some places chaotic structure (see Figure 2.10) is generated from the very simple, non-random formula which May derived from population dynamics. The accepted name for a non-random process such as this is a deterministic process.

Figure 2.10: 24 pin print-out of the Feigenbaum diagram

An example derivation of the equation can be obtained by considering the spread of a virus through a group of people. A virus has been used for the demonstration because viruses are among the simplest living organisms. In this case the population is defined as the number of people infected by the virus at any one time, and we are trying to find a method of predicting it at any stage in the future.

The accepted method of deriving a formula is to start from a very basic equation and then elaborate on each part of it in turn using other equations and basic facts. The equation we shall use as a starting point is shown below. It allows us to calculate the number of infected people at the end of any one week. The equation basically states that the number of people infected at the end of the current week (pnew) is equal to the number of people infected at the end of the previous week (p) plus the number of people who have been infected since then (n).

pnew=p+n

This equation is technically correct, but to calculate the population for any given week it is necessary to know the number of people infected during that week. This means that the equation can only be used to work out populations of weeks gone by, and cannot predict the future.

However, the need to know the number of newly infected people (n) can be eliminated because it can be shown that n is directly proportional to the number of people already infected at the end of the previous week (p). This is due to that fact that as the number of people infected increases so does the likelihood of an uninfected person coming into contact with an infected one, which in turn increases the chances of the healthy person being infected.

n α p

Of course, infected people cannot be re-infected so, assuming that all unwell people are contagious, the number of newly infected people must also be proportional to the number of healthy people. If, for example, everyone in the group was infected nobody else would be able to become infected, so n would be zero. Assuming that there are no in-betweens the number of healthy people must be equal to the total number of people in the group less the ones who are infected.

To eliminate the need to know the actual number of people in the group (and to show that this is unimportant to the equation) we can use percentages instead. In mathematics the numbers from 0 to 1 inclusive are generally used to represent percentages from 0% to 100%, where 0.5 is 50% for example (Appendix B gives more details on this and other pure maths topics). From now on we shall assume that p and pnew represent the percentage of the group infected at different times, and that n represents the percentage of the group who become infected during the week. The percentage of healthy people in the group is therefore 100%-p, as shown in Figure 2.5.

Figure 2.5: Pictorial representation of the method used to calculate the percentage of healthy people in the group

This can be represented mathematically as (1-p), remembering that 1 represents 100%:

n α (1-p)

Because n is directly proportional to both p and (1-p) we can say that n is proportional to p*(1-p). An equation combining these facts is:

n=p*(1-p)

You may be able to guess that if this definition of n were substituted into the equation for p shown above, p would steadily rise to 100% and then stay there. If this were the case in the real world, and the group of people were the earth’s population then the existence of the human race would be in doubt. In reality, however, certain people may not meet as many people as others, and a meeting does not mean definite infection. Also people may recover from the virus, either using the body’s immune system, or with the aid of medical treatment (depending on the nature of the virus). For this reason a constant, c, must be included in the equation. This represents the success rate at which the virus spreads and persists in the human body, and will hereafter be referred to as the contamination constant.

Because n is proportional to c we now have the following definition of n:

n=c*p*(1-p)

The result of substituting this definition into the equation for p is the May equation:

pnew=p+c*p*(1-p)

The beauty of this equation is its universality. As already discussed the use of percentages makes it independent of the number of people in the group. Since the number of people in any particular group is certain to change, due to births and deaths, this is of fundamental importance.

The introduction of the constant, c, also makes the equation more versatile by allowing it to be customised for any combination of virus and host. This means that the hosts do not have to be people, they could just as easily be animals or even computer disks! For a different virus/host combination a different value of c is all that is required.

Applying the Equation

Now that we have our equation we can use it to predict the likely number of infected people at the end of any week. Let us create a theoretical starting situation (call it week 0) in which 30% of the group have the virus (p=0.3) and the virus has a contamination constant of 1.9 (c=1.9). The calculations to determine the number of people infected after the first four weeks are shown below.

Time Calculation People infected

Week 0

p = 0.3

30%

Week 1

p = 0.3 + 1.9 * 0.3 * (1 - 0.3) = 0.7

70%

Week 2

p = 0.7 + 1.9 * 0.7 * (1 - 0.7) = 1.1

110%

Week 3

p = 1.1 + 1.9 * 1.1 * (1 - 1.1) = 0.89

89%

The results are calculated by feeding the previous week’s population into the next week’s equation. This equation is therefore said to be dependant on mathematical feedback, the flow of which is shown in Figure 2.6. Feedback occurs in many real world situations. For example audio feedback is induced in a PA system if the microphone is brought too close to the amplifier.

Figure 2.6: Mathematical feedback in the May equation

The Sierpiński triangle also relies on feedback. Each time the rules are applied their starting point is taken from the previous iteration. Like the Sierpiński triangle, our population simulation is an iterative process, but in this case an iteration is defined as one application of the formula. Such processes are notoriously difficult to predict, as there is no way of telling what a particular set of initial conditions will lead to without working through the whole process. Unfortunately many natural things like the weather are subject to feedback, and this is why powerful super-computers are needed to calculate forecasts. The more formal name used to describe such a process is 'non-linear dynamic system', and it is this definition which will be found in the more mathematically minded chaos works.

Figure 2.7a: Time series graph for the May equation where p=1.9

From the time series graph of these results (see Figure 2.7a) it can be seen that the population oscillates back and forth a few times before coming to rest at 100%. This behaviour is typical of a feedback dependant process. It begins in an unstable state and gradually returns to its equilibrium position, where it is stable. A real life example of this effect in action can be demonstrated by disturbing a simple pendulum, it will oscillate with a decreasing frequency until it hangs vertically in equilibrium. The graph also reveals a flaw in our formula - during the first few weeks the number of people infected actually exceeds 100%. This is not actually an error because values above 100% can be taken to represent over-population, where there are more occurrences of the virus than there are people.

As we already know, populations below 100% occur when some people remain uninfected by the virus, in other words when the virus is in a state of under-population. During a period of under-population the virus generally thrives, spreading itself to as many healthy people as possible. But when the people become over-populated many of the viruses cannot sustain reasonable living conditions and therefore die. This results in the population always being attracted to the best possible population value, 100%, and for this reason p=1 is said to be the attractor of the May equation (this is only true for values of c within a certain range).

Figure 2.7b: Time series graph for the May equation where c=2.3

As previously mentioned a different value of c can be used to represent a different virus. Figure 2.7b shows graphically the population of the virus with c=2.3. Again, there is some oscillation in the first few weeks before the virus settles down, but unlike the previous example this virus settles down into a two point oscillation rather than a single value.

Figure 2.7c: Time series graph for the May equation where c=2.5

Figure 2.7d: Time series graph for the May equation where c=2.9

Further values for c give even more interesting oscillations. For instance when c=2.5 the population settles down to oscillate between four values (see Figure 2.7c). By the time viruses such as that with a contamination constant of 2.9 are reached the population graphs have degenerated into chaos (see Figure 2.7d), with p jumping between many different values with no apparent rhythm.

A simple graph plotting program can be used to experiment with various viruses (i.e. values of c). Such a program is shown in Listing 2.2.

``  Unresolved directive in chapter2.adoc - include::../../src/amiga/2.2[]``
Listing 2.2: A simple program to plot a May time series graph for any value of the contamination constant, c.

This program requests a value for the contamination constant,c, and then performs the May equation for this value over a period of 80 weeks, drawing a time series graph (with time plotted horizontally and population vertically) as it goes. Remember that, in Amiga BASIC, you must select the program’s output window by clicking on it before you can type in the constant at the prompt. To keep the Listing short no axes are plotted, but this does not matter because the actual values are unimportant, the program is simply used to demonstrate different behaviour patterns.

The main body of the program is the FOR…​NEXT loop containing the easily recognisable May equation, and the LINE command which draws a line from one point to the next. Note that because the range of the week variable (1 to 80) is small in relation to the horizontal screen resolution (0 to 639) used by LINE, the x position passed to LINE must be multiplied by eight so that the whole screen is used.

Similarly the population variable,p, has to be made larger due to its very small range (0 to 1.3), by multiplying it by 140. However, due to the inverted nature of the Amiga’s vertical axis the amplified value of p must be taken away from the maximum vertical position (200) to ensure that the graph is plotted the right way up. More details on the Amiga’s co-ordinate axes are given in the graphs section of Appendix B.

Experimentation with this program should establish the following facts:

• p always takes a few weeks to settle into a pattern

• In most cases the higher the value of c the more values p oscillates between after settling down

• p always oscillates between an even number of values, except in the chaotic regions, and when c is one of a certain set of values (try c=2.83)

• Values of c above 3.0 give meaningless results, as the attractor for this part of the diagram is minus infinity.

Using this program it would be very time consuming, and difficult, to find the exact value of c at which the number of values oscillated between curiously goes from 1 to 2. A better way to discover this and other things is to combine the graphs of all the possible values of c in one graphical representation. The composite graph that results from doing this is referred to as the Feigenbaum diagram, the generation of which is a relatively simple affair, now that we have a program to calculate successive values of p.

To draw the Feigenbaum diagram it is necessary to compress the time series graph for each value of c so that it fits into one vertical column of pixels on the screen (i.e. so that it is one dimensional). All these compressed graphs are then drawn horizontally across the screen, so that they join together to form a two dimensional full-screen map of p against c. This process is shown pictorially in Figure 2.8.

Figure 2.8: Relationship between May time series graphs and the Feigenbaum diagram

The first stage of creating the compressed graph can be achieved by adapting the time series program to the form shown in Listing 2.3. There are two important changes here. Firstly the virus population is now calculated for 100 weeks, but the first 50 of these calculations are ignored to ensure that p is only plotted after stabilisation. Because the system relies of feedback the first 50 calculations must be performed, even though they are not used.

Figure 2.9: One dimensional graph of the four-point oscillation associated with c=2.5 (the pixels have been enlarged to aid visibility).

Secondly all points on the graph are now plotted at the same horizontal position, and are not connected by lines, making oscillatory patterns easily identifiable. A sample run of the program with c=2.5 is shown in Figure 2.9, with the four points of oscillation easily visible. You could try experimenting with this program using the notable values of c mentioned above to gain an understanding of the program’s function.

``  Unresolved directive in chapter2.adoc - include::../../src/amiga/2.3[]``
Listing 2.3: The program for plotting compressed graphs

If graphs of this sort are now plotted side by side for values of c between 1.8 and 3.0 inclusive (values below 1.8 are uninteresting, as they stabilise at a single value) we get the Feigenbaum diagram shown in Figure 2.10. Listing 2.3 can be altered for this purpose by replacing the human input for c by one produced by a BASIC FOR…​NEXT loop. Listing 2.4 has had this modification included so that it now produces the whole Feigenbaum diagram.

``  Unresolved directive in chapter2.adoc - include::../../src/amiga/2.4[]``
Listing 2.4: Program to draw the complete Feigenbaum diagram

The Feigenbaum diagram is not one of the most aesthetically pleasing fractals, but it is a stunning illustration of how a simple, non random, iterative process can produce a finely structured image in some places, and total chaos in others.

One result of much scientific study of the diagram is that conventions have been created to describe some of its most important features. The section of seemingly random pixels to the right of the diagram is referred to as the chaotic region. Note that this area is not actually random, as the results are the same each time the program is executed.

The point where a single line splits into two is called a bifurcation and the splitting, or bifurcating, which occurs at these points is referred to as period doubling, as the number of equilibrium states (the period) doubles. The period is the number of possible values that p oscillates between after stabilising, for example the section between 2.44 and 2.54 is period four. By the very nature of period doubling almost all sections past c=1.95 have an even number of possible states, but there are however some 'windows' of order occurring in the chaotic regions which are of period three. Table 2.1 shows the period for various sections of the Feigenbaum diagram.

Range of c (approx) Period

0.00 to 1.95

1

1.95 to 2.44

2

2.44 to 2.54

4

2.54 to 2.56

8

2.83 to 2.84

3

Figure 2.11: Mathematically enlarged section of the Feigenbaum diagram’s main window (24 pin)

Figure 2.11 shows a section of the main window in the chaotic region enlarged to about seven times its original size. This reveals some interesting features. The most striking of these is the very noticeable self-similarity. Inside the window a miniature Feigenbaum can be found, complete with characteristic period doubling and window-ridden chaotic region. This section also highlights the incredibly random nature of the chaotic region. Comparing this to a screen full of pixels randomly placed using Amiga BASIC’s RND function reveals a striking similarity.

Enhancements to the Feigenbaum Program

Musical fractals are rare, but the Feigenbaum diagram is one of the few that lends itself to the addition of sound. By playing a tone proportional in pitch to the value of the population, p, whenever a point is plotted the sound can be used to give an audio representation of how the Feigenbaum diagram degenerates into chaos. At the left hand side of the diagram the period is 1, so the tone is constant. After bifurcation the tone oscillates between two pitches, and by the end of the diagram chaos has broken out, creating a bubbling sound of pseudo-random tones. Adding the following command after the PSET line will create the sound:

``  SOUND p*1000+200,1``

A piece of information lost in the Feigenbaum program is the number of times each point has been plotted, which would be useful when looking at the number of equilibrium states in different parts of the Feigenbaum diagram. For example, in the period one section of the diagram each point is plotted 50 times, whereas in the period two region every point is plotted 25 times, and so on until in the chaotic region most points are not plotted more that once.

The number of iterations for each point could be conveyed either by using colour or by plotting the diagram in three dimensions, where the third dimension represents the number of times the point is plotted. The three dimensional 'Feigenbaum landscape' method is discussed in Chapter 7, as it obviously requires a knowledge of three dimensional drawing techniques. The colour method is somewhat simpler and the program given in Listing 2.4 can easily be adapted to use colour. The default screen set-up used by Amiga BASIC only allows two bit planes, limiting the number of colours to four. It is therefore necessary, if we are to make full use of the Amiga’s colour capabilities, to open a new screen with four bit planes (giving a total of 16 colours) before drawing the diagram. This can be done by adding the following two lines at the head of Listing 2.4:

``````  SCREEN 1,640,200,4,2
WINDOW 2,"Colour Feigenbaum Diagram",(0,0)-(617,180),15,1``````

In simple terms these two lines open a high resolution, non-interlaced, screen with a 16 colour capacity) and then open a window on that screen to which all subsequent Amiga BASIC output will be sent. Comprehensive details of the SCREEN and WINDOW commands can be found in appendix A and in the Amiga BASIC reference manual.

The colour plotting algorithm can now be incorporated into the program by replacing the PSET line with the following program fragment:

``````  xp=(c-1.8)*520              'Calculate x and...
yp=200-p*140                '...y position of point
oldcolour=POINT(xp,yp)      'Determine old colour
COLOR (oldcolour+1) MOD 16  'Set new colour
PSET(xp,yp)                 'Plot the point``````

Note the way that the x and y position of the point are calculated and stored (in `xp` and `yp`) to save having to calculate them a second time later on. The actual colour in which a point is plotted is determined by finding the old colour of the point at the plotting position (found using the POINT function) and adding 1. Because we are limited to 16 colours, modulus 16 is taken of the colour value in order to keep the colour within the relevant 0 to 15 range (for more information about MOD see Appendix B).

At first this process may sound like the start of a complicated magic trick, but it is actually the simplest way of changing the colour of a pixel each time it is plotted so that the colour of the point shows how many times it has been visited.

Conclusion

In this chapter it has been shown, using the Feigenbaum diagram, that a simple mathematical equation, with no random elements can give rise to what appears to be a totally random result, and that a random set of rules can give rise to a pattern of very rich structure, such as the Sierpiński triangle.

This two way relationship between order and chaos is more than anecdotal, it represents a whole new way of perceiving natural processes occurring around us. For although our simple virus population program has inaccuracies it still illustrates the important point that a natural process such as population growth and decline can be described in one simple equation with nothing more than some common sense and a little help from the laws of statistics. Is it possible, then, that even the most complex natural processes can be reduced to simple equations, thus allowing them to be predicted with relative ease? This important question is dealt with in the next chapter.

3. Weather, Chemistry and Strange Attractors

Weather forecasters have a bad reputation. People find it incredible that a forecaster can proudly announce "there will be no hurricane" merely hours before the country is hit by some of the worst storms in living memory. In Britain forecasters are particularly unfortunate because our island location means that they have to deal with a wide variety of weather, from blizzards to heat waves.

The seeds of modern forecasting methods were sewn in the 50s and 60s, with the advent of computers. Meteorologists had long believed that they could represent the world’s weather using a complex mathematical model, and now hoped that computers would be able to perform the multitudinous calculations needed to run such models. With the advantages of accurate weather forecasting being obvious to even the most unscientific of people much time and money was put into meteorology research. This was spent on various projects, ranging from collecting weather data to building specialised weather forecasting computers. For decades the most powerful civilian computers were always to be found at weather forecasting centres.

One beneficiary of the forecasting boom was Edward Lorenz, a research programmer at the Massachusetts Institute of Technology. Being a keen mathematician and an experienced meteorologist Lorenz was quite at home writing programs on his rather elementary computer (by today’s standards) to model simple weather systems. In his first model Lorenz simplified the weather down to a fluid dynamics system, sea and air being fluids, and chose 12 relevant equations each representing a different element of the weather, and each interacting with the others.

The model was very simple, some thought too simple, but the output from it had some interesting similarities with the behaviour of real weather. Everything that happened in the Lorenz model followed on from what had happened previously. It was by no means random, but the exact same series of events never happened twice. Lorenz managed to interest many colleagues in the model, and despite the fact that all it did was continuously churn out numbers on a Teletype (he had no graphics facilities) they would enjoy visiting Lorenz’s office to see what the simulated weather was up to.

Because the program was rather slow Lorenz decided to try and simplify his model further, as far as possible without destroying its realistic unpredictability. He eventually trimmed the model down to the three non-linear differential equations:

dx/dt = a*(y-x)

dy/dt = b*x-y-x*z

dz/dt = x*y-c*z

where:

x, y and z are variables representing different weather aspects;

t is time;

a, b and c are constants.

The way that these equations were derived places their detailed explanation beyond the scope of this book, but knowledge of how to use them is important for later examples. As with the May equation described in the previous chapter, the non-linearity of the equations is due to the mathematical feedback created after successive iterations. The feedback in the Lorenz model is a good deal more complex than in the population simulation because there are three equations and three variables, as shown in Figure 3.1.

Figure 3.1: Feedback in the Lorenz model

The Lorenz equations are referred to as differential equations because, unlike normal ones, they are used to find the rate of change of variables, rather than the absolute values. For example dx/dt represents the rate of change of x, or more explicitly the (small) change in x that occurs in the (small) time, t. The d in front of the variables stipulates that the changes in them must be small - technically the d should be written as δ (the Greek letter delta). The reason that the change must be small is that the equations were not designed for large changes, to perform a large change you have to carry out the equivalent group of small ones.

The actual values of x, y and z are calculated by repeatedly adding the change in the variable (e.g. dx) calculated using the formula, to the previous value of the variable (e.g. x). Lorenz’s original programs used these equations to simply produce a list of numbers, similar to those shown below.

 iteration:0 x=1.000 y=1.000 z=1.000 iteration:1 x=1.000 y=1.260 z=0.983 iteration:2 x=1.026 y=1.518 z=0.970 iteration:3 x=1.075 y=1.780 z=0.959 iteration:4 x=1.146 y=2.053 z=0.953 iteration:5 x=1.236 y=2.342 z=0.951

These figures meant a lot if you could decipher them, but general trends and repetition were hard to spot. However, Lorenz invented a primitive kind of graphics on his system to make the results more accessible. He chose one of the variables and instead of plotting it numerically he scaled it so as to be in the range of 0 to 79 and then plotted an 'a' character on the current printer row, that number of spaces across the page. Over a period of a few days a rather long time series graph would be produced showing the strange non-periodic behaviour of the model. Similar output is shown in Figure 3.2.

Figure 3.2: Simulation of Lorenz’s first time series graphs

Graphics are now almost a standard feature of modern microcomputers, and an acceptable line graph of the information can easily be created on the Amiga with the aid of a simple program, such as the one given in Listing 3.1.

``  Unresolved directive in chapter3.adoc - include::../../src/amiga/3.1[]``
Listing 3.1: A program to produce a time series graph from the Lorenz equations

This program can be broken down into two main sections. The first is the initialisation of all the constants, variables and graphics facilities. The second is the loop in which successive values of x, y and z are calculated, and where z is plotted vertically against time. It is important that all three variables are calculated, even though only one is displayed, because all the variables interact.

There are several experimental alterations that can be made to this program, the most obvious being the values of the constants a, b and c, and the initial values of x, y and z. You could also change the LINE call so that x or y are plotted instead of (or as well as) z. The detail of the plot may also be altered by using the <F102>detail<F> variable, which determines how many 'dummy' iterations are executed for each point plotted on the screen - the larger the number the less the detail (but the more time you can fit on the screen).

As previously discussed, the dummy calculations must be processed even though they are not plotted, as large changes in the variables (and in time) cannot be made due to the nature of the equations. Note that if you do change this variable or the constants you will probably have to change the scaling factor used to calculate the y parameter of the LINE call. However, this is not all that easy, as it is impossible to calculate the maximum and minimum values of x,y and z without actually calculating the result of the three equations for every iteration.

Figure 3.3: Time series graph produced by Listing 3.1

The kind of output that you can expect is shown in Figure 3.3. No axes are drawn as actual values are irrelevant without a full description of the equations' derivation.

The graphs produced by this program are not all that aesthetically pleasing, but they do allow the non-periodic repetition of similar events to be spotted with relative ease. An example of such an event is the vertical oscillation which gets larger and larger before collapsing and beginning again, with a slight difference.

The Butterfly Effect

The is no doubt that Lorenz’s carefully planned experimental models were very enlightening, but probably his most important discovery was accidental. Wishing to hurriedly continue an earlier run of his model Lorenz used initial values of x, y and z taken from an earlier printout. He assumed the fact that the values on the printout were rounded off to three decimal places (from his computer’s more accurate internal values) to save space would not make much difference to his readings.

To simulate the original run of the model we can alter the relevant lines in Listing 3.1 so that the initial values of x, y and z are set up as follows (it is important to type the numbers in full):

```x=8.164710103228572
y=8.967625255594612
z=25.50035258277785```

These numbers are of the same precision as those used internally in the program by Amiga BASIC, and are said to be shown to 16 significant figures. If the program is now executed (with the `detail` variable set to 5) a graph similar to that shown in Figure 3.4a will be produced, which is not particularly unusual. However, if the program is executed with the initial values rounded off to four significant figures completely different results can be seen.

Figure 3.4a: Time series graph using initial conditions accurate to seven significant figures

The initial values can be rounded off by changing the relevant lines to read as follows:

```x=8.165
y=8.968
z=25.50```

Rounding these numbers off is effectively what Lorenz did by typing in the old printed out values. The output from these less precise initial conditions is shown in Figure 3.4b.

Figure 3.4b: Time series graph using initial conditions accurate to three significant figures

Comparing the two lines reveals that although they follow the same path for the first few iterations they gradually diverge into totally different patterns. Figure 3.4c illustrates this more concisely by super-imposing the two lines on the same graph. The effect witnessed here is what is technically known as sensitive dependence on initial conditions, meaning that the equations are so sensitive to the initial conditions supplied to them that even the smallest difference in the initial values quickly grows into a large difference.

Figure 3.4c: Figures 3.4a and 3.4b superimposed

Lorenz immediately realised what his results meant. They showed that his model, and hence the weather, had a sensitive dependence on initial conditions. Ultimately this meant that weather prediction was practically impossible, because even the smallest error from a thermometer or other instrument would send a forecaster’s model off down an increasingly wrong track. Other meteorologists took a more optimistic view, believing that they could alter future weather by making the relevant small changes such as heating small areas of land.

They were of course right, but Lorenz pointed out they could not control the weather with any precision, as the effects of the small changes they made would be impossible to predict because other, undetected, small changes would have equally large results. Some argued that a complex system such as the weather was so sensitive that the path of a hurricane could be reversed by a butterfly taking off from the other side of the globe. Hence the name 'butterfly effect' was coined as the common term for sensitive dependence on initial condition.

A more down-to-earth example of the butterfly effect is the story of an airline company going bankrupt. A man leaves his house on foot one morning, with the purpose of going to the travel agents to book his holiday. After walking a few yards down the street he looks down and finds that his shoelace has come undone. Naturally he bends down to tie the lace, but in doing so he loses his balance and falls into the road. As he stumbles to his feet a passing bus collides with him and he is subsequently taken to the local hospital. Because of this the man is unable to book his flight before the deadline, meaning that the airline company is one person short of the quota needed to make the flight financially viable. The flight is cancelled, and ticket holders refunded, but with the absence of the income from that flight the already troubled airline company gets into financial difficulties and goes bankrupt. It seems incredible, you never read of companies going bankrupt because of undone shoelaces, but it could happen, and it would be almost impossible to predict.

As the airline example shows, weather is not the only thing that would suffer because of the butterfly effect, world money markets, politics and the lives of individuals can be drastically affected by even the smallest changes. The airline discussed above could have been a major company, and its downfall could have pushed an already gloomy economy below a critical level and toppled the government and lead to civil war.

The butterfly effect raises important questions about the wisdom of organisations around the world who expend considerable resources on weather prediction. However, in the short term the errors are small, so forecasting a day or two in advance is just about possible with a fast computer, accurate algorithms and a bit of luck, but past that predicting the weather becomes more of an art that a precise science.

At first it seems strange that something as simple as the Lorenz model is so unpredictable, when more complex periodic systems such as the ebb and flow of the tides, and planetary motions are so easy to forecast. The reason that periodic systems are much easier to predict is that their repetitions are easy to observe and describe mathematically. Almost all non-linear dynamic systems (including the Lorenz model) are non-periodic, meaning that they are prone to the butterfly effect, and hence unpredictability and chaos.

The Lorenz Attractor

The Lorenz model has three variables, three directions in which it can change, and is therefore said to be three dimensional. This means that instead of simply plotting a time series graph of z against time we could plot x against y against z on three dimensional axes, like those shown in Figure 3.5.

Figure 3.5: Typical three dimensional axes

The shape produced is known as the Lorenz attractor, an infinitely long complex spiral which never intersects itself. Plotting it in all three dimensions is rather involved and unnecessary, as a good representation can be achieved by plotting z against x in two dimensions.

Figure 3.6: The Lorenz attractor for a=10, b=28, c=8/3

Figure 3.6 shows the graph that results from this process, but to fully appreciate the line’s behaviour it is necessary to see it being drawn. Listing 3.2 can be used for this purpose.

``  Unresolved directive in chapter3.adoc - include::../../src/amiga/3.2[]``
Listing 3.2: Plotting z against x in two dimensions to produce the Lorenz attractor

The Lorenz attractor takes its name from the mathematical definition of such an object. An attractor is a set of points, in an imaginary mathematical space known as phase space. Each point on the attractor represents the state of the system at a particular time, and the line showing the sequence of these points is known as the trajectory. By looking at the trajectory of the attractor the behaviour of the system can be determined. The idea of attractors has been around for some time, but non-periodic attractors (called strange attractors) such as the one created from the Lorenz model are relatively new.

As this program shows, the trajectory behaves unpredictably, switching from lobe to lobe almost at random, but it does produce a very complex structure. One of the curiosities of the line is that it never intersects itself, although on the two dimensional screen of the computer it looks like it does. Because the Lorenz attractor is non-periodic it cannot intersect itself, if it did it would rejoin an old trajectory and periodically repeat. However, if the line never crosses and goes on forever it must be true that an infinitely long line can exist in a finite three dimensional space. As with many chaos ideas this obviously contradicts the laws of physics, but in fact it is possible if the space in which the line is enclosed is what is known as folding space (see the books listed in the Bibliography for details).

The Lorenz attractor is similar to many natural systems which exist in real space. For instance the turbulence in a fluid and the trajectory of an object orbiting a two planet system have similar properties. A recently developed executive toy, consisting of a metal pendulum oscillating over a bed of magnets, also exhibits similar behaviour.

The butterfly effect raises an interesting question: what is the real Lorenz attractor? Lorenz drew his first diagrams using output from a program using numbers correct to six decimal places. The Amiga BASIC variables which we used had a precision of 15 decimal places, so our program will produce a different output to that achieved by Lorenz. So who is right? This is a good question, but although the values may differ the general characteristics are the same.

The Rössler Attractor

Using three non-linear differential equations derived from chemistry it is possible to draw another strange attractor, the Rössler attractor. This exhibits all the same features as the Lorenz attractor, but is a different shape (see Figure 3.7). The equations are as follows:

dx/dt = -(y + z)

dy/dt = x + y * a

dz/dt = b + z * (x-c)

Figure 3.7: The Rössler attractor

Listing 3.3 can be used to draw the attractor. Note that in most respects it is the same as the Lorenz program (Listing 3.2), except that the three equations are different and a pseudo-3D effect is created by plotting y+2*z (instead of just z) against x.

Listing 3.3: Program to draw the Rössler attractor
``  Unresolved directive in chapter3.adoc - include::../../src/amiga/3.3[]``

Again the constants and initial variables can be altered to produce a different Rössler attractor.

Enhancements to the Attractor Programs

It has been mentioned in this chapter that the trajectory of a strange attractor never intersects itself. This is obviously true because the systems on which strange attractors are based are non-periodic (non-repeating), and for the trajectory to go through the same point twice would obviously be a repetition. Because the attractors are three dimensional but the Amiga’s screen is only two dimensional it often looks as though intersections are occurring, where lines are actually passing behind or in front of one another.

The simplest way of verifying the non-periodic nature of the attractor is to use colour to represent the third dimension. If each point was plotted in a colour proportional to the value of y in the program given in Listing 3.2 then, ideally, two adjacent pixels would be touching only if they were both of the same colour. However, the range of values for y is considerably larger than the range of colours on the Amiga making this method unusable.

An alternative, and very impressive, way of getting round this problem is to write a program to produce an isometric 3D attractor which rotates around the z axis. The rotation gives the impression that the attractor is truly three dimensional, and just by watching it spin it can be seen that the line doesn’t cross. Such a program is beyond the scope of this book, but a less ambitious one allowing static 3D versions of the attractor to be drawn at any angle of rotation would allow the same observations to be made. Three dimensional drawing methods are introduced in Chapter 7.

There are many more strange attractors which are just as interesting as those discussed here. Information on other attractors can be found in some of the books listed in the bibliography near the end of this book.

4. The Mandelbrot Set

If you asked a scientist or computer user, of any nationality, what single word comes into their mind when chaos or fractals are mentioned they would most probably say 'Mandelbrot'. The colourful Mandelbrot set is by far the most famous image from the world of chaos, and has now become the emblem of the science. The set’s fame has been brought about by its propagation throughout the world by the media, who picked up on this particular branch of modern science because it produces aesthetically pleasing pictures, in colour, which can be appreciated by everyone. The set has now found its way into all manner of places, from T-shirts and record covers to newspaper articles about the economy. Unfortunately, however, the origins and meaning of the set are often distorted or even lost during the journey to the popular press. Even computer magazines only print programs to draw pretty pictures, rather than explain the underlying concepts. This chapter has been written to help rectify the situation, using various programs to demonstrate the theory.

Benoit Mandelbrot could not have anticipated that his work would have such wide popular appeal when he saw the set for the first time in early 1980. Mandelbrot’s first plots were ill-defined monochrome printouts, interesting not for their appearance, but for their relevance to pure mathematics research. This research concerned the behaviour of iterative processes involving complex numbers. Something which often appears contradictory in chaos literature is that the process for plotting the Mandelbrot set is said to be very simple (just one equation), but complex numbers make the process seem a good deal more involved. In actual fact complex numbers are really just simple mathematical tools, and not a crucial part of the process at all.

It is interesting to note that the Mandelbrot set, unlike most of the fractals discussed so far, originated from purely mathematical beginnings, not from a study of nature. However, the set does exhibit the interesting features associated with other fractals, most notably the infinite complexity which results in the most spectacular self-similarity seen in this book.

This chapter describes two methods for drawing and experimenting with the Mandelbrot set. The first uses a simple scenario of a circle, a line, and some elementary maths. The second method follows Mandelbrot’s original technique of using complex numbers. If you feel your maths is not up to the relative challenge of complex numbers you will probably find it easier to follow the circle method, but a subsequent glance at complex numbers may prove fruitful. Note that both methods result in the same program.

The Circle Method

Like the strange attractors discussed in the previous chapter the Mandelbrot set is based on a deterministic non-linear process, but whereas the results of the Lorenz and Rössler equations are drawn directly onto the screen as a graph the relationship between the Mandelbrot equation and the resulting set is more complicated.

Figure 4.1: The Mandelbrot circle

The first part of the process takes place in a circle such as that shown in Figure 4.1, with centre (0,0) and radius 2. Inside this circle there is a single point, at an initial position (a,b), to which the following non-linear equations can be applied (for details of indices and the ^ notation see Appendix B):

p~new = p^2 - q^2 + a

qnew = 2*p*q + b

In BASIC this is most suitably written:

`pnew = p*p - q*q + a`

`qnew = 2*p*q + b`

Figure 4.2 helps to explain further these two equations by showing the feedback that occurs between them on successive iterations.

Figure 4.2: Feedback in the Mandelbrot equations

So far we have only specified that a is the initial horizontal position of the point and b is the initial vertical position. No mention has been made of _ p_, pnew, q or qnew. Naturally these four variables will be subject to change as the equation is applied, and the values of _ Pnew_ and qnew will automatically be worked out because they are the results of the equation. That still leaves p and q unaccounted for, and although these will inherit values from pnew and qnew during subsequent applications of the formula they require some initial values to get the process started. For the time being we shall say that p and q are both set to 0 before the process is invoked, so that the first position calculated is the initial position, (a,b).

Now that we have our set of initial values, and the full equation, we can perform a 'dry run' of the process for one point, as shown below. In this example our starting point will be (-1,0.5), i.e. a=-1, b=0.5.

Initially: a = -1 b=0.5 p=0 q=0

Application of the equation gives:

pnew = _p*p - q*q + a q~new = 2*p*q + b

`= 0*0 - 0*0 + (-1)			= 2*0*0 + 0.5`
`= -1				= 0.5`

So the first position is (a,b) as expected.

p now becomes the latest value of pnew, and q becomes qnew, this can be represented as:

<$table4>p = pnew (= -1) q = qnew (= 0.5) The equation can now be re-applied…​ pnew = p*p - q*q + a qnew = 2*p*q + b `= (-1)*(-1) - 0.5*0.5 + (-1) = 2*(-1)*0.5 + 0.5` `= -0.25 = -0.5` p = pnew q = qnew …​and re-applied…​ pnew = p*p - q*q + a qnew = 2*p*q + b `= (-0.25)*-(0.25) - (-0.5)*(-0.5)+(-1) = 2*(-0.25)*(-0.5)+0.5` `= -1.1875 = 0.75` p = pnew q = qnew Only three iterations have been performed here to save space, but had more iterations been performed the results in Table 4.1 would have been observed.  Iteration number p q 0 0.0000 0.0000 1 -1.0000 0.5000 2 -0.2500 -0.5000 3 -1.1875 0.7500 4 -0.1523 -1.2813 5 -2.6184 0.8904 etc…​ These results are reminiscent of those produced by the virus equation in Chapter 2, in that they are very complex, and appear to be almost random. Of course the results are not random as they come from a deterministic process, which if repeated will yield the same values. Figure 4.3: The path of the point with initial position (-1,0.5) Returning to the wider situation we can now plot the path of the point in the circle, which results in the image shown in Figure 4.3. This output can easily be achieved on the Amiga by calculating the values of p and q (hence the position of the point) for successive iterations, and then drawing straight lines between these positions. A program to do this is shown in Listing 4.1. Note that the iteration count is displayed in the top left of the output window to give an impression of just how many calculations are taking place. Listing 4.1: Plotting the path of the point in the circle `` Unresolved directive in chapter4.adoc - include::../../src/amiga/4.1[]`` The program automatically stops if the line breaches the circle, but can also be halted by hitting the right Amiga key and the period (.) key simultaneously. Checking whether the point is outside the circle is done mathematically, using the technique shown in Figure 4.4. Figure 4.4: Method of testing whether a given point is in the circle The diagram shows that if (p,q) and (0,0) are taken as being two corners of a right angled triangle, the distance between these two corners is equal to the length of the triangle’s longest side, the hypotenuse. The length of the hypotenuse (r) can be calculated from the length of the other two sides using Pythagorus' theorem, the formula for which is shown below. r2=p2+q^2 If the point (p,q) is inside the circle it is impossible for the hypotenuse to be greater than 2, so: r⇐2 squaring both sides gives: r^2⇐4 and after substitution into the Pythagorus equation: p2+q2⇐4 It is this final form of the equation which is used in the WHILE line to check whether the program should end. Appendix B provides a more detailed description of Pythagorus' theorem. Because the range of values for p and q is so small the plotting co-ordinates must be significantly amplified in order to make full use of the screen. This is done, as usual, in the PSET and LINE lines. Also familiar is the reversal of the q position, required because the Amiga’s screen co-ordinates are numbered from top to bottom rather than in the conventional bottom to top way (see the graphs section of Appendix B for details). From Figure 4.3 it can be seen that the route which the point takes is not especially structured. This is true of most paths plotted by the program, but sometimes patterns do seem to emerge, for example the pentagram in Figure 4.5a which is the path of the point with initial position (-0.52,-0.55). Figure 4.5a: Behaviour of the point with initial position (-0.52,-0.55) This is all very interesting, and you may want to spend some time experimenting with various initial positions of the point, but how does this relate to the Mandelbrot set? Well, each set of initial values (a and b) can be put into one of two separate categories depending on the path of the corresponding point, generated by applying the Mandelbrot equation. The two categories of points are: 1. Those whose paths rapidly accelerate to infinity 2. Those whose paths never leave the confines of the circle It would obviously take an infinite amount of time to thoroughly test a point in order to determine whether it flies off to infinity or not, but it can be shown that if the point breaches the boundary of the radius two circle it will eventually reach infinity. Using this reasoning we can say that the initial point (-1.5,0.5) demonstrated above is an example of a category one point. An example of a category two point is (0.35,0.35), whose path is shown in Figure 4.5b. Figure 4.5b: Behaviour of the point with initial position (0.35,0.35) This was tested using the program above, and had still not exited the circle after 20,000 iterations. Unfortunately there is no simple way to prove that a line will not reach infinity, so the brute force method of applying the equation many times must be used. The point at which we terminate the program and assume the point will never leave the circle is arbitrary, but iteration ceilings of 32 are generally typical (details of how to determine this number are given later). Figure 4.5c: Behaviour of the point with initial position (0.3505,0.35) To underline the behaviour of the equation Figure 4.5c shows another type one line, this time produced from the initial point (0.3505, 0.35). It is interesting to note that just the small change (a is 0.0005 units larger than in Figure 4.5b) produces an eminently different result. Such sensitivity to small changes in the initial conditions is typical of non-linear systems, as discussed in the previous chapter. The categorisation of points according to their behaviour in the circle is the key to creating the Mandelbrot set. What we do is to take every point on another plane, where -2⇐p⇐2 and -2⇐q⇐2, put them through the equation and then categorise them. If all the type two points are coloured black and all others are left white we get the image shown in Figure 4.6, commonly known as the Mandelbrot set. Figure 4.6: Mandelbrot set produced using the 24 pin printing techniques introduced in Appendix D Figure 4.7: Relationship between the Mandelbrot circle and the set Figure 4.7 has been included to help visualise the process and shows the relationship between points of the Mandelbrot set and their associated circle planes. A program to produce a plot similar to the one shown in Figure 4.6 is given in Listing 4.2. `` Unresolved directive in chapter4.adoc - include::../../src/amiga/4.2[]`` Listing 4.2: A program to produce the Mandelbrot set Note that not all the points on the 4x4 plane can be tested because, as on any plane, there are an infinite number, meaning that they would take an infinite amount of time to process and they would not fit on the screen. Instead we plot a small portion of regularly spaced samples, allowing a reasonable approximation of the Mandelbrot set to be built up. An approximation is the best that can be achieved because of the limitations imposed by the Amiga’s screen. The actual number of points calculated is determined by the distance between them. In Listing 4.2 this is the STEP size of the FOR…​NEXT loops. This distance is optimised here in order to give the maximum detail in the smallest amount of time, so naturally varies depending on what screen resolution is being used. In simple terms we are running a simple, two dimensional, strange attractor with a series of different initial conditions. However, the ceiling of 32 iterations for each point means that mathematical feedback has a good chance to work without the problems introduced by the butterfly effect. Adding Colour Technically the Mandelbrot set is defined as the map of all the category two points on the plane (i.e. all the points shown in Figure 4.6), but most people are more familiar with the version surrounded by colourful curved lines. These lines are known as contour lines, and are analogous to those found on a map, except that instead of representing height they represent the ease with which the points left the circle. Listing 4.2 can easily be adapted to display these contours. Because the colour of each point represents the ease with which it left the circle, the colour can obviously be deduced from the number of iterations performed before the point escaped - the fewer iterations it took the easier it was. May of the commands used to plot colour Mandelbrot sets were also required when plot colour Feigenbaum diagrams in chapter two. It is therefore recommended that you refer back to that chapter if you are unsure of how to use colour (or open new screens and windows) in Amiga BASIC programs. In the Amiga’s high resolution non-interlaced mode we can determine the colour by taking MOD 16 of the number of iterations required to free each point from the circle. This gives a range of colours between 0 and 15. Obviously colours are repeated, meaning that no one colour can be assigned to any one iteration number, but the colour cycle ensures that no two adjacent contours are the same. The iteration ceiling of 32 ensures that all members of actual Mandelbrot set are left black in line with conventions (32 MOD 16 = 0). For more information on the MOD operator see Appendix B. To add colour the following two lines must be inserted at the head of the program in order to open a screen capable of displaying 16, rather than the usual four, colours: `````` SCREEN 1,640,200,4,2 WINDOW 2,"Colour Mandelbrot Set",(0,0)-(617,180),15,1`````` The IF…​END IF section of Listing 4.2 should be removed, and the following routine substituted in its place: `````` COLOR iteration MOD 16 `Set colour PSET(320+a*100,100-b*50) `Plot point`````` If, however, your Mandelbrot is destined for a monochrome monitor or a printer you will only have a palette of two colours to choose from. In this case the most effective method of displaying the set is to plot a point in the foreground colour if the number of iterations needed to remove that point from the circle is odd, but otherwise to leave it black. An alternative IF…​END IF replacement is shown below (note that the SCREEN and WINDOW lines are not required for monochrome plots). Because the maximum iteration number is even (32), the actual set is automatically left black. `````` IF iteration/2 <> INT(iteration/2) `If iteration number is odd... PSET(320+a*50,100-b*50) `...plot the point END IF`````` Figure 4.8a: Nine pin low-resolution screen dump of Mandelbrot set with colour contours Figure 4.8b: Nine pin low-resolution screen dump of Mandelbrot set with monochrome contours Standard nine pin dot-matrix screen dumps showing the output of both versions of Listing 4.2 are given in Figures 4.8a and 4.8b. Figure 4.8b most clearly shows the boundaries of the contour lines, and as this book is set in only two colours this monochrome convention will be used in all subsequent Mandelbrot related figures. Naturally Figure 4.8a would have benefited from being dumped to a colour printer, and colour monitor users will probably want to replace the monochrome plotting method used in the examples with the 16-colour equivalent described above. The majority of Amiga BASIC Mandelbrot programs can take the best part of 2 hours to produce the whole contoured set, although exact times depend on the nature of the program and the screen resolution. Hopefully you can now appreciate why they take so long, remember that for each point on the Mandelbrot set the program must determine and categorise the behaviour of the corresponding line in the circle. If we assume that on average 16 iterations are required to categorise each point, and that the size of the relevant screen area is 400x200 pixels (low resolution) then altogether 400x200x16 = 1,280,000 iterations must take place. That’s 1,280,000 applications of the two equations, the circle breaching test and all the associated program instructions. For a high resolution interlaced plot (400x400 pixels) the number of iterations is around the 2560,000 mark. This is a clear demonstration of why computers are now so relevant to chaos research. Of course, not all points require this number of iterations, and if you observe the contour drawing version of the program in action you can easily spot the changes in calculation speed. This is most noticeable in and around the actual Mandelbrot set, where more iterations must be performed per point than anywhere else. The Complex Number Method As mentioned in the first part of this chapter one of the Mandelbrot set’s most outstanding features is that such an intricate shape can be created from such simple rules. The circle and line method is not a particularly good way of showing this, as it is rather involved. By contrast the use of complex numbers can allow the Mandelbrot to be created from the single equation: znew = z^2 + c This is a very simple equation with mathematical feedback, much like the virus population one described in Chapter 2. In plain English this equation states that z is equal to the previous value of z multiplied by itself, plus a constant, c. In BASIC the equation might be entered as: z = z*z + c If we repeatedly apply this to normal numbers such as z=0.5 and c=1 we find that z tends to zero or infinity in a very ordered fashion. However, if complex numbers are used more interesting results can be observed. Before proceeding further a basic understanding of complex numbers may be necessary for some readers. Complex Numbers The mere mention of complex numbers can strike fear into many non-mathematicians. The word complex in their name is unfortunate because, rather than meaning complicated, this actually reflects the fact that such numbers consist of many parts. Complex numbers are nothing more than a mathematical tool, and are not necessarily a key part in the Mandelbrot process, as the circle method shows. A complex number actually comprises only two parts, called the real and imaginary parts. A typical complex number is shown below, at first it may look like an algebraic equation, but it actually tells us that it is a complex number consisting of four real parts and two imaginary ones: `4+2i` Other examples of complex numbers are: ```6-4i 6 real and -4 imaginary parts -2+2i -2 real parts and 2 imaginary ones 5+i 5 real and 1 imaginary part (i 1i) 0.745+0.113i 0.745 real and 0.113 imaginary parts``` The i not only allows us to distinguish the imaginary part from the real part, but also represents the square root of -1 (see Appendix B for details of indices and roots). If you try to calculate the square root of -1 on a pocket calculator using the usual method for square-rooting numbers you will probably be presented with an error message. This is because the roots of negative numbers are impossible to calculate. Despite this the root of -1 is in fact quite useful in mathematics so mathematicians are not satisfied with the error message on their calculators and instead have assigned the letter i, for imaginary, to it. Several clever things can be done with this number. For example if it is multiplied by itself the answer -1 is obtained, just as multiplying the square root of 4 by itself gives 4. This helps us to manipulate complex numbers, as real numbers can be produced from the imaginary parts. ```SQR(4) * SQR(4) = 4 SQR(-1) * SQR(-1) = -1 i * i = -1 because i<186>SQR(-1)``` Adding complex numbers is a simple case of gathering together the real and imaginary parts into one number, e.g. ```4+5i + 2+3i = 4+2 + 5i+3i = 6+8i 2-i + 3+3i = 2+3 + -i+3i = 5+2i``` Subtraction is similar, except that compatible parts are taken away from each other instead of being gathered together, e.g. ```10+15i - 5+5i = 10-5 + 15i-5i = 5+10i 6+3i - 4-2i = 2+5i (two minuses give a plus, 3--2 = 3+2 = 5)``` Multiplication in the Mandelbrot equation is limited to multiplying a number by itself, often referred to as squaring the number. In general, using the variables a and b to represent the complex number’s real and imaginary parts respectively, we can use the equation below to determine the result of squaring a number: `(a+bi)^2 = a^2 + (2*a*b)i - b^2` An example of this formula in action is shown below for the calculation of (4+2i)^2. ```a=4, b=2, therefore the answer is 4^2 + 2*4*2*i - 2^2 = 16 + 16i - 4 = 12 + 16i``` because Amiga BASIC, like most other computer languages, cannot deal directly with complex numbers this simple formula will prove invaluable in Mandelbrot related programs. Applying Complex Numbers to the Equation now that we can add, subtract and square complex numbers we can set about using them in the Mandelbrot equation shown below. The only unknown is which values to put in for c and z. `_z~new~_ = z^2 + c` In the process of plotting the Mandelbrot set, this equation is used to test various complex numbers between -2-2i and 2+2i, to see if they do one of two things, either they become very large very quickly, or they stay small. A number is tested by placing it into the equation as the constant, c, with z initially set to zero (0+0i). An example of a number (-1+0.5i) under the first few iterations of the test is shown below. ```_z~new~_ = z^2 + c = 0^2 - 1+0.5i = -1+0.5i z= _z~new~_ _z~new~_ = z^2 + c = (-1+0.5i)^2 - 1+0.5i = (-1)^2 + (2*(-1)*0.5)i - (0.5)^2 - 1+0.5i = 1 - i - 0.25 - 1+0.5i = -0.25-0.5i z = _z~new~_ _z~new~_ = z^2 +c = (-0.25-0.5i)^2 - 1+0.5i = (-0.25)^2 + (2*(-0.25)*(-0.5))i - (-0.5)^2 - 1+0.5i = 0.0625 + 0.25i - 0.25 - 1+0.5i = -1.875+0.75i``` The number is categorised into one of the behavioural groups described in the last paragraph by taking the modulus of the resulting complex number after each stage of the calculation. The modulus of a complex number is a simple way of describing its size using normal numbers, if the number is represented as a+bi, the modulus is the square root of (a2+b2). It can be shown that if the modulus of the complex number, z, becomes greater than two then it is destined to become larger very quickly. If the number’s modulus stays below 2 for an arbitrary number of iterations, say 32, it is assumed that the number is destined to stay small. This method of categorisation is applied to a regularly spaced selection of complex numbers between -2-2i and 2+2i, and all or these which have been assumed to be of the 'staying small' type are plotted on a plane of real against imaginary parts, such as that shown in Figure 4.9. Figure 4.9: The 4x4 complex plane Note that because BASIC cannot manipulate complex numbers directly they have to be explicitly represented as two separate variable parts as follows: `z = p + qi` `_z~new~_ = _p~new~_ + _q~new~_i` `c = a + bi` Incorporating this, and the rest of the process, into a program gives Listing 4.3, which when executed produces the Mandelbrot set. `` Unresolved directive in chapter4.adoc - include::../../src/amiga/4.3[]`` Listing 4.3: Program to produce the Mandelbrot set This listing is exactly the same as the one derived from the circle method (Listing 4.2), proving that both methods lead to the same result. Even the WHILE line is the same, because the equation for the distance of a point in the circle from the origin is the same as the one for the modulus of a complex number. There is not even a speed difference between the programs. However modern C compilers have extensive complex number handling routines as standard, making this alternative more attractive from a performance and clarity viewpoint. Manipulating the Mandelbrot Set The Mandelbrot set is an incomprehensibly complex object, so a single plot cannot possibly contain enough detail to show all the intricacies of the set. For this reason it is useful to be able to mathematically manipulate the set using the Amiga, something which can easily be done in Amiga BASIC. Zooming In Most programs written to display the Mandelbrot set allow the user to zoom in and pan around the set, so as to examine certain parts in more detail. This is easy to do using either method, but from now on, for convenience, the process will be discussed in the circle context only. Figure 4.10: 24 pin print-out of the contoured Mandelbrot set Assume that we wanted to zoom in on the section shown in Figure 4.10. By intelligent guess-work, or by measuring the diagram, it can be shown that the area in question lies in the region of the plane where 0.25⇐a⇐0.5 and 0.5⇐b⇐0.75. This is all the information we need to magnify this section, and only three minor alterations need to be made to Listing 4.2 in order to do it. The first is to make the program only test points in the selected region, rather than in the whole plane, in order to do this the two FOR…​NEXT loop initialisation lines should be altered to read as follows: ```FOR a=0.25 TO 0.5 STEP 0.01 FOR b=0.5 TO 0.75 STEP 0.02``` This has the desired effect of only showing the relevant area, but the points are plotted in the same place as before, rather than filling up the whole screen. The movement of the points is done by editing the line which draws <$ob>the points. The relevant replacement for the PSET line is:

`      PSET(a*1600-240,200-(b*800-400))`

Running the program now will display the right section at the relevant size, but the pixels are very spaced out. This is because the STEPs of the FOR…​NEXT loops were optimised for the full, four by four, Mandelbrot set. Because the height and width of our section is only one sixteenth the size of the total Mandelbrot the STEP sizes must be reduced to this proportion. The revised FOR…​NEXT initialisation lines are:

```FOR a=0.25 To 0.5 STEP 0.000625
FOR b=0.5 To 0.75 STEP 0.00125```
Figure 4.11: Mathematical enlargement of section from figure 4.10 (24 pin). amin=0.25, amax=0.5, bmin=0.5, bmax=0.75.

Figure 4.11 shows the kind of output that you can expect from the altered program. A further enlargement is shown in Figure 4.12. This can be created using the general purpose Mandelbrot program, given as Listing 4.4, in which the values of the <F102>amin, amax, bmin<F> and <F102>bmax<F> variables may be altered in order to display any part of the set.

Figure 4.12: Mathematical enlargement of section from figure 4.11 (24 pin). amin=0.34, amax=0.38, bmin=0.63, bmax=0.67

The section shown in Figure 4.12 represents a 100x magnification of Figure 4.1. If any higher magnification factors are used the maximum number of iterations carried out before assuming the number will never leave the circle (<F102>max_iteration<F> in Listing 4.4) will have to be increased from 32. This is because, as more magnification is performed, more contour lines become visible, so more calculations must be done to resolve points on the contours from points which actually belong to the 'official' Mandelbrot set. When the magnification level is particularly high, and the numbers involved are therefore smaller, it is necessary to use an accurate number rounding algorithm instead of the simple INT function provided by Amiga BASIC. An example Amiga BASIC function definition for rounding numbers is provided in appendix A, in GFA BASIC the built-in ROUND() function is accurate enough.

There is no simple formula for calculating the maximum number of iterations for a given part of the set, because the contours are of varying widths. As a rule of thumb the maximum iteration number is good if the border of the actual set is jagged and free of smooth lines. Naturally the lower this number is the better, as the set will take less time to generate.

There is one other important point to make about magnifying parts of any fractal, which is that the aspect ratio (the ratio of height to width) should always be preserved. This is the reason why the sets plotted by Listings 4.2 to 4.4 only take up the middle section of the screen. Extending it to fill the rest of the space would cause it to be artificially elongated. Such tampering with the aspect ratio causes inaccurate plots to be produced, making self-similarity hard to distinguish.

Self-similarity is much in evidence in Figure 4.11, which is reassuring proof that the Mandelbrot set is a fractal, in case you didn’t already know! The tiny replicas of the Mandelbrot shape are thought to be found even when the set is infinitely magnified, at a lecture Professor Mandelbrot exhibited a picture showing a section of the set magnified by a factor of 10^23 (100,000,000,000,000,000,000,000 times its original size), and still the original Mandelbrot image could be seen.

Most Mandelbrot exploration programs allow the user to select the area to magnify using the mouse, and indeed a program for this purpose is included on this book’s support disk. To save repetition, the techniques used in such a program are not explained here. However a fully documented example for similar exploration of Julia sets is included in the next chapter.

``  Unresolved directive in chapter4.adoc - include::../../src/amiga/4.4[]``
Listing 4.4: General purpose colour Mandelbrot program

In his book, The Fractal Geometry of Nature (see bibliography at the back of the book for details), Mandelbrot discusses the likeness between the boundary of the Mandelbrot set and the coastline (land/sea boundary) of an island. Comparing the set (Figure 4.6) with a map of Great Britain does reveal that although the Mandelbrot is symmetrical and contains near-perfect circles, there are similarities. For instance both have thin cracks running into them and thin ligaments coming out, known to geographers as rivers and peninsulas respectively. By far the most stunning feature, however, is the similarity between the jagged nature of the Mandelbrot boundary and the coastline.

We have already seen that after multiple magnifications the Mandelbrot set maintains the same complexity, and theoretically the same is true of the British coastline. If you were to parachute to the ground from a considerable height the amount of 'crinkliness' visible along the coastline would be constant during the whole descent. This is true because more detail is visible as the coastline gets nearer, but at the same time the field of view is reduced. Even when the parachutist lands the level of visible crinkliness will be the same (if he lands near enough to the coastline to see it) because each tiny piece of rock has an outline of equal detail (in relation to its size) to that of the British coastline. Because of this feature coastlines are classed as fractals and researchers have even determined a fractal dimension, between one and two, for the British coastline. This dimension indicates that the coastline is something between a line (one dimensional) and a plane (two dimensional).

Other Ways of Displaying the Mandelbrot Set

So far the Mandelbrot set has only been plotted on a two dimensional plane, with contour lines showing the ease with which points left the circle. However, this is not necessarily the best way to observe the set, and other methods are available. The first of these is to plot the contours in a more map-like style, using thin monochrome lines rather than wide bands of colour.

Figure 4.13: Geographical map showing variation in land height (relief)

A typical geographical map section is shown in Figure 4.13, here the lines are used to represent height, with every point on a particular line being at the same height. This allows the height of any particular area to be found, for example it can easily be seen that the point A is at a height of over 40 metres above sea level.

Figure 4.14: Anaemic Mandelbrot set

A similar diagram of the Mandelbrot set is shown in Figure 4.14, in this case the lines represent the ease with which the points on the plane left the circle. An 'anaemic' Mandelbrot set of this type can be created using the program given in Listing 4.5. This is basically a re-working of Listing 4.2 which, instead of colouring groups of points depending on the ease at which they left the circle, plots lines to show the boundaries between these groups.

``  Unresolved directive in chapter4.adoc - include::../../src/amiga/4.5[]``
Listing 4.5: Program to produce the anaemic Mandelbrot set

The boundaries between these groups are detected using a fairly simple algorithm. In the simplest terms a point is plotted if the number of iterations required to free it from the circle is different to that of the points adjacent to it. Because the program traverses the screen from left to right, bottom to top, any point being plotted (providing that it is not right at the edge of the plane) will have the point adjacent to it on the left and the point below it already calculated. To save time these two pre-calculated points are the only adjacent points tested against the point being investigated. However, the iteration numbers associated with these adjacent points were not stored in our original program, because such information was not needed, and it would have required an array with a 80,000 integer capacity.

In Listing 4.5 a 400 element array is set up at the start of the program to store every value in a particular column. This array is updated whenever a new point is calculated, just before the end of the FOR…​NEXT loop. This ensures that for any point the iteration numbers of the two adjacent points are stored in the array, and that as soon as an array element is no longer needed it is overwritten with the latest calculation. This may seem rather complicated at first but a look at Listing 4.5, together with its action, should rectify the situation.

Enhancements to the Mandelbrot Programs

Because it is difficult with colours alone to ensure that each contour around the Mandelbrot set a unique colour it would be useful to have another way of representing these contours. A popular method for doing this is to draw an isometric three dimensional Mandelbrot, known as a Mandelbrot landscape, where points are elevated from the plane by differing amounts depending on the ease with which they left the circle. Because isometric drawing is an essential element of this technique 3D Mandelbrots are discussed in the chapter on fractal landscapes (Chapter 7).

Internal Structure

In the figures shown in this chapter the area inside the Mandelbrot boundary has always been coloured black, but it is actually possible to colour the inside of the set with contours similar to those found outside. Although these internal contours are derived from the Mandelbrot process they have little importance to the study of the set, and are simply included to add visual effect. I do not intend to discuss internal structure in detail, but Figure 4.15 has been included to give an in idea of how it manifests itself in the Mandelbrot set.

Figure 4.15: Mandelbrot set with internal structure

5. Julia Sets

By looking at the Julia set plots included in this chapter it is clear that there are many similarities between these and the more familiar Mandelbrot set. This is because the same equation is employed to draw both fractals. It is merely the way in which it is applied that make them different. This close link between the two fractals means that Julia sets provide just as much scope for experimentation as their more famous counterpart.

<$ha>The Julia Process Both methods used for drawing the Mandelbrot set, shown in the last chapter, are suitable for drawing Julia sets. However, the circle method will be used here in order to avoid the intricacies of complex numbers. In Chapter 3 it was suggested that the Lorenz attractor was based on three equations containing three variables, x, y and z, so the system could be treated as being three dimensional. In Chapter 4 we found that to plot the path of the point (p,q) with initial position (a,b) in the Mandelbrot circle the following pair of equations was required: pnew = p^2 - q^2 + a qnew = 2*p*q + b This system clearly has four variables (pnew and qnew are not counted as they become p and q on subsequent iterations). Following the theory of the Lorenz equations it is clear that the Mandelbrot set can be treated as a four dimensional object. While it is challenging enough to represent three dimensional objects on a flat screen, displaying a four dimensional object is practically impossible. To overcome this problem only two variables from the Mandelbrot equation are used to plot the set. These are a and b which are used to position each point - p and q are used in the calculations but are later discarded. <$fc>Figure 5.1: Comparison of Julia and Mandelbrot processes

Figure 5.1 shows the relationship between the line in the circle for a Mandelbrot set, and next to it the relationship between a similar line in a circle and a Julia set. In the case of the Julia set the initial position of the point in the circle (a,b) is kept constant for the whole set, and the variable initial values of p and q are used to determine the position of the point to be plotted. Apart from this subtle change the rest of the process is the same for both fractals, including the test for escape from the circle, the method of determining the colour of contours and the method used to magnify a section of the set.

You may have noticed that Julia sets keep being mentioned in plural form. This is because there are an infinite variety of possible sets which can be produced using the Julia process, each with a unique initial position constant (a,b). In most chaos literature this constant is expressed in the equivalent complex number form of a+bi.

The Amiga BASIC program in Listing 5.1 draws the Julia set with constants a=-1.16 and b=-0.25 (-1.16-0.25i).

<$ob>DEFDBL a,b,p,q,x,y,pnew,qnew <$ob>DEFINT iteration <$ob>' <$ob>REM Set constants (easily editable for any set)<$ob>a=-1.16<$ob>b=-.25<$ob>'<$ob>REM Plot set<$ob>FOR x=-2 TO 2 STEP .01<$ob> FOR y=-2 TO 2 STEP .02<$ob> p=x<$ob> q=y<$ob> iteration=0<$ob> WHILE p*p+q*q<4 AND iteration<33<$ob> pnew=p*p-q*q+a<$ob> qnew=2*p*q+b<$ob> p=pnew<$ob> q=qnew<$ob> iteration=iteration+1<$ob> WEND<$ob> IF iteration=33 THEN<$ob> PSET(320+x*100,100-y*50)<$ob> END IF<$ob> NEXT y<$ob>NEXT x <$fc>Listing 5.1: An Amiga BASIC program which can be used to draw any Julia set

In many respects this program is similar to the equivalent Mandelbrot example given in Listing 4.2. Because the Julia plane extends from -2 to 2 in both directions, as in the Mandelbrot set, even the scaling is the same. The main differences are that a and b are set as constants at the start of the program and p and q, instead of being set to zero at the start of the circle process, are given values (x and y) derived from the nested FOR…​NEXT loops. This program can be made to draw any Julia set simply by changing the values of the constants a and b, some interesting values are suggested in Table 5.1.

<$tc>Table 5.1: Some suggested Julia set constants <$table3>a b Equivalent complex number <$table3>-1.16 -0.25 -1.16-0.25i <$table3>0.32 0.04 0.32+0.04i <$table3>-1.25 -0.01 -1.25-0.01i <$table3>0.00 -1.00 -i

Contours can be added to Julia sets using the same methods used in the previous chapter for the Mandelbrot set, where the contours are created by colouring points depending on the ease with which they left the circle. To save you referring back to Chapter 4, the possible replacements for the IF…​END IF section of the program are as follows:

For 16 colours (low resolution only):

<$ob> COLOR iteration MOD 16 !Set colour <$ob> PSET(320+x*100,100-y*50) !Plot point

For monochrome contours (high or low resolution):

<$ob> IF iteration/2<>INT(iteration/2) THEN !If iteration umber is odd…​ <$ob> PSET(320+x*100,100-y*50) !…​plot the point <$ob> END IF Like the Mandelbrot set, Julia sets can also be drawn as 3D landscapes (see Figure 5.4), or with internal structure. <$hb>Derivation from the Complex Number Method

Again, the method is similar to that given for the Mandelbrot set. The equation shown below still holds but c is kept constant the whole time, rather than being different for each point on the complex plane.

znew = z^2 + c

Because c is constant the position of the point is determined by the initial value of z which, instead of always being zero (0+0i), is different for each point on the plane. The Julia set drawn is dependant on the complex constant, c, whose real and imaginary parts are equivalent to the circle method’s a and b variables respectively.

<$ha>A GFA BASIC Julia Set Explorer In general I have attempted to write the programs for this book in as compact a form as possible. However, there are programs available which allow you to explore fractals with all the convenience offered by graphical user interfaces, disk filing operations and printout facilities. This type of program tends to be very long and hence is unsuitable for a book primarily written as an introduction to the theory of chaos. However, as an insight into how such programs are created, a Intuition-based example to generate and explore Julia sets is presented in this chapter. The program is written in GFA BASIC, as this is the more suited to this type of task than Amiga BASIC. This example will prove especially useful if you want to develop your own fractal exploration programs, based on processes given in this book or elsewhere. The full Julia program is shown in Listing 5.2. This is provided mainly for educational purposes and to allow important techniques to be gleaned. Due to the length of this listing it would be irresponsible for me to recommend typing it in. If you are interested in seeing the program in operation I suggest you purchase the support disk instead as this contains the entire program. Planning Before actually tackling a project of this sort a clear set of objectives and a plan for developing the program must be devised. The objectives for the Julia program are that it should: <$bullet>Be able to plot (quickly) any Julia set requested by the user

<$bullet>Allow zooming (of any magnification) into any section of the set <$bullet>Give options for filing sets and sections to and from disk

<$bullet>Display parameters in a form compatible with other literature <$bullet>Allow for the filing of such parameters for later use

<$bullet>Give the user the option to dump images to a printer <$bullet>Be implemented in an intuitive graphical environment.

This list basically brings together the Julia set facilities discussed above and the features that would be expected from a commercial art package or word processor. One of the most important points here is the parameter saving option as this would allow (for example) a Julia landscape program to draw three dimensional versions of sets already plotted and saved in two dimensions.

The intuitive graphical environment mentioned in the list must be Intuition, as it is resident in the Amiga operating system and is well supported by GFA BASIC. The use of Intuition in the Julia program comprises a pull down menu bar, alert requesters, simple windowing and mouse input.

A discussion of the menu bar (given later) provides a more detailed breakdown of the options available because, in a program of this sort, all the options must be selected either directly or indirectly from the menu bar.

Universal Programming Techniques

There is no escaping the fact that Julia set generators must perform many time-consuming calculations. For this reason a number of steps have been taken to ensure that the Julia program is as fast as possible. The most general is the application of structured programming techniques, including the use of local variables wherever possible. Unlike global variables, these are local to a particular procedure and are lost when that procedure ends, which means that at any time there are fewer variables for BASIC to keep track of than there would be if all the variables were global. Also, when dealing with whole numbers, integer variables have been employed instead of the floating point ones selected by default. Integer variables are signified by having their names suffixed by a per cent character (%) and are processed significantly faster than their floating point counterparts. Of course DEFINT could be used to declare variables as integers, but the suffix notation makes it possible to see instantly what type a particular variable is.

I have used comments and long procedure and variable names throughout the program to make it easier to understand but, when running under the interpreter, these will slow the program’s execution. Obviously the comments may be omitted, but shortening names can introduce errors if sufficient care is not taken. If you intend to compile the program the long names and comments will not matter because these are automatically stripped out when creating the executable file.

<$hb>The Menu Bar Because each menu option is consecutively numbered it is imperative that all possible options are considered before incorporating the menu into a program. The conventional way of organising menu bars is for the file, or project menu to come first, followed by a series of program-specific menus. The file menu always contains file manipulation options, some of which may be program-specific, followed, after a blank line, by an option to quit the program. These conventions have been adhered to in the Julia program, for which the full menu structure is shown in Figure 5.5. The first program-specific menu is the picture menu, whose options are fairly self-explanatory. The colours menu is used to select which of the three methods will be used to display the set, so the three options are mutually exclusive. The action of each option is: Mono: The set is drawn with monochrome contours (default) Set Only: Only members of the actual set are drawn, in black Colour: The set is drawn with colour contours. A tick is placed next to the menu item representing the selected colouring method. To save accidental loss of pictures, a tick is also used in the file menu (next to the save option) to inform the user of whether or not the set in memory has been saved. GFA BASIC stipulates that the menu structure must be stored in an array with each menu item in a separate element of the array, placed consecutively in the order in which they are to appear. An empty string ("") is used to mark the end of each menu, and a pair of empty strings mark the end of the last menu (the end of the menu structure). The information for the Julia menus is placed in an array called <F102>menu$<F> in the <F102>initialise<F> procedure (see Listing 5.2). This array could have been filled using more compact code, but by using the longhand method it is easy to determine how the menus are arranged and what number each menu item has, merely by glancing at the program.

A further array, called <F102>status%<F>, is used to store information about the status of menu items - whether they are disabled for example. This information must be stored in this way because the built-in GFA BASIC menu routines lose such information when the menu bar is switched off (as it frequently is in this program). The status of a menu item is set using a command of the form:

<$ob>status%(item_number)=flag where <F102>item_number<F> is the number of the menu item and <F102>flag<F> is one of the following values: <$table2>80 Enabled, no tick

<$table2>336 Enabled and ticked <$table2>0 Light (disabled)

The default status for some items is set just after the initialisation of the menu structure. Of course the contents of the array are not passed to the GFA BASIC menu handler automatically. This job is done by the <F102>menu_status<F> procedure which uses a version of the MENU command to set the menu item status. This procedure must be called whenever the menu is switched on after being off, or whenever an item’s status is altered while the menu is turned on.

<$hb>Menu Handling The menu bar is set into action near the end of the <F102>initialise<F> procedure. Here the menu is turned on, given the latest status information using a call to <F102>menu_status<F>, and instructed to call the <F102>branch<F> procedure upon a menu item being selected. After the initialisation procedure is complete an endless DO…​LOOP loop is entered which continually checks for menu activity using the SLEEP command. On selection of an item (i.e. when the user has pulled a menu down and clicked on one of the options) the <F102>branch<F> procedure is automatically called by the GFA BASIC menu handler. When called the <F102>branch<F> procedure immediately switches off the menu bar to stop confusion occurring due to options being selected while other options are being carried out. The procedure then branches off to one of 11 different procedures depending on the item selected (the item’s number is returned by the GFA BASIC menu handler in element zero of the built-in <F102>menu()<F> array). On returning from the selected procedure the menu is switched back on and its status updated by calling the <F102>menu_status<F> procedure. The <F102>branch<F> procedure then returns control to the menu testing loop at the head of the program. The various procedures called from <F102>branch<F> perform all the program’s main operations and are discussed below. <$hc>Desk_About

This very simple procedure is called to display an alert requester proclaiming the name and source of the program whenever the about this program option is selected from the file menu. The procedure first initialises the local integer variable, <F102>button%<F>, before calling <F102>bar_text<F> to display a description of the action selected in the main window’s title bar. After this and a quick piece of house-keeping the procedure’s main function, the creation of the alert requester, is carried out. This is performed using the GFA BASIC ALERT command for which the syntax is:

<$ob>ALERT dummy,message$,default,button$,button% where: <F102>dummy<F> is a dummy variable, ignored in Amiga versions of GFA BASIC up to at least 3.5 <F102>message$<F> holds the message to be displayed, in which the vertical bar character (|) is used to mark the beginning of a new line).

<F102>default<F> holds the number of the default button, which will be displayed with a thickened outline and automatically selected if the <Return> key is pressed while the alert requester is on the screen. It is possible to have no default button by placing a zero here.

<F102>button$<F> holds the text for one, two or three buttons, separated by vertical bars (|). <F102>button%<F> is a variable in which the number of the chosen button is returned after it has been selected. Note that the buttons are numbered consecutively as they appear in <F102>button$<F>, starting from 1.

The ALERT command has been explained in detail because it is widely used in the Julia set program for conveying messages to the user. It is particularly useful in a program of this sort because the screen area obscured by the alert requester is automatically restored when the requester is removed.

The structure of this procedure is similar to all procedures called as a result of a menu item being selected, with the declaration of local variables and the calling of <F102>bar_text<F> prior to performing the selected actions.

<$HC>Dummy Calls to <F102>dummy<F> are included in the ON GOSUB line of the <F102>branch<F> procedure to fill in the gaps where no procedure needs to be called (e.g. where unselectable dotted lines are included as menu options). As careless typing could cause <F102>dummy<F> to be called instead of the correct procedure, a procedure definition for <F102>dummy<F> has been included which displays an alert requester warning of the error. File Operations Creating, or zooming into, a Julia set takes a relatively long time and it is therefore useful to be able to save the image on the screen to disk. <$fc>Figure 5.3: Julia set with a=0.32, b=0.043 (24 pin)

<$fc>Figure 5.4: Julia landscape produced using techniques from Chapter 7 <$fc>Figure 5.5: The menu structure of the Julia set exploration program

The actual format used for the Julia files will vary depending on what you intend to do with them, obviously a standard format such as IFF would be ideal for importing into other programs but this would make it hard to incorporate the constants for the set in the file. The reward for including the constants is that it is then possible for the program to tell exactly what part of what Julia set has been loaded.

The <F102>file_load<F> and <F102>file_save<F> procedures provided in Listing 5.1 incorporate only the necessary user interface routines, but they have been written in such a way that it is a simple task to include your own file handling code.

<$fc>Figure 5.6: Julia set with a=-0.75, b=-0.11 produced by Listing 5.2 <$hc>File_load

If the status of menu option 3 is set to 16 (if the Julia set currently in memory has not been saved) confirmation of the menu selection is sought using an alert requester. This step is necessary as loading a new file overwrites the set previously in memory. If the user confirms the intention to load a new file a file selector is displayed using the FILESELECT command, the syntax for which is:

<$ob>FILESELECT title$,buttontext$,searchpath$,selected$where: <F102>title$<F> holds a message to be displayed at the head of the file selector window.

<F102>buttontext$<F> provides the text which is to be placed on the confirmatory button. <F102>searchpath$<F> contains the path to the directory (drawer) in which the file is to be found and also provides a suggested file to load. A typical searchpath string might be <F102>"DF0:/CHAOS/Untitled"<F> which would cause the file selector to display all the files in the ram disk’s CHAOS directory and would provided the name Untitled as a file suggestion.

<F102>selected$<F> is a variable in which the full name (including the path) of the selected file is returned. An empty string is returned if the cancel button is hit. The relevant picture loading routine should be included at the point indicated in the listing so that once a file has been successfully selected the procedure goes on to open it for input and load the picture. After successfully loading the file the zoom in and view parameter menu options are enabled, and a tick is placed next to the save option (obviously if the file has been loaded from disk it must already have been saved). If the file does not exist an alert requester is generated to inform the user of the problem. <$hc>File_save

This procedure is very similar to <F102>file_load<F> but many of the checks are not performed as we know that the file will exist. Again, the actual file handling code is missing but the basic structure is in place. After the file has been successfully saved a tick is placed next to the save menu item (item 3) by setting the 3rd element of the <F102>status<F> array to 1. Note that there is no disk full or write protection check, and any file of the same name as the one being saved could be overwritten without warning.

<$HC>File_delete It is sometimes necessary when saving Julia plots to delete old files to make room for the new one. This can leave the user in a bit of a predicament if a delete menu option is not included. The operation of this procedure is very simple, relying on the FILESELECT command discussed above for file selection and using the KILL command to delete the selected file. <$HC>File_Quit

This procedure is called when the user indicates the desire to leave the program, by selecting the quit option from the file menu. Because of the severity of this action a confirmatory alert requester is displayed before returning to the editing screen (returning to the Workbench screen if the program has been compiled).

<$HC>Picture_julia Apart from loading a file from disk the only way to get a Julia set on the screen is by selecting the plot Julia set option from the picture menu. <F102>Picture_julia<F> is the procedure called when this option is selected. If there is an unsaved set on the screen when the option is chosen the program produces a confirmatory alert requester before going on to set the <F102>xmin, xmax, ymin, ymax<F> and <F102>max_iteration%<F> variables to their default values for a full set plot. A window is then displayed in the bottom left of the screen in which the user must enter the two parts of the constant required to draw the set. The window is opened with the <F102>open_window<F> procedure and closed after use by a call to <F102>close_window<F>. The save, zoom and view parameters menu items are then enabled, as they now represent viable options, and the <F102>plot_set<F> routine is called to draw the set. <F102>Plot_set, open_window<F> and <F102>close_window<F> are described later. <$HC>Picture_zoom

In many respects this procedure is similar to <F102>picture_julia<F> because it initiates the drawing of a Julia set using the <F102>plot_set<F> routine. The differences are that the constant need not be entered because it is inherited from the set already on the screen, but the dimensions of the set and the iteration ceiling are selected to determine the section of the set instead. The user selects the area to enlarge simply by moving to the relevant area, holding the left mouse button down and dragging an 'elastic' rectangle out from the point until the desired area is enclosed (pressing the right hand button will cancel the process at this stage).

This action is simple from the user’s point of view but, behind the scenes, there is a lot of activity. As well as animating the rectangle in such a way that it does not distort the set already on display the program also keeps the aspect ratio of the rectangle the same as that of the screen. Therefore the magnified section of the set does not have to be artificially elongated when plotted. Once the rectangular area has been selected a series of calculations is performed to determine the values of the <F102>xmin, xmax, ymin<F> and <F102>ymax<F> variables for the new plot and a value for <F102>max_iteration%<F> (the iteration ceiling) is determined. The value of <F102>max_iteration%<F> is increased proportionally to the level of magnification because the more detailed the plot the more calculations are required to distinguish between different types of points.

The rule used to calculate the iteration ceiling was arrived at empirically because the complex nature of Julia sets prevents us from ever knowing accurately the amount of detail likely to be visible in the enlarged section. Also the result of <F102>(max_iteration% MOD 16)<F> must always be zero so that the set is black in the centre, in line with the accepted conventions.

<$HC>Picture_parameter As well as displaying the constant for the set currently on the screen (and its position on the Julia plane) this procedure also displays the iteration ceiling, in an editable text field. This latter option is very useful as it allows the user to alter the level of detail resolved, either to speed up plotting or to by-pass the somewhat inaccurate method of calculation used when zooming into a set. The procedure uses <F102>open_window<F> and <F102>close_window<F> (both discussed shortly) to manipulate the window in which the information is displayed. Two notable GFA BASIC functions are used in <F102>picture_parameter<F>: PRINT USING which is employed to format the floating point numbers, and FORM INPUT which places a string on the screen, allows the user to edit it and then returns the result. Note that some conversion between strings and numeric variables (using STR$ and VAL) must be performed here because the variable being edited, <F102>max_iteration%<F>, is numeric and FORM INPUT will only accept strings.

<$hc>Picture_print This procedure is called when print is selected from the picture menu and, after seeking confirmation from the user, it dumps the current screen image to a printer using the HARDCOPY command. HARDCOPY is a built-in GFA BASIC command which prints in whatever graphics format is selected in your Workbench preferences file. <$hc>Colour

The <F102>colour<F> procedure deals with all three options of the colours menu, and is used to select the way in which the set is displayed. The structure of the procedure differs slightly from that of others called from the menu handler, in that no message is placed in the menu bar because the speed of the procedure is such that the user would not have time to read it. After initialising the local variables a FOR…​NEXT loop is entered in which the selected option has a tick placed next to it by amending the relevant <F102>status%<F> entry. The other two options have their status set to zero (unticked). This loop ensures that only one option can be ticked at a time. Note that this procedure in no way alters the way in which the set is drawn, it merely changes the status array to allow the plotting procedure to know which option should be used the next time it is called.

Procedures not Called Directly from the Menu Handler

Routines to perform common operations, such as printing messages in the menu bar, have been written as procedures so that they can be easily used in many different program locations. Such routines are described here.

<$HC>Plot_set This can be thought of as a universal Julia set plotter as it can draw any section of any set with any iteration ceiling, using any of the three main colouring methods. The procedure fetches the information about the set from the global variables, of which <F102>a<F> and <F102>b<F> describe the constant, <F102>xmin, xmax, ymin<F> and <F102>ymax<F> give the position in the set of the section to be drawn and <F102>max_iteration%<F> holds the iteration ceiling. The colouring method is determined by interrogating the <F102>status%<F> entry for each of the three colours menu options to identify which one is ticked. Note that two important variables, <F102>xincrease<F> and <F102>yincrease<F>, are calculated before plotting begins. Theses variable holds the size of a screen pixel in terms of the 4x4 Julia set plane and are used when converting the position of a point on the screen to the corresponding point on a Julia set. <F102>Xincrease<F> and <F102>yincrease<F> are a global variables because, as well as being used in <F102>plot_set<F>, they are used in <F102>picture_zoom<F> to convert the co-ordinates of the rectangle drawn by the user into the position of the co-ordinates of the Julia set section. <$hc>Open_window and Close_window

Together these provide all the necessary functions for manipulating the parameter/input window. It is relatively easy to open a window in GFA BASIC using the OPENW command, however, since a similar size window is required for two different parts of the program it seemed sensible to place the window opening and closing code in this pair of procedures.

6. Imitating Nature - Plants, Shrubs and Trees

Using chaos to draw abstract patterns such as the Mandelbrot set is certainly interesting, but many researchers are now turning their attention to the application of chaos techniques in mathematical generation of images depicting scenes from nature. It is relatively easy for an artist to copy such an image from the real world but generating a landscape, plant, animal or other natural object from scratch using nothing more than a set of rules and some mathematics is considerably more difficult. This challenging research not only produces pleasing pictures but also gives an insight into the processes which shape the world around us.

Ideally it should not be possible to tell the difference between a computer generated scene and a similar real life one. Due to the limited graphics facilities this ideal cannot be fully realised on the Amiga, but a useful insight into the relevant techniques and algorithms is given in this and the two following chapters concerned with imitating nature. This particular chapter deals with the re-creation of a variety of plants.

Before continuing it is imperative to point out that the algorithms provided in this chapter cannot be implemented in Amiga BASIC without considerable difficulty because Amiga BASIC does not support recursion (the calling of a sub-routine from within itself). The examples in this chapter are therefore given in GFA BASIC. Note, however, that it is relatively easy to convert the programs to Hisoft BASIC or to C.

What is a Plant?

Although every plant species, and every individual within that species, is different it is easy to distinguish between plants and other objects just by looking at them. The way humans do this is to recognise a particular set of characteristics peculiar to plants. If the object we are trying to identify exhibits most of the characteristics from this set we can be fairly sure that it is a plant, even if we have never seen a plant of the same species before. To draw plants which are recognisable as such we need to know what these characteristics are, and how they can be re-created. This is not something specific to computer generated plants, artists also work in this way.

By briefly observing a variety of common plants a short list of features can be constructed which describes the silhouette of a typical plant. The list reads as follows:

<$bullet>A single stem protrudes from the ground with many branches <$bullet>Most branches split into several smaller sub-branches until the smallest branches (e.g. twigs in trees) are reached

<$bullet>Each plant has a large number of branch ends <$bullet>The nature of the branching is consistent throughout the plant.

Any biologist could easily discredit this definition, by quoting numerous examples which contradict it, but it serves as a good visual guide for the majority of plants. Note that flowers, fruit and leaves are not mentioned so as to keep the list valid for any time of the year and for as many plant types as possible.

Some of the characteristics typifying plants have been mentioned elsewhere in this book when describing fractals, the most obvious example being the branching of the Feigenbaum diagram. Because the branching is consistent throughout the length of the fractal the Feigenbaum diagram is said to be self-similar. Plants can also be shown to be self-similar using the same argument. A fern, for example, is composed of many miniature fern shapes, each closely resembling the plant of which it is part. Similarly the tiny veins on the underside of a leaf give a good approximation of the structure of the parent tree because, as stated above, the nature of the branching is the same throughout the tree.

Also like our computer generated fractals, the complex structure of a tree is created from a very simple process. Plants must have a large leaf-to-size ratio in order to trap as much light as possible for energy-gathering during photosynthesis. The best way for a plant to achieve this is to make every branch split into many sub-branches to create as many branch ends, hence as many leaves, as possible. This simple process is actually very successful in creating a very large surface area on which light can be trapped. To demonstrate the difference in surface area between a tree and a similar geometrical shape imagine painting a Christmas tree. You would need much more paint to cover its intricate structure of branches and needles than you would to cover a smooth faced cone of the same size.

In having such a large surface area trees are like the Mandelbrot and Julia sets, whose boundaries use up more space than expected for the number of dimensions in which they exist. By continuing this train of thought it is possible to give real plants, which exist in three dimensions, a fractal dimension of just over three. Because of this, and their self-similarity, plants can be classed as fractals, and by thinking of them in this context they can easily be drawn using similar methods to those employed in the fractal generators of earlier chapters. The exact set of rules, or algorithm, for drawing trees on a computer can be determined simply by observing the branching structure of plants in detail. Fractal plants, unlike previous fractals in this book, are created using a rule-based iterative process, rather than a mathematical one. This means that the process is easy to understand, but incorporating it into a program can be a challenge due to the reluctance of computers to deal in anything other than numbers.

Describing a Plant

The basic branching algorithm is relatively simple, but before considering this in depth we must devise a method of describing the structure of a plant which the Amiga can deal with. Not surprisingly GFA BASIC doesn’t have built-in data structures for storing descriptions of trees! For demonstration purposes we can try to build up a description of the grass shown in Figure 6.1. This plant is referred to as a grass because it is simple in structure and only branches in one direction, making the principles easier to understand.

<$fc>Figure 6.1: Example of a simple grass The first thing that can be done to make the grass easier for the computer to handle is to split it up into straight segments of uniform length (see Figure 6.2). By doing this the branch lengths are said to be quantised, with the smallest possible branch being one segment long, and the length of every branch being a multiple of one segment length. Angles are also quantised for convenience, in steps of 45 degrees. Curved branches can be built up by putting several angled segments together. For more realistic, more detailed, plants shorter segments and smaller angles should be used, but because longer descriptions are needed to store such plants the values given above will be used in this initial discussion. <$fc>Figure 6.2: Segmented version of Figure 6.1 which the computer can deal with

Simple notation can now be used to describe such plants using just three characters. The string shown below for the grass in Figure 6.2 will be used as the basis for the following discussion.

"1[1]1[11]"

This is a normal ASCII text string, not the most efficient method of storage as there are only three possible characters, but it is easy to deal with in GFA BASIC and makes future expansion simpler.

Each branch segment is represented in the string by the one character (1). The square brackets are collectively used to describe the tree’s branches, where an open bracket ([) represents a 45 degree clockwise split from the current position (this usually represents the start of a branch) and a closed bracket (]) represents the end of a branch. Whole branches can easily be identified in such strings, as they are like miniature trees, with an equal number of open and closed brackets surrounding them. For example the [1] in the above description represents the first 45 degree branch.

Branches can also be given sub-branches. For example the longer branch ([11]) could be changed to [1[1]1], meaning that it had a single segment branching off at 45 degrees halfway along its length. The grass with this altered branch is shown in Figure 6.3.

<$fc>Figure 6.3: The grass with description "1[1]1[1[1]]" Curved branches may also be described using this notation. The 180 degree curve in Figure 6.4 is described using the string: "1[1[1[1[1]]]]" <$fc>Figure 6.4: 180 degree curve in steps of 45 degrees

This is a particularly good example of how an open bracket does not necessarily represent the beginning of a new branch. The versatility of the notation means that it is possible to describe any plant which only branches in one direction, no matter how complex it is.

Turning Descriptions into Drawings

In order to write a BASIC program to turn such strings into pictures on the Amiga’s screen we must define a method of interpreting each character into an action that the program can perform. To do this it is necessary to formulate a precise course of action for each of the three characters mentioned above. These actions are as follows:

<$hang1>1 Segment: Draw a uniform length line from the current position, at an angle equal to that held in the current angle count. The end of the new branch becomes the new graphics drawing position. <$hang1>[ 45 degree clockwise split: Add 45 to current angle count.

<$hang1>] End of branch: Return to the beginning of the current branch (i.e. where the matching open bracket occurred) and decrease the angle count by 45 degrees. The terms angle count and drawing position are best described using the turtle analogy on which the Logo programming language is based. Instead of the conventional approach of plotting points and lines at the graphics cursor, we imagine that a graphics turtle is being used. In this context a turtle is similar to a cursor in that it can be placed at any pixel position on the screen, but in addition to a position the turtle also has an angular direction associated with it (the angle count). The direction is usually specified in degrees with zero being vertical and positive angles representing clockwise rotation from the vertical. Negative values are used to indicate anti-clockwise angles. Note that a line drawn from the current turtle position does not need to have its destination specified, as this is automatically calculated from the length of the line and the angle count. Figure 6.5 shows the path that a logo-type turtle would take when plotting the plant described by the example string introduced above. Experienced GFA BASIC users will know that the language has a suite of built-in commands to facilitate Logo-style drawing. I have, however, refrained from using such features here in order to make conversion to other languages easier. <$fc>Figure 6.5: The path followed by a Logo turtle plotting the grass in Figure 6.3

The three actions have been incorporated into the <F102>draw_plant<F> procedure in Listing 6.1, which is used to draw the grass from Figure 6.3. Note that GFA BASIC does not support turtle graphics so the program uses simple trigonometry to convert between the actions and the DRAW TO command. Note that <F102>draw_plant<F> merely turns the string supplied in <F102>plant$<F> into a drawing, it doesn’t actually generate plants. <$ob>COLOR 1<$ob>plant$="1[1]1[1[1]1]" !Put sample string into plant$<$ob>unit_length=10 !Set up constants for unit length…​<$ob>unit_angle=PI/4 !…​and unit angle<$ob>'<$ob>GOSUB draw_plant !Call procedure to draw the plant <$ob>'<$ob>PROCEDURE draw_plant<$ob> angle_count=0 !Initialise angle count<$ob> pointer=0 !Set up initial string pointer position<$ob> GOSUB draw_branch(320,180)<$ob>RETURN<$ob>'<$ob>PROCEDURE draw_branch(x,y)<$ob> INC pointer !Increment pointer position and…​<$ob> character$=MID$(plant$,pointer,1) !…​note character in that position<$ob> '<$ob> REPEAT<$ob> SELECT character$<$ob> CASE "1" !If the character is a 1, draw a segment<$ob> PLOT x,y<$ob> x=x+SIN(angle_count)*unit_length*2<$ob> y=y-COS(angle_count)*unit_length<$ob> DRAW TO x,y<$ob> '<$ob> CASE "[" !If it’s a [ then branch off<$ob> ADD angle_count,unit_angle<$ob> GOSUB draw_branch(x,y)<$ob> ENDSELECT<$ob> INC pointer !Note next character<$ob> character$=MID$(plant$,pointer,1)<$ob> UNTIL character$="]" OR character$=""<$ob> '<$ob> SUB angle_count,unit_angle !If it’s a ] (end of branch) then<$ob>RETURN !decrease angle count and return <$fc>Listing 6.1: The <F102>draw_plant<F> procedure and sample calling code

The first section of the program initialises all global variables before calling the procedure to plot the plant. The variables are defined as follows:

<$hang1><F102>plant$<F> Holds the string describing the grass, in plant generation programs later in this chapter the string will be created by the program.

<$hang1><F102>unit_length<F> Is the uniform segment length in low resolution pixels. This and <F102>unit_angle<F> are variable so as to enable easy experimentation. <$hang1><F102>unit_angle<F> Determines the angle at which branches split off. In this example it is 45 degrees, but as GFA BASIC only accepts radians it must be set to the equivalent radian value of PI/4 (see Appendix B for details of conversions).

<$hang1><F102>angle_count<F> Is used by the <F102>draw_branch<F> procedure to store the current angle count, described above. <$hang1><F102>pointer<F> Is used to store the current position in the string, i.e where to look for the next character to interpret. It is initially set to 0 because the first thing that the <F102>draw_branch<F> procedure does is increment it (see below).

After initialisation the <F102>draw_branch<F> procedure is called with the arguments 160 and 180. These are the x and y co-ordinates of the current drawing position that initially define the position of the bottom of the plant’s stem.

The second section of the program is the actual definition of the <F102>draw_branch<F> procedure. It is not immediately obvious how this operates, but all the techniques used are perfectly legal and can easily be understood if analysed logically. The first task of draw_branch is to increment the string pointer and copy the character at that new position in <F102>plant$<F> into the <F102>character$<F> variable. A REPEAT…​UNTIL loop is then entered which checks for the open square bracket ([) and one (1) characters.

If the loop detects that the character is a 1 a segment of unit length is drawn at the current angle. This results in the end of the new line becoming the current drawing position. The next character is then extracted from the string, and if it is not a closed square bracket the loop then repeats.

If the character is not a 1, but an open square bracket, the <F102>unit_angle<F> value is added to the current angle count and the <F102>draw_branch<F> procedure is then called with the latest x and y arguments. If this happens the REPEAT…​UNTIL loop does not continue until the <F102>draw_branch<F> call returns.

If a closed bracket is taken from the string, or the end of the string is reached (<F102>character$=""<F>), the REPEAT…​UNTIL loop is terminated, <F102>unit_angle<F> is subtracted from the angle count and the procedure returns to where it was called. In GFA BASIC the only time a procedure can return cleanly is at its end. This is obviously true for <F102>draw_branch<F>. Usually the procedure will return to the open bracket clause of the IF…​ENDIF construct. However, in the case where the end of the string has been reached the procedure will automatically return to where it was first called, at the head of the program, resulting in the termination of the procedure and a return to the editing screen. The method used here of calling a procedure from inside itself is known as recursion and, although a very powerful tool, it is sometimes hard to follow. Its use not only means that the same procedure can be used to draw all manner of branches, from twig to whole plant, but also means that we do not have to waste time and memory keeping track of where branching points occurred on the screen. Each time the <F102>draw_branch<F> procedure invokes itself the previous group of arguments are stored by BASIC on a kind of internal stack. This means that as the procedure returns to where it was called the old values of <F102>x<F> and <F102>y<F> are taken off the stack and automatically become the current drawing position co-ordinates. It is easy to experiment with different grasses by just altering the description string at the start of the program. Note that any valid string may be used but plants represented by particularly long strings may not fit on the screen. In this case you should reduce the <F102>unit_length<F> and <F102>unit_angle<F> constants as necessary. Some of the examples given later suggest suitable values for these constants which vary from the customary 10 pixel and 45 degree assignments. The most important thing to ensure, when trying to create a natural looking plant, is that each open bracket has a matching closed bracket. String Generation Interesting though the plant drawing program is, it will not actually generate plant descriptions. If you want a detailed piece of grass you have to type in a long string, which does not necessarily look very realistic. It is a fairly trivial matter to write a procedure to generate a plant’s string description now that the necessary notational conventions have been established. Like other fractals, plants are generated using a simple, structure enriching, iterative process. The process used here is initially quite simple, but can be enhanced fairly easily. Basically every iteration sees each segment being replaced by a larger, more complex branch, thus making the plant larger while keeping the same level of relative complexity. This is achieved in the program by searching through the whole plant description (<F102>plant$<F>), replacing all the 1 characters with a more complex user-defined string, set to determine the type of plant.

Listing 6.2 contains the <F102>replace_chars<F> procedure to accomplish this task, along with the old plant drawing routine. In this example all 1s are replaced by the contents of <F102>one$<F>, initially set to "11[1[1[1]]]", although any valid plant-type structure may be used. The plant description with which the program starts is the simplest possible structure, "1". It is important to note that there is no easy way to insert characters into a string, so what the program actually does is copy the contents of <F102>plant$<F> into the variable <F102>newplant$<F>, character by character. However, if a one is found then the replacement string is copied across instead. Once the replacements are complete the contents of <F102>newplant$<F> are copied back to <F102>plant$<F> and <F102>newplant$<F> is cleared ready for the next iteration.

<$ob>COLOR 1<$ob>plant$="1"<$ob>one$="11[1[1[1]]]"<$ob>unit_length=2 !Set up unit length…​<$ob>unit_angle=PI/16 !…​and unit angle<$ob>'<$ob>REPEAT<$ob> CLS<$ob> GOSUB draw_plant !Draw plant so far <$ob> GOSUB replace_chars !Make plant grow<$ob>UNTIL MOUSEK=2<$ob>'<$ob>PROCEDURE draw_plant<$ob> angle_count=0 !Initialise angle count<$ob> pointer=0 !Set up initial string pointer position<$ob> GOSUB draw_branch(320,180)<$ob>RETURN<$ob>'<$ob>PROCEDURE draw_branch(x,y)<$ob> INC pointer !Increment pointer position and…​<$ob> character$=MID$(plant$,pointer,1) !…​note character in that position<$ob> '<$ob> REPEAT<$ob> SELECT character$<$ob> CASE "1" !If the character is a 1, draw a segment<$ob> PLOT x,y<$ob> x=x+SIN(angle_count)*unit_length*2<$ob> y=y-COS(angle_count)*unit_length<$ob> DRAW TO x,y<$ob> '<$ob> CASE "[" !If it’s a [ then branch off<$ob> ADD angle_count,unit_angle<$ob> GOSUB draw_branch(x,y)<$ob> ENDSELECT<$ob> INC pointer !Note next character<$ob> character$=MID$(plant$,pointer,1)<$ob> UNTIL character$="]" OR character$=""<$ob> '<$ob> SUB angle_count,unit_angle !If it’s a ] (end of branch) then<$ob>RETURN !decrease angle count and return<$ob>'<$ob>PROCEDURE replace_chars<$ob> pointer=0 !Set initial pointer position<$ob> REPEAT<$ob> INC pointer !Increment pointer position<$ob> character$=MID$(plant$,pointer,1) !Extract character<$ob> SELECT character$<$ob> CASE "1" !If it’s a 1 then replace it with one$<$ob> newplant$=newplant$+one$<$ob> DEFAULT<$ob> newplant$=newplant$+character$!Otherwise leave it unchanged<$ob> ENDSELECT<$ob> UNTIL pointer=LEN(plant$)<$ob> plant$=newplant$!Replace old plant$ with new one<$ob> newplant$="" !Clear newplant$just in case<$ob>RETURN <$fc>Listing 6.2: The <F102>draw_plant<F> and <F102>replace_chars<F> procedures being used together to generate and plot a grass If you type in and run the program a tiny vertical sprig of grass will appear at the bottom of the screen, which will then begin to grow. Figure 6.6 shows the appearance of the grass on successive iterations. Note that as the plant, and hence its string description, gets more complex the iterations take considerably longer to perform. If you have a compiler you may wish to use it here, although the calculations are still quicker than those required to plot the Mandelbrot set. <$fc>Figure 6.6: The appearance of the grass with replacement string "11[1[1[1]]]" during the first six iterations

There are no limits imposed by the plant algorithms regarding the number of iterations performed but complexity is limited because GFA BASIC can only cope with strings of less than 32,768 characters. If the plant’s description goes beyond this limit an error message will be displayed and the program will terminate. This is not as great a limitation as it seems because such complex plants could only be accommodated by the screen with a major loss of detail.

The program produces some quite interesting, sometimes unexpected, results making experimentation very worthwhile. Table 6.1 shows some notable values for <F102>unit_angle, unit_length, plant$<F> and <F102>one$<F>. When using your own values be sure to have matching brackets, and also try to plan ahead so that the tree gets larger as well as getting more complex. For example a replacement for one of "11[1]" will cause the replaced branch to get longer whereas "1[1]" will just change it into a two branches, each similar in length to the original.

<$tc>Table 6.1: Variable assignments for a variety of grasses <$table4> <F102>plant$one$ unit_length unit_angle<F>

1 1 1[11[[1]]] 3 PI/12

1 1 [1]1[[1[[1]]1[[1]]]] 4 PI/16

1 1 1[1[[1[1[[1]]]]]] 4 PI/16

<$ha>Trees Only grasses have been drawn thus far because their branches split off in a single direction. Now that the basic algorithm has been established it is fairly easy to alter it for the bi-directional type of branching found in trees and shrubs. All that needs to be done is to add two new characters, the curly brackets ({ and }), to the set recognised by the <F102>draw_branch<F> plant visualisation procedure. The curly brackets are used to enclose descriptions of branches which split off anti-clockwise. The relevant actions for each curly bracket, shown below, are very similar to those associated with the equivalent square brackets. <$hang1>{ Anti-clockwise split: Subtract <F102>unit_angle<F> from current angle count.

<$hang1>} End of branch: Return the graphics cursor to the beginning of the current branch (i.e. where the matching open bracket occurred) and increase the angle count by <F102>unit_angle<F>. Listing 6.3 shows the adapted definition of the <F102>draw_branch<F> procedure required for bi-directional branching. The program code used to check for and act on curly brackets is very similar to that for square brackets. There is no need to edit the definition of <F102>replace_chars<F> because brackets are not replaced under any circumstances. <$ob>PROCEDURE draw_branch(x,y)<$ob> INC pointer !Increment pointer position and…​<$ob> character$=MID$(plant$,pointer,1) !…​note character in that position<$ob> '<$ob> REPEAT<$ob> SELECT character$<$ob> CASE "1" !If the character is a 1, draw a segment<$ob> PLOT x,y<$ob> x=x+SIN(angle_count)*unit_length*2<$ob> y=y-COS(angle_count)*unit_length<$ob> DRAW TO x,y<$ob> '<$ob> CASE "[" !If it’s a [ then branch off clockwise<$ob> ADD angle_count,unit_angle<$ob> GOSUB draw_branch(x,y)<$ob> CASE "{" !If it’s a { then branch off anti-clockwise<$ob> SUB angle_count,unit_angle<$ob> GOSUB draw_branch(x,y)<$ob> ENDSELECT<$ob> INC pointer !Note next character<$ob> character$=MID$(plant$,pointer,1)<$ob> UNTIL character$="]" OR character$="}" OR character$=""<$ob> '<$ob> SELECT character$<$ob> CASE "]" !If it’s a ] (end of branch) then<$ob> SUB angle_count,unit_angle !decrease the angle count..<$ob> CASE "}" !If it’s a } then increase angle…​<$ob> ADD angle_count,unit_angle<$ob> ENDSELECT<$ob>RETURN !…​and return from the procedure <$fc>Listing 6.3: Adapted definition of draw_branch for bi-directional branching If you haven’t altered the definition of <F102>plant$<F> and <F102>one$<F> at the start of the program you will see exactly the same output as before. However, the whole program can now cope with curly brackets and if the variables are set to some of the values from Table 6.2 trees with bi-directional branching will be generated. Note that it is perfectly natural for anti-clockwise branches to have clockwise sub-branches and vice-versa. Some example trees are shown in Figures 6.7 and 6.8. <$fc>Figure 6.7: Tree produced using the fourth set of values from Table 6.2

<$fc>Figure 6.8: Tree produced using the third set of values from Table 6.2 <$tc>Table 6.2: Variable assignments for a variety of trees

<$table4><F102>plant$ one$unit_length unit_angle <F> <$table4>1 11{1{1}[1]}[1[1{1}]] 3 PI/8 <$table4>1 11[1[[1[1[[1]]]]]]{1{{1{1{{1}}}}}} 4 PI/16 <$table4>1 11[1{1}1]{11[1]1} 3 PI/8 <$table4>1 1{{11}}[[11]]1[{{11}}[[11]]1[{{1}}[[1]]1]] 8 PI/10 When is a Plant not a Plant? The plant generation program developed in this chapter can also be used to create a wide variety of non-biological structures. For example, by replacing the variable assignment lines with those shown below it is possible to produce a very complex snowflake structure (see Figure 6.9). <$fc>Figure 6.9: Snowflake generated using the plant program

Naturally the character replacement technique is exactly the same as that used for the plants above, and even the replacement string looks familiar, but in this case the plant initially has six branches, each one at 60 degrees to its neighbour. In this example <F102>one$<F> specifies that each segment should be replaced by a slightly longer segment with a bifurcation at the end. Note that because the strings include curly brackets they can only be implemented in the bi-directional plant drawing program. <$ob>plant$="{{{11}}}{{11}}{11}[11][[11]]11" <$ob>one$="11[1]{1}" <$ob>unit_length=1 <$ob>unit_angle=PI/3 In addition the origin of the plant must be moved from the bottom of the screen to the middle, because the flake grows in all directions. The plant’s origin is changed by altering the line which calls the <F102>draw_branch<F> procedure so that it reads: <$ob> GOSUB draw_branch(320,100)

which causes a different initial drawing position to be sent to the procedure.

It is accepted that every snowflake should be unique in structure. This is hard to achieve on a computer but the structure of the one in the plant program can be drastically altered by simply changing the contents of <F102>one$<F>. As an example, "11[1]{1}1" would give trifurcations at the end of each branch. The C-Curve If you are already familiar with chaos theory you may be able to identify some other, non-branching applications of the character replacement routine used in the plant generation program. Some of the earliest fractals were based on the idea of replacing straight line sections of shapes with more complex line structures, something which the plant program excels at. These were discovered first because good approximations of them can easily be drawn by hand, without the need for excessive number crunching. <$fc>Figure 6.10: The C-curve after 11 iterations

One such fractal is the C-curve, shown in Figure 6.10. The basic algorithm for drawing it is to start with a single line and then repeatedly replace every straight line with two sides of a right angled triangle. The effect of the first few iterations is shown in Figure 6.11.

<$fc>Figure 6.11: The replacement process used to create the C-curve After several iterations the curve becomes very complex in structure, resembling an elaborate letter C, hence the name. The relevant modifications to the program are: <$ob>plant$="[[1]]" <$ob>one$="{1[[1{" <$ob>unit_length=3 <$ob>unit_angle=PI/4 <$ob> GOSUB draw_branch(200,160)

These values have been specially formulated to give the most detailed and largest C-curve possible. In this case the complexity of the fractal is limited by the number of recursive procedure calls that GFA BASIC can cope with, rather than the length of the description string. The recursion problem becomes particularly acute in the C-curve because the replacement string (<F102>one$<F>) contains unmatched open brackets. This means that the <F102>draw_branch<F> procedure has to call itself many more times than usual (once for every bracket), so memory is rapidly used up to store the ever-growing stack of arguments from previous calls. Because the problem originates from excessive recursion the Memory full error box will appear while drawing is taking place, rather than during string processing. This is how this type of error can be distinguished from one caused by <F102>plant$<F> becoming too long. Note that for this particular fractal (and the Koch curve below) a non-recursive drawing procedure would be preferable as it would facilitate the creation of more detailed images.

The Koch Curve

The Koch curve is a kind of aesthetically pleasing version of the C-curve, which originates from an equilateral triangle (see Figure 6.12).

<$fc>Figure 6.12: The replacement process used to create the Koch curve The middle third of each side of the triangle is replaced by two sides of a smaller equilateral triangle, leaving a 12-sided shape. The same replacement process is then repeated on every side of the new shape, and on every side of the subsequent shape. After an infinite number of iterations the shape becomes the Koch curve (Figure 6.13), often referred to as the Koch flake due to its snowflake-like structure. <$fc>Figure 6.13: Approximation of the Koch curve

The replacement lines needed to produce the curve using the plant program are:

<$ob>plant$="1[[1[[1]]]]" <$ob>one$="1{1[[1{1" <$ob>unit_length=200 <$ob>unit_angle=PI/3

<$ob>GOSUB draw_branch(240,200) To make the Koch curve appear to change, rather than grow, the unit_length variable must initially be large and then be reduced to a third of its value after each iteration. This is done by altering the first REPEAT…​UNTIL loop in the plant generation program so that it reads as follows: <$ob>REPEAT <$ob> CLS <$ob> GOSUB draw_plant !Draw plant so far <$ob> GOSUB replace_chars !Make plant grow <$ob> unit_length=unit_length/3 !Reduce size by a third <$ob>UNTIL mousek=2 Of course, the curve created here is only an approximation of the real Koch curve due to the limitations of the computer, which make it impossible to perform an infinite number of iterations. Like many other fractals in this book the Koch flake has an infinite length. Although this is impossible to prove practically it can be done using elementary maths. Assuming that the length of each side of the initial triangle is one metre, its total perimeter is three metres, as shown in Figure 6.14. After one iteration the length of each side becomes 11/3 (4/3) metres long due to the introduction of the triangular deviation, giving a perimeter of 3*(4/3) = 4 metres (an increase of one third). <$fc>Figure 6.14: The size of the equilateral triangle from which the Koch curve is produced

On the next iteration each side is only a third of a metre long, and is only extended by a ninth of a metre to 4/9 metre. But because there are 12 sides the overall length increases by a third to 12*(4/9) = 5 metres. On each iteration the length of the perimeter increases by a third (represented mathematically as multiplying by 4/3), so on successive iterations the perimeter gets larger at an increasing rate. For example, after two iterations the perimeter would be:

<$table2>3 * 4/3 * 4/3 = 5.33m After four iterations the perimeter would be: <$table2>3 * 4/3 * 4/3 * 4/3 * 4/3 = 9.48m

As a general rule the perimeter length can be expressed as:

length = 3 * (4/3)^n

where n is the number of iterations. Using this expression it is possible to see that the length quickly becomes very large. After only 45 iterations (n=45) it is already over one million metres. Substituting infinity for the iteration number shows that at this stage the perimeter is itself infinite in length. However, the area of the shape stays small implying that it is possible to store an infinitely long line in a finite space (less than one square metre in this example). To reflect the way that the Koch curve is more like a two dimensional shape than a one dimensional line (in its ability to fill space) it is given a fractal dimension of 1.2818.

The infinite complexity of the Koch curve is very difficult to comprehend. Magnifying one section reveals that it is made up of similar, equally detailed sections, and further magnification reveals even smaller and more detailed sections. Unlike normal geometrical shapes it would be impossible to magnify the curve so as to expose the straight lines from which it is totally composed, a concept familiar from earlier chapters.

In classical mathematics the slope (or gradient) of a curve at any point can normally be found by mathematically 'magnifying' it until it approximates to a straight line, for which the gradient can easily be determined. However, because the Koch curve changes direction at every point and never approximates to a straight line it is impossible to find a single value for the gradient at any given point. This is just one of many reasons why classical mathematicians have always viewed the Koch curve (and similar fractals) with a noticeable amount of disdain.

Further Experimentation

Growing artificial plants is a very interesting pastime, so I have attempted to write the accompanying program in such a way as to make experimentation in this area as easy as possible. The way the fractal grows can be altered simply by changing the few variables initialised at the start of the program. Many other natural fractal patterns such as lightning and frost on a window can be created using the methods and programs discussed in this chapter. A variety of more specific enhancements are outlined below.

Branches with Variable Thickness

One of the most noticeable simplifications of the current model is that the thickness of each branch is in no way related to the load it has to bear. Although grasses have an acceptable appearance, trees look as though their trunks are as thin as twigs. The GFA BASIC DEFLINE command allows us to set the line width, which is then used in all future calls to DRAW, CIRCLE and the other line drawing commands. By reducing the line width in this way whenever a new branch begins (whenever an open bracket is reached) the branches realistically appear to get thinner towards their ends.

More Complex Complexity

The rule set used throughout this chapter only allows for one type of segment character, "1". The addition of a second type of segment, which has a different replacement string, means that still more complex and diverse branching effects can be achieved using the existing algorithm. The character we shall use for this other type of segment is the zero (0). The changes required to implement the new segment type are fairly minor. The first is near the beginning of the <F102>draw_branch<F> procedure. This line:

<$ob>IF character$="1"

should be changed to:

<$ob>IF character$="1" OR character$="0" so that the program draws a segment if a 1, or a zero, is found. The <F102>replace_chars<F> procedure should then be edited to read: <$ob>PROCEDURE replace_chars <$ob> pointer=0 !Set initial pointer position <$ob> REPEAT <$ob> INC pointer !Increment pointer position <$ob> character$=MID$(plant$,pointer,1) !Extract character <$ob> IF character$="1" !If it’s a 1 then replace it with one$ <$ob> newplant$=newplant$+one$ <$ob> ENDIF <$ob> IF character$="0" !If it’s a 0 then replace it with zero$ <$ob> newplant$=newplant$+zero$ <$ob> ENDIF <$ob> IF character$="{" OR character$="}" OR character$="[" OR character$="]" <$ob> newplant$=newplant$+character$ !Otherwise leave it unchanged <$ob> ENDIF <$ob> UNTIL pointer=LEN(plant$) <$ob> plant$=newplant$ !Replace old plant$with new one <$ob> newplant$="" !Clear newplant$ just in case <$ob>RETURN Changes are needed here because it can no longer be assumed that if the character taken from the string is not a 1 it will not need replacing - each character must be explicitly checked for. This modified program will still draw all the fractals discussed earlier in the chapter, but just as all occurrences of "1" are replaced by the contents of <F102>one$<F>, all zeros are replaced by the contents of <F102>zero$<F>. The example variable assignments below will allow you to pilot the new program. <$ob>plant$="10" <$ob>one$="10[01]{10}" <$ob>zero$="1[1{01}0]{11}" <$ob>unit_length=5 <$ob>unit_angle=PI/6 Note that any number of different segment types may be added to the program in the same way. If there are more than two possible types of segment it is advisable to use an array to store the relevant replacement strings, failure to do this will make the <F102>replace_chars<F> procedure very much longer and will result in slower execution. Using C The tree drawing routine is a very good example of a program that should have been written in C. The messy global variables and constants would not be a problem in C as values can easily be returned at any point in a routine and true constants are well supported. Although GFA BASIC allows functions to return values, its functions are limited to one line, which is too short in this case. C also makes available more efficient, user-definable data structures in which plants could be stored. 7. Imitating Nature - Fractal Landscapes Probably the slowest natural process discussed in this book is landscape dynamics. Although normally appearing stationary for many thousands of years, the Earth’s landscape is actually undergoing constant change due to erosion, volcanic activity and more gradual land movements. The lack of speed makes the process slightly easier to predict than something like the weather, but no easier to simulate because it is still inherently chaotic. This chapter describes some of the simplest methods for creating natural looking landscapes and proceeds on to discuss the application of landscape drawing techniques to fractals such as the Feigenbaum diagram and the Mandelbrot set. Landscape generation programs produce aesthetically pleasing results and provide an interesting introduction to the processes which create terrestrial landscapes, making experimentation very worthwhile. Note that the programs provided in this chapter are all written in GFA BASIC. This is because when drawing solid looking landscapes it is necessary to draw filled polygons with 50 corners. Unlike GFA’s POLYFILL command the equivalent Amiga BASIC command (AREAFILL) only allows filled polygons with 20 corners, so is not suitable. Note that it is possible to alter the landscape drawing routine used in this chapter so that three filled polygons, or a single unfilled polygon, may be used instead. The algorithms discussed for generating landscape surfaces are, in any case, very simple and can be converted to run in any language. Something common to all landscapes is that they are three dimensional, they stretch out horizontally in two directions and also vary in height. To draw such objects it is obviously necessary to find a technique which will allow plotting in three dimensions. The most popular way of doing this is called isometric drawing. Isometric Drawing On a two dimensional plane, such as that used to display most graphs, there are two directions of movement, represented by two axes at right angles to one another, as shown in Figure 7.1. <$fc>Figure 7.1: Typical two dimensional plane

The two directions of movement are vertical (represented by the y axis) and horizontal (represented by the x axis). In three dimensional space, however, there is an additional direction of movement, and hence three axes (x, y and z). As the screen of the Amiga is only two dimensional it is of course impossible to display three dimensional objects, but take a look at Figure 7.2.

<$fc>Figure 7.2: Isometrically drawn three dimensional space and cube This shows the three axes used to represent three dimensional space, with a cube placed where they meet. Obviously this diagram is not really three dimensional, but the illusion of three dimensions is produced sufficiently to convey the fact that the object shown is a 3D cube. This effect is achieved by plotting the x and z axis as normal but drawing the y axis at 30 degrees to the x axis, and treating all corresponding movements in the y direction in the same way. The technique is called isometric drawing and is frequently used by draughtsmen when producing illustrations of their designs, and can also be found in some computer games and demos. Drawing objects such as spheres and cubes isometrically can be a fairly complicated process as a method must be devised for storing the structure of objects, but plotting landscapes is somewhat easier due to landscapes simply being uneven planes. <$fc>Figure 7.3: Mathematical splitting of y into horizontal and vertical components

To perform isometric drawing in GFA BASIC an algorithm is required to convert a three dimensional (x,y,z) co-ordinate to the two dimensional form accepted by the BASIC PLOT and DRAW functions. This algorithm involves splitting the isometrically drawn y position into its constituent horizontal and vertical components, as shown in Figure 7.3, and adding these components to the x and z positions for the point. By merging y with x and z in this way we have effectively eliminated y and are now left with a two dimensional co-ordinate of the form (x,z) which can be used directly with PLOT, POLYLINE or any other graphics function available from BASIC.

If the components of y are to be added to x and z to produce a composite co-ordinate it is necessary to determine the size of the components. This is done by calculating the ratio between them using elementary trigonometry, and since it is a ratio, rather than an absolute value, it is the same for any value of y. The mathematical tangent function can be used to find this universal ratio because the tangent of an angle is defined as the ratio of the length of the angle’s opposite and adjacent sides in a right angled triangle. From Figure 7.3 it is clear that the ratio of the z component to the x component is the tangent of 30 degrees. The value of this tangent can be found using the tan button of a pocket calculator, which should return a result of approximately 0.6. The mathematical interpretation of this process is shown below (angles, right angled triangles and the tangent function are discussed in detail in Appendix B).

<$table5>tan 30 = z component = 6 (or 0.6 or 6:10) <$table5> x component 10

Now that the ratio of the z component to the x component is known it is possible to plot points on the screen isometrically. The method for doing this is best conveyed by the example plane drawing program in Listing 7.1. Here a pair of FOR…​NEXT loops are used to completely traverse a 48x48 point plane, on which each point is plotted isometricaly. The (x,y) position of the point passed to the PLOT command is calculated as follows:

<$table4>Position passed Offset x and z parts y components <$table4>to PLOT <$table4>horizontal = x*5 + 1.7*y <$table4>vertical = 148 - y

There are several important points to note in these equations. Firstly the z position of every point is zero, so this is not included in the calculation of the y position. Secondly the horizontal component of y is equal to 1.7*y and the vertical component is equal to y (or 1*y), so the ratio is correct for 30 degrees (1/1.7 = 0.6). Also note that the vertical position passed to PLOT is inverted by subtracting it from 148 - necessary because of the Amiga’s upside-down vertical axis. The x position is multiplied by 5 to ensure that the maximum screen area is utilised, although in this example the fact that z is equal to zero means that there are large spaces at the top and bottom of the screen.

<$ob>FOR y=47 TO 0 STEP -1<$ob> FOR x=0 TO 47<$ob> PLOT x*10+y*3.4,148-y<$ob> NEXT x<$ob>NEXT y <$fc>Listing 7.1: The plane drawing program

<$fc>Figure 7.4a: Isometric 3D plane produced by Listing 7.1 The plane produced by Listing 7.1 (in high resolution mode) is shown in Figure 7.4a. A landscape could easily be represented here by giving each point a height determined by some particular method, but since each point is plotted as a single pixel, with no connections to adjacent pixels the points would quickly become confused making the image meaningless. <$fc>Figure 7.4b: Testcard landscape drawn using single pixels

Figure 7.4b shows the single pixel plotting method being used to draw a landscape 'testcard'. A more suitable plotting method involves joining together all points that have the same y position in order to create a series of 48 lines depicting the plane.

<$fc>Figure 7.4c: Testcard consisting of lines As Figure 7.4c shows this is a great improvement, but the plane appears transparent in places because lines which are meant to be behind or below other parts of the plane can still be seen. The process of removing unwanted lines to make an object appear solid is called 'hidden line removal', the result of which is shown in Figure 7.4d. <$fc>Figure 7.4d: Testcard drawn using lines and hidden line removal

The best way to remove hidden lines in our case is by using the GFA BASIC POLYFILL command to draw each line and fill the gap between it and the last line drawn, in the background colour. Figure 7.5 shows the way in which a polygon is constructed from the two most recent lines so that the gap between them is totally enclosed and can be automatically filled in.

<$fc>Figure 7.5 Application of filled polygons to remove hidden lines The programming needed to implement this hidden line removal system is beyond the scope of this book, and unrelated to chaos theory, so an 'off-the-shelf' procedure has been included in Listing 7.2 which draws a plane with hidden lines removed, from the data describing the height of each point stored in the 48x48 element <F102>z%<F> array. The size and position of this plane, and the relationship between <F102>z%<F> co-ordinates and their screen positions are shown in Figures 7.6 and 7.7. This procedure definition will need to be included in all programs given later in this chapter which contain calls to landscape. <$fc>Figure 7.6: Size of the plane produced by the landscape procedure (in low-resolution pixel units)

<$fc>Figure 7.7: Relationship between elements of the <F102>z%<F> array and the plane drawn by landscape <$ob>PROCEDURE landscape<$ob> CLS<$ob> DEFFILL 1,0<$ob> BOUNDARY 1<$ob> COLOR 1<$ob> '<$ob> FOR y=47 TO 0 STEP -1 <$ob> FOR x=0 TO 47<$ob> poly_x%(x)=x*10<$ob> poly_y%(x)=-z%(x,y)<$ob> NEXT x<$ob> IF y<47<$ob> FOR x=0 TO 47<$ob> poly_x%(95-x)=x*10+3.4<$ob> poly_y%(95-x)=-z%(x,y+1)-1<$ob> NEXT x<$ob> POLYFILL 96,poly_x%(),poly_y%() OFFSET y*3.4,(148-y)<$ob> REM Enhancements may be added here (see text)<$ob> ELSE<$ob> POLYLINE 48,poly_x%(),poly_y%() OFFSET y*3.4,(148-y)<$ob> ENDIF<$ob> NEXT y<$ob>RETURN <$fc>Listing 7.2: The <F102>landscape<F> procedure, which uses hidden line removal to render a landscape from the contents of the <F102>z%<F> array Using the Landscape Procedure To use this procedure it is necessary to dimension two 96 element temporary storage arrays - <F102>poly_x%<F> and<F102> poly_y%<F> and also the 48 by 48 element <F102>z%<F> array. The height of each of the 2304 (48*48) points on the plane should then be placed in the correct elements of the <F102>z%<F> array before <F102>landscape<F> is called to plot the plane. As an example the program shown in Listing 7.3 would be used to set the height of the point in the middle of the plane to be 20 (low resolution) pixels high and draw the resulting landscape. Running this short program will determine whether the <F102>landscape<F> procedure has been entered correctly, and also demonstrates the minimum program construction needed to use <F102>landscape<F>. <$ob>REM Set up arrays<$ob>DIM poly_x%(95)<$ob>DIM poly_y%(95)<$ob>DIM z%(47,47)<$ob>'<$ob>REM Set height of point<$ob>z%(23,23)=20<$ob>'<$ob>REM Draw landscape<$ob>GOSUB landscape <$fc>Listing 7.3: Setting the height of a single point

Naturally it would be impractical to set the height of all 2,304 points in this way, so a nested pair of FOR…​NEXT loops are generally used to fill the array instead. An example of a program using such a method, to draw the testcard, is shown in Listing 7.4. Note that the testcard is actually a three dimensional cosine curve on which the height of each point is determined according to its distance from the centre of the plane, and is in no way related to fractal techniques. Fractal landscapes are discussed in the next section.

<$ob>REM Set up arrays<$ob>DIM z%(47,47)<$ob>DIM poly_x%(95)<$ob>DIM poly_y%(95)<$ob>'<$ob>REM Set point heights in z% array<$ob>PRINT "Calculating…​"<$ob>'<$ob>FOR x=0 TO 47<$ob> FOR y=0 TO 47<$ob> z%(x,y)=COS(SQR(24-x)2)+((24-y)2/2)*20<$ob> NEXT y<$ob>NEXT x<$ob>'<$ob>REM Draw landscape<$ob>GOSUB landscape

<$fc>Listing 7.4: Using a FOR…​NEXT loop to set the height of all points on the plane <$hb>Improving the Landscape Procedure

The planes produced by <F102>landscape<F> can be made to look more solid by joining points along the direction of the y axis to make the plane look like it is composed of tiny quadrilaterals, as shown in Figure 7.8. This effect is most suitably used in the Amiga’s high resolution interlaced mode, where the small pixels mean that the landscape looks less cluttered.

<$fc>Figure 7.8: Testcard landscape composed of quadrilateral segments The <F102>landscape<F> procedure can be adapted for this purpose by adding the following lines to the procedure definition given in Listing 7.2, at the point marked for enhancements. FOR x=0 TO 47 DRAW poly_x%(x)+y*3.4,poly_y%(x)+148-y TO poly_x%(95-x)+y*3.4,poly_y%(95-x)+148-yNEXT x Other improvements which could be included in this section of the routine, such as the addition of shading, are discussed at the end of this chapter. Generating Pseudo-Natural Landscapes Now that we have a method of plotting a landscape we can begin to use chaos techniques to determine the height of each point on the plane, in order to give the landscape a more natural appearance. There are numerous methods available for producing landscapes with features similar to those found in the real world, and more realistic ones are being produced all the time. Upon looking at a natural landscape it is evident that the process used to set the height of each point on the plane would have to be random, and a first attempt at a landscape program may look like Listing 7.5. <$ob>DIM z%(47,47)<$ob>DIM poly_x%(95)<$ob>DIM poly_y%(95)<$ob>'<$ob>PRINT "Calculating…​"<$ob>'<$ob>FOR x=0 TO 47<$ob> FOR y=0 TO 47<$ob> z%(x,y)=RANDOM(21)-10 !Set height of point<$ob> NEXT y<$ob>NEXT x<$ob>'<$ob>GOSUB landscape

<$fc>Listing 7.5: A totally random landscape program This sets the height of each point randomly within the range of -10 to 10, and certainly does not produce a natural landscape, as Figure 7.9 demonstrates. <$fc>Figure 7.9: Totally random landscape

The problem is that the degree of randomness is too great, causing incredibly steep slopes and sudden points of height to be produced. On a real landscape such steep slopes would be made gentler by a gradual process of weathering. It is true that most landscape generators rely on some kind of random process, but to generate a realistic landscape it is necessary to control the randomness in some way so that it does not get too out of hand, as in the Sierpiński triangle. Two generators based on controlled randomness are presented below.

Inheritance

The main problem with Listing 7.5 is that there is no continuity between adjacent points, meaning that huge unnatural pillars of land can stand out of the ground. Using a process of inheritance the height of each point is always related in some way to the height of at least one neighbouring point. A simple example of this can be demonstrated in two dimensions using the program in Listing 7.6. <$ob>y=0 !Set initial y position<$ob>PLOT 20,100-y !Plot the position<$ob>FOR x=21 TO 620 STEP 2<$ob> y=y+RANDOM(5)-2 !Calculate new y value<$ob> DRAW TO x,100-y !Draw line to new point<$ob>NEXT x

<$fc>Listing 7.6: Relating the height of a point to its neighbour using random inheritance In this example the first point (at the left of the screen) is set to be at a height of zero pixels. A FOR…​NEXT loop then begins to draw the rest of the points, whose heights are calculated by adding a random number, between -2 and 2, to the height of the previous point. The result, shown in Figure 7.10, is reminiscent of the cross section of a rather mountainous landscape. <$fc>Figure 7.10: Random inheritance in two dimensions

The roughness of the terrain can easily be altered by changing the range of the random number added to the previous point, for instance the following line will give gentler slopes:

<$ob> y=y+RANDOM(3)-1 By flattening points below a certain height, referred to as the sea level, it is possible to create the effect of lakes, rivers and seas between pieces of land, as shown in Figure 7.11. <$fc>Figure 7.11: Cross-section of a landscape with seas and lakes

Replacing the line which draws between points with the IF…​ENDIF construct below will draw all points with y positions below zero as if their y positions were zero. In this example the sea level can be said to be zero. Note that the value of y is not actually altered, as this would upset the inheritance process.

<$ob> IF y>0 !If y is above sea level <$ob> DRAW TO x,100-y !Draw line to new point <$ob> ELSE <$ob> DRAW TO x,100 !Draw line at sea level <$ob> ENDIF This method of flattening sub-zero sections of the plane to create areas of water can be used in all of the natural landscape generators detailed in this chapter. Enhancing the inheritance process for three dimensions is relatively easy, since the plane plotting program has already been provided. However, in three dimensions, the inheritance pattern is slightly different because each point has up to two adjacent points from which to inherit values, as shown in Figure 7.12. <$fc>Figure 7.12: Inheritance on a plane in three dimensions

It can be seen from this diagram that the first point calculated, at (0,0), is completely random because there are no adjacent points from which heights can be inherited. Other points along the left and lower sides of the plane have only one adjacent point whose height has been calculated, so they inherit their height from a single point only, as in the two dimensional example in Listing 7.6. All other points on the plane have two processed neighbours, so the heights of the two adjacent points are combined, by adding them together and dividing the result by 2. The resulting value is used to determine the height of the new point. Listing 7.7 uses this method to fill the <F102>z%<F> array with values and will plot the result.

<$ob>DIM z%(47,47)<$ob>DIM poly_x%(95)<$ob>DIM poly_y%(95)<$ob>'<$ob>PRINT "Calculating…​"<$ob>'<$ob>FOR x=0 TO 47<$ob> FOR y=0 TO 47<$ob> IF x=0 AND y=0 !(0,0) - no inheritance<$ob> z%(x,y)=RANDOM(5)<$ob> ENDIF<$ob> IF x=0 AND y>0 !Left edge of plane - one neighbour<$ob> z%(x,y)=z%(x,y-1)+RANDOM(5)-2<$ob> ENDIF<$ob> IF x>0 AND y=0 !Nearside edge - one neighbour<$ob> z%(x,y)=z%(x-1,y)+RANDOM(5)-2<$ob> ENDIF<$ob> IF x>0 AND y>0 !Rest of plane - two neighbours<$ob> z%(x,y)=(z%(x,y-1)+z%(x-1,y))/2+RANDOM(5)-2<$ob> ENDIF<$ob> NEXT y<$ob>NEXT x<$ob>'<$ob>GOSUB landscape

<$fc>Listing 7.7: Using up to two neighbours to determine the height of the new point As the program shows, the very artificial inheritance method of determining the height of the points produces surprisingly realistic landscapes, an example of which is shown in Figure 7.13. <$fc>Figure 7.13: Pseudo-natural landscape created by Listing 7.7

A more refined type of inheritance can be created by basing the change in height between one point and its neighbour on the height difference between the neighbour and its other neighbour. This creates slightly more natural landscapes with gentler slopes.

Faulting

In the real world landscapes are formed as the result of a variety of processes, including large movements of land, one type of which forms the basis of the faulting method of landscape generation. Geologists have found that the Earth’s crust can be thought of as a relatively thin layer of solid land floating on a sea of molten rock, just like ice floating on water. Because of pressure points and other phenomena occurring inside the Earth it is never possible for the crust to be in one piece. Instead it is made up of enormous plates of land, which overlap and slide underneath each other where they meet. When the edges of two such plates push together, but do not slide over each other, a large force is exerted between them.

Initially this force goes unnoticed, but it eventually builds up to such an extent that one of the plates must suddenly rise or fall in order to relieve the pressure. This sudden movement of land is called a fault, shown pictorially in Figure 7.14. After hundreds of faults occurring over millions of years spectacular mountainous landscapes can be created.

<$fc>Figure 7.14: A geographical fault It is relatively easy to perform faulting across the plane produced by the landscape procedure, and the results can be very realistic. After lengthy experimentation I have found that approximately 100 faults produce very realistic landscapes. The program in Listing 7.8 will randomly produce one 100 straight faults, after every 10 of which the landscape plane is displayed. <$ob>DIM z%(47,47)<$ob>DIM poly_x%(95)<$ob>DIM poly_y%(95)<$ob>DIM side(1),x(1),y(1)<$ob>'<$ob>FOR iteration=1 TO 100<$ob> '<$ob> REM Select sides to draw line between<$ob> side(0)=RANDOM(4)<$ob> REPEAT<$ob> side(1)=RANDOM(4)<$ob> UNTIL side(0)<>side(1) !Make sure that they’re different<$ob> '<$ob> REM Pick a random position on each side to draw between<$ob> FOR f=0 TO 1<$ob> IF ODD(side(f))<$ob> x(f)=RANDOM(48)<$ob> y(f)=((side(f)-1)/2)*48<$ob> ELSE<$ob> x(f)=(side(f)/2)*48<$ob> y(f)=RANDOM(48)<$ob> ENDIF<$ob> NEXT f<$ob> '<$ob> REPEAT<$ob> change%=RANDOM(5)-2 !Select random change in height<$ob> UNTIL change%<>0 !for -2 to 2 (but not 0)<$ob> '<$ob> REM Plot fault on screen<$ob> CLS !Clear the screen<$ob> PRINT AT(1,8);"Iteration:";iteration<$ob> BOX -1,-1,48,48 !Draw the plane boundary<$ob> DRAW x(0),y(0) TO x(1),y(1) !Draw the fault line<$ob> DEFFILL 1,2,8 !Set fill pattern to solid black<$ob> FILL RANDOM(2)*47,RANDOM(2)*47 !Fill one side of line<$ob> '<$ob> REM Alter z% array from screen data<$ob> FOR x=0 TO 47<$ob> FOR y=0 TO 47<$ob> IF POINT(x,y)>0 !If colour of point is black…​<$ob> z%(x,y)=z%(x,y)+change% !…​alter the height by change%<$ob> ENDIF<$ob> NEXT y<$ob> NEXT x<$ob> '<$ob> REM Plot landscape every ten iterations<$ob> IF INT(iteration/10)=iteration/10<$ob> GOSUB landscape<$ob> ENDIF<$ob> '<$ob>NEXT iteration

<$fc>Listing 7.8: Program to create a landscape using random faulting The method used to generate faults here can best be understood by following the sequence on the screen, as this is where all the processing occurs. First a box is drawn on the screen around the area x=0, y=0 to x=47, y=47 - the enclosed area represents the plane. A straight line is then randomly drawn from one side of the box to another. The area on one, randomly selected, side of the line is then filled in using the FILL command. A nested pair of FOR…​NEXT loops containing a call to the POINT function are then used to determine the colour of each pixel in the area representing the plane. All those which are filled in then have their corresponding point on the 3D landscape plane raised (or lowered) by a uniformly random amount, the rest of the points being left as they were. Figures 7.15a to 7.15d show the faulting program in action. <$fc>Figure 7.15a - d: The process after 10, 40, 75 and 100 iterations

Determining which points lie on which side of the line using the FILL and POINT commands is something of a cheat to avoid mathematics, but it is very effective and also means that a range of differently shaped faults can be accommodated. For example, a lunar surface can be generated by randomly placing randomly sized circles inside the box to create crater effects.

Using Landscape Techniques to Plot Fractals in 3D

Earlier in the book it was suggested that some fractals, such as the Feigenbaum diagram and Mandelbrot set may benefit from being displayed in three dimensional landscape form. As we now have the <F102>landscape<F> procedure, which will produce a landscape from the contents of the two dimensional <F102>z%<F> array, such images are easy to create. The same processes and calculations as those used in the original fractals must be followed of course, but instead of plotting to the screen values are placed in the <F102>z%<F> array and plotted later by calling <F102>landscape<F>.

The Feigenbaum Diagram

In Chapter 2, where the Feigenbaum process was discussed, it was mentioned that colour could be used to represent the number of times a point was plotted, but colours did not present this very clearly. A more effective way of showing this information is to draw the Feigenbaum diagram in three dimensions, with the height of each point representing the number of times a point was plotted.

This is very easy to perform using the <F102>landscape<F> procedure, as the <F102>z%<F> element for each point can be incremented each time it is visited, eventually giving a height proportional to the number of times the point would have been plotted. The program to draw the 3D diagram is given as Listing 7.9 and the result is shown in Figure 7.16.

<$fc>Figure 7.16: Feigenbaum landscape It is easy to see here that the period one section is the highest part of the diagram (raised 50 pixels above the plane) but at bifurcation the height of the points is halved (to 25). In the chaotic region at the far end of the plane the points vary in height, but are generally very low (one or two pixels high). <$ob>DIM z%(47,47)<$ob>DIM poly_x%(95)<$ob>DIM poly_y%(95)<$ob>'<$ob>PRINT "Calculating…​"<$ob>'<$ob>FOR c=1.8 TO 3 STEP 0.025<$ob> y=INT((c-1.8)40)<$ob> p=0.3<$ob> FOR iteration=0 TO 100<$ob> IF iteration>25<$ob> x=p*34<$ob> INC z%(x,y) !Raise point by one pixel<$ob> ENDIF<$ob> p=p+c*p(1-p)<$ob> NEXT iteration<$ob>NEXT c<$ob>'<$ob>GOSUB landscape

<$fc>Listing 7.9: Program to produce a Feigenbaum landscape similar toFigure 7.16 Much of this listing is familiar, the only difference being that instead of scaling the co-ordinates up to the size of the screen they are scaled to fit on the 48x48 point plane. The corresponding fact that only 48 values of c need to be tested makes the 3D method considerably faster than plotting to the screen. Because of the variety of height on the plane it is quite difficult to recognise the Feigenbaum diagram when it is drawn in this way, as the foreground obstructs most of the diagram. To see the structure of the diagram more easily points on the plane can be lowered in proportion to the number of times they were visited, rather than raised. This can easily be achieved by changing the <F102>INC z%(x,y)<F> line to read: <$ob> DEC z%(x,y) !Depress point by one pixel

The Mandelbrot Set

As explained in Chapter 4, the coloured contours around the Mandelbrot set are used to indicate the number of iterations required to free the corresponding point from the circle. The drawback of coloured contours is that the Amiga’s standard high-resolution palette of 16 colours makes it impossible for each iteration number to have a unique colour. By plotting the set in three dimensions this problem is overcome because each iteration number can have a unique height instead (provided that the iteration ceiling does not exceed the vertical resolution of the screen). Listing 7.10 is a 3D conversion of the multi-purpose Mandelbrot plotter given in Listing 4.4, it will draw any section of the set (specified using the variables at the beginning of the program) in three dimensions. The result of running the program with the default values of <F102>xmin=-2, xmax=2, ymin=-2, ymax=2 <F>and <F102>max_iteration=33<F> is shown in Figure 7.17.

<$fc>Figure 7.17: Mandelbrot landscape <$ob>DIM z%(47,47)<$ob>DIM poly_x%(95)<$ob>DIM poly_y%(95)<$ob>'<$ob>COLOR 1<$ob>BOX 0,0,47,47<$ob>'<$ob>REM User editable constants<$ob>amin=-2 !Lowest value of a<$ob>amax=2 !Highest value of a<$ob>bmin=-2 !Lowest value of b<$ob>bmax=2 !Highest value of b<$ob>max_iteration=33 !Iteration ceiling<$ob>'<$ob>REM Calculate ranges, offsets and multipliers<$ob>a_range=amax-amin<$ob>b_range=bmax-bmin<$ob>a_mult=47/a_range<$ob>b_mult=47/b_range<$ob>a_offset=(0-amin)*a_mult<$ob>b_offset=(0-bmin)*b_mult<$ob>'<$ob>FOR a=amin TO amax STEP (a_range/47)<$ob> FOR b=bmin TO bmax STEP (b_range/47)<$ob> p=0<$ob> q=0<$ob> iteration=0<$ob> REPEAT<$ob> pnew=p*p-q*q+a<$ob> qnew=2*p*q+b<$ob> p=pnew<$ob> q=qnew<$ob> INC iteration<$ob> UNTIL p*p+q*q>=4 OR iteration=max_iteration<$ob> xp=ROUND(a*a_mult+a_offset) !Calculate pixel position<$ob> yp=ROUND(b*b_mult+b_offset) !<$ob> z%(xp,yp)=iteration !Place height in z% array<$ob> IF ODD(iteration)<$ob> PLOT xp,yp !Plot to screen<$ob> ENDIF<$ob> NEXT b<$ob>NEXT a<$ob>'<$ob>REM Plot landscape<$ob>GOSUB landscape

<$>Listing 7.10: 3D GFA BASIC version of the Mandelbrot plotter given in Listing 4.4 As with the Feigenbaum diagram, the program is much faster than the equivalent two dimensional one because it is scaled to plot to the 48x48 plane rather than the whole screen. In order to trace the progress of the calculations a tiny replica of the plane is produced while the array is being filled, before the call to <F102>landscape<F> is made. Unlike the faulting method of landscape generation this replica is not used during the creation of the 3D Mandelbrot set. <$fc>Figure 7.18: Depressed Mandelbrot landscape

As Figure 7.18 shows, the set can be inverted by altering the line which sets the height of each point in the <F102>z%<F> array <F102>(Z%(Xp,Yp)=Iteration)<F> to read:

<$ob> Z%(Xp,Yp)=-Iteration The method and scaling detailed here can also be used to plot three dimensional Julia sets. Further Ideas for Fractal Landscapes As indicated earlier in this chapter there are many techniques available for generating natural looking landscapes, but only a small selection have been discussed here. Details of more complicated, and more realistic, landscape generators can be found in the books listed in Chapter 9. Even if you find the naturalness of the landscapes given here acceptable it is still possible to improve their appearance in a variety of ways, described below. The techniques used to plot the Mandelbrot set and Feigenbaum diagram may also be enhanced to produce better quality output, and both Martin fractals (discussed in the next chapter) and Julia sets may be drawn in landscape form using the methods described above. Shading Colour can be put to good use when drawing both pseudo-natural and fractal landscapes. It is particularly suited to landscapes which are composed of tiny quadrilaterals, because each can be filled in with a colour determined by some shading method. The easiest of all colouring methods involves filling areas of water in blue and areas of land with green. However, the simplicity of this solution is such that the black divisions between the segments of landscape still need to be drawn to separate foreground hills from those further back. Ideally the segments of land should be coloured differently depending on their height, or maybe depending on the angle they make with the vertical and the direction in which they face. This is one of the few occasions when using a colour monitor can actually increase the apparent resolution because each colour segment of land can be as small as one pixel, unlike the monochrome segments which need to be larger to convey their position and inclination properly. Of course the best (but slowest) shading method would be ray tracing, the intricacies of which are best left to a specialised book on 3D graphics. Such a book would also explain rotation in three dimensions which would allow landscapes to be viewed from different angles, making more of their structure visible. Although experimentation with shading is to be recommended, three dimensional techniques of this complexity often require significant mathematical and programming skill. Trees and Bushes A interesting addition to a pseudo-natural landscape image is foliage. The plant generation techniques discussed in the previous chapter can be applied to create a variety of trees and bushes on landscapes. The starting and replacement strings can be made to vary from plant to plant using the same type of inheritance as that discussed above. Each plant could have random description strings generated based on those of neighbouring plants, causing realistic grouping of similar plants. Facts about possible growing positions such as height, gradient and proximity to water could also be used when determining what type of plants the computer should place where. Adding plants can be very enlightening, but a difficult task if approached in the wrong way. However, the following guide-lines should make it relatively easy to incorporate plants into a landscape generator: <$bullet>It is more convenient to place plants so that their trunks emanate only from the points where four segments meet, because the elevation of these points above the plane can quickly be found by interrogating the <F102>z%<F> array.

<$bullet>The plant drawing routine does not need to be altered because the plant must still be drawn with a vertical trunk, and a two dimensional rendering is still perfectly acceptable. <$bullet>The addition of plants to monochrome landscapes is not advisable because the plants tend to get lost among the lines that make up the landscape.

8. Imitating Nature - Cell Culture

Unlike the two previous chapters on imitating nature this one does not set out to copy a natural object using a set of rules, but rather shows how a purely mathematically based process can generate images which resemble natural structures. The title is actually very much of a generalisation as the images shown here could be said to resemble a range of natural objects, the only similarity being that they look like the type of thing you would expect to see through the microscope of a biology laboratory.

The Martin Process

Of course the common factor linking these images is the process used to draw them, developed by Barry Martin of Birmingham’s Aston University. The process is very similar to those covered at the beginning of the book in that a pair of variables are repeatedly transformed by two non-linear equations and the results plotted on a two dimensional plane, as shown in Figure 8.1.

<$fc>Figure 8.1: Relationship between the two equations and a Martin fractal This process is particularly similar to the one used to draw the Lorenz attractor except that only two variables are used, rather than three, and the points are simply plotted rather than being joined together. Also, the Martin equations (shown below) are considerably simpler, relying only on a few elementary functions. xnew=y-SGN(x)*SQR(ABS(b*x-c)) ynew=a-x Here x and y are the two variables and a, b and c are constants similar to those found in the Lorenz equations. The important parts of these equations are the SGN, ABS and SQR functions, whose BASIC interpretations are shown below. SGN(x) is used to determine the sign of a number (whether it is positive or negative). The function returns 1 if x>0 (positive), -1 if x<0 (negative) and 0 if x=0. SQR(x) returns the square root of a number (<214>x). So, for example, SQR(9) would return 3, because 3*3 = 9. More information on square roots and indices can be found in Appendix B. ABS(x) returns the absolute value of x, in other words the distance of x from 0. This is sometimes useful for checking if a variable is within a certain range. For example <F102>IF ABS(x)⇐3<F> could be used to test whether x is between -3 and +3 inclusive. The combination of these functions and the feedback between x and y (shown in Figure 8.2) on successive applications of the equations makes the Martin process unpredictable and immensely complex. By plotting x against y on the Amiga’s screen, after each iteration, two dimensional Martin fractals can be created. These demonstrate just how chaotic and varied the output from these two simple equations can be. <$fc>Figure 8.2: Feedback in the Martin equations

An Amiga BASIC program to perform the Martin process and simultaneously plot the resulting fractal is shown in Listing 8.1. A GFA version of this program is provided later in the chapter.

<$ob>DEFDBL a-c,x,y,xnew,ynew <$ob>REM Request constants<$ob>INPUT "a=",a<$ob>INPUT "b=",b<$ob>INPUT "c=",c<$ob>CLS<$ob>'<$ob>REM Set initial values for variables<$ob>x=0<$ob>y=0<$ob>i=0<$ob>'<$ob>REM Main loop draws fractal and checks for mouse button press<$ob> WHILE MOUSE(0)=0 'Perform Martin process until mouse button is <$ob>pressed<$ob> '<$ob> REM Calculate new values of x and y<$ob> xnew=y-SGN(x)SQR(ABS(b*x-c))<$ob> ynew=a-x<$ob> x=xnew<$ob> y=ynew<$ob> '<$ob> PSET(300+x.8,100-y*.4) 'Plot point<$ob> i=i+1 'Increment iteration count<$ob> WEND <$fc>Listing 8.1: An Amiga BASIC program to perform the Martin process and simultaneously plot the result

Upon running the program you will be prompted to enter the three constants, remember that it is necessary to click in the output window before Amiga BASIC will allow you to respond to these prompts. The values entered for these determine which of the many possible fractals will be drawn. To pilot the program and to demonstrate a typical Martin, if there is such a thing, the following constants can be used:

<$ob>a=45 <$ob>b=2 <$ob>c=-300 After entering the above values, in response to the relevant prompts, you will soon (during the first 5,000 iterations) see a small group of near-perfect circles appear, connected together by pairs of lines, as shown in Figure 8.3a. <$fc>Figure 8.3a: Martin fractal where a=45, b=2, c=-300 after 5,000 iterations

At this stage it looks more like a network of underground passages than a biological construction. After a short wait (22,000 iterations complete) all of the exits from this network are blocked off and a collection of tiny lines appears in the space outside, see Figure 8.3b.

<$fc>Figure 8.3b: Martin fractal where a=45, b=2, c=-300 after 26,000 iterations A period of inactivity then follows in which the only sign of change is the creation of a circle in the centre of the fractal. This inactivity is suddenly ended by a 'burst' of growth, which occurs after 85,000 iterations have been completed. This growth is much more sustained than before, continuing until around the 100,000th iteration. The scene after this burst is shown in Figure 8.3c and clearly highlights the new, more natural looking, growth dwarfing the original network. <$fc>Figure 8.3c The Martin fractal after 70,000 iterations

Further bursts which occur later thicken up existing lines and fill in gaps, eventually producing the picture shown in Figure 8.3d. The main features of this fractal, such as the high level of complexity and the bursting phenomenon are typical products of this program.

<$fc>Figure 8.3d: The Martin fractal after 1,000,000 iterations Like many other fractals it is very interesting to watch this one grow. By slowing the program down it can be seen that the image is built by pixels being dropped, as if from an invisible satellite in irregular orbit around the centre of the screen. Using this analogy the periods of inactivity seen between the bursts in the above example can be said to be caused by the satellite dropping pixels on top of ones already on the screen, thus making no visible change. This infers that the satellite has predictable periodic behaviour, but this can never be the case with a non-linear dynamic system. In fact, during periods of inactivity, the satellite is following a subtly different path, too similar to be noticed in the crude approximation produced on the screen. However, the same type of positive feedback that causes the butterfly effect amplifies this small deviation and eventually it suddenly becomes very large, causing the satellite to assume a much different path, noticeable on the screen as a sudden burst of new pixels. Return of the Butterfly The notion of non-periodicy and the butterfly effect in a system such as this will already be well known from Chapter 3, but the bursting phenomenon is far more visual here than in the Lorenz attractor. These kinds of erratic flurries of activity are also familiar to share dealers, stock brokers, weather forecasters and other people whose professions are concerned with the world of non-linear dynamics. In these areas of the real world such bursts manifest themselves as sharp swings in share prices, sudden drops in temperature and rapid increases in oil prices. It is highly probable that such sudden unexpected events are created by previously unnoticed small changes getting larger by ever increasing amounts. The now legendary stock market crash of 1987 was apparently due to dealers over-reacting to a small devaluing of stocks, sending shares into a downward spiral. The fall was further magnified by the newly installed computer dealing systems which sold shares instantaneously, without the protective natural damping provided by slower human reactions. Nobody could have noticed the start of this spiral in advance because it would have been too small, meaning that such events can only be spotted when their magnitude is such that they are hard to prevent. A similar manifestation of the butterfly effect could cause a system, apparently continuing normally, to suddenly collapse without warning. This is an interesting suggestion and has precipitated much debate in recent years. The fact that the weather and the Earth’s eco-system are vulnerable to this type of positive feedback makes the possibility of a sudden end to the human race quite feasible. Further Experimentation with the Martin Program The number of different images that can be produced with Listing 8.1 is truly daunting, so a selection of some of the most interesting constants are shown in Table 8.1. Note that although most of these examples will fit on the screen using the current scaling, some of them may benefit from being enlarged. This can easily be done by changing the co-ordinate multipliers in the PSET line. For example the following line will produce output of twice the size: <$ob>PSET(300+x*1.6,100-y*0.8)

<$tc>Table 8.1: Some interesting constants for the Martin program <$table3>a b c <$table3>68 75 83 <$table3>90 30 10 <$table3>10 -10 100 <$table3>-200 -4 -80 <$table3>-137 17 -4 <$table3>10 100 -10 <$table3>-137 17 4 The position of the image on the screen may also need to be altered. This can be done by editing the two large numbers in this line, as it is these which collectively determine the offset. An example of a Martin enlargement is shown in Figure 8.4. <$fc>Figure 8.4: Mathematical enlargement of Figure 8.3d

Extensions to the Martin Program

There is very little point in adding an expeditious user interface to the Martin program because very little input is required. However, there are a few simple enhancements that can be included in other areas of the program.

Colour

If you have a colour monitor you may want to add colour to the Martin fractals in order to gain a fuller understanding of the process which creates them. To do this there must be something that the colours can represent, as discussed in earlier chapters. Here they can be used to show the number of times a particular pixel position has been visited. As well as making the fractals more aesthetically pleasing this also makes it possible to see which points are being plotted at any given time, even during a period of inactivity.

The simplest, but probably slowest, way to facilitate colour is to work out the position of the next pixel to be plotted, get the colour of the screen at this position using the POINT function, put it in a variable, increment it and then plot the pixel in the colour denoted by the variable. By doing this the colour of a point on the screen will change every time a pixel is plotted on it. Using this method has the advantage that everything is stored on the screen, so no extra variables or arrays need to be declared.

Exactly the same technique was used in Chapter 2 to add colour to the Feigenbaum diagram. It can by implemented in the Martin program by making the alterations shown below . Note that, because of the 16 colour limit, the variable containing the colour value is of modulus 16, so after a point has been visited 16 times it returns to being white and the colour cycle starts again.

The following lines should be included at the head of the program, to open a window on a high resolution 16 colour screen:

<$ob>SCREEN 1,640,200,4,2 <$ob>WINDOW 2,"Colour Martin Fractal",(0,0)-(617,180),15,1

The PSET line which plots the points which comprise the Martin fractal should be replaced by the program fragment shown below:

<$ob> xp=300+x*0.8 !Calculate x and…​ <$ob> yp=100-y*0.4 !…​y position of point <$ob> oldcolour=POINT(xp,yp) !Determine old colour <$ob> COLOR (oldcolour+1) MOD 16 !Set new colour <$ob> PSET(Xp,Yp) !Plot the point Sound If you have recovered from the musical Feigenbaum diagram in Chapter 2 you may want to add sound to the Martin program. This can be very useful for alerting you to pixels being plotted on the screen which would otherwise have gone unnoticed. The most efficient way to create such a pixel alarm is to start a tone whenever a new pixel is plotted, and halt the tone when one is plotted on an already occupied space. Again, the POINT function can be used to check the intended plotting position to see whether it is free. To add the alarm the routine below should be included in the program given in Listing 8.1, in place of the PSET line. Note that, in an effort to increase program speed, no pixel is plotted if the chosen space is already full. <$ob> xp=300+x*0.8 !Calculate x and…​ <$ob> yp=100-y*0.4 !…​y positions <$ob> IF POINT(xp,yp)=0 !If position is empty…​ <$ob> SOUND 500,1 !..beep and…​ <$ob> PSET(xp,yp) !…​plot point <$ob> END IF !Otherwise do nothing This technique allows even the smallest of additions to the fractal to be noticed, and bursts of activity can easily be distinguished due to their characteristically continuous tones. This will relieve you of the boredom of having to sit and stare at the screen in order to be sure of seeing all of the important changes. Making the pitch of sound produced proportional to the x or y variable may be a way of further enhancing the program. This would allow the listener to judge where the point was being drawn and would also give a second way of conveying the nature of chaotic movement. In the ultimate Martin sound system the Amiga’s stereo effects could be used to express the horizontal position of each point as it was drawn. <$hb>The GFA BASIC Solution

The number and variety of images which can be generated from the simple Martin equation is truly astounding, making experimentation very worthwhile. The GFA BASIC program given in Listing 8.2 is essentially the same as the Amiga BASIC one in listing 8.1, but contains several special facilities to aid experimentation.

9. The Future

Although the initial chaos boom has now subsided, the outlook for applied chaos research looks promising. There are now many dedicated chaos researchers and it is very probable that the number of prominent scientists interested in chaos will grow dramatically in the future as the older generation of scientists is replaced by one more willing to use computers.

<$ha>Can the Future of Chaos be Predicted? In any science there is normally a time lag between leading edge research and the state of the amateur scientist’s contribution, due to professional researchers having superior resources. This is true in chaos, and by browsing through recent publications such as those listed in the bibliography it is possible to predict the nature of future home computer based research. Because the advances in computer technology are so fast the chaos time lag is shorter than in most sciences, at around 10 years. To provide an insight into the home computer chaos research of tomorrow a few of the more interesting areas of today’s professional activity are discussed below. <$hb>Fractal Maths in Data Compression Applications

The advantages of data compression are obvious, money can be saved in telephone bills when sending files via modem and for any given file less storage space is required on a computer’s storage device. After discovering that something as complex as the Mandelbrot set can be described by one simple equation (see Figure 9.1a) scientists speculated that all images may be able to be treated like fractals and compressed so as to be described by a small equation, as shown in Figure 9.1b. Theoretically such compression is possible, but with current technology the complex fractal mathematics involved in compressing a single television quality image in software takes such a long time that it is of no practical use.

<$fc>Figure 9.1a: Compression and decompression of the Mandelbrot set to and from its equation <$fc>Figure 9.1b: Mathematical compression and decompression of a photograph

Recently, however, an American team co-led by prominent chaos writer Michael Barnsley claims to have developed an IBM PC-based image compression system which can compress and decompress moving VGA graphics (640x480 pixels, 16 colours) in real time. The system can run in real time because most of the processing is carried out by a purpose-built expansion card, which has its own high speed microprocessor. The incredible speed is matched by the reduction in data size, which allows almost one minute of moving video images to be squeezed onto one 800k floppy disk. The extent to which the system relies on state-of-the-art hardware is reflected in the price of $25,000, which obviously puts it out of the reach of most PC users. Since a computer image is constructed from digital data it should be possible to extend this technique to other forms of computer data, including text and program code. Performing such data manipulation on the Amiga would be almost impossible, considering that the decompression of the Mandelbrot set, from its equation, can take several hours (depending on the resolution and language used). <$hb>Telephones with Minds of Their Own

The vast telecommunications network that now stretches across the world has also been attracting the attention of chaos researchers. It has been suggested that tiny bursts of static interference in telephone lines can be amplified by positive feedback at telephone exchanges and other installations causing them to grow into large, audible sounds. Arthur C. Clarke first proposed this idea many years ago, and it has now been embraced by science-fiction fans who can imagine the world telephone network acting like a complex biological organism, dialling up people and sending them intelligible messages and faxes at will. In reality feedback on this scale does not happen but stories of people’s telephones doing strange things, normally explained with reference to ghosts and the underworld, could well be caused by this typical manifestation of the butterfly effect. This type of chaos is not unique to telephone lines though, it could also occur in video and computer networks.

<$hb>The New Art The recent publication of a range of books dedicated to fractal generated plants, animals, landscapes and abstract patterns has helped to establish fractal graphics as an important art form in its own right. This is one area of current interest which can be explored with humble home computers like the Amiga, as described in Chapters 6, 7 and 8. Adding picture saving routines to the programs given in this book make it possible to load fractal images into standard Amiga art packages and create a wide range of original compositions, an example of which is shown in Figure 9.2. Chaos can also be applied to one of the other popular programming areas, ray tracing. If you have ever experimented with ray tracing shading algorithms you will know that creating the illusion of a naturally rough surface, such as a carpet or road, is no easy task. An interesting technique for overcoming the problems normally associated with such surfaces is to use fractal landscape generation techniques, shown in Chapter 7, to create fairly flat landscapes on surfaces designated as being rough. The light incident on such surfaces will then appear to be reflected naturally. Such sophisticated ray tracing on the Amiga would be very time consuming, but it is not impossible. <$ha>What Use is Chaos?

As discussed in Chapter 4, the computer press have been keen to demonstrate the use of chaos for drawing pictures but are reticent about the relevant theory. Even when the implications are discussed the attitude adopted is generally negative, with emphasis being placed on chaos hampering weather predictions and possibly being associated with the end of civilisation. However, chaos is like fire, it is unpredictable and destructive. Yet if its power is tapped it can be used in a variety of useful applications. Although fractal maths has been used in image compression, few applications have yet been written which actually use chaos. An example of how powerful a tool chaos could be is clearly visible in nature, in the anti-viral immune system of the human body.

When a virus enters the body the immune system produces a random variety of anti-bodies and sends them into the blood stream. The progress of each type of anti-body is then monitored and after a period of time the identity of the successful one is fed back to the immune system. The body then produces the thousands of clones of the successful anti-body which are necessary to kill off the virus. This use of randomness and feedback (i.e. chaos) provides the body with a system which can cope with an enormous range of viruses, even new strains to which the body is not accustomed. This method has significant advantages over the equivalent deterministic process, which would need to know the details of every virus that currently exists and all those which will exist in the future.

The evolutionary process of natural selection also uses chaos, to preserve life on earth. New organisms are produced randomly through mutations and the success of each organism is fed back into the system in line with the 'survival of the fittest' theory. By moving forward on a broad, ever changing front nature can ensure that life still exists even in the bleakest of conditions, when the dinosaurs became extinct for example.

<$hb>Derivation of Pi Using the Monte Carlo Method The power of controlled randomness can also be demonstrated by considering the challenge of finding the value for the commonly used mathematical constant, pi (<80>). Interrogating a pocket calculator, or GFA BASIC, will give a pre-determined value for pi in the region of 3.141592654 but pi is irrational, meaning that an infinite number of decimal places are required to convey its value accurately. Many people have experimented with pi over the centuries and many different ways of deriving its value have been devised, the most accurate of which recently yielded an answer with 10 billion decimal places. Not surprisingly, most of these methods are deterministic but one, known as the Monte Carlo method, can be based on something as random as the scattering of raindrops on the ground. Before executing this method it is necessary to find a large area of unsheltered flat land, on which the pattern shown in Figure 9.3 should be drawn as accurately as possible. Before describing the next stages of the process some peculiarities of the pattern must be considered. As the figure shows, it consists of a circle, of radius two metres, inside a square of side four metres. Remembering that the area of a rectangle is its height multiplied by its length and that for a circle it is pi multiplied by the square of the radius: Area of square = 4*4 = 16 Area of circle = PI*(2)^2 = 4*PI The ratio of the area of the circle to the area of square is therefore: 4*PI:16 = 4*PI = PI or PI/4 `16 4` Note that the actual size of the pattern on the ground in not important, as long as the ratio of the area of the circle to the square is constant at pi/4. However, the larger the area the more accurate the results will be for reasons discussed later. After inscribing the pattern on the ground it is necessary to wait for a light rain shower to occur. After the shower the total number of raindrops which landed in the entire square area (the number of drops captured) should be recorded. Of these, the number which fell inside the circle (the number of hits) should also be noted. The hit ratio, the ratio of hits to raindrops captured, should then be calculated. Because the raindrops can be assumed to have fallen totally randomly we can say that they were distributed evenly across the square area. Since the ratio of the area of the circle to the area of the square was calculated above as pi/4 it follows that the ratio of the number of hits to the number of raindrops captured also equals pi/4. Therefore, by multiplying the hit ratio by four, pi can be determined. In practice this method would prove very time consuming, as a large area would be needed to determine pi with any accuracy meaning that a large number of raindrops would need to be counted. However, the outcome of the process can be speeded up by modelling it in BASIC, this also eliminates the problems associated with waiting for a light shower and ensuring that the raindrops do not evaporate before the counting is complete. <$fc>Figure 9.3: Pattern used when determining the value of pi

A program to simulate the process and calculate the result after each drop has landed is shown in Listing 9.1. This uses a pattern of the same size as that shown in Figure 9.3. This size is convenient because the circle is the same as the one used to produce the Mandelbrot set, so exactly the same equation can be used to determine whether a given drop fell inside or outside the circle.

<$ob>DEFDBL result <$ob>DEFLNG throws,hits <$ob>throws=0 'Set initial throw count…​<$ob>hits=0 '…​ and hit count TO zero<$ob>'<$ob>WHILE 1 'Repeat indefinitely<$ob> '<$ob> x=RND*4-2 'Select random x position<$ob> y=RND*4-2 'Select random y position<$ob> '<$ob> throws=throws+1 'Increment the throw counter<$ob> '<$ob> IF x2+y2<4 THEN 'If point lies inside circle…​<$ob> hits=hits+1 '…​increment hit counter<$ob> END IF<$ob> '<$ob> result=4*hits/throws 'Calculate 4*ratio (i.e. PI)<$ob> LOCATE 1,1<$ob> PRINT "Throws:",throws<$ob> PRINT "PI:",result 'Print result<$ob>WEND <$fc>Listing 9.1: Simulation of determining PI using raindrops

After initialisation of the raindrop and hit counts, the program enters a continuous cycle in which a point on the board is selected at random within the range -2<x<2 and -2<y<2 and the raindrop count (held in <F102>rain_drop<F>) is incremented. The point is then tested and, if it fell inside the circle (x2+y2<4), the <F102>hits<F> variable is also incremented. After each drop has landed (i.e. after each iteration) the ratio of <F102>hit<F> to <F102>rain_drop<F> is calculated and multiplied by 4 to give the value of pi, which is printed on the screen. After several thousand iterations the first few digits of pi have been determined.

This is one of the slowest processes for deriving pi but it does show how randomness can be put to a useful purpose. If civilisation knew no other way of calculating pi this method would obviously be very important. The model produced by Amiga BASIC is flawed, however, because the random number generator is not truly random causing accuracy to be limited. Also the random nature of the algorithm causes the result to drift annoyingly away from the real value of pi during some stages of processing, although the laws of statistics lead us to believe that the accuracy of the result does generally increase with the number of iterations performed.

The way that the ratio drifts randomly around pi/4 is shown in Figure 9.4, which clearly demonstrates that pi/4 is the attractor of the process. The graph is similar to the deterministic time series shown in Chapter 2 (Figure 2.7a), where the population was attracted to 100%, but the oscillations around the pi/4 attractor are much less ordered due to their random origins.

<$fc>Figure 9.4: Graph showing the hit ratio (drawn as a solid line) drifting around the real value of pi/4 (dotted line) Chaos techniques such as those used in the immune system would be of particular use in artificial intelligence applications where, at present, all possible scenarios must be considered when developing a set of rules to solve a certain problem. It is even conceivable that future high level computer languages might have built in chaos commands. <$ha>Life After Chaos and Fractals on Your Amiga

I have covered most of the major facets of chaos theory in this book, but the constantly expanding nature of the science means that there are still many more experiments which may be performed on the Amiga. Even during the closing stages of this book’s production new theories and ideas that could have been included were being brought to my attention. A large number of further experiments can be based on those in this book, or on the purpose-written support disk (see front of book). Some ideas for such experiments are discussed in the section below.

<$hb>Two Headed Mandelbrots and Similar Feigenbaums There is an interesting piece of information that can be found about the Feigenbaum diagram which forms a good basis for a challenging programming problem. By measuring the length of each section of the diagram, results similar to those given in Table 9.1 may be obtained. <$tc>Table 9.1: Length of different sections of the Feigenbaum diagram

<$table3>Period Length/mm Ratio <$table3>1 690 <$table3> 4.6 <$table3>2 149 <$table3> 4.8 <$table3>4 31 <$table3> 4.4 <$table3>8 7

From these it is possible to work out the ratio between the length of one section and the next. Note that the size of the diagram being measured is irrelevant but the bigger the diagram the more accurate the results. If the measuring stage of the experiment is performed accurately enough all the ratios should be of exactly the same value, approximately 4.7 (the mean average of those in the table is 4.6). Mitchell Feigenbaum was the first to discover this common ratio, he managed to calculate it as being approximately 4.669201660910299097.

This value is only approximate because, like the pi constant (<80>) discussed above, it is an irrational number, i.e. it has an infinite number of decimal places. To save the trouble of writing out all the decimal places in full mathematicians have labelled it with the Greek letter delta (<100>).

Even though the results in Table 9.1 were taken from an enormous print-out of the Feigenbaum diagram the value of delta is still quite inaccurate. A more sensible way of calculating the ratios is by using a BASIC program to determine the exact value of c associated with each bifurcation, and then use these values to calculate delta. A program of this sort would take a long time to execute, but would make a very worthwhile overnight project.

<$fc>Figure 9.5: 'Spiky' Mandelbrot set After several years of home computer activity related to the Mandelbrot set it could be assumed that there is nothing further to explore. However, even with this veteran fractal there are many more experiments to conduct. Some variations on the set, such as those shown in Figures 9.5 and 9.6, can be generated by altering the Mandelbrot programs without any preconceived plan. In general, however, it is better to have at least some idea of what effect alterations are likely to have. Some suggestions for simple variations on the Mandelbrot set are described below. <$fc>Figure 9.6: Mandelbrot 'web' (based on Listing 4.5)

The first variation on the original Mandelbrot set is caused by altering the lines in the Mandelbrot program which set the initial values of p and q, so that these variables are something other than zero. The changing of these initial values is similar to the Julia process described in Chapter five, but in that process a and b are treated as constants while the initial values of p and q vary according to the position of the pixel being plotted. In the altered Mandelbrot process, however, a and b vary as normal and the initial values of p and q are constants as normal.The result of such a subtle change to p and q can be achieved by replacing the original variable assignment lines in a Mandelbrot program, such as Listing 4.4 with:

<$ob>p=0.6 <$ob>q=0.2

The output from the resulting program (see Figure 9.7) looks like a cross between a Julia and a Mandelbrot set, reflecting the hybrid process. There are a vast number of different sets that can be creating using this method, due to the many possible combinations of p and q. However, an important thing to remember when experimenting in this way is that both variables should always be in the range of -2 to 2.

<$fc>Figure 9.7: Mandelbrot set with amin=-2,amax=2, bmin=-2, bmax=2 and initial values for p and q of 0.6 and 0.2 respectively Another variation on the set can be achieved by altering the Mandelbrot equation which, recalling Chapter 4, was as follows for the complex number method: znew = z^2 + c This could be expanded to: pnew = p^2 - q^2 + a qnew = 2*p*q + b The second Mandelbrot set variation is produced by changing the index of z from 2 to 3. The new equation, with its resulting expansion, is as follows: znew = z^3 + c pnew = p^3 - 3*p*q^2 + a qnew = 3*p^2*q - q^3 + b The two equations for pnew and qnew can easily be incorporated into a Mandelbrot program such as the one given in Listing 4.4, and will produce output similar to that shown in Figure 9.8. Although the more complicated calculations slow the program down this is still an area worthy of experimentation, try raising z to the power of 4, 5 or even a floating point number (mathematical skills are required to achieve this). <$fc>Figure 9.8: 'Two-headed' Mandelbrot set with amin=-2, amax=2, bmin=-2, bmax=2

<$fc>Figure 9.2: Original fractal composition produced using Deluxe Paint and programs from this book The increased popularity of chaos in the last few years has resulted in many new works becoming available, which collectively describe all the major aspects of chaos theory. If you are interested in creating fractals which have not been covered in this book, or if you want to find out more about chaos history, you will find additional information in the bibliography provided later in the book. It contains brief descriptions of the most notable titles, suggested as a logical progression from this book. 10. Appendix A <$cn>Maximising the Amiga’s Features in Chaos Applications

Throughout this book the example programs have been kept as short as possible to save room and to keep to a minimum the amount of typing, in an attempt to make the theory and principles as clear as possible. This has been done partly by gathering together the most widely used procedure definitions and techniques and placing them in this appendix, thereby avoiding unnecessary repetition.

This appendix provides suggestions of how the Amiga’s powerful features can be used to the full in chaos programs. Note that the routines shown here are not specific to drawing fractals, and you may well want to apply them to other applications.

<$ha>Graphics The high colour, high resolution, graphics facilities provided by the Amiga’s custom chips are undoubtedly one of the computer’s most notable features. The programs in this book can easily be extended in order to take advantage of all the available graphics facilities, allowing the creating of truly stunning fractal images. All of the examples in this book have been implemented in the non-interlaced high resolution mode, giving a resolution of 640x200 pixels and a maximum of 16 colours. This mode was chosen because it is suitable for all types of monitor. However, experienced Amiga programmers will know that higher resolution modes, and more colours are available. If you don’t mind having fewer colours, the 640x400 pixel resolution provided by the interlaced high resolution mode is ideal for chaos work, particularly if it is to printed out. On the other extreme, there is the non-interlaced low resolution mode, providing up to 32 colours, but only up to a resolution of 320x200. Unlike most other computers, it is not necessary to actually change screen modes on the Amiga, instead a new screen of the required mode can be opened, which will co-exist with screens which are already open. An interesting side effect of this feature is that we can open, for example, a high resolution monochrome and low resolution colour screen and plot to both simultaneously. This has the advantage that two different plots may be created from a single set of calculations, and thanks to the Amiga’s powerful graphics chips plotting to two screens instead of one doesn’t noticeably slow the program down. A practical example of this technique is provided later in this appendix. The decision regarding which mode to use for fractal work is difficult, and varies from fractal to fractal. The main considerations to take into account are summed up below: <$bullet>In lower resolutions some pixel-based fractals (such as the Mandelbrot set) are drawn much faster, as it is not necessary to calculate the behaviour of as many points.

<$bullet>Drawing times of fractals such as the Lorenz attractor and Martin images are unrelated to the screen mode being used. <$bullet>Carefully chosen colours can make an image appear to be of a higher resolution than it actually is.

<$bullet>Monochrome screen dumps of colour images look terrible, it is much better to aim for the maximum resolution if you intend to produce screen dumps. <$hb>Amiga BASIC Graphics Facilities

Amiga BASIC provides the essential foundations for utilising the Amiga’s graphics facilities. Before plotting in a certain mode it is first necessary to open a screen of the relevant type and then to open a suitably sized window to which BASIC output can be sent.

As an example, a 32 colour low resolution (320x200) screen, and accompanying window, can be opened using the following lines:

<$ob>SCREEN 1,320,200,5,1 <$ob>WINDOW 2,"32 Colour Low Resolution Mode",(0,0)-(297,180),15,1

A 16 colour high resolution, interlaced, screen (640x400) can be facilitated by the lines below:

<$ob>SCREEN 2,640,400,4,4 <$ob>WINDOW 3,"High Resolution 16 Colour Mode",(0,0)-(617,380),15,2

A detailed discussion of the SCREEN and WINDOW commands can be found in Amiga BASIC manual.

Since the example programs provided throughout this book assume a screen resolution of 640x200 it is necessary to change the scaling used to transform the small numbers produced by the fractal equations to large numbers suitable for plotting on the screen. To illustrate this, and the dual screen plotting methods mentioned earlier, a dual mode Julia set plotter is provided in Listing A.1.

<$ob>REM Open low resolution screen<$ob>SCREEN 1,320,200,5,1<$ob>WINDOW 2,"32 Colour Julia Set",(0,0)-(297,180),15,1 <$ob>'<$ob>REM Open high resolution screen<$ob>SCREEN 2,640,400,4,4<$ob>WINDOW 3,"High Resolution 16 Colour Julia Set",(0,0)-(617,380),15,2<$ob>'<$ob>REM Set constants<$ob>a=-1.16<$ob>b=-.25<$ob>'<$ob>REM Plot set<$ob>FOR x=-2 TO 2 STEP .01<$ob> FOR y=-2 TO 2 STEP .01<$ob> p=x<$ob> q=y<$ob> iteration=0 <$ob> ' <$ob> REM Perform Julia calculations<$ob> WHILE p*p+q*q<4 AND iteration<32<$ob> pnew=p*p-q*q+a<$ob> qnew=2*p*q+b<$ob> p=pnew<$ob> q=qnew<$ob> iteration=iteration+1<$ob> WEND <$ob> '<$ob> REM Plot high-res point<$ob> WINDOW OUTPUT 3 'Select high resolution window<$ob> COLOR iteration MOD 16<$ob> PSET(320+x*100,200-y*100) <$ob> ' <$ob> REM Plot low-res point<$ob> WINDOW OUTPUT 2 'Select low resolution window<$ob> COLOR iteration MOD 32<$ob> PSET(160+x*50,100-y*50) <$ob> '<$ob> NEXT y<$ob>NEXT x <$fc>Listing A.1: A dual screen Julia set plotter based on Listing 5.1 The first part of this program comprises the two program fragments given above, these are used to open the two screens and windows. After this the constants for the set are initialised, by default these are a=-1.16 and b=-0.25 but other values can easily be substituted in their place. The main calculation loop then follows, this is similar to the loop in Listing 5.1 except for the inclusion of two point plotting routines instead of one and a halving in the STEP size of the FOR…​NEXT loop which controls the value y. The STEP size has been reduced as more detailed calculation must be performed to make maximum use of the high resolution screen. Note that in programs such as this it is always necessary to adapt the detail of calculation to cater for the highest resolution being used, therefore slowing the program down to the speed of the higher resolution. The high resolution plotting routine is very similar to the one originally used in chapter five but notice that window 3, the window on the high resolution screen, must be selected as the current output window before the point is plotted. The scaling used in the PSET line here is also different from that found in Listing 5.1, the y position of the point is multiplied by twice as much in this version to compensate for the doubling of the vertical resolution. This is necessary because the screen size and aspect ratio are different for the high resolution interlaced screen than for the non-interlaced screen used originally. The low resolution plotting routine contains similar changes, the relevant window must be selected and the scaling is unique to the low resolution (320x200 pixel) screen. In addition the colour is of MOD 32, rather than MOD 16, reflecting the increased number of colours available in the lower resolution mode. <$hb>GFA BASIC Screen Modes

GFA BASIC allows both of the modes discussed above for Amiga BASIC, the relevant program code for using these modes is given below. Note that in GFA BASIC there is no need to open a window before commencing plotting, in addition the vertical resolution of the screens is 28% larger than in Amiga BASIC. Of course the changes in scaling associated with altering the screen mode, discussed above, are also required in GFA BASIC.

For 640x512 pixels with 16 colours:

<$ob>OPENS 1,0,0,640,512,4,32772 For 320x256 with 32 colours: <$ob>OPENS 1,0,0,320,256,5,0

Also, a low resolution 320x200 pixel mode, with 64 colours, can be made possible using the lines below:

<$ob>OPENS 1,0,0,320,256,6,128 In addition to these modes, and the other standard Amiga modes, GFA BASIC has provision for a special screen mode called Hold And Modify (HAM), this gives a low resolution of 320x200 pixels, but an incredible 4096 colours, more than can actually be used in most fractals. <$hb>Colour Cycling

Colour cycling is the name given to the effect often used to animate fractals, it basically relies on rapid changes of the hardware colour palette. This effect is a favourite of PC users, who have access to the public domain Fractint program. There is no mathematical value to be gained from colour cycling, it is simply included for aesthetics.

Listing A.2 demonstrates the GFA BASIC techniques used in colour switching, it opens a 32 colour low resolution screen and draws 64 concentric circles, one in each of the foreground colours. Colour cycling then begins, creating a ripple effect emanating from the common centre of circles.

<$ob>OPENS 1,0,0,320,256,6,128 !Open 64 colour screen<$ob>TITLES #1,"Colour Cycling With 64 Colours"<$ob>'<$ob>FOR c=63 DOWNTO 1 !Draw 64…​ <$ob> COLOR c !…​uniquely coloured…​<$ob> PCIRCLE 160,128,c*2 !…​concentric circles<$ob>NEXT c<$ob>'<$ob>DO<$ob> FOR mix=4095 DOWNTO 0<$ob> FOR c=1 TO 63 !Change the shade of colour c …​ <$ob> SETCOLOR c,(mix+c) MOD 4096 !…​to one of 4096 shades in turn…​<$ob> NEXT c !…​thus creating the illusion of movement<$ob> NEXT mix<$ob>LOOP !Repeat cycle forever <$fc>Listing A.2: A simple demonstration of colour cycling in GFA BASIC

<$ha>Printing Fractals If you have a suitable printer connected to your Amiga and the relevant options have been set in the Workbench Preferences program a screen dump of the current window can be created from GFA BASIC by including the line below at the end of any fractal program: <$ob>HARDCOPY

<$ha>Useful Routines The programs in this book have generally been kept as short as possible to make experimentation with them less of a chore. The few routines discussed below have been specially selected to aid augmentation of the programs for ease of use. For example, in GFA BASIC the routine which waits for a mouse button to be pressed allows a fractal image to be viewed after completion without the End of program alert requester getting in the way. All of the routines given here should be appended to the program in which they are used and called from wherever needed in the recommended way. <$hb>Waiting for a Mouse Button Press

In Amiga BASIC the sub-program definition below provides a simple and compact way of waiting for a press of the left mouse button.

<$ob>SUB WAITMOUSE STATIC<$ob>WHILE MOUSE(0)=0<$ob>WEND<$ob>END SUB

The line below would be used to call the program:

<$ob>WAITMOUSE In GFA BASIC the built-in MOUSEK variable reports the mouse button status. The procedure below utilises this in order to wait for a press of the mouse button. <$ob>PROCEDURE waitmouse<$ob> REPEAT<$ob> UNTIL MOUSEK>0<$ob>RETURN The procedure is most suitably called using: <$ob>@waitmouse

<$hb>Rounding Numbers to the Nearest Integer The INT function provided by Amiga BASIC’s always rounds numbers down, for example INT(2.99) would give an answer of two, even though the nearest value is obviously three. In normal uses INT is acceptable but in fractal programming, where attention to mathematical detail is of great importance, a more accurate rounding algorithm is required. The program fragment below provides this in the form of a compact function. DEF FNround(float) = SGN(float)(INT(ABS(float))+INT(2(ABS(float)-INT(ABS(float))))) Note that GFA BASIC provides the ROUND() function for rounding numbers to the nearest integer, so a routine such as this is not necessary. 11. Appendix B <$cn>Mathematics in BASIC

This book has been written to appeal to a wide audience of Amiga users, from those who are accomplished chaos mathematicians, but are unable to incorporate the theory into Amiga programs, to programmers who are not familiar with the relevant mathematics. This appendix is for reference by those in the second category, and explains the elementary concepts referred to in some of the earlier chapters. If further assistance is required in this area then a good GCSE or O Level mathematics text is recommended, as I have attempted to keep the majority of this to that level.

Topics are discussed roughly in order of importance, with particular attention being paid to their BASIC implementations.

Simple Algebraic Principles and Calculations

Algebra encompasses a wide range of fundamental mathematics principles, but probably its most well know aspect is the use of letters to represent numbers. This allows written statements to be expressed in mathematical form. The letters are known as variables, and have similar properties to the variables found in programming languages, except that they may only be one character long. Algebraic variables may be used in all the same ways as normal numbers. For example the algebraic expression shown below involves addition and multiplication.

y + y + y = 3 * y

The 3*y section would normally be written simply as 3y, but the BASIC notation adhered to in this book means that the longer equivalent must be used. The equation states that any number (here represented by y) added to itself three times equals 3 multiplied by the number. To use this equation values must be substituted in place of y, as an example the equation when y is four is shown below. Note that in BASIC a number can be assigned to a variable using a command of the form <F102>y=4<F>.

4 + 4 + 4 = 3 * 4

This equation holds for any number, and since it is an equation both sides always equate to the same value. Although this example is of little practical use there are many useful equations expressed in a similar form, such as Pythagorus' theorem.

Variables whose values are always the same are known as constants. By convention different variables are used for different number types. For example x, y, z and i, j, k are usually used to represent variables whereas a, b, c, k and Greek letters such as <80> are employed to represent constants. When used alone, z is traditionally used to represent complex numbers, whose constituent parts are usually represented by a for the real part and b for the imaginary one. Apart from having different names there is no distinction between constants and variables in most dialects of BASIC but in C it is possible to declare variables as being static, meaning that they are treated as constants, which are processed faster.

z = a + bi

Operator Priorities

In equations with more than two terms it is often important to ensure that the priority of the operators is taken into account when working out the answer. If this is not done then unexpected answers may be produced when carrying out calculations on the computer which, like a calculator, always uses the correct priorities.

For example, the two calculations below yield the same answer, despite the fact that the numbers are in a different order.

<$table2>1 + 2 * 3 = 1 + 6 = 7 NOT 3 * 3 = 9 <$table2>2 * 3 + 1 = 6 + 1 = 7 NOT 2 * 4 = 8

The two calculations equate to the same value because multiplication (and division) has a higher priority than addition (and subtraction), so all multiplying is done first. If the addition needs to be done first brackets must be placed around the addition section in order to increase its priority. The result of the addition of brackets to the last example is:

(1 + 2) * 3 = 3 * 3 = 9

2 * (3 + 1) = 2 * 4 = 8

<$HB>Ratios A ratio is a way of expressing the relationship between two numbers (or two sets of numbers) without quoting their absolute values. For instance the scale of a model aircraft might be said to be '1 to 48', often written as 1:48. This scale is the ratio of the size of the model to the size of real aircraft, and in this case it states that for a single unit length on the model the equivalent real length is 48 unit lengths. Because the ratio applies to the whole plane it is possible to work out the dimensions of any part of the real aircraft using measurements from the model. For instance if the wing-span of the model is 20cm the wing-span of the real aircraft would be 48*20cm = 960cm or 9.6m. The colon character (:) is not generally used in maths so ratios are expressed as fractions instead, for example 1:48 would be written as 1/48. This makes manipulation of ratios easier because to scale a length up we just divide the model’s length by the ratio, and when scaling down the lengths are multiplied by the ratio. For example the wing-span of the full size aircraft could be calculated as follows (recalling that the model’s wing-span is 20cm): 20cm / (1/48) = 960cm If the fuselage length of the real plane was known to be 816cm the following calculation would be used to calculate the length of the model: 816cm * (1/48) = 17cm Of course 1/48 could instead be expressed as the decimal number 0.02083. Fractional ratio notation is of particular use when a ratio doesn’t contain a 1. For example if the scale of the aircraft above was 3:48 the model’s wing-span could not simply be multiplied by 48 as it was for the 1:48 scale. <$HB>Fractions as Percentages

Per cent literally means per 100, so percentages are usually given as integers between 0 and 100, for example 93%. However, in maths, and in programming, the per cent character (%) is meaningless in this sense so fractions derived from a ratio are used instead:

<$table3>Percentage Ratio Decimal <$table3>93% = 93/100 = 0.93

Using fractions makes percentage calculations much easier, for example to find 93% of 500 we can simply multiply the decimal equivalent of the percentage by 500:

0.93 * 500 = 465

Fractions can be used to represent any percentage and the list below demonstrates most of the common forms.

<$table2>2% = 0.02 <$table2>25% = 0.25 <$table2>50% = 0.50 (usually written as 0.5, without redundant zero) <$table2>100% = 1 (decimal part is totally redundant here) <$table2>137% = 1.37 <$table2>12.5% = 0.125 (decimal part may be extended indefinitely)

<$ha>Graphs and Co-ordinate Geometry The only way of creating pictures of any kind in mathematics is by using graphs, and fractals are no exception. A typical graph is a two-dimensional plane on which points and lines can be plotted, like the one shown in Figure B.1. Each point on the plane has a unique two part co-ordinate associated with it, consisting of the horizontal and vertical positions of the point. These positions are given numerically with reference to a pair of graticuled lines, called axes, which run through the graph at right angles to one another. The horizontal axis is generally called the x axis and gives the x position of points on the graph. Similarly the y position of points is read from the vertical y axis. Thus the example point in Figure B.1 can be said to be at the position where x=2, y=3. For neatness the x and y values are normally put together in a single co-ordinate packet of the form (x,y). In the current example this would be written as (2,3). <$fc>Figure B.1: Typical two dimensional graph

A typical application of a two dimensional graph is the display of time series charts, where a quantity (on the vertical axis) is plotted against time (on the horizontal axis). A selection of time series graphs can be found in Chapters 2 and 3.

The graph used in the examples above was a standard two dimensional one, but one dimensional graphs with a single axis and three dimensional ones with three axes (x, y and z) are also used in this book. These all work on exactly the same principle. In fact as many dimensions as necessary can be used as long as there is one axis to represent each dimension.

Note that x and y are the normal variable names for the horizontal and vertical positions of a point but any variables can be used. To find out which variable represents which direction you can simply look at the axes, as these are labelled with the variable names.

The Amiga’s screen can be thought of as a graph, with the axes shown in Figure B.2. Graphs can be drawn on the screen by specifying point and line positions using commands like PLOT and DRAW. Unfortunately, however, the range of the screen axes is fixed, so graphs need to scaled up and down in order to make maximum use of the screen area. Also, the y axis is inverted making direct graph plotting rather difficult. A bonus of the Amiga’s screen (in low and medium resolution modes) is that in addition to a position co-ordinate each point also has a colour associated with it. This allows more information to be displayed in a smaller space, for example it can be used to represent the third dimension of the Lorenz attractor.

<$fc>Figure B.2: The Amiga’s screen as a graph Indices Commonly known as powers, indices are a very useful mathematical tool, used in many applications. An example of index notation is shown below. The number shown in normal sized type is called the base and the smaller superscript number is known as the index. 24 The equivalent BASIC notation (used in this book), where the base and index are separated by a ^ character, is: 2^4 We say that this example shows 2 to the power of 4, or that 2 is raised to the power of 4. In BASIC the ^ character can effectively be pronounced to the power of. Certain indices have particular names associated with them. For example a number raised to the power of 2 is said to be squared and a number raised to the power of 3 is said to be cubed, e.g. 5^3 is pronounced five cubed. It is also common to say the a number is the square of another number. For example 16 is the square of 4, because 4^2 equals 16. Similarly 125 might be called the cube of 5 (details of how to calculate these values are given below). In its simplest form index notation is basically a quick way of expressing calculations which would otherwise be quite tiresome to write. The index denotes how many times the base is multiplied by itself, so the actual value of 2^4 can be calculated as shown below. <$table3>2^4 = 2 * 2 * 2 * 2 = 16

Similarly:

<$table3>10^2 = 10 * 10 = 100 <$table3>x^4 = x * x * x * x

<$HB>Square Roots The opposite action to raising a number to a power is called rooting a number. The only type of rooting used in this book is square rooting, which has the opposite effect to raising a number to the power of 2. An example of this is shown below. Note that the BASIC SQR syntax is used to denote square rooting in this book. SQR(64) = 8 (because 8^2 = 8*8 = 64) The MOD Operator The conventional use of the BASIC MOD operator is for finding the remainder after integer division has been performed. However, it is also very useful for keeping numbers within certain ranges without giving adjacent numbers the same value. This technique is used, for example, in Mandelbrot and Julia programs when converting iteration numbers of unknown range to the 16 colour range of the Amiga’s non-interlaced high resolution screen mode. An everyday use of MOD can be found in the display on a 24 hour digital watch. The part of the display showing hours is incremented every hour, on the hour, but when the end of the 23rd hour occurs (at 23:59:59) the hour value is set back to zero, from where the process begins again. The hour value is said to be modulus 24 because it has 24 possible values (0 to 23 inclusive). This modulus solution must be used here in order to prevent parts of the time from becoming too large to handle. The only alternative is to stop the hour value when it gets to 23, which would be equally unworkable. <$fc>Figure B.3: 24-hour clock face showing the action of MOD 24

The 24 hour clock face in Figure B.3 graphically demonstrates a practical example of MOD at work. The hour hand of the clock can rotate an infinite number of times without ever stopping or leaving the range of 0 to 23, but no two adjacent hours ever have the same value.

<$fc>Figure B.4: Four-hour clock face demonstrating MOD 4 The clock face analogy can be used to demonstrate all instances of modulus arithmetic. For example Figure B.4 shows the application of MOD 4. This clock has four possible values, 0 to 3 inclusive. Note that the range always begins at 0, an offset would need to be added to the result if a different base was required. Table B.1 demonstrates the use of the four hour clock by showing the value of a number, x, against (x MOD 4). <$tc>Table B.1: The action of MOD 4

<$table2>x (x MOD 4) <$table2>0 0

<$table2>1 1 <$table2>2 2

<$table2>3 3 <$table2>4 0

<$table2>5 1 <$table2>6 2

<$table2>7 3 <$table2>8 0

<$table2>9 1 <$table2>etc. etc.

The following BASIC program can be used to demonstrate the MOD operator (the program will work in both Amiga and GFA BASIC). On each cycle of the WHILE…​WEND loop a variable, x, is incremented and then printed on the screen next to the corresponding value of (x MOD 4).

<$ob>x=0 <$ob>WHILE 1 <$ob> PRINT x,(x MOD 4) <$ob> INC x <$ob>WEND <$ha>Angles - Degrees and Radians

An angle describes the size of a turn. Angles are usually expressed in degrees, where a full circle consists of 360 degrees, a semi-circle has 180 degrees and so on. However, in Amiga and GFA BASIC the more mathematical radian equivalent is used, where angles are normally expressed in terms of the constant, PI. PI is a common mathematical constant known to BASIC, and is equal to approximately 3.1416.

<$fc>Figure B.5: Common degree to angle conversions Figure B.5 shows the equivalent degree and radian values for some common angles. In the case of more general values the two equations below can be used to convert between radians and degrees. <$ob>degree_value = 180 * (radian_value/PI)

<$ob>radian_value = degree_value * (PI/180) Peculiarities of Right-angled Triangles A right angled triangle is one which contains a right-angled (90 degree) corner. An example of a triangle of this type, with the right-angle and hypotenuse marked, is shown in Figure B.6. The hypotenuse is just the mathematical name for the longest side, which can always be found directly opposite the right-angle. <$fc>Figure B.6: A typical right-angled triangle

<$HB>Pythagorus' Theorem This simple law states that the square of the longest side of a right-angled triangle (the hypotenuse) is equal to the sum of the squares of the two other sides. Taking the triangle in Figure B.6 as an example we can form the following equation: c^2 = a^2 + b^2 This only gives the square of c. To find the actual length both sides of the equation must be square rooted: c = SQR(a^2 + b^2) <$HB>Simple Trigonometry

There are three useful equations which allow the length of a right-angled triangle’s sides to be calculated from one of its angles and vice-versa. These are shown in Figures B.7a to B.7c.

<$fc>Figure B.7: Application of the three common trigonometrical ratios Note that opposite is used to refer to the length of the side opposite the angle, A, and adjacent is the length of the non-hypotenuse side which touches the angle. SIN, COS and TAN are mathematical functions, which are built into BASIC. The relevant syntax for these functions is: <$ob>result = SIN(angle)

<$ob>result = COS(angle) <$ob>result = TAN(angle)

12. Appendix C

<$cn>Listing Conversion As mentioned at the beginning of this book, it is now possible to investigate Chaos theory using a range of popular personal computers and computer languages. The Amiga is only one such computer, and Amiga BASIC is only one of the many available languages. This appendix is provided as a guide to converting the Amiga example programs given in this book to other languages, including BBC BASIC, the language used on BBC micro and Archimedes computers. The conversion notes have been written specifically to include all programming techniques used in this book without digressing into more advanced, unnecessary, material. For each computer language syntax comparison tables are given. These may not on their own allow a program to be converted but they should, in conjunction with the text and example programs, steer you toward a suitable conversion. Because Amiga BASIC specific commands have been avoided wherever possible and because the majority of the programs consist mainly of standard maths and loops, conversion to other BASICs is relatively easy. <$ha>Acorn BBC and Archimedes Computers

Unlike many other computers, all BBC machines have a very good built-in language, BBC BASIC, which is rich in structure and exploits the computers to the full.

There are two distinct classes of BBC computers, the original eight-bit machines (Model As, Model Bs and Masters) and the newer 32-bit reduced instruction set computer (RISC) machines - the Archimedes. Both computers use syntax compatible versions of BBC BASIC, but the original BBC design is naturally much slower than the new one.

Table C.2 is provided as a rough guide to converting the listings given throughout this book to BBC BASIC. Note that BBC BASIC version five, built into Archimedes machines, provides enhanced programming structures, such as IF…​THEN…​[ELSE]…​ENDIF, and extra commands and functions including mouse pointer position functions. All programs written in early versions of BBC BASIC will, however, work on the Archimedes.

Amiga BASIC BBC BASIC

REM REM *CALL name PROCname FOR…​NEXT FOR…​.NEXT IF…​THEN…​[ELSE]…​END IF IF…​THEN…​[ELSE] INPUT INPUT LINE [(a,b)]-(x,y) DRAW x,y MOUSE(0) Mouse not standard on all BBCs *PALETTE VDU 19…​ PRINT PRINT PSET(x,y) PLOT 69,x,y *SCREEN MODE *SOUND frequency,duration SOUND channel,volume,frequency,duration *SUB name STATIC…​END SUB DEFPROCname…​ENDPROC *WHILE…​WEND REPEAT…​UNTIL

*Indicates nearest, not exact, replacement

<$tc>Table C.2: Amiga BASIC conversion table for BBC BASIC. A BBC BASIC version of the Sierpinski triangle program is given in listing C.3 below. This is written so as to be compatible with all versions of BBC BASIC, although it will not run on the BBC model A since it uses screen mode 0. DIM x%(2), y%(2) 'The following assignments should be changed 'if a screen mode other than 0 is to be used ScreenWidth% = 640 ScreenHeight% = 256 MODE 0 x%(0) = 0 x%(1) = ScreenWidth% / 2 x%(2) = ScreenWidth% y%(0) = ScreenHeight% y%(1) = 0 y%(2) = ScreenHeight% vertex% = RND(3)-1 px% = x%(vertex) py% = y%(vertex) REPEAT vertex% = RND(3)-1 px% = px% + (x%(vertex%) - px%) / 2 py% = py% + (y%(vertex%) - py%) / 2 PSET (px%, py%) UNTIL INKEY$<>""

<$fc>Listing C.3: A BBC BASIC program which produces the Sierpinski triangle If you are using an Archimedes rather an 8-bit BBC, you will obvious want to maximise your computer’s graphics capabilities by running the program in the highest resolution screen mode possible. With a standard monitor this is mode 17, which provides a resolution of 1056x256 pixels. To alter the program to run in mode 17, change the MODE command accordingly and edit the lines which initialise the screen size variables to read as follows: ScreenWidth% = 1056 ScreenHeight% = 256 <$ha>Using Other Languages on the Amiga

If you are using an Amiga, but Amiga BASIC is not your preferred programming language you will undoubtedly want to know how to convert the example listings in this book to a form suitable for your chosen language.

It would be impossible to discuss all the differences between the popular Amiga languages in this appendix, so instead only the most important differences for fractal work are mentioned. A table to aid BASIC program conversion is shown below.

Amiga BASIC/Hisoft BASIC GFA BASIC

' or REM ! or REM *CALL name GOSUB name FOR…​NEXT FOR…​.NEXT IF…​THEN…​[ELSE]…​END IF IF…​[THEN]…​[ELSE]…​ENDIF *IF…​THEN…​[ELSE]…​END IF SELECT…​CASE…​DEFAULT…​ENDSELECT INPUT [FORM] INPUT LINE [(a,b)]-(x,y) DRAW [a,b] TO x,y MOUSE(0) MOUSEK *PALETTE SETCOLOR PRINT PRINT PSET(x,y) PLOT x,y *SCREEN OPENS *SOUND frequency,duration SOUND frequency,duration *SUB name STATIC…​END SUB PROCEDURE name…​RETURN *WHILE…​WEND REPEAT…​UNTIL WHILE…​WEND WHILE…​WEND *WINDOW OPENW

*Indicates nearest, not exact, replacement

<$tc>Table C.3: Syntax comparison of Amiga BASIC and GFA BASIC. All Amiga BASIC programs should work without alteration in Hisoft BASIC. <$hb>Hisoft BASIC

Hisoft BASIC is intended to be used for the same type of applications as Amiga and GFA BASIC and therefore contains many of the same commands. Hisoft BASIC syntax is refreshingly similar to that used in Amiga BASIC, making program conversion easy, but a host of powerful extra commands are also provided. Hisoft BASIC is well structured, allowing constants and other special variables to be defined, which is very useful in chaos programming. When converting the programs in this book to Hisoft format you should use such special definitions wherever possible to improve speed and clarity.

<$hb>AMOS BASIC and Blitz BASIC The key point to remember when converting programs to AMOS is that all non-integer variables must be suffixed by a hash character (#). Because AMOS is geared toward game creation, where floating point numbers are rare, it treats all variables as integers unless instructed otherwise in this way. In most other respects AMOS has the same syntax as Amiga BASIC except for a few notable exceptions, the most prominent of which is the use of line numbers and GOSUBs rather than named procedures. The many extra non-standard commands provided by AMOS make some tasks, such as graphics manipulation, much simpler. This discussion also holds for Blitz BASIC, as it is intended for the same programming tasks as AMOS. Note, however, that Blitz BASIC is compiled and is the fastest BASIC implementation for the Amiga. <$hb>GFA BASIC

The GFA BASIC interpreter is the ideal language for chaos programming as it is well structured and very fast (faster even than compiled Hisoft BASIC). The fact that a compiler is not necessary to achieve this speed make it very easy to make small changes to programs and see their effect without delay. The easy access to the Amiga’s user interface, discussed in chapter five, make GFA BASIC ideal for development of easy to use fractal 'explorers'. If you have a special interest is writing chaos programs in GFA BASIC then I recommend the Atari edition of this book in which all the examples are written in the ST version of GFA BASIC, which is syntax compatible with the Amiga version.

<$hb> C Although specific applications for C have been cited throughout this book, no example code has been given. Originally I had intended to present all the example programs in both C and BASIC, but it was decided that this would only make the principles of chaos theory harder to grasp. Also, it was assumed that most C programmers would be able to read BASIC and should therefore be able to convert the programs with relative ease. C has the advantage of being generally quicker than BASIC, especially interpreted BASIC, and is also more easily ported between different computers. 13. Appendix D <$cn>Use of Amiga Peripherals

There is now a wide variety of expansion options and peripherals available for the Amiga range of computers, but to enable all Amiga owners to use this book the major programming examples have been written around the basic Amiga 500. This system has 512k of RAM and an internal floppy drive with a minimum formatted capacity of about 800k. However, fractal generation will benefit from some peripherals, this chapter introduces the techniques necessary to adapt existing programs to make the most of such peripherals.

<$ha>Monitors The listings contained in this book can easily be adapted to run in any Amiga screen mode, therefore making the most of any monitor. Techniques for using the various modes are given in appendix A. <$ha>Printers

If the fractals you create are intended to reach a wider audience than fellow Amiga owners, a printer will become a necessity. Figures D.1a to D.1c show the difference between the output of the three main monochrome printing methods.

<$fc>Figure D.1a: Low-resolution non-interlaced screen dump to 9 or 24 pin printer (320x200) <$fc>Figure D.1b: High-resolution interlaced screen dump to 9 or 24 pin printer (640x400)

<$fc>Figure D.1c: Very high resolution image produced on a 24 pin printer by Listing D.1 (1416x1416 pixels) <$hb>Nine Pin

The most common Amiga printers are of the monochrome nine pin dot-matrix type, and with the help of a good screen dump routine they can quickly produce an accurate copy of the screen display, with the same resolution as the screen. Colour may also be represented on such printers using shades of grey, although the effect is not very good.

<$hb>Twenty-four Pin Twenty-four pin dot-matrix printers are now becoming far more numerous. They can produce output of around 180 dpi (dots per inch) which is far superior to nine pin printers. However, when dumping screen images the resolution of the dump is of course limited by the resolution of the screen, and since nine pin printers can easily manage 640x512 pixel plots a 24 pin printer would appear totally unnecessary. However, higher resolution plots can be produced by by-passing the screen altogether and writing directly to the printer instead. Unfortunately this technique can only really be used to draw fractals which are created by placing pixels in rows, like the Feigenbaum diagram and Mandelbrot/Julia sets. When plotting such fractals near laser quality output of 1416x1896 pixels is possible on an A4 sheet of paper. The program given in Listing D.1 demonstrates the 24 pin printing method for the Mandelbrot set. After requesting information on the section of the set that you want to plot, the program will print that section with a resolution of 1416x1416 pixels. It sounds like a quick process, but remember that there are 1414x1416 = 2,005,056 pixels (compared to 64,000 or at most 256,000 when plotted to the screen) and each pixel must be tested using one of the time consuming methods discussed in Chapter four. Also, the iteration ceiling must be fairly high in order to differentiate between pixels which belong to the actual Mandelbrot set and those that do not. Many calculations are necessary and plots taking over 24 hours to produce are not uncommon when running the program in interpreted BASIC. At least the Amiga’s multitasking allows you to get on with something else while Mandelbrot calculation is underway. This program has been written specifically for the Epson LQ-500, but it should work with any compatible 24 pin printer. However, it must be stressed that I have not tested the program on any other printers. <$ob>' 1416 x 1416 Mandelbrot Direct To Printer (V1.30A)<$ob>' (c)1991 Conrad Bessant <$ob>'<$ob>DEFINT bit,byte,yl,x,t,value<$ob>DEFDBL a,b,c,i,j,imin,jmin<$ob>'<$ob>PRINT "Please enter required parameters:"<$ob>PRINT<$ob>INPUT "Minimum value of a (amin):",imin<$ob>INPUT "Maximum value of a (amax):",imax<$ob>INPUT "Minimum value of b (bmin):",jmin<$ob>PRINT "Please wait…​"<$ob>s=imax-imin<$ob>n=1416/s<$ob>'<$ob>OPEN "PAR:" FOR OUTPUT AS #2 'Open connection to parallel printer <$ob>PRINT #2,CHR$(27);CHR$(51);CHR$(24); 'Set line feed size<$ob>FOR yl=0 TO 1415 STEP 24<$ob> PRINT #2,<$ob> PRINT #2,CHR$(27);CHR$(42);CHR$(39);CHR$(136);CHR$(5); 'Tell printer that 1416 <$ob> columns of 180 DPI graphics are coming<$ob> FOR x=0 TO 1415<$ob> FOR byte=1 TO 3<$ob> value=0<$ob> FOR bit=0 TO 7<$ob> y=1416-(yl+byte*8-bit-1)<$ob> i=x/n+imin<$ob> j=y/n+jmin<$ob> a=0<$ob> b=0<$ob> t=0<$ob> WHILE a*a+b*b<4 AND t<120<$ob> t=t+1<$ob> c=a<$ob> a=i+(c*c-b*b)<$ob> b=j+2*c*b<$ob> WEND<$ob> IF t MOD 2=0 THEN<$ob> value=value+(2^bit)<$ob> END IF <$ob> NEXT bit<$ob> PRINT #2,CHR$(value);<$ob> NEXT byte<$ob> NEXT x<$ob>NEXT yl<$ob>PRINT #2,CHR$(10)<$ob>PRINT #2,CHR$(12)<$ob>CLOSE #2 'Close printer connection <$fc>Listing D.1: Amiga BASIC program to output the Mandelbrot set to an Epson LQ-500 24 pin printer It would be tedious to explain this program in detail, but the key difference between this and the similar screen oriented ones given in Chapter 4 is the order in which the pixels are plotted. Because the screen is effectively a random access device it is possible to plot pixels in the most convenient way, usually in vertical columns drawn from left to right. Conversely, the printer expects data to be sent in a precise format and as a result the order in which the pixels are plotted has to be tailored to suit the printer. The way in which the image in built up on the LQ-500 is shown in Figure D.2. In addition, the iteration ceiling of 120 is much larger than those used earlier in the book. This is because higher resolution plots need a proportionately higher level of calculation detail. Because the printer is fairly unintelligent we need to send the pixel information in low level bits, with control codes governing the way the data is interpreted. However, all information sent to the printer using the LPRINT command is automatically piped through an AmigaDOS converter, which takes out printer-specific control codes. Since the Mandelbrot program relies entirely on control codes this is a serious problem. However, by opening a direct connection to the printer with the <F102>OPEN "PAR:" FOR OUTPUT AS #2<F> command and subsequently using PRINT #2 to write the data, the image arrives at the printer in it’s raw, unconverted, form. The time consuming calculations required to plot a Mandelbrot set at this resolution make conversion of the Amiga BASIC program to the faster GFA BASIC an obvious enhancement. Because the program is more printer-specific than language-specific such conversion is not difficult. In fact the only major change required is the replacement of the line used to open the direct connection to the printer, which should be altered to read as shown below: <$ob>OPEN "O",#2,"PAR:"

As mentioned in the introductory chapter, GFA BASIC does not allow the use of the apostrophe character (') for end of line comments, this should be bourne in mind when converting the Mandelbrot program. Remember also that END IF should be written as one word, ENDIF, in GFA BASIC.

<$fc>Figure D.2: Structure of a 180dpi picture on the Epson LQ-500 <$hb>Colour

Obviously a colour printer would be desirable for producing some fractals, and the resolution of such printers can normally cope with Amiga screen dumps up to 640x400 pixels in size. However, colour printer drivers for the Amiga are scarce and you may be forced to write your own.

<$hb>Plotters Due to the constraints imposed by their design, plotters can realistically be used only for drawing fractals consisting mainly of lines. The fact that plotters physically draw with real pens means that they can produce very high resolution Lorenz attractors, fractal plants and landscapes if accessed directly (without plotting to the screen). Other fractals, such as the Mandelbrot set, are very difficult to produce on a plotter because the drawing of single dots is slow and screen dump programs for plotters are rare. <$ha>Extra RAM

Fractal generation is a relatively memory efficient process because the programs are short and most of the information is stored in the screen memory area reserved by AmigaDOS at all times. Extra memory is only really useful for animating fractals. On a two megabyte machine, for example, it would be possible to create 36 three dimensional pictures (frames) of the Lorenz attractor, each one rotated 10 degrees in relation to the previous one. After being drawn on the screen (not necessarily with any great speed) each frame could be stored in RAM. Once all the frames had been stored they could then be copied onto the screen in quick succession to create the illusion of a rotating three dimensional Lorenz attractor. Having tried this technique on the Acorn Archimedes I can testify that the end result is very effective and dispels any doubt about the Lorenz line intersects itself. GFA BASIC is particularly suited to this task, as it’s matrix manipulation facilities allow effortless manipulation of the three dimensional attractor.

From Chapter 6 you may recall that extra memory is also useful when drawing very detailed fractal plants.

<\$cn>Bibliography

James Gleick, Chaos: Making a New Science (Viking Press, 1987)

This is the most popular introduction to chaos, because of its agreeable style and accessibility. It describes a wide variety of fractals and chaotic processes, backed up by a host of real world chaos examples. The addition of historical information and interviews with people such as Mitchell Feigenbaum help create a useful insight into the modern scientific community, which is very enlightening if you have even considered becoming part of it. Unfortunately, few mathematical details are given, making it difficult to recreate the experiments discussed.

Karl-Heinz Becker and Michael Dorfler, Dynamical Systems and Fractals: Computer experiments in Pascal (Cambridge, 1989)

In stark contrast to Chaos this book is composed mainly of technical information such as equations and programs, with little historical background or real world examples to relate to. The programs are written in Pascal for the Apple Macintosh, but are explained well enough to make conversion to other computers relatively easy (as long as you can read Pascal).

It contains several unexpected topics, such as a chapter entitled Journey to the Land of Infinite Structures in which the authors take the reader on a flying trip around the Mandelbrot set. However, the book does not suffer from these interludes or from the fact that it has been translated from German and the resulting mix of hard information and some lighter moments make this the ideal companion for amateur scientists who wish the carry out their own chaos experiments.

Heinz-Otto Peitgen and Peter H. Richter, The Beauty of Fractals, (Springer, 1985)

The wealth of fractal pictures (many in colour) adorning The Beauty of Fractals has apparently made it the first ever coffee table mathematics book. However, there is some detailed theory buried among the pictures and it has to be said that such high quality fractals would be hard to find anywhere else (although some of these are reproduced in Chaos). This is really the sort of book to get out of the library, as it doesn’t contain much in the way of reference material and is quite expensive. However, if your coffee table is looking slightly bare then it may be worth getting your own copy.

Benoit Mandelbrot, The Fractal Geometry of Nature (Freeman, 1982)

Mandelbrot washes down the sometimes difficult, and now famous, theory with his personal recollections, comprising details of the equipment he was using and what he was trying to achieve when he formulated certain theories. Some of his original output is also included, most notably the first ever printout of the Mandelbrot set. Being one of the first books on chaos, people are inclined to dismiss this as being nothing more than a historic artifact. However, much of the work is still very relevant. In fact a great deal of the modern chaos literature is based on theories derived from this and Mandelbrot’s other book, Fractals: Form, Chance, and Dimension (Freeman, 1977).

Predrag Cvitanovic, Universality in Chaos (Adam Hilger, 1989)

This is a collection of notable scientific papers on chaos published over the last 20 years. Such papers provide the information and inspiration for almost all chaos books, and are where the now famous discoveries were first published. Pieces by various authors are featured including Edward Lorenz, Mitchell Feigenbaum and Robert May. Although the papers are only fully comprehensible to science graduates useful information for chaos programming on the Amiga can be gleaned.

1. Sulivan, The Big Red Book of C (Sigma Press, 1983)

In order to utilise the full power of the Amiga, programs should really be developed using a C compiler, the structure and speed of the language make it ideal for chaos work and if Chaos, Fractals and the Amiga has urged you to make the step up to C you really need a suitable tutorial text from which to learn it. Although I know of only a few such books written specifically for the Amiga, there are many more general works, of which The Big Red Book represents a good value example. All the main aspects of the language are described but supplementary information, specific to the Amiga, would be required before embarking on any major programming project.

If you wish to continue chaos programming using GFA or Amiga BASIC a number of books are available to help you master the general programming techniques. Amiga magazines, particularly the more serious ones, usually carry adverts and reviews of new titles in this field.

It is important to note that this is only a very small selection of relevant literature and that much more material is available in the form of scientific papers, articles and other books. Such material is often highly mathematical, so I would not advise moving up to these books right away. However, most of the chaos books above have their own unique bibliographies which will give more specialised titles, hopefully leading you gently in the direction in which you want to go.