Calc - Java Calculator for cell-phones and MIDP devices

Programming examples

In the program listings, the most basic commands are written as just the command name, such as "ENTER", "*", and numbers such as "0.5". Commands that are found deeper down in the menu hierarchy will be written prefixed with some parts of the menu choices needed to get to the command. This is done to remind the user how to find the command, it is not meant as a complete description of how to find the command. For example, if the program listing contains "stack/RCL st# 1", what this actually means is "main/special/stack/more/RCL st#/<0-3>/<1>".

Programming hint: If a program needs parameters pushed onto the stack before running, it is useful to push some dummy arguments on the stack before starting to record the program, to be able to follow the programming with actual data.

Index of programming examples

  1. Drawing a function
  2. Quadratic equation
  3. Solving an equation for different variables
  4. DDD.MMSS conversions
  5. Biorhythms
  6. Advanced curve fitting example
  7. Spline drawing
  8. Plotting a circle
  9. Plotting a sphere
  10. Mandelbrot fractals
  11. Linear Programming
  12. Plotting curves

Drawing a function

Many people have requested a more thorough explanation of how to do such a simple thing as drawing a function. To recap the long and perhaps tedious documentation, a funtion is a program that uses the lowest element on the stack (x) as input, calculates the function value based on this input, and leaves the resulting function value in the lowest element on the stack as output.

We will proceed to program and draw two different programs of different complexity. The first program represents one of the simplest functions possible:

    f(x) = x

Let's see how we can make a program to calculate this function. When the program is called by the graph drawing operation, the input value x will be present as the lowest element on the stack. Since the input value of this particular function equals the output value, we don't have to do anything to it to it to calculate the function value. Also, since the output value of the function should be located in the lowest stack element when the program finishes, the program does not have to do anything either. Therefore, the operations you need to execute to program this function are as follows:
prog/new "f1"
prog/finish

To draw the function, we first need the boundaries of the draw area pushed on the stack. We will use the area -10<x<10 and -10<y<10. After entering the limits, we may execute the draw operation:
-10
10
-10
10
prog/draw/y=f(x)/f1

This looks rather boring. However, the function is not completely useless. If we use other drawing modes (with the same limits as above), we get some nice results:
prog/draw/r=f(θ)/f1
prog/draw/z=f(z)/f1

Now for a more complex example. The function we want to draw now is the following:

    f(x) = 2x - 3x - 1

One of the first problems we have to deal with is that we need 'x' twice to calculate this function, and there is only one 'x' on the stack when the program starts. We could use "mem/STO" to store the x in a memory location, and "mem/RCL" to get it back whenever we need it, but a more common method is just to dupliacte the number on the stack, using ENTER. One possible solution to program the function follows. Before programming, clear the stack and enter a dummy input value on the stack, e.g. '1', to be able to follow the calculation:
stack/clear              ; clear the stack
1                        ; enter dummy input value (be wild: try '2' instead)
                         ; (text after the ';' is just comments)
prog/new "f2"            ; now we start programming our new function
ENTER                    ; copy x. Now there are two x'es on the stack
math/simple/x           ; calculate x, using one of the x'es
2
*                        ; now we have 2x
stack/x<->y              ; bring down the other x
3
*                        ; now we have 2x and 3x on the stack
-                        ; subtract the two, now we have 2x-3x
1
-                        ; subtract 1, now the stack contains f(x)
prog/finish              ; program finished

Let's see how the program works. Clear the stack, and enter the number 3. Then press "prog/run/f2". The '3' disappears and is replaced by '8'. We have now calculated f(3)=8.

To draw the function, follow the steps above, but select program f2 instead of f1 for drawing.

Quadratic equation

This program solves a second-order polynomial of the type

    ax + bx + c = 0

The method used is accurate also if b >> 4ac, as described here. Before running the program, push the coefficients a, b, and c to the stack, in that order. On exit, the stack contains the original coefficients plus the two roots (real or complex).

Program listing:
prog/new "Quad"          ; start recording a new program
stack/RCL st# 1          ; b
math/simple/x
stack/RCL st# 3          ; a
stack/RCL st# 2          ; c
*
4
*
-
math/simple/sqrt         ; sqrt(b-4ac)
stack/RCL st# 2          ; b (as the stack moves the index changes)
prog/util/sgn
*                        ; sgn(b)*sqrt(b-4ac)
stack/RCL st# 2          ; b
+
-2
/                        ; q = -(b+sgn(b)*sqrt(b-4ac))/2
stack/RCL st# 1          ; c
stack/RCL st# 1          ; q
/                        ; first root, c/q
stack/x<->y              ; q again
stack/RCL st# 4          ; c
/                        ; second root, q/a
prog/finish              ; program finished

Solving an equation for different variables

Assume you have an equation having different variables, for example the equation describing an accelerated motion from a starting height h0, with accelleration a, and initial velocity v0. You want to know the height h at time t. Then the equation looks like

    h = h0 + v0t  + at

To solve this equation for different variables, we make use of some registers, one for each variable, the memory monitor mode and the equation solver, and even function plotting.

Since the built in solve command solves equations for 0, we transform the above equation into

    0 = h0 + v0t + at - h

We use the following memory registers for the different variables, just in order of their appearance in the equation:

M1=h0, M2=v0, M3=t, M4=a, M5=h

The register M0 will have a special role here: It will the determine the variable we want to solve the equation for, e.g. set M0=3 for solving for M3=t.

For following this example, switch on the memory monitor mode ("mode/monitor/mem/<4-7>/<6>") and enter some values:

M0=3, M1=100, M2=10, M4=-9.80665, M5=0

(The 9.80665 can be entered using "special/conv/const/astro/gn".)

Now we enter the equation as a program:
prog/new "Speed"         ; start a new program
mem/RCL 0                ; find out for which variable to solve
prog/mem/STO[x]          ; store value on stack in this variable
clear                    ; remove var number+value
clear
mem/RCL 1                ; start equation: recall h_0
mem/RCL 2                ; recall v_0
mem/RCL 3                ; recall t
*                        ; v_0 * t
+                        ; h_0 + v_0 * t
mem/RCL 4                ; recall a
mem/RCL 3                ; recall t
math/simple/x           ; t2
*                        ; a * t2
2
/                        ; a * t2 / 2
+                        ; h_0 + v_0 * t + a * t2 / 2
mem/RCL 5                ; h
-                        ; h_0 + v_0 * t + a * t2 / 2 - h
prog/finish

Now we can use the solver to find out which value for t matches the given values for the other variables. To start the solver, we need to specify a lower bound and an upper bound where to search for the zero of the equation, eg:
0                        ; lower bound
100                      ; upper bound
prog/more/solve/Speed    ; call solver for our equation

The result shows that after about 5.65 seconds we will experience a pretty hard landing on the ground. :)

If we wanted to know what initial speed is needed so that a ball starting at a height of 100m is at 50m after 1 second, we would enter:

M0=2 (solve for M2=v0), M1=100; M3=1, M4=-9.80665, M5=50
0                        ; lower bound
100                      ; upper bound
prog/more/solve/Speed    ; call solver

This gives us "nan" as result, since our guess of the bounds was inappropiate. To find out which guess we could use, plot the error of the equation:
0
100
0
100
prog/draw/y=f(x)/Speed

The output looks like there is a zero for negative values of v0. Hence try
-100
0
prog/more/solve/Speed

This gives the result: approx -45.1 So we have to hit the ball quite hard ... :)

One suggestion: The hardest part is to remember which variable is in which register. The memory mode displaying the entered/solved values can be quite helpful for guessing which is which, but probabely you will need a pice of paper here.

Submitted by Goetz Schwandtner

DDD.MMSS conversions

The built-in conversions to and from the DH.MS date-time representation is quite impractical if you are calculating map or stellar coordinates and want to convert between decimal degrees and degrees, minutes and seconds. However, a few small programs using the already built-in conversions can help you out:

Program listing:
prog/new "->DMS"         ; start recording first program
ENTER                    ; make copy of input argument
math/misc/int/trunc      ; integer part
stack/x<>y
math/misc/int/frac       ; fractional part
conv/time/->DH.MS
+
prog/finish              ; first program finished


prog/new "->D"           ; start recording second program
ENTER                    ; make copy of input argument
math/misc/int/trunc      ; integer part
stack/x<>y
math/misc/int/frac       ; fractional part
conv/time/->H
+
prog/finish              ; second program finished


prog/new "DMS+"          ; start recording third program
ENTER                    ; make copy of first argument
math/misc/int/frac       ; fractional part
stack/x<>y
math/misc/int/trunc      ; integer part
stack/x<>st# 2           ; get second argument
ENTER                    ; make copy of second argument
math/misc/int/frac       ; fractional part
stack/x<>y
math/misc/int/trunc      ; integer part
stack/x<>st# 2           ; get first fractional part
conv/time/DH.MS+         ; sum everything...
+
+
prog/finish              ; second program finished

Note: If you are able to enter the sign "" in the program name, it will magically turn into the sign "→" in the menu. In that case you should call your programs "D.MS", "D", and "D.MS+" for clarity.

Biorhythms

The theory of biorhythms is speculative, but slightly entertaining. To program it, we use the fact that real numbers are drawn as a pink line, imaginary parts of complex numbers are drawn as a yellow line, and points where the imaginary and real parts of a complex number coincide, are drawn with a white color. In addition we use the fact that if a program result alternates between values, many graphs can be drawn together. This program is only intended to be "draw"n, if you "run" the program, you will get alternating results.

Program listing:
prog/new "BioRt"         ; start recording a new program
conv/time/now            ; get current time
mem/RCL 0                ; get birthday
+/-
conv/time/DH.MS+         ; subtract dates
conv/time/->H            ; hours since birthday
24
/                        ; days since birthday
+                        ; add x offset (on stack at start of program)
trig/pi
ENTER
+
*                        ; multiply by 2*pi to convert to radians
ENTER
ENTER                    ; 3 copies now on stack
23
/
trig/sin                 ; physical cycle is 23 days
stack/x<->y
28
/
trig/sin                 ; emotional cycle is 28 days
trig/coord/r->cplx       ; make complex number, re=emotional, im=physical
stack/x<->y
33
/
trig/sin                 ; intellectual cycle is 33 days
ENTER
trig/coord/r->cplx       ; make complex number with re=im
mem/RCL 1                ; memory location 1 alternates between 0 and 1
+/-
1
+
mem/STO 1                ; alternating number calculated and stored
prog/util/select         ; select y or z depending on alternating number
prog/finish              ; program finished

Preparations: To try the program, store your birthday in memory location 0, store 0 in memory location 1, make sure calculator is in RAD mode, and push suitable graph limits on the stack, for instance -14, 14, -1.2, 1.2, before finally executing "prog/draw/y=f(x)/BioRt". This will draw the biorhythms for 14 days before and after present. The pink graph shows your emotional cycle, the yellow graph shows your physical cycle, and the white graph shows your intellectual cycle, with zero on the x axis representing current time.

Advanced curve fitting example

Suppose we have a set of data observations and we want to find curve that fits the points:
   x   y
 -------
   0   0
   5   6
  10  10
  20  17
  30  21
  50  27

Looking at the points, the curve that seems to start at the origin, grow linearly at first, and gradually lose its steepness towards the right. One function that behaves this way is:

    y = a ln(bx + 1)

Now we want to find out if such a curve fits the observed data well. The statistics module can fit a curve to the function y = a ln(x) + b. If we add a parameter c to the x values before summing the statistics, we can fit a curve to the function y = a ln(x + c) + b. This curve corresponds to the function above iff. ln(c) = -b/a. Consequently, programming the entering of the above data points with modified x values, and calculating the resulting ln(c) + b/a, we can find our curve using the "solve" operation.

Program listing:
prog/new "fit"           ; start recording a new program
mem/STO 0                ; save input value, c
stat/clear
0 ENTER
0
mem/RCL 0
+
stat/SUM+                ; x1+c,y1
6 ENTER
5
mem/RCL 0
+
stat/SUM+                ; x2+c,y2
10 ENTER
10
mem/RCL 0
+
stat/SUM+                ; x3+c,y3
17 ENTER
20
mem/RCL 0
+
stat/SUM+                ; x4+c,y4
21 ENTER
30
mem/RCL 0
+
stat/SUM+                ; x5+c,y5
27 ENTER
50
mem/RCL 0
+
stat/SUM+                ; x6+c,y6
mem/RCL 0
math/pow/ln              ; ln(c)
stat/result/alnx+b/a,b   ; calculate a,b
/
+                        ; output value, ln(c)+b/a
prog/finish              ; program finished

Running the program with a few test values shows that c=1 and c=20 is located on either side of the solution, and entering those limits before executing "prog/solve/fit", yields the result c=10.76295. Looking at the curve with "stat/result/alnx+b/draw" shows a pretty good fit (although translated 10.76 along the x axis), and the correlation coefficient "stat/result/alnx+b/r" returns 0.9996, which indicates that the data is modeled well by the curve. The coefficient "b" of the original function can be calculated from the coefficients of our current function as eb/a, and the coefficient "a" is the same.

By letting the input value of a program modify a set of statistics, we have used the solve function to perform a curve fitting out of the ordinary. The resulting function is as follows:

    y = 15.724 ln(0.092911x + 1)

Spline drawing

Given a number of 2D control points, a smooth parametric curve consisting of segments of third degree polynomials x=f(t) and y=g(t) can be fitted to the control points. This is called a spline. The calculations are simplified by the use of matrices.

Program listing:
prog/new "Splin"         ; start recording a new program
mem/RCL 0                ; matrix of control points
stack/x<->y              ; t
math/matrix/size
clear                    ; rows = number of control points
3
-                        ; rows-3 = number of curve segments
*                        ; t*(rows-3) = current curve segment
enter
math/misc/int/frac       ; t within this curve segment
mem/STO 2
clear
math/matrix/split        ; extracting control points for this curve segment...
stack/x<->y
clear
4
math/matrix/split
clear                    ; G = [ G0  G1  G2  G3 ]^T  (control points)
mem/RCL 2
3
math/pow/y^x
mem/RCL 2
math/simple/x
mem/RCL 2
1
math/matrix/concat
math/matrix/concat
math/matrix/concat       ; T = [ t^3  t  t  1 ]
mem/RCL 1                ; M (spline basis matrix)
*
stack/x<->y
*                        ; [ x  y ] = T*M*G
-1
math/matrix/split        ; x, y
stack/x<->y
trig/coord/r->cplx       ; z = x+yi
prog/finish              ; program finished

Preparations: Create an Nx2 matrix containing the control points and store them in memory location 0. Create a 4x4 spline basis matrix and store it in memory location 1. Push suitable graph limits on the stack before executing "prog/draw/z=f(t)/Splin".

Example set of 2D control points:
     /  2 22 \       ...
     |  7 19 |     | 21  6 |
     | 14 18 |     | 23  0 |
     | 13 22 |     | 26  7 |
     |  2 11 |     | 30 18 |
     |  4  2 |     | 29 22 |
     | 11  2 |     | 25  9 |
     | 18  9 |     | 27  0 |
     | 21  9 |     | 31  5 |
     | 22  9 |     | 34  9 |
     | 19  9 |     | 37  9 |
     | 16  6 |     | 38  9 |
     | 16  2 |     | 35  9 |
     | 20  2 |     | 32  6 |
     | 22  5 |     | 32  2 |
     | 22  7 |     | 36  2 |
     | 23 10 |     | 38  3 |
       ...         \ 39  4 /

Examples of spline basis matrices:
 Catmull-Rom spline:     B-spline:

     / -1  3 -3  1 \         / -1  3 -3  1 \
 1/2*|  2 -5  4 -1 |     1/6*|  3 -6  3  0 |
     | -1  0  1  0 |         | -3  0  3  0 |
     \  0  2  0  0 /         \  1  4  1  0 /

A Catmull-Rom spline interpolates the control points but is only first order continuous. A B-spline does not interpolate the control points (it most often winds between the points), but it is second order continuous, and therefore often smoother.

Plotting a circle

There are several possibilites to plot a circle or ellipse, using different plotting commands.
  1. Polar graph, circle centered at (0,0):
    prog/new/"a" 
    clear ; To remove input (phi) 
    1 ; Radius 
    prog/finish 
     
    -1 
    1 
    -1 
    1 
    prog/draw/r=f(phi)/"a" 
    
  2. Parametric curve, in this case an ellipse centered at (5,5):
    prog/new/"b" 
    2 
    trig/pi 
    * 
    * ; 2*pi*t 
    ENTER 
    trig/sin 
    3 ; y radius 
    * 
    5 ; y center 
    + 
    stack/x<>y 
    trig/cos 
    4 ; x radius 
    * 
    5 ; x center 
    + 
    trig/coord/r->cplx 
    prog/finish 
     
    0 
    10 
    0 
    10 
    prog/draw/z=f(t)/"b" 
    
  3. Positive and negative halves combined into one complex number. This exploits the fact that prog/draw/y=f(x) draws the imaginary part of a complex number in a different colour.
    prog/new/"c" 
    math/x 
    +/- 
    1 
    + 
    math/sqrt(x) 
    ENTER 
    +/- 
    trig/cord/r->cplx 
    prog/finish 
     
    -1 
    1 
    -1 
    1 
    prog/draw/y=f(x)/"c" 
    
  4. Function result alternating between positive and negative halves. This exploits the fact that points on the graph appear semi-randomly, so you might draw any number of curves simultaneously by alternating between them.
    0 
    mem/STO 0 
     
    prog/new/"d" 
    math/x 
    +/- 
    1 
    + 
    math/sqrt(x) 
    ENTER 
    +/- 
    mem/RCL 0 
    +/- 
    1 
    + 
    mem/STO 0 
    prog/util/select 
    prog/finish 
     
    -1 
    1 
    -1 
    1 
    prog/draw/y=f(x)/"d" 
    
  5. Advanced example using the f(z)=z plot.
  6. Plotting a sphere

    A simple example for demonstrating how to plot 3D graphs with the f(z)=z plot, i.e. functions with a tuple (x,y) as input and a real value f(x,y) as output, is a unit sphere. In the f(z)=z plot, the angle of the complex number f(z) determines the color for each point, while abs(f(z)) determines the intensity. If the function f(z) has only real function values, then the positive values will be shades of green and the negative values shades of magenta. It is a good idea to reduce the values plotted to the interval [0,1], since the colors in the intervals [j,j+1] for any natural number j are identical and it is not easy to determine the point with the maximal value then. So an appropriate scaling is necessary which depends on the problem.
    We use the function f(z)=sqrt(max(0,(1-abs(z)))), which displays a nice unit sphere:

    Program listing:
    prog/new "Spher"
    trig/cplx/abs   ; be sure to have a complex number on stack
                    ; entry x before!
    math/simple/x^2
    1
    -
    +/-
    0
    prog/util/max   
    math/simple/√x
    prog/finish
    
    Plot this with the following bounds and command:
    -1 1 -1 1 prog/draw/f(z)=z/"Spher"
    

    Mandelbrot fractals

    Another nice example for the f(z)=z plotting mode is plotting a Madelbrot fractal. It also demonstrates the simple usage of the DSE loop command. Remember to turn on mode/monitor/prog when programming this, it will make it easier to spot and fix programming errors. (Consult the manual to learn more)

    Program listing:
     
    prog/new/"mbrot" 
    mem/STO 0 
    16 
    mem/STO 1 
    clear 
    prog/flow/LBL 0 
    x  
    mem/RCL 0  
    +  
    prog/flow/DSE 1 
    prog/flow/GTO 0 
    prog/util/abs  
    math/pow/log_2  
    math/pow/log_2  
    ENTER  
    16  
    /  
    trig/tanh  
    trig/coord/p->r  
    trig/coord/r->cplx  
    prog/finish
    
    The calculations at the end of the program approximates the number of iterations from log(log(|z|)), then calculates a new complex number based on this, to give the resulting image some nice colors. One problem now is that the somewhat ad-hoc color calculation towards the end only works for around 16 iterations, so if you increase the number of iterations, it might not look as good. Of course, now that we have control flow, we can do more sophisticated things such as actually setting the color to black (0) when the mandelbrot loop has not diverged. I leave this as an exercise. You can draw it using something like this:
    -2 1 -1.5 1.5 prog/draw/z=f(z)/"mbrot"
    

    Linear Programming

    Assume you are given a system of linear inequalities in two variables, e.g. you have two kinds of resources at hand and want to produce some products from them, but you have to maintain the capacity constraints of your factory machinery. You want to know what the maximum profit is. I leave it to your imagination to come to this abstract form of the problem.
    Maximize x+y
    subject to inequalities: x+2y ≤ 70, 2x+y ≤ 50, x≥0, y≥0
    
    To solve this, we first transfer all inequalities to the form a1*x+a2*y≤b, where we can change ≥ to ≤ by multiplying with -1. We can then then collect all these inequalities into a matrix A, each line being one inequaliy: A*(x,y) ≤ b This leads to the following standard form: max c*(x,y), s.t. A*(x,y)≤b, with c and b vectors and A a matrix over the reals.
        / 1  2 \       / 70 \
    A = | 2  1 | , b = | 50 | , c = ( 1, 1 )
        | -1 0 |       |  0 |
        \ 0 -1 /       \  0 /
    
    We will now plot the set of solutions to the inequalities using the f(z)=z plot function as a green area, where bright green stands for big values of x+y and dark green for small values of x+y, the function we want to maximize. First enter the matrixes above and store A in memory 4, b in memory 5 and c in memory 6.

    Program listing:
    prog/new/"Simpl"                ; shorthand for "Simplex"  :) 
    trig/complex/cplx->r            ; split stack x+yi into x and y
    matrix/stack                    ; transform into vector [x,y]
    ENTER                           ; copy for later usage
    mem/RCL 4                       ; A
    stack/x<->y
    *                               ; A*[x,y] computed
    mem/RCL 5                       ; b
    -                               ; transform to A*[x,y]-b <= 0
    matrix/a_max                    ; check if all components/lines ...
    stack/x<->y
    clear
    prog/flow/x<=0                  ; ... satisfy <= 0
    prog/flow/GTO 1                 ; if yes, goto calculation of c*x
    clear                           ; case: some ineqality not satisfied
    clear
    0                               ; ... hence function result 0
    prog/flow/STOP
    prog/flow/LBL 1                 ; case: all ineq. satisfied
    clear
    mem/RCL 6                       ; c
    stack/x<->y
    *                               ; calculate c*x
    0.025
    *                               ; scale to 0.025*c*x    -> (1)
    prog/finish
    
    Then you may plot the simplex by
    -5 40 -5 40 prog/plot/f(z)=z/"Simpl"
    
    Pick the corner of the green area which has the brightest color -- that will be your solution. More information about solving linear programs can be found in the literature, including algorithms for cases with more than 2 variables.

    Note: The scaling (1) is used to reduce the values plotted to the interval [0,1], since the colors in the intervals [j,j+1] for any natural number j are identical and it is not easy to determine the point with the maximal value then. So an appropriate scaling is necessary which depends on the problem.

    Plotting curves

    Some curves are the a set of zeros of a two dimensional function, like x+3y-2x+3y-1=0. We show how to plot this equation using the f(z)=z plot. But instead of plotting f(x,y)=0, we use the following form for our example:
    | x+3y-2x+3y-1 | < ε
    where we choose ε (epsilon) to be some suitable value. We will use ε=0.1 here. Ok, first type in the function (it is useful to put something like 3+4i on the stack first, so you can have a look at what happens if running the program for x=3 and y=4).

    Program listing:
    prog/new "curve"  
    monitor/prog/4       ; not needed, but useful if you want to see what you typed :) 
    trig/complex/cplx->r ; transform x+iy input to x and y coords 
    ENTER                ; now we start to compute x^2+3y^2-2x+3y-1  
    2 
    * 
    stack/x<->y 
    x^2 
    +                    ; we have x+2x on stack 
    stack/x<->y          ; now the stuff involving y ... 
    ENTER 
    3 
    * 
    stack/x<->y 
    x^2 
    3 
    * 
    +                    ; now we have 3y+3y 
    + 
    -1 
    +                    ; now we have x+2x+3y+3y-1 
    prog/util/abs 
    0.1                  ; epsilon 
    prog/flow/x>y?       ; compare if | ... | < 0.1 
    prog/flow/GTO 1      ; it WAS smaller 0.1 
    clear 
    clear 
    0                    ; in case | ... | ≥ 0.1, function output will be 0 
    prog/flow/STOP
    prog/flow/LBL 1  
    clear 
    clear 
    1                    ; in case | ... | < 0.1, function output will be 1 
    prog/finish
    
    Now we plot a nice curve (ellipse, that is) by
     
    -3 1 -2 1 prog/draw/z=f(z)/"curve" 
    
    Depending on your device's speed you will sooner or later see the curve, being refined from a rather bulky version in the beginning. Remember that you can redraw the coordinate system by pressing 0 in graph mode. If you want to have a finer version, you may decrese epsilon to a smaller value. The smaller the value of epsilon is, the better the curve looks and the more accurate the plot is. But ... if it is too small, you will have to wait ages, before something starts appearing, maybe forever.

    More

    If you think you've made a clever program, send it to me (roarl at users.sourceforge.net), and I might publish it here.