Fast forward 30^{+} years and I now own two computers that were incomparably cheaper to purchase than the Apple III. Each of them is equipped with a fast "Pentium" processor and several gigabytes of RAM and churns out what I used to get after days of runtime within minutes or hours instead. The modern computers are also capable of showing whopping 2^{24} = 16,777,216 so-called "24-bit" colors instead of only 2^{ 1} or two "1-bit" colors and coloring was no longer simple.

According to the American Academy of Ophthalmology, good human eyes can distinguish "only" up to 10 million colors out of the 16 million^{+} colors, so I expect that the 24-bit colors would stay with us for a long time without a substantial change regardless of the progress in computer technologies. ..... J. Sekno, 2015.

and a canvas comprising 3,200 × 1,850 pixels. To demonstrate our coloring schemes, we briefly review a few basic ideas used in the main website and plot a few Julia fractals generated by (1.1) on the canvas in a primitive way.

So, we first choose a value of the parameter p, say,

(1.2) p = (-0.77146, -0.10119),

which is a complex number belonging to the largest circular atom of period2 in the Mandelbrot set near its boundary. As explained in the main website, there is no mystery in choosing a number that looks like that. Throughout this website, p means the parameter in (1.2) and "Julia fractals" means Julia fractals of p generated by (1.1).

We also view the complex plane as the set of all initial values z_{0} of the orbits z_{n} generated by (1.1) with (1.2). By "orbits" or "orbits of z_{0}," we will always mean these orbits. To every complex number z_{0} in the complex plane, there corresponds the orbit of z_{0} that may diverge to ∞, converge to a cycle or behave unpredictably.

In addition, we superimpose the canvas on the rectangular region in the complex plane bounded by xmin = -1.65, xmax = 1.65, ymin = -0.9539 and ymax = 0.9539, where 0.9539 ≈ 1.65 × 1850/3200. As explained in our main website, we identify each pixel (i, j) (which is a tiny rectangle) in the canvas with a unique representative complex number z_{0} belonging to the pixel. Thus, we write a pixel not only as (i, j) using the pixel coordinates of the canvas but also as z_{0} which is the complex number identified with (i, j).

by (1.1) and the triangle inequality, hence, the computer paints the pixel white and stops iterating (1.1) to avoid wasting the computing time. Observe that m = 1 is the least iteration index such that m < M and |z_{m}| > 2, which determines that the orbit of z_{0} diverges to ∞.

This gives us another way of carrying out the divergence scheme. For every pixel (i, j) in the canvas, er complex number z_{0}, iterate (1.1) at most M times and compute the nonnegative integer

(1.3) d(i, j) = the least iteration index m < M satisfying |z_{m}| > θ,

if such m exists, and d(i, j) = M, otherwise. (1.3) defines the array of the least iteration indices over the canvas, which we call the iteration array for short. Note that Figure 2 is now given by the instruction that is far clearer and simpler than the awkward if-statement:

If d(i, j) is even, color the pixel (i, j) black, and if not, color the pixel (i, j) white.

If we add one more line

If d(i, j) = M, color the pixel (i, j) red

at the end, we get the first image of Figure 3 as output. The red figure and its boundary represent (computer approximations of) the filled-in Julia set and the Julia set of p = (-0.77146,-0.10119), respectively.

If d(i, j) = M so the pixel (i, j) belongs to the filled-in Julia set, let

(1.4) c(i, j) = the least iteration index m < M satisfying |z_{m+k} - z_{m}| < ε,

if such m exists, and c(i, j) = M, otherwise. If d(i, j) < M so the orbit of z_{0} diverges to ∞ at the pixel (i, j) then set c(i, j) = M.

The nonnegative integers c(i, j) define the iteration array for the convergence scheme over the canvas. Along with the iteration array d(i, j) for the divergence scheme, it contains all the information needed for plotting the Julia fractals with the 24-bit colors. Note that the Julia set (which is the boundary of the filled-in Julia set) is approximated by the set of all pixels (i, j) on the canvas such that c(i, j) = M and d(i, j) = M, which is a region of interest in chaos theory.

The goal of this website is to plot the specific image, which is shown below, using the iteration arrays. Along the way, we will cover most of the important coloring techniques and at the same time see why we use such a large number M as the maximum number of iterations. Our canvas size is also large as it gives us, among other things, an option of getting a high resolution printout of the image.

Because of the large numbers, it took our Pentium machine 13 hours to compute the iteration arrays, even though we use the fact that the Julia fractals by the Mandelbrot equation are symmetric about the origin of the complex plane and halve the computing time. So, the first thing we do after the computation is done is to save the iteration arrays in the hard disk. The computer can convert them to a visual image such as the one shown below in a matter of seconds.

This is because the color cube is structured like Rubik's cube comprising 256 × 256 × 256 small subcubes (instead of 3 × 3 × 3 subcubes) and each color is a small subcube of the color cube rather than a point in the 3D space. It is similar to a canvas comprising pixels that are rectangles instead of points. Thus, the color cube is a solid cube rather than a discrete set of points with integer components.

Even though colors are not points, there are some analogies between the color cube and the Cartesian 3D space with x-, y-, z-axes. For example, each color in the color cube is uniquely represented as an ordered triple (R, G, B) of integers between 0 and 255, called bytes and it makes the total number of colors in the color cube 256^{ 3} = 2^{24} = 16,777,216. Also, each color can be written uniquely as the sum of the aforementioned shades of the primary colors. For example, we have

Even though a color is not a point, we'll still treat it like a "point" when we talk about it. For example, black = (0, 0, 0) is the "origin" of the color coordinate system and white = (255, 255, 255) is the "vertex" of the color cube opposite from the origin (black). In general, a color is darker if it is closer to black and lighter if it is closer to white. In fact, we can darken a color by a (componentwise) scalar multiplication (followed by rounding the components to the nearest bytes) of the color by a scalar strictly between 0 and 1. For example, we have:

Recall that M = 40000 is the maximum number of iterations of the Mandelbrot equation (1.1) used to compute the iteration arrays c(i, j) and d(i, j) over the canvas comprising 3,200 × 1,850 pixels (i, j). For each pixel (i, j) in the canvas, m = c(i, j) is the least iteration index of an orbit z_{n} generated by the convergence scheme and m = d(i, j) by the divergence scheme. Hence, both are integers belonging to the closed interval [0, M] of real numbers, i.e., both are integer-valued functions mapping the canvas as the set of pixels into the interval [0, M].

Our goal is to paint a colorful Julia fractal on the canvas, whose framework has been done as in Figure 2, using the iteration arrays and the 24-bit colors. Compared to computing the iteration arrays that took our Pentium machine 13 hours, the painting can be done in a matter of seconds, but like painting a newly built house, it can improve or ruin the entire effort. So we proceed carefully.

The diagram shown below illustrates the overview of our painting process and points out that painting on the canvas is the composition of three functions m = c(i, j) or m = d(i, j) followed by a function t = f(m) followed by a line segment function. Since we have already found two of the component functions in § 1 and § 2, all we need at this stage is to find a function t = f(m) that maps an arbitrary subinterval [m_{1}, m_{2}] of [0, M] into the unit interval [0, 1]. As we shall see, once it is done, we can paint what is hidden in the iteration arrays on the whole canvas in full color.

is used instead of (3.2), the color C changes with m from C_{2} to C_{1} ; see the second graph of Figure 6 with the special case q = 1. If q = 2 in (3.2) and (3.3), the corresponding graphs have parabolic curves instead of the slanted lines in Figure 6.

FOR i = 0 TO 3199
FOR j = 0 TO 1849
m = c(i, j) power = 1 IF m ≤ 26500 THEN ColorUp(i, j, false, 0, 26500, power, altgreen, black)
ELSE ColorPixel(i, j, black) ENDIF-ELSE
ENDFORj
ENDFORi

Note that the complement (background) of the filled-in Julia set is painted by the earlier program in Example 1. If we change power = 1 in the fourth line by power = 4, we get the second image with bigger green "eyes."

The new program uses the number 26,500 in place of 80 in the program in Example 1. It shows that the convergence of the orbits is much slower than the divergence to ∞ and explains why we need a large maximum number of iterations like M = 40000. We can actually use two different values for M, one for the divergence scheme and the other for the convergence scheme. But again, we keep the programs as simple as possible on this website so the important ideas don't get lost.

Periodic Wave Functions: Recall that the ColorUp and ColorDn subroutines are based on the first and second graphs of Figure 6, respectively; hence, a repetition of ColorUp followed by ColorDn provides us with the effect of coloring pixels through a "periodic wave function." The following subroutine is based on the idea:

SUBROUTINE ColorWave(i, j, dvg,m_{1}, m_{2}, freq, power, C_{1}, C_{2})
// All variables except freq are the same as in the subroutine ColorUp
// freq is an integer and designates the frequency of the wave cycles
hp = round((m_{2} - m_{1})/(2*freq)) // hp stands for half the period of a wave
IF dvg THEN m = d(i, j) ELSE m = c(i, j) ENDIF-ELSE
m_{0} = m_{1} counter = 1 WHILE m_{0} < m_{2} IF m_{0} ≤ m < (m_{0} + hp) THEN
CASE (counter MOD 2) = 1:
ColorUp(i, j, dvg,m_{0}, m_{0} + hp, power, C_{1}, C_{2}) ENDCASE
CASE (counter MOD 2) = 0:
ColorDn(i, j, dvg,m_{0}, m_{0} + hp, power, C_{1}, C_{2}) ENDCASE
ENDIF
m_{0} = m_{0} + hp counter = counter + 1 ENDWHILE
ENDSUBROUTINE

If we combine this program with the program for the second image of Example 1, which paints the complement (background) of the filled-in Julia set, we get the first image of Figure 10 as the output image. Note that the periodic wave function with frequency 36 and power 1.8 starts at m_{1} = 3410 and not m_{1} = 26500. We find all these numbers by trial-and-error computer experiments.

If we increase the frequency from 36 to 40 in ColorWave, we get the second image of Figure 10.

As we have seen, the ColorUp and ColorDn subroutines are based on the very simple functions in Figure 6 but they generate powerful subroutines for fractal plotting. Here's an additional subroutine:

for all i and j provides us with the color change through the "triangular wave function" given by the first image shown below:

A large part of the first image in Figure 1 shown at the outset uses a similar color change with a greater frequency of the cycles. Furthermore, the periodic functions don't have to be continuous. For example, we can easily modify the last subroutine to do the coloring shown below:

This is a type of coloring used to paint a part of a fractal like: