Chip 2001 November/Chip_2001-11_cd1.bin/sharewar/chaospro/cpro302.exe/LKM.CCL

Text File | 2001-07-15 | 116.9 KB | 5,977 lines

comment {

Coloring schemes developed by Kerry Mitchell

Compilation dated 15jul2001

Includes:
  Basic
  Basic Plus
  Emboss
  Gaussian Integer
  Statistics
  Polar Curves
  Astroid
  5 Point Star
  Conic Sections
  Conic Lite
  2 Sections
  Conic Lite Old
  Bubbles
  Range
  Passthru
  Distance to a Point
  Pythagorean Triple
  Spiral
  Grid
  Single Truchet
  Double Truchet
  Rose Curve String Art
  Astroid String Art
  Crosshatch
}


basic { // Kerry Mitchell 20sep98
//
// The standard coloring parameters:
// polar angle (includes all 4 quadrants
// instead of only 2 from atan),
// magnitude of z, real(z), imag(z),
// and iteration number.
//
real Var_x;
real Var_y;
parameter int colorby;
real t;

    void init(void)
    {
    }
    void loop(void)
    {
    }
    void final(void)
    {
        Var_x=real(z);
        Var_y=imag(z);
        if ((colorby==1))
        {// polar angle
            t=atan2(z);
            t=t/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==2))
        {// magnitude
            index=cabs(z);
        }
        else if ((colorby==3))
        {// real(z)
            index=Var_x;
        }
        else if ((colorby==4))
        {// imag(z)
            index=Var_y;
        }
        else
        {
                     // iteration
            index=0.01*numiter;
        }
    }
    void description(void)
    {
        this.title="Basic";

        colorby.caption = "color by";
        colorby.default = 0;
        colorby.enum = "iteration\npolar angle\nmagnitude\nreal(z)\nimag(z)";

    }
}


basic-plus { // Kerry Mitchell 26feb99
//
// The standard coloring parameters:
// polar angle, magnitude, real, imag
// for z, pixel, z-pixel, & z/pixel
// and iteration number.
//
// Replaces the "Basic" coloring formula
//
complex var;
real t;
parameter int vartype;
parameter int colorby;

    void init(void)
    {
        var=(0.0,0.0);
        t=0.0;
    }
    void loop(void)
    {
    }
    void final(void)
    {
//
// establish variable used for coloring
//
        if ((vartype==1))
        {// pixel
            var=pixel;
        }
        else if ((vartype==2))
        {// z-pixel
            var=z-pixel;
        }
        else if ((vartype==3))
        {// z/pixel
            var=z/pixel;
        }
        else
        {
                    // z
            var=z;
        }
//
// pick coloring characteristic
//
        if ((colorby==1))
        {// polar angle
            t=atan2(var);
            t=t/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==2))
        {// magnitude
            index=cabs(var);
        }
        else if ((colorby==3))
        {// real
            index=real(var);
        }
        else if ((colorby==4))
        {// imag
            index=imag(var);
        }
        else
        {
                     // iteration
            index=0.01*numiter;
        }
    }
    void description(void)
    {
        this.title="Basic Plus";

        vartype.caption="variable";
        vartype.default=0;
        vartype.enum="z\npixel\nz-pixel\nz/pixel";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="iteration\npolar angle\nmagnitude\nreal part\nimag part";

    }
}


emboss { // Kerry Mitchell 13sep98
//
// For use with "Emboss" formulas and
// gradient.  Uses 3 colors:
// #index values of 0.2, 0.5, and 0.8.
//

    void final(void)
    {
        if ((real(z)<imag(z)))
        {
            index=0.2;
        }
        else if ((imag(z)<real(z)))
        {
            index=0.8;
        }
        else
        {

            index=0.5;
        }
    }
    void description(void)
    {
        this.title="Emboss";
    }
}


comment { ; narrative copyright Kerry Mitchell 20sep98

Gaussian Integer Coloring

Gaussian integers are complex numbers such that both the real and
imaginary parts are integers.  Examples are:  (0,0), (-2,3), (17,-5),
and (1000000,123456789).  The gaussian scheme is concerned with how
the orbit behaves relative to a Gaussian integer.

To find the Gaussian integer which the orbit most closely approaches,
the built-in function round() is used.  Round(z) returns a complex
number whose components are the rounded components of z.  This is a
Gaussian integer.  The distance from z to round(z) is simply the
magnitude of z - round(z).  The "minimum distance" method tracks this
distance and records its smallest value.  The "iteration @ min" colors
by the iteration number when this minimum is reached, and the "angle @
min" methods colors by the polar angle of z when the minimum was reached.
The corresponding "maximum" methods ("maximum distance", "iteration @
max", and "angle @ max") work in a similar fashion for the maximum value
of z - round(z).  The final method, "average distance", simply colors by
the mean of the distance over the course of the orbit.  Note that the
angle used here is the actual polar angle, with a range of 360 degrees,
instead of the angle returned by the built-in "atan" function, which
only has a range of 180 degrees.

These methods essentially overlay a 1x1 square grid onto the complex
plane and ask how close (or far) the orbit comes to a node in the grid.
Two options allow varying the size of the grid; these are called the
"normalization" methods.  The first is to normalize z by the #pixel.
This has the effect of asking how close the orbit comes to Gaussian
integer multiples of the #pixel value (or, z/#pixel = has integer real
and imaginary parts).  The other normalization method employs a user-
specified complex factor.  Then, the coloring is performed according to
how close the orbit comes to Gaussian integer multiples of this factor.
For example, if the factor was set to (1.2,5.1), then the methods would be
coloring by the orbits approach to nodes of a rectangular grid that was
1.2 units long by 5.1 units high.  This might be useful with Julia set
images, as the factor can be set to the Julia parameter.  The default,
"none" normalization, uses a factor of 1.
}


gauss { // Kerry Mitchell 20sep1998, last updated 28apr2001
//
// Colors by orbit's relationship
// to Gaussian Integers. See comment
// block for more information.
//
// Updated 02apr99 to include options for finding
// the integer:  added trunc, floor, and ceil.
// Thanks to Marcelo Anelli for the idea.
//
// Updated 28apr2001 to include normalizing by some function of z.
// Added to new 'color by' choices:  max/min ratio and
// min/mean/max angle.  Added "randomize" feature to randomly
// change z before finding the Gaussian integer.  This will
// change the regular grid-like distribution of points.
//
real r;
real rmin;
real rmax;
real rave;
real total;
real t;
int iter;
int itermin;
int itermax;
complex zmin;
complex zmax;
parameter int norm;
complex normfac;
parameter complex fac;
real logfac;
parameter real logseed;
complex temp2;
parameter bool randomize;
parameter complex randomsize;
parameter int inttype;
complex temp;
complex remain;
parameter int colorby;

    void init(void)
    {
        r=0.0;
        rmin=1.0e12;
        rmax=0.0;
        rave=0.0;
        total=0.0;
        t=0.0;
        iter=0;
        itermin=0;
        itermax=0;
        zmin=(0.0,0.0);
        zmax=(0.0,0.0);
        if ((norm==1))
        {// pixel normalization
            normfac=pixel;
        }
        else if ((norm==2))
        {// factor normalization
            normfac=fac;
        }
        else if ((norm==3))
        {// f(z) normalization
            normfac=normfunc(z);
        }
        else
        {
                               // no normalization
            normfac=(1.0,0.0);
        }
        logfac=logseed;
    }
    void loop(void)
    {
        iter=iter+1;
        temp2=z;
        if ((randomize==true))
        {
            logfac=4*logfac*(1-logfac);
            temp2=temp2*(1-randomsize*logfac);
        }
        if ((inttype==1))
        {// trunc
            temp=trunc(temp2/normfac);
        }
        else if ((inttype==2))
        {// floor
            temp=floor(temp2/normfac);
        }
        else if ((inttype==3))
        {// ceil
            temp=ceil(temp2/normfac);
        }
        else
        {
                               // round
            temp=round(temp2/normfac);
        }
        remain=temp2-temp*normfac;
        r=cabs(remain);
        total=total+r;
        rave=total/iter;
        if ((r<rmin))
        {
            rmin=r;
            zmin=temp2;
            itermin=iter;
        }
        if ((r>rmax))
        {
            rmax=r;
            zmax=temp2;
            itermax=iter;
        }
    }
    void final(void)
    {
        if ((colorby==1))
        {// iteration @ min
            index=0.01*itermin;
        }
        else if ((colorby==2))
        {// angle @ min
            t=atan2(zmin);
            t=t/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==3))
        {// maximum distance
            index=rmax;
        }
        else if ((colorby==4))
        {// iteration @ max
            index=0.01*itermax;
        }
        else if ((colorby==5))
        {// angle @ max
            t=atan2(zmax);
            t=t/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==6))
        {// average distance
            index=rave;
        }
        else if ((colorby==7))
        {// min/mean/max angle
            zmax=(rave-rmin)+flip(rmax-rave);
            t=atan2(zmax);
            t=t/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==8))
        {// max/min ratio
            index=rmax/(rmin+1.e-12);
        }
        else
        {
                       // minimum distance
            index=rmin;
        }
    }
    void description(void)
    {
        this.title="Gaussian Integer";

        inttype.caption="integer type";
        inttype.default=0;
        inttype.enum="round(z)\ntrunc(z)\nfloor(z)\nceil(z)";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="minimum distance\niteration @ min\nangle @ min\nmaximum distance\niteration @ max\nangle @ max\naverage distance\nmin/mean/max angle\nmax/min ratio";


        norm.caption="normalization";
        norm.default=0;
        norm.enum="none\npixel\nfactor\nf(z)";


        fac.caption="normalizing factor";
        fac.default=(2.0,1.0);


        randomize.caption="randomize?";
        randomize.default=false;
        randomize.hint="Applies a random factor to z every iteration before finding the Gaussian integer.";


        randomsize.caption="random size";
        randomsize.default=(0.1,0);
        randomsize.hint="Size of random factor, if 'randomize?' is checked.";


        logseed.caption="random seed";
        logseed.default=0.1;
        logseed.min=0.0;
        logseed.max=1.0;
        logseed.hint="Randomize seed, between 0 and 1.";


        normfunc.caption="normalizing function";
        normfunc.default="ident" ;
        normfunc.hint="For 'f(z)' normalization.";

    }
}


comment { ; copyright Kerry Mitchell 04nov98

Statistics

This started out with an article written by Stephen Ferguson on an
analysis of fractal dimension, which he based on an algorithm by Holger
Jaenisch.  The formulas here modify Stephen's analysis and add other
standard statistical measures.  Since the measures are typically used
with bounded datasets, this coloring method may be more applicable as
an "inside" scheme, but there's nothing stopping a user from employing
it as an "outside" scheme as well.

The statistical measures implemented are:  minimum, maximum, range,
mean, standard deviation, coefficient of variation, and fractal
dimension.  All are defined only for real variables, so there are
4 choices for reducing the complex #z to a real number:  real(#z),
imag(#z), the magnitude of #z, and imag(#z)/real(#z).  The last
method is an attempt to capture the polar angle of #z, without the
discontinuities involved in actually using the angle arg(#z).

Once the choice of variable has been made, the first 3 measures are
simple enough to compute.  Simply monitor the orbit, updating the
minimum and maximum as they change.  Once the iteration has ceased,
the range is just (maximum - minimum).  The mean is computed by keeping
a running sum of the variable, then dividing it by the number of
iterations.

The standard deviation is a measure of the spread of the data, and is
defined in terms of the sum of the squared differences between each
datum and the mean.  This can also be computed by keeping a running
sum of the square of the variable.  This sum is used with the sum for
the mean to determine the standard deviation.  The coefficient of
variation is a normalized standard deviation:  it's the ratio of the
standard deviation to the mean.  Finally, the "fractal dimension"
computed here is not the true fractal dimension, but an approximation
to it.  It's the standard deviation normalized by the range.

}


statistics { // Kerry Mitchell 03nov98
//
// Colors according to various
// statistical properties of the
// iterate.  Primarily an "inside"
// coloring method.
//
int iter;
real r;
real rmin;
real rmax;
real sumsum;
real sum;
parameter int vartype;
parameter int stattype;

    void init(void)
    {
        iter=0;
        r=0.0;
        rmin=1e12;
        rmax=-rmin;
        sumsum=0.0;
        sum=0.0;
    }
    void loop(void)
    {
        iter=iter+1;
        if ((vartype==0))
        {// real(z)
            r=real(z);
        }
        else if ((vartype==1))
        {// imag(z)
            r=imag(z);
        }
        else if ((vartype==2))
        {// magnitude(z)
            r=cabs(z);
        }
        else
        {
                          // imag(z)/real(z)
            r=imag(z)/real(z);
        }
        sum=sum+r;
        sumsum=sumsum+r*r;
        if ((r<rmin))
        {
            rmin=r;
        }
        if ((r>rmax))
        {
            rmax=r;
        }
    }
    void final(void)
    {
        if ((stattype==0))
        {// minimum
            index=rmin;
        }
        else if ((stattype==1))
        {// maximum
            index=rmax;
        }
        else if ((stattype==2))
        {// range
            index=rmax-rmin;
        }
        else if ((stattype==3))
        {// mean
            index=sum/iter;
        }
        else if ((stattype==4))
        {// standard deviation
            index=sqrt(sumsum-sum*sum/iter);
        }
        else if ((stattype==5))
        {// coefficient of variation
            index=sqrt(sumsum-sum*sum/iter)/sum*iter;
        }
        else
        {
                          // fractal dimension
            index=sqrt(sumsum-sum*sum/iter)/(rmax-rmin);
        }
    }
    void description(void)
    {
        this.title="Statistics";

        vartype.caption="variable";
        vartype.default=2;
        vartype.enum="real(z)\nimag(z)\nmagnitude\nimag/real";
        vartype.hint="variable for which statistics will be calculated";


        stattype.caption="statistic";
        stattype.default=6;
        stattype.enum="minimum\nmaximum\nrange\nmean\nstd. deviation\nvariation\nfractal dimension";
        stattype.hint="statistic that will be calculated";

    }
}


comment { ; copyright Kerry Mitchell 15nov98

Polar Curves

Typically, points in a plane are thought of in terms of their x- and
y-coordinates, that is, how far away (and on which side) the point is
from the horizontal x-axis and the vertical y-axis.  Another way of
looking at the point is with polar coordinates, which specify the
distance of the point from the origin (r) and its direction (t).  The
two methods are equivalent:

x = r*cos(t), y=r*sin(t), or
r^2 = x^2 + y^2, tan(t) = y/x.

Polar curves are curves that specify r as a function of t, instead of
y as a function of x.  The curve used in this coloring method is:

r = [a * f(b*t)]^n + r0,

where f is one of UltraFractal's builtin functions (e.g., sin, cos, exp,
etc.), a is an amplitude scaling factor, and b is a frequency factor.
The exponent n is useful for making the curve wider or thinner, and r0
is a expansion or contraction constant.

Some special curves can be generated using this function.  Spirals can
be made by using the "ident" function.  Here, the exponent n controls
how tightly wound the spiral is.  However, only one revolution of the
spiral will be drawn, as t is limited to the range 0 to 2*pi.

"Rose" curves are made by using either sin or cos functions.  The
parameter a controls the size of the curve.  The number of "petals"
is b, so long as b is a positive integer.  Increasing n from 1 will
make the petals thinner; decreasing it toward 0 will make them thicker.
Leave r0=0 for the standard rose curve, where the petals all join at
the origin.

Since the sin and cos functions generate negative values, the "rose"
curves will have some regions of negative r.  How this is handled
depends on r0 and "negrflag", the negative r flag.  Setting negrflag to
1 will make the routine ignore negative r values.  This, with r0 is set
to 0, will cause the rose curve to have "b" number of petals, all of
them the same size.  Setting negrflag=2 will make the routine consider
negative r's the same as positive r's.  Thus, the rose curve will have
2*b petals.  Increasing r0 from 0 will make r positive more often than
negative.  This will also increase the number of petals from b to 2*b,
but half of the petals will be small and half will be large.

The best way to see the effects of the parameter choices is to use the
"draw curves" setting of the "color by" parameter.  Here, the image
isn't a fractal, but rather the polar curve determined by a, b, n, r0,
the chosen function, and the negative r flag.

The other "color by" settings use polar curves to color fractals.  There
are 2 basic rendering methods:  how close the orbit comes to the polar
curve, and how often the orbit is inside the polar curve.  The minimum
approach method also has a few variants:  iteration number when the
minimum was reached, and angle of #z when the minimum was reached.

}


polar-curves { // Kerry Mitchell 15nov98
//
// colors by the relationship between
// the orbit and a user-specified polar
// curve (such as a spiral or a rose curve)
//
real err;
real errmin;
real r;
real twopi;
real t;
int iter;
int itermin;
int incount;
complex zmin;
parameter real a;
parameter real b;
parameter real n;
parameter real r0;
parameter int rneg;
parameter int colorby;

    void init(void)
    {
        err=0.0;
        errmin=1.0e12;
        r=0.0;
        twopi=2.0*pi;
        t=0.0;
        iter=0;
        itermin=0;
        incount=0;
        zmin=(0.0,0.0);
    }
    void loop(void)
    {
        iter=iter+1;
        t=imag(log(z));
        if ((t<0.0))
        {
            t=t+twopi;
        }
        r=a*real(fn1(b*t));
        if ((r>=0.0))
        {
            r=r^n;
        }
        else
        {

            r=-((-r)^n);
        }
        r=r+r0;
        if ((rneg==0))
        {
            err=r-cabs(z);
        }
        else
        {

            err=|r|-|z|;
        }
        if ((err<0.0))
        {
            incount=incount+1;
            err=-err;
        }
        if ((err<errmin))
        {
            errmin=err;
            zmin=z;
            itermin=iter;
        }
    }
    void final(void)
    {
        if ((colorby==0))
        {// minimum distance
            index=cabs(log(errmin));
        }
        else if ((colorby==1))
        {// iteration @ min
            index=0.01*itermin;
        }
        else if ((colorby==2))
        {// angle @ min
            t=atan2(zmin);
            t=t/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==3))
        {// in fraction
            index=incount/iter;
        }
        else
        {
                          // draw curves
            t=imag(log(pixel));
            if ((t<0.0))
            {
                t=t+twopi;
            }
            r=a*real(fn1(b*t));
            if ((r>=0.0))
            {
                r=r^n;
            }
            else
            {

                r=-((-r)^n);
            }
            r=r+r0;
            if ((rneg==0))
            {
                err=r-cabs(pixel);
            }
            else
            {

                err=|r|-|pixel|;
            }
            index=cabs(log(cabs(err)));
        }
    }
    void description(void)
    {
        this.title="Polar Curves";

        a.caption="amplitude";
        a.default=1.0;


        b.caption="frequency";
        b.default=1.0;


        n.caption="exponent";
        n.default=1.0;


        r0.caption="baseline r";
        r0.default=0.0;


        rneg.caption="r<0 flag";
        rneg.default=0;
        rneg.enum="ignore r<0\ntreat as r>0";
        rneg.hint="1 to ignore negative r values, 2 to treat as positive";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="minimum distance\niteration @ min\nangle @ min\nin fraction\ndraw curves";

    }
}


comment { ; copyright Kerry Mitchell 15nov98

Astroid

The astroid is a figure from analytic geometry, resembling a four-
pointed star with concave sides.  Its defining equation is:

x^(2/3) + y^(2/3) = a^(2/3)

where a determines the size of the figure, similar to the radius of
a circle.  This equation can be generalized by changing the exponent
of 2/3 to any value n.  If n is between 0 and 1, the figure resembles
the standard astroid.  The sides go from being straight lines for n=1
to lying atop of the coordinates axes as n approaches 0.  For n>1, the
figure becomes convex, and is a circle for n=2.  As n approaches
infinity, the figure approaches a square.  The astroid is further
generalized by allowing it to be placed somewhere other than at the
center of the complex plane.  The figure's orientation and location
in the plane are determined by a "center" and "rotation angle" parameters.

This generalized astroid is the basis of this coloring scheme.  As with
other plane figures, the astroid can be compared to the Mandelbrot and
Julia orbits.  The image can be colored by how close the orbit comes
to the astroid or how often the orbit lands inside the astroid.  A final
choice of the "colorby" parameter draws the astroid only, so the effects
of the parameter choices can be seen.

}


astroid { // Kerry Mitchell 15nov98
//
// colors according to the relationship
// between the orbit and an astroid
// (4-pointed star figure)
//
real err;
real errmin;
real Var_x;
real Var_y;
real t;
real aton;
parameter real a;
parameter real n;
int iter;
int itermin;
int incount;
complex zmin;
parameter real rotangle;
complex rot;
complex temp;
parameter complex Parm_center;
parameter int colorby;

    void init(void)
    {
        err=0.0;
        errmin=1.0e12;
        Var_x=0.0;
        Var_y=0.0;
        t=0.0;
        aton=a^n;
        iter=0;
        itermin=0;
        incount=0;
        zmin=(0.0,0.0);
        t=rotangle/180.0*pi;
        rot=exp(flip(t));
    }
    void loop(void)
    {
        iter=iter+1;
        temp=rot*(z-Parm_center);
        Var_x=cabs(real(temp));
        Var_y=cabs(imag(temp));
        err=Var_x^n+Var_y^n-aton;
        if ((err<0.0))
        {
            incount=incount+1;
            err=-err;
        }
        if ((err<errmin))
        {
            errmin=err;
            zmin=z;
            itermin=iter;
        }
    }
    void final(void)
    {
        if ((colorby==0))
        {// minimum distance
            index=cabs(log(errmin));
        }
        else if ((colorby==1))
        {// iteration @ min
            index=0.01*itermin;
        }
        else if ((colorby==2))
        {// angle @ min
            t=atan2(zmin);
            t=t/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==3))
        {// in fraction
            index=incount/iter;
        }
        else
        {
                          // draw curve
            temp=rot*(pixel-Parm_center);
            Var_x=cabs(real(temp));
            Var_y=cabs(imag(temp));
            err=Var_x^n+Var_y^n-aton;
            index=cabs(log(cabs(err)));
        }
    }
    void description(void)
    {
        this.title="Astroid";

        a.caption="size";
        a.default=1.0;


        n.caption="exponent";
        n.default=0.67;


        rotangle.caption="rotation, degrees";
        rotangle.default=0.0;


        Parm_center.caption="center";
        Parm_center.default=(0.0,0.0);


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="minimum distance\niteration @ min\nangle @ min\nin fraction\ndraw curve";

    }
}


comment { ; narrative copyright Kerry Mitchell 25nov98

Seeing Stars

Most fractal images involve circles in some respect:  either stopping
the iteration when the orbit moves outside of a given circle, or
coloring by how close the orbit comes to a certain circle, or some
variation.  This formula uses a 5 point star instead of a circle.
The fractal can be colored by how close the orbit comes to a star, or
how often the orbit lands inside the star.

The star is represented by the 5 outer points.  These are equally spaced
on a circle.  The center and size of the circle are user-specified, as
is the rotation angle of the star.

How z at any iteration relates to the star (inside, outside, how close)
is determined by looking at each of the 10 sides, one at a time.  Each
side can be represented by a line A*x + B*y + C = 0, where A, B, and C
come from the coordinates of the 2 outer points that are joined to make
the side.  Given the numbers A, B, and C, the quantity

        q = A*real(z) + B*imag(z) + C

is computed.  If q is positive, then z is on one side of the line, and
if q is negative, then z is on the other side of the line.  If |q| is
very small, that means that z is very close to that side of the star.
Taking the signs of q for all 10 sides will determine if z is inside or
outside of the star.  Using the smallest value of |q| gives the distance
from z to the star.

To see this in action, use the "draw star" setting of the "color by"
parameter.  The outline of a star will be drawn, with the position,
orientation and size that you choose.  This setting along with the
"minimum distance" setting, uses a logarithmic transfer function from
distance to #index.  This helps highlight the star without needing to
find out if the point is exactly on the star or not.  The "in fraction"
will count how many times the orbit lands inside the star, then set
#index to the ratio of inside hits to total iterations for that pixel.

}


5star { // Kerry Mitchell 25nov98
//
// colors by relationship to
// 5-point star
//
real err;
real errmin;
real Var_x;
real Var_y;
real xcen;
parameter complex Parm_center;
real ycen;
real phi;
parameter real phiangle;
real twopi;
real temp;
real t;
real t0;
real t1;
real t2;
real t3;
real t4;
real t5;
real t6;
real t7;
real t8;
real t9;
real x0;
parameter real r;
real y0;
real x1;
real y1;
real x2;
real y2;
real x3;
real y3;
real x4;
real y4;
int flag;
int incount;
int iter;
int itermin;
complex zmin;
parameter int colorby;

    void init(void)
    {
        err=0.0;
        errmin=1.0e20;
        Var_x=0.0;
        Var_y=0.0;
        xcen=real(Parm_center);
        ycen=imag(Parm_center);
        phi=phiangle/180.0*pi;
        twopi=2.0*pi;
        temp=0.1*twopi;
        t=0.0;
        t0=phi;
        t1=t0+temp;
        t2=t1+temp;
        t3=t2+temp;
        t4=t3+temp;
        t5=t4+temp;
        t6=t5+temp;
        t7=t6+temp;
        t8=t7+temp;
        t9=t8+temp;
        x0=r*cos(t0)+xcen;
        y0=r*sin(t0)+ycen;
        x1=r*cos(t2)+xcen;
        y1=r*sin(t2)+ycen;
        x2=r*cos(t4)+xcen;
        y2=r*sin(t4)+ycen;
        x3=r*cos(t6)+xcen;
        y3=r*sin(t6)+ycen;
        x4=r*cos(t8)+xcen;
        y4=r*sin(t8)+ycen;
        flag=0;
        incount=0;
        iter=0;
        itermin=0;
        zmin=(0.0,0.0);
    }
    void loop(void)
    {
        iter=iter+1;
        Var_x=real(z);
        Var_y=imag(z);
        t=imag(log(z-Parm_center));
        if ((t<0))
        {
            t=t+twopi;
        }
  //
  //  compute how close iterate is to each side of star
  //  and on which side
  //
        if (((t>t0)&&(t<=t1)))
        {
            err=Var_x*(y0-y2)+Var_y*(x2-x0)-x2*y0+x0*y2;
            if ((err<0))
            {
                flag=1;
            }
        }
        else if (((t>t1)&&(t<=t2)))
        {
            err=Var_x*(y1-y4)+Var_y*(x4-x1)-x4*y1+x1*y4;
            if ((err<0))
            {
                flag=1;
            }
        }
        else if (((t>t2)&&(t<=t3)))
        {
            err=Var_x*(y1-y3)+Var_y*(x3-x1)-x3*y1+x1*y3;
            if ((err<0))
            {
                flag=1;
            }
        }
        else if (((t>t3)&&(t<=t4)))
        {
            err=Var_x*(y2-y0)+Var_y*(x0-x2)-x0*y2+x2*y0;
            if ((err<0))
            {
                flag=1;
            }
        }
        else if (((t>t4)&&(t<=t5)))
        {
            err=Var_x*(y2-y4)+Var_y*(x4-x2)-x4*y2+x2*y4;
            if ((err<0))
            {
                flag=1;
            }
        }
        else if (((t>t5)&&(t<=t6)))
        {
            err=Var_x*(y3-y1)+Var_y*(x1-x3)-x1*y3+x3*y1;
            if ((err<0))
            {
                flag=1;
            }
        }
        else if (((t>t6)&&(t<=t7)))
        {
            err=Var_x*(y3-y0)+Var_y*(x0-x3)-x0*y3+x3*y0;
            if ((err<0))
            {
                flag=1;
            }
        }
        else if (((t>t7)&&(t<=t8)))
        {
            err=Var_x*(y4-y2)+Var_y*(x2-x4)-x2*y4+x4*y2;
            if ((err<0))
            {
                flag=1;
            }
        }
        else if (((t>t8)&&(t<=t9)))
        {
            err=Var_x*(y4-y1)+Var_y*(x1-x4)-x1*y4+x4*y1;
            if ((err<0))
            {
                flag=1;
            }
        }
        else
        {

            err=Var_x*(y3-y0)+Var_y*(x0-x3)-x0*y3+x3*y0;
            if ((err<0))
            {
                flag=1;
            }
        }
        err=cabs(err);
        if ((err<errmin))
        {
            errmin=err;
            itermin=iter;
            zmin=z;
        }
        if ((flag==1))
        {
            incount=incount+1;
        }
    }
    void final(void)
    {
        if ((colorby==0))
        {// minimum distance
            index=cabs(log(errmin));
        }
        else if ((colorby==1))
        {// iteration @ min
            index=0.01*itermin;
        }
        else if ((colorby==2))
        {// angle @ min
            t=atan2(zmin);
            t=t/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==3))
        {// in fraction
            index=incount/iter;
        }
        else
        {
                          // draw star
            Var_x=real(pixel);
            Var_y=imag(pixel);
            t=imag(log(pixel-Parm_center));
            if ((t<0))
            {
                t=t+twopi;
            }
            if (((t>t0)&&(t<=t1)))
            {
                err=Var_x*(y0-y2)+Var_y*(x2-x0)-x2*y0+x0*y2;
            }
            else if (((t>t1)&&(t<=t2)))
            {
                err=Var_x*(y1-y4)+Var_y*(x4-x1)-x4*y1+x1*y4;
            }
            else if (((t>t2)&&(t<=t3)))
            {
                err=Var_x*(y1-y3)+Var_y*(x3-x1)-x3*y1+x1*y3;
            }
            else if (((t>t3)&&(t<=t4)))
            {
                err=Var_x*(y2-y0)+Var_y*(x0-x2)-x0*y2+x2*y0;
            }
            else if (((t>t4)&&(t<=t5)))
            {
                err=Var_x*(y2-y4)+Var_y*(x4-x2)-x4*y2+x2*y4;
            }
            else if (((t>t5)&&(t<=t6)))
            {
                err=Var_x*(y3-y1)+Var_y*(x1-x3)-x1*y3+x3*y1;
            }
            else if (((t>t6)&&(t<=t7)))
            {
                err=Var_x*(y3-y0)+Var_y*(x0-x3)-x0*y3+x3*y0;
            }
            else if (((t>t7)&&(t<=t8)))
            {
                err=Var_x*(y4-y2)+Var_y*(x2-x4)-x2*y4+x4*y2;
            }
            else if (((t>t8)&&(t<=t9)))
            {
                err=Var_x*(y4-y1)+Var_y*(x1-x4)-x1*y4+x4*y1;
            }
            else
            {

                err=Var_x*(y3-y0)+Var_y*(x0-x3)-x0*y3+x3*y0;
            }
            err=cabs(err);
            index=cabs(log(cabs(err)));
        }
    }
    void description(void)
    {
        this.title="5 Point Star";

        Parm_center.caption="center of star";
        Parm_center.default=(0.0,0.0);


        r.caption="size of star";
        r.default=0.25;
        r.min=0.0;


        phiangle.caption="rotation, degrees";
        phiangle.default=0.0;
        phiangle.min=0.0;
        phiangle.max=36.0;
        phiangle.hint="only use angles between 0 and 36 degrees";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="minimum distance\niteration @ min\nangle @ min\nin fraction\ndraw star";

    }
}


comment { ; copyright Kerry Mitchell 20dec98

Conic Sections

Conic sections are sections of cones. Specifically, take a double-ended
cone, like 2 funnels placed tip-to-tip. Then, form the intersection of
the (double) cone with a plane. The intersection, usually one or two
curves, is a conic section. Conic sections can be a point, one line,
two lines, a parabola, an ellipse, a circle, or a hyperbola, depending
on the orientation of the plane relative to the cone. Analytically, they
can all be expressed by the same formula:

        Ax^2 + Bx + Cy^2 + Dy + Exy + F = 0,

where the parameters A through F determine the shape of the section, and
x and y are the 2 spatial coordinates. For example, the line y=x can be
represented as

        x - y = 0, or
        A = 0, B = 1, C = 0, D = -1, E = 0, F = 0.

A circle centered at (1,0) with a radius of 2 would have the equation

        (x - 1)^2 + y^2 = 4, or
        x^2 - 2x + y^2 - 3 = 0, giving
        A = 1, B = -2, C = 1, D = 0, E = 0, F = -3.

How can these shapes be used to render fractals?  Firstly, the "draw
section" setting simply draws the section determined by the six parameters
A - F, to give the user an idea of how the parameter choices affect the
section generated.

These coloring formulas record how the orbit interacts with the given
section.  In "Conic Sections", all six section parameters are input to
determine the curve.  Then, the image can be colored according to the
distance of the orbit's approach to the section or how often the orbit
landed inside the section.  Here's how the "colorby" parameters work in
"Conic Sections":

        "minimum distance":  closest approach to section
        "iteration @ min":  iteration number at minimum approach
        "angle @ min":  polar angle of z at minimum approach
        "maximum distance":  furthest approach from section
        "iteration @ max":  iteration number at maximum approach
        "angle @ max":  polar angle of z at maximum approach
        "in fraction":  how many times orbit was inside section, as a
fraction of the total number of iterations
        "in/out ratio": ration of number of times orbit was inside the
section to the number of times the orbit was outside of the section
        "lsb binary":  builds up a binary index, adding a bit to the
right for each iteration, "1" if the orbit was inside the section
and "0" if it was outside
        "msb binary":  builds up a binary index, adding a bit to the
left for each iteration, "1" if the orbit was inside the section
and "0" if it was outside
        "draw section":  simply draws the chosen conic section

While using all six section parameters is very powerful, it is less
than user-friendly.  Hence "Conic Lite".  Here, the user chooses the
type of section, and only enters the relevant geometric characteristics.
The choices are:

        "line":  enter angle of line to horizontal, and point through
which line passes
        "circle":  enter center and radius
        "ellipse"  enter semi-major and semi-minor axes and center
        "horiz. hyperbola":  enter semi-major and semi-minor axes and
center
        "vert. hyperbola":  enter semi-major and semi-minor axes and
center

The coloring options have also been reduced, to:

        "minimum distance":  closest approach to section
        "iteration @ min":  iteration number at minimum approach
        "angle @ min":  polar angle of z at minimum approach
        "draw section":  simply draws the chosen conic section

The "2 Sections" coloring method combines the results from 2 conic
sections into one.  The choices of sections are the same as in "Conic
Lite".  At each iteration, the distance from the first curve (call it
x) and the distance from the second curve (call it y) are computed.
From these, an overall distance metric is determined:

        "x^2 + y^2":  standard magnitude
        "|x| + |y|":  manhattan metric
        "|x*y|":  hyperbolic
        "|x-y|":  umm, call it "|x-y|"

When the chosen metric is at its smallest value, the x and y values
are combined into a complex number, zerrmin.  Also, the iterate z is
stored as zmin.

The coloring options ("color by" parameter) are:

        "2 curve min":  magnitude of zerrmin
        "2 curve angle":  polar angle of zerrmin
        "z angle @ min":  polar angle of zmin
        "z mag @ min":  magnitude of zmin
        "iteration # @ min":  iteration number at metric minimum
        "draw section":  simply draws the 2 sections

The standard forms of the conic sections are given below. With a bit of
algebraic twiddling, they can be transformed into the general form, for
use in the coloring formulas.

point (h,k) x=h, y=k; or (x-h)^2 + (y-k)^2 = 0 [circle of radius 0
centered at (h,k)].

vertical line through (h,k): x=h

non-vertical line with slope m, through (h,k): y-k = m*(x-h)

parabola, opening up or down, with vertex at (h,k): y-k = 4*p*(x-h)^2
[p determines width; +/up, -/down]

parabola, opening left or right, with vertex at (h,k): x-h = 4*p*(y-k)^2
[p determines width; +/right, -/left]

ellipse centered at (h,k), semimajor axis alpha, semiminor axis beta:
(x-h)^2/alpha^2 + (y-k)^2/beta^2 = 1

circle centered at (h,k), with radius r: (x-h)^2 + (y-k)^2 = r^2

hyperbola centered at (h,k), semimajor axis alpha, semiminor axis beta,
opening left/right: (x-h)^2/alpha^2 - (y-k)^2/beta^2 = 1

hyperbola centered at (h,k), semimajor axis alpha, semiminor axis beta,
opening up/down: (y-k)^2/beta^2 - (x-h)^2/alpha^2 = 1

coordinate rotation, from (u,v) to (x,y), through an angle theta:
u = x*cos(theta) + y*sin(theta)
v = -x*sin(theta) + y*cos(theta)
[rotating sections is how to generate non-zero E parameters]

}


conic-sections { // Kerry Mitchell 20dec98
//
// colors by relationship between orbit
// and fully-specified conic section
//
real cerr;
real cerrmin;
real cerrmax;
real Var_x;
real Var_y;
real err;
real t;
real total;
real lilbit;
real lilbase;
real bigbit;
real bigbase;
real bit;
int iter;
int incount;
int itermin;
int itermax;
complex zmin;
complex zmax;
parameter real aa;
parameter real bb;
parameter real cc;
parameter real dd;
parameter real ee;
parameter real ff;
parameter int colorby;

    void init(void)
    {
        cerr=0.0;
        cerrmin=1.0e12;
        cerrmax=0.0;
        Var_x=0.0;
        Var_y=0.0;
        err=0.0;
        t=0.0;
        total=0.0;
        lilbit=0.0;
        lilbase=0.5;
        bigbit=0.0;
        bigbase=1.0;
        bit=0.0;
        iter=0;
        incount=0;
        itermin=0;
        itermax=0;
        zmin=(0.0,0.0);
        zmax=(0.0,0.0);
    }
    void loop(void)
    {
        iter=iter+1;
        Var_x=real(z);
        Var_y=imag(z);
        err=(aa*Var_x+bb)*Var_x+(cc*Var_y+dd)*Var_y+ee*Var_x*Var_y+ff;
        bit=0.0;
        if ((err<0.0))
        {
            incount=incount+1;
            bit=1.0;
        }
        lilbit=lilbit+bit*lilbase;
        lilbase=lilbase*0.5;
        bigbit=bigbit+bit*bigbase;
        bigbase=bigbase*2.0;
        cerr=cabs(err);
        if ((cerr<cerrmin))
        {
            cerrmin=cerr;
            zmin=z;
            itermin=iter;
        }
        if ((cerr>cerrmax))
        {
            cerrmax=cerr;
            zmax=z;
            itermax=iter;
        }
        total=total+cerr;
    }
    void final(void)
    {
        if ((colorby==0))
        {// minimum distance
            index=0.5*cabs(log(cerrmin));
        }
        else if ((colorby==1))
        {// iteration @ min
            index=0.01*itermin;
        }
        else if ((colorby==2))
        {// angle @ min
            t=atan2(zmin)/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==3))
        {// maximum distance
            index=cerrmax;
        }
        else if ((colorby==4))
        {// iteration @ max
            index=0.01*itermax;
        }
        else if ((colorby==5))
        {// angle @ max
            t=atan2(zmax)/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==6))
        {// in fraction
            index=incount/iter;
        }
        else if ((colorby==7))
        {// in/out ratio
            index=incount/(iter-incount);
        }
        else if ((colorby==8))
        {// lsb binary
            index=lilbit;
        }
        else if ((colorby==9))
        {// msb binary
            index=0.5*bigbit/bigbase;
        }
        else
        {
                          // draw section
            Var_x=real(pixel);
            Var_y=imag(pixel);
            err=(aa*Var_x+bb)*Var_x+(cc*Var_y+dd)*Var_y+ee*Var_x*Var_y+ff;
            index=0.5*cabs(log(cabs(err)));
        }
    }
    void description(void)
    {
        this.title="Conic Sections";

        aa.caption="a";
        aa.default=0.0;
        aa.hint="coefficient of x^2 in conic section equation";


        bb.caption="b";
        bb.default=1.0;
        bb.hint="coefficient of x in conic section equation";


        cc.caption="c";
        cc.default=0.0;
        cc.hint="coefficient of y^2 in conic section equation";


        dd.caption="d";
        dd.default=-1.0;
        dd.hint="coefficient of y in conic section equation";


        ee.caption="e";
        ee.default=0.0;
        ee.hint="coefficient of x*y in conic section equation";


        ff.caption="f";
        ff.default=0.0;
        ff.hint="constant term in conic section equation";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="minimum distance\niteration @ min\nangle @ min\nmaximum distance\niteration @ max\nangle @ max\nin fraction\nin/out ratio\nlsb binary\nmsb binary\ndraw section";
        colorby.hint="see lkm-pub.ucl for explanations";

    }
}


conic-lite { // Kerry Mitchell 20dec98
//
// colors by relationship between orbit
// and conic section.  tastes great, less
// choices, less filling, easier to use.
//
real cerr;
real cerrmin;
real Var_x;
real Var_y;
real t;
int iter;
int itermin;
real h;
parameter complex Parm_center;
real k;
real aa;
real bb;
real cc;
real dd;
real ff;
real a2;
real b2;
complex zmin;
parameter int type;
parameter real radius;
parameter real xaxis;
parameter real yaxis;
parameter real theta;
parameter int colorby;
complex err;

    void init(void)
    {
        cerr=0.0;
        cerrmin=1.0e12;
        Var_x=0.0;
        Var_y=0.0;
        t=0.0;
        iter=0;
        itermin=0;
        h=real(Parm_center);
        k=imag(Parm_center);
        aa=0.0;
        bb=0.0;
        cc=0.0;
        dd=0.0;
        ff=0.0;
        a2=0.0;
        b2=0.0;
        zmin=(0.0,0.0);
//
// set up constants depending on chosen type of section
//
        if ((type==1))
        {// circle
            aa=1.0;
            bb=-2.0*h;
            cc=1.0;
            dd=-2.0*k;
            ff=sqr(h)+sqr(k)-sqr(radius);
        }
        else if ((type==2))
        {// ellipse
            a2=sqr(xaxis);
            b2=sqr(yaxis);
            aa=b2;
            bb=-2.0*h*b2;
            cc=a2;
            dd=-2.0*k*a2;
            ff=b2*sqr(h)+a2*sqr(k)-a2*b2;
        }
        else if ((type==3))
        {// horizontal hyperbola
            a2=sqr(xaxis);
            b2=sqr(yaxis);
            aa=b2;
            bb=-2.0*h*b2;
            cc=-a2;
            dd=2.0*k*a2;
            ff=b2*sqr(h)-a2*sqr(k)-a2*b2;
        }
        else if ((type==4))
        {// vertical hyperbola
            a2=sqr(xaxis);
            b2=sqr(yaxis);
            aa=-b2;
            bb=2.0*h*b2;
            cc=a2;
            dd=-2.0*k*a2;
            ff=-b2*sqr(h)+a2*sqr(k)-a2*b2;
        }
        else
        {
                          // line
            t=theta/180*pi;
            aa=0.0;
            bb=sin(t);
            cc=0.0;
            dd=-cos(t);
            ff=-(bb*h+dd*k);
        }
    }
    void loop(void)
    {
        iter=iter+1;
        Var_x=real(z);
        Var_y=imag(z);
        cerr=cabs((aa*Var_x+bb)*Var_x+(cc*Var_y+dd)*Var_y+ff);
        if ((cerr<cerrmin))
        {
            cerrmin=cerr;
            zmin=z;
            itermin=iter;
        }
    }
    void final(void)
    {
        if ((colorby==0))
        {// minimum distance
            index=0.5*cabs(log(cerrmin));
        }
        else if ((colorby==1))
        {// iteration @ min
            index=0.01*itermin;
        }
        else if ((colorby==2))
        {// angle @ min
            t=atan2(zmin)/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else
        {
                          // draw section
            Var_x=real(pixel);
            Var_y=imag(pixel);
            err=(aa*Var_x+bb)*Var_x+(cc*Var_y+dd)*Var_y+ff;
            index=0.5*cabs(log(cabs(err)));
        }
    }
    void description(void)
    {
        this.title="Conic Lite";

        type.caption="section type";
        type.default=1;
        type.enum="line\ncircle\nellipse\nhoriz. hyperbola\nvert. hyperbola";


        Parm_center.caption="center";
        Parm_center.default=(0.0,0.0);


        theta.caption="angle of line";
        theta.default=45.0;
        theta.hint="angle to horizontal, degrees";


        radius.caption="radius";
        radius.default=1.0;
        radius.min=0.0;


        xaxis.caption="x axis";
        xaxis.default=1.5;
        xaxis.min=0.0;
        xaxis.hint="for ellipses and hyperbolas";


        yaxis.caption="y axis";
        yaxis.default=0.75;
        yaxis.min=0.0;
        yaxis.hint="for ellipses and hyperbolas";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="minimum distance\niteration @ min\nangle @ min\ndraw section";

    }
}


2-sections { // Kerry Mitchell 20dec98
//
// colors by the relationship between
// the orbit and 2 conic sections
//
real Var_x;
real Var_y;
real t;
real err1;
real err2;
real cerr;
real cerrmin;
int iter;
int itermin;
real h1;
parameter complex center1;
real k1;
real a1;
real b1;
real c1;
real d1;
real f1;
real h2;
parameter complex center2;
real k2;
real a2;
real b2;
real c2;
real d2;
real f2;
real tempa;
real tempb;
complex zmin;
complex zerrmin;
complex temperr;
parameter int type1;
parameter real radius1;
parameter real xaxis1;
parameter real yaxis1;
parameter real theta1;
parameter int type2;
parameter real radius2;
parameter real xaxis2;
parameter real yaxis2;
parameter real theta2;
parameter int metrictype;
parameter int colorby;

    void init(void)
    {
        Var_x=0.0;
        Var_y=0.0;
        t=0.0;
        err1=0.0;
        err2=0.0;
        cerr=0.0;
        cerrmin=1.0e12;
        iter=0;
        itermin=0;
        h1=real(center1);
        k1=imag(center1);
        a1=0.0;
        b1=0.0;
        c1=0.0;
        d1=0.0;
        f1=0.0;
        h2=real(center2);
        k2=imag(center2);
        a2=0.0;
        b2=0.0;
        c2=0.0;
        d2=0.0;
        f2=0.0;
        tempa=0.0;
        tempb=0.0;
        zmin=(0.0,0.0);
        zerrmin=(0.0,0.0);
        temperr=(0.0,0.0);
//
// conic section parameters for first curve
//
        if ((type1==1))
        {// circle
            a1=1.0;
            b1=-2.0*h1;
            c1=1.0;
            d1=-2.0*k1;
            f1=sqr(h1)+sqr(k1)-sqr(radius1);
        }
        else if ((type1==2))
        {// ellipse
            tempa=sqr(xaxis1);
            tempb=sqr(yaxis1);
            a1=tempb;
            b1=-2.0*h1*tempb;
            c1=tempa;
            d1=-2.0*k1*tempa;
            f1=tempb*sqr(h1)+tempa*sqr(k1)-tempa*tempb;
        }
        else if ((type1==3))
        {// horizontal hyperbola
            tempa=sqr(xaxis1);
            tempb=sqr(yaxis1);
            a1=tempb;
            b1=-2.0*h1*tempb;
            c1=-tempa;
            d1=2.0*k1*tempa;
            f1=tempb*sqr(h1)-tempa*sqr(k1)-tempa*tempb;
        }
        else if ((type1==4))
        {// vertical hyperbola
            tempa=sqr(xaxis1);
            tempb=sqr(yaxis1);
            a1=-tempb;
            b1=2.0*h1*tempb;
            c1=tempa;
            d1=-2.0*k1*tempa;
            f1=-tempb*sqr(h1)+tempa*sqr(k1)-tempa*tempb;
        }
        else
        {
                           // line
            t=theta1/180*pi;
            a1=0.0;
            b1=sin(t);
            c1=0.0;
            d1=-cos(t);
            f1=-(b1*h1+d1*k1);
        }
//
// conic section parameters for second curve
//
        if ((type2==1))
        {// circle
            a2=1.0;
            b2=-2.0*h2;
            c2=1.0;
            d2=-2.0*k2;
            f2=sqr(h2)+sqr(k2)-sqr(radius2);
        }
        else if ((type2==2))
        {// ellipse
            tempa=sqr(xaxis2);
            tempb=sqr(yaxis2);
            a2=tempb;
            b2=-2.0*h2*tempb;
            c2=tempa;
            d2=-2.0*k2*tempa;
            f2=tempb*sqr(h2)+tempa*sqr(k2)-tempa*tempb;
        }
        else if ((type2==3))
        {// horizontal hyperbola
            tempa=sqr(xaxis2);
            tempb=sqr(yaxis2);
            a2=tempb;
            b2=-2.0*h2*tempb;
            c2=-tempa;
            d2=2.0*k2*tempa;
            f2=tempb*sqr(h2)-tempa*sqr(k2)-tempa*tempb;
        }
        else if ((type2==4))
        {// vertical hyperbola
            tempa=sqr(xaxis2);
            tempb=sqr(yaxis2);
            a2=-tempb;
            b2=2.0*h2*tempb;
            c2=tempa;
            d2=-2.0*k2*tempa;
            f2=-tempb*sqr(h2)+tempa*sqr(k2)-tempa*tempb;
        }
        else
        {
                           // line
            t=theta2/180*pi;
            a2=0.0;
            b2=sin(t);
            c2=0.0;
            d2=-cos(t);
            f2=-(b2*h2+d2*k2);
        }
    }
    void loop(void)
    {
        iter=iter+1;
//
// see how far iterate is from curves
//
        Var_x=real(z);
        Var_y=imag(z);
        err1=(a1*Var_x+b1)*Var_x+(c1*Var_y+d1)*Var_y+f1;
        err2=(a2*Var_x+b2)*Var_x+(c2*Var_y+d2)*Var_y+f2;
        temperr=err1+flip(err2);
//
// create distance metric and compare to current minimum
//
        if ((metrictype==0))
        {// magnitude
            cerr=cabs(temperr);
        }
        else if ((metrictype==1))
        {// manhattan metric
            cerr=cabs(err1)+cabs(err2);
        }
        else if ((metrictype==2))
        {// hyperbolic
            cerr=cabs(err1*err2);
        }
        else
        {
                          // |x-y|
            cerr=cabs(err1-err2);
        }
        if ((cerr<cerrmin))
        {
            cerrmin=cerr;
            itermin=iter;
            zmin=z;
            zerrmin=temperr;
        }
    }
    void final(void)
    {
//
// set index according to colorby parameter
//
        if ((colorby==0))
        {// minimum magnitude
            index=0.5*cabs(log(cerrmin));
        }
        else if ((colorby==1))
        {// angle @ min
            t=atan2(zerrmin)/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==2))
        {// z angle @ min
            t=atan2(zmin)/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==3))
        {// z mag @ min
            index=cabs(zmin);
        }
        else if ((colorby==4))
        {// iteration # @ min
            index=0.01*itermin;
        }
        else
        {
                           // draw sections
            Var_x=real(pixel);
            Var_y=imag(pixel);
            err1=cabs((a1*Var_x+b1)*Var_x+(c1*Var_y+d1)*Var_y+f1);
            err2=cabs((a2*Var_x+b2)*Var_x+(c2*Var_y+d2)*Var_y+f2);
            if ((err1<err2))
            {
                cerr=err1;
            }
            else
            {

                cerr=err2;
            }
            index=0.5*cabs(log(cerr));
        }
    }
    void description(void)
    {
        this.title="2 Sections";

        type1.caption="section 1 type";
        type1.default=0;
        type1.enum="line\ncircle\nellipse\nhoriz. hyperbola\nvert. hyperbola";


        center1.caption="section 1 center";
        center1.default=(0.0,0.0);


        theta1.caption="section 1 angle";
        theta1.default=45.0;
        theta1.hint="angle of line to horizontal, degrees";


        radius1.caption="section 1 radius";
        radius1.default=1.0;
        radius1.min=0.0;


        xaxis1.caption="section 1 x axis";
        xaxis1.default=1.0;
        xaxis1.min=0.0;
        xaxis1.hint="for ellipses and hyperbolas";


        yaxis1.caption="section 1 y axis";
        yaxis1.default=1.0;
        yaxis1.min=0.0;
        yaxis1.hint="for ellipses and hyperbolas";


        type2.caption="section 2 type";
        type2.default=1;
        type2.enum="line\ncircle\nellipse\nhoriz. hyperbola\nvert. hyperbola";


        center2.caption="section 2 center";
        center2.default=(0.0,0.0);


        theta2.caption="section 2 angle";
        theta2.default=0.0;
        theta2.hint="angle of line to horizontal, degrees";


        radius2.caption="section 2 radius";
        radius2.default=1.0;
        radius2.min=0.0;


        xaxis2.caption="section 2 x axis";
        xaxis2.default=1.0;
        xaxis2.min=0.0;
        xaxis2.hint="for ellipses and hyperbolas";


        yaxis2.caption="section 2 y axis";
        yaxis2.default=1.0;
        yaxis2.min=0.0;
        yaxis2.hint="for ellipses and hyperbolas";


        metrictype.caption="metric";
        metrictype.default=0;
        metrictype.enum="x^2 + y^2\n|x| + |y|\n|x*y|\n|x-y|";
        metrictype.hint="how the minimum iterate is determined";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="2 curve min\n2 curve angle\nz angle @ min\nz mag @ min\niteration # @ min\ndraw section";
        colorby.hint="see lkm-pub.ucl file for explanations";

    }
}


conic-lite-old { // Kerry Mitchell
//
// *** this has been superceded by "conic lite",
// *** but is offered for compatibility
//
// colors by relationship between orbit
// and conic section.  tastes great, less
// choices, less filling, easier to use.
//
real cerr;
real cerrmin;
real Var_x;
real Var_y;
real t;
int iter;
int itermin;
real h;
parameter complex Parm_center;
real k;
real aa;
real bb;
real cc;
real dd;
real ff;
real a2;
real b2;
complex zmin;
parameter int type;
parameter real radius;
parameter real xaxis;
parameter real yaxis;
parameter real theta;
parameter int colorby;
complex err;

    void init(void)
    {
        cerr=0.0;
        cerrmin=1.0e12;
        Var_x=0.0;
        Var_y=0.0;
        t=0.0;
        iter=0;
        itermin=0;
        h=real(Parm_center);
        k=imag(Parm_center);
        aa=0.0;
        bb=0.0;
        cc=0.0;
        dd=0.0;
        ff=0.0;
        a2=0.0;
        b2=0.0;
        zmin=(0.0,0.0);
//
// set up constants depending on chosen type of section
//
        if ((type==1))
        {// circle
            aa=1.0;
            bb=-2.0*h;
            cc=1.0;
            dd=-2.0*k;
            ff=sqr(h)+sqr(k)-sqr(radius);
        }
        else if ((type==2))
        {// ellipse
            a2=sqr(xaxis);
            b2=sqr(yaxis);
            aa=b2;
            bb=-2.0*h*b2;
            cc=a2;
            dd=-2.0*k*a2;
            ff=b2*sqr(h)+a2*sqr(k)-a2*b2;
        }
        else if ((type==3))
        {// horizontal hyperbola
            a2=sqr(xaxis);
            b2=sqr(yaxis);
            aa=b2;
            bb=-2.0*h*b2;
            cc=-a2;
            dd=2.0*k*a2;
            ff=b2*sqr(h)-a2*sqr(k)-a2*b2;
        }
        else if ((type==4))
        {// vertical hyperbola
            a2=sqr(xaxis);
            b2=sqr(yaxis);
            aa=-b2;
            bb=2.0*h*b2;
            cc=a2;
            dd=-2.0*k*a2;
            ff=-b2*sqr(h)+a2*sqr(k)-a2*b2;
        }
        else
        {
                          // line
            t=theta/180*pi;
            aa=0.0;
            bb=sin(t);
            cc=0.0;
            dd=-cos(t);
            ff=-(bb*h+dd*k);
        }
    }
    void loop(void)
    {
        iter=iter+1;
        Var_x=real(z);
        Var_y=imag(z);
        cerr=cabs((aa*Var_x+bb)*Var_x+(cc*Var_y+dd)*Var_y+ff);
        if ((cerr<cerrmin))
        {
            cerrmin=cerr;
            zmin=z;
            itermin=iter;
        }
    }
    void final(void)
    {
        if ((colorby==0))
        {// minimum distance
            index=cerrmin;
        }
        else if ((colorby==1))
        {// iteration @ min
            index=0.01*itermin;
        }
        else if ((colorby==2))
        {// angle @ min
            t=atan2(zmin)/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else
        {
                          // draw section
            Var_x=real(pixel);
            Var_y=imag(pixel);
            err=(aa*Var_x+bb)*Var_x+(cc*Var_y+dd)*Var_y+ff;
            index=0.5*cabs(log(cabs(err)));
        }
    }
    void description(void)
    {
        this.title="Old Conic Lite";

        type.caption="section type";
        type.default=1;
        type.enum="line\ncircle\nellipse\nhoriz. hyperbola\nvert. hyperbola";


        Parm_center.caption="center";
        Parm_center.default=(0.0,0.0);


        theta.caption="angle of line";
        theta.default=45.0;
        theta.hint="angle to horizontal, degrees";


        radius.caption="radius";
        radius.default=1.0;
        radius.min=0.0;


        xaxis.caption="x axis";
        xaxis.default=1.5;
        xaxis.min=0.0;
        xaxis.hint="for ellipses and hyperbolas";


        yaxis.caption="y axis";
        yaxis.default=0.75;
        yaxis.min=0.0;
        yaxis.hint="for ellipses and hyperbolas";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="minimum distance\niteration @ min\nangle @ min\ndraw section";

    }
}


comment { ; narrative copyright Kerry Mitchell 06feb99

The Bubble Method

The bubble method is an extension of Fractint's bof60 scheme.  In
bof60, the interior of the fractal is colored by how close the iterate
comes to the origin.  In the bubble method, a specific value is set as
the threshold.  At each duration, the magnitude of the iterate is compared
to the threshold.  If the current magnitude is smaller, it becomes the new
threshold level.  The effect is to cover the fractal with "bubbles",
circles the radius of the threshold.  Below a certain threshold value,
the image is a dust of small, disconnected bubbles.  At a particular
theshold value, which varies with the parameter c, the circles all
touch.  Beyond this, the circles squish into each other, like mounds
of soap bubbles.

When a new threshold level is set, the iterate and the iteration number
are stored.  At bailout, the final new values are available for use in
the coloring.  The index value is a weighted average of the final threshold,
the iteration number when the last threshold was set, and the polar angle
of the last threshold's iterate.  The weights of each contribution can be
individually set using the "mag. weight", "iteration weight", and "angle
weight" parameters.  Using only the magnitude to color results in a series
of concentric circles, and using only the angle results in a radial color
pattern.  Combining equal amounts of magnitude and angle gives a swirl
effect.  Changing the angle weight to a negative value will reverse the
direction of the swirl.  A judicious choice of color density will eliminate
the branch cuts that can be seen with the swirls.

}


bubbles { // Kerry Mitchell 06feb98
//
// colors by the "bubbles" variation of
// Fractint's bof60 and bof61 methods
//
real weighttotal;
parameter real weightr;
parameter real weightt;
parameter real weighti;
real wr;
real wt;
real wi;
real twopi;
real rz;
real rmin;
parameter real threshold;
real tmin;
complex zmin;
bool minset;
real iter;
real itermin;

    void init(void)
    {
        weighttotal=weightr+cabs(weightt)+weighti;
        if ((weighttotal==0.0))
        {
            weighttotal=1.0;
        }
        wr=weightr/weighttotal;
        wt=cabs(weightt)/weighttotal;
        wi=weighti/weighttotal;
        twopi=2.0*pi;
        rz=0.0;
        rmin=threshold;
        tmin=0.0;
        zmin=(0.0,0.0);
        minset=false;
        iter=0.0;
        itermin=0.0;
    }
    void loop(void)
    {
        iter=iter+1;
        rz=cabs(z);
        if ((rz<rmin))
        {
            minset=true;
            rmin=rz;
            zmin=z;
            itermin=iter;
        }
    }
    void final(void)
    {
        if ((minset==true))
        {
            rmin=rmin/threshold;
            tmin=atan2(zmin);
            if ((tmin<0.0))
            {
                tmin=tmin+twopi;
            }
            tmin=tmin/twopi;
            if ((weightt<0.0))
            {
                tmin=1.0-tmin;
            }
            iter=itermin/iter;
            index=cabs(wr*rmin+wt*tmin+wi*iter);
        }
        else
        {

            index=0.0;
            solid=true;
        }
    }
    void description(void)
    {
        this.title="Bubbles";

        threshold.caption="bubble radius";
        threshold.default=1.0;
        threshold.min=0.0;


        weightr.caption="mag. weight";
        weightr.default=1.0;
        weightr.hint="relative weight of magnitude component";


        weightt.caption="angle weight";
        weightt.default=1.0;
        weightt.hint="relative weight of polar angle";


        weighti.caption="iteration weight";
        weighti.default=1.0;
        weighti.hint="relative weight of iteration count";

    }
}


comment { ; narrative copyright Kerry Mitchell 21feb99

Range

In the range coloring schemes, the pixels are only colored when the
iterate (or a component of it) falls within a specified range.
(Otherwise, the #solid flag is set to true.)  The range is given in
terms of the center and width.  For example, setting the center to
1.0 and the range to 0.1 would yield a range of 0.95 to 1.05.

The simplest case is shown in the "Range Lite" formula.  Here, the
orbit is monitored until the magnitude of the iterate falls within
a specified range.  When this occurs, that iterate is saved.  Upon
bailout (either escape or maximum iterations), the point can be colored
according to either the magnitude or the polar angle of the iterate.
Either the first iterate to fall into the range or the last iterate
can be used.  The difference is in the overlapping of the rings of
color.  Using the first iterate gives the appearance of the "earlier"
rings (from lower iteration numbers) stacked on top of the "later"
rings.  Using the laster iterate reverses the order of the stacking.
Depending on the range chosen, the results can be circular rings,
pinched loops (figure "8"), or open loops.  The shapes are definitely
non-fractal, but are assembled in such a way to recover the underlying
fractal structure.

This method opens itself up to many variations, which can be explored
using the full "Range" coloring.  In essence, some variable is monitored
to see when it falls in the specified range.  When it does, another
variable is used for coloring.  There are several options for the
"range variable", the variable that is checked against the range:
the magnitude of z, its real or imaginary part, and its polar angle.
Plus, any of the standard functions can be used on z, and the real
or imaginary parts used.  A final choice, "heart", is a combination
of the magnitude and polar angle, that covers the image with heart
shapes.  The same choices are available for the "coloring variable",
the variable used to set #index.

Since #index takes on values from 0.0 to 1.0, the coloring variable
may need adjusting before being used as the #index.  If the coloring
variable is the same as the range variable, then the range parameters
are used to scale the variable into #index.  Otherwise, the coloring
variable can be easily scaled into the #index variable.  Specifically,

#index = slope * scale_function(coloring_variable) + offset.

For example, if real(z) is used as the coloring variable, then the
tanh() function can scale real(z) to -1.0 to 1.0.  Then, using slope
of 0.5 and an offset of 0.5 will give the final range of 0.0 to 1.0
for the #index.

Finally, the "color by" parameter determines which iterate or iterates
are used in the coloring.  Using the "first" or "last" settings will
vary the apparant stacking of the color rings.  The "average" setting
will cause all the coloring variables for which the range variables
fall into the range to be averaged together, and that mean used for
the coloring.  If the range width is set fairly small, then there won't
be much difference between "average" and "first" or "last".  But there
can be a quite significant change if the range width is large.  The
last choice, "fraction in range", simply uses the ratio of number of
times the range variable falls within the range, to the total number of
iterations for that pixel.

}


range-lite { // Kerry Mitchell 21feb99
//
// Simplified version of "Range" coloring.
// Colors by magnitude or angle when magnitude
// of z falls in prescribed range.
//
int rangeiter;
real rangemin;
parameter real rangecenter;
parameter real rangewidth;
real rangemax;
real rz;
real twopi;
real r1;
real rn;
real t1;
real tn;
parameter int colorby;

    void init(void)
    {
        rangeiter=0;
        rangemin=rangecenter-0.5*rangewidth;
        rangemax=rangemin+rangewidth;
        rz=0.0;
        twopi=2.0*pi;
        r1=0.0;
        rn=0.0;
        t1=0.0;
        tn=0.0;
    }
    void loop(void)
    {
        rz=cabs(z);
        if (((rz>=rangemin)&&(rz<=rangemax)))
        {
            rangeiter=rangeiter+1;
            rn=rz;
            tn=atan2(z);
            if ((rangeiter==1))
            {
                r1=rn;
                t1=tn;
            }
        }
    }
    void final(void)
    {
        if ((rangeiter==0))
        {
            solid=true;
        }
        else
        {

            if ((colorby==0))
            {// first mag
                r1=(r1-rangemin)/(rangemax-rangemin);
                index=r1;
            }
            else if ((colorby==1))
            {// last mag
                rn=(rn-rangemin)/(rangemax-rangemin);
                index=rn;
            }
            else if ((colorby==2))
            {// first angle
                if ((t1<0.0))
                {
                    t1=t1+twopi;
                }
                t1=t1/twopi;
                index=t1;
            }
            else if ((colorby==3))
            {// last angle
                if ((tn<0.0))
                {
                    tn=tn+twopi;
                }
                tn=tn/twopi;
                index=tn;
            }
        }
    }
    void description(void)
    {
        this.title="Range Lite";

        rangecenter.caption="range center";
        rangecenter.default=1.0;


        rangewidth.caption="range width";
        rangewidth.default=0.1;


        colorby.caption="color by";
        colorby.default=3;
        colorby.enum="first mag\nlast mag\nfirst angle\nlast angle";

    }
}


range { // Kerry Mitchell 21feb99
//
// Colors by 2nd number when 1st
// number falls in prescribed range.
// Entirely too many choices for 1st
// and 2nd numbers.  :-)
//
int rangeiter;
real rangemin;
parameter real rangecenter;
parameter real rangewidth;
real rangemax;
real rvar;
real cvar;
real cvar1;
real cvarn;
real cvarave;
real twopi;
parameter int rangevar;
parameter int colorvar;
parameter real slope;
parameter real offset;
parameter int colorby;

    void init(void)
    {
        rangeiter=0;
        rangemin=rangecenter-0.5*rangewidth;
        rangemax=rangemin+rangewidth;
        rvar=0.0;
        cvar=0.0;
        cvar1=0.0;
        cvarn=0.0;
        cvarave=0.0;
        twopi=2.0*pi;
    }
    void loop(void)
    {
        if ((rangevar==0))
        {
            rvar=cabs(z);
        }
        else if ((rangevar==1))
        {
            rvar=real(z);
        }
        else if ((rangevar==2))
        {
            rvar=imag(z);
        }
        else if ((rangevar==3))
        {
            rvar=atan2(z);
        }
        else if ((rangevar==4))
        {
            rvar=real(rangefunc(z));
        }
        else if ((rangevar==5))
        {
            rvar=imag(rangefunc(z));
        }
        else if ((rangevar==6))
        {
            rvar=cabs(atan2(z))/pi-cabs(z);
        }
        if (((rvar>=rangemin)&&(rvar<=rangemax)))
        {
            rangeiter=rangeiter+1;
            if ((colorvar==0))
            {
                cvar=cabs(z);
            }
            else if ((colorvar==1))
            {
                cvar=real(z);
            }
            else if ((colorvar==2))
            {
                cvar=imag(z);
            }
            else if ((colorvar==3))
            {
                cvar=atan2(z);
                if ((cvar<0.0))
                {
                    cvar=cvar+twopi;
                }
                cvar=cvar/twopi;
            }
            else if ((colorvar==4))
            {
                cvar=real(colorfunc(z));
            }
            else if ((colorvar==5))
            {
                cvar=imag(colorfunc(z));
            }
            else if ((colorvar==6))
            {
                cvar=cabs(atan2(z))/pi-cabs(z);
            }
            if ((colorvar==rangevar))
            {
                cvarn=(cvar-rangemin)/(rangemax-rangemin);
            }
            else
            {

                cvarn=slope*real(scalefunc(cvar))+offset;
            }
            if ((rangeiter==1))
            {
                cvar1=cvarn;
            }
            cvarave=cvarave+cvarn;
        }
    }
    void final(void)
    {
        if ((rangeiter==0))
        {
            solid=true;
        }
        else
        {

            if ((colorby==0))
            {// first
                index=cvar1;
            }
            else if ((colorby==1))
            {// average
                index=cvarave/rangeiter;
            }
            else if ((colorby==2))
            {// last
                index=cvarn;
            }
            else if ((colorby==3))
            {// % in range
                index=rangeiter/numiter;
            }
        }
    }
    void description(void)
    {
        this.title="Range";

        rangecenter.caption="range center";
        rangecenter.default=1.0;


        rangewidth.caption="range width";
        rangewidth.default=0.1;


        rangevar.caption="range variable";
        rangevar.default=0;
        rangevar.enum="magnitude\nreal(z)\nimag(z)\npolar angle\nreal f(z)\nimag f(z)\nheart";


        colorvar.caption="coloring variable";
        colorvar.default=0;
        colorvar.enum="magnitude\nreal(z)\nimag(z)\npolar angle\nreal f(z)\nimag f(z)\nheart";


        slope.caption="slope";
        slope.default=1.0;
        slope.hint="for adjusting coloring variable into 0.0 - 1.0 range for #index; index=slope*fn(var))+offset";


        offset.caption="offset";
        offset.default=0.0;
        offset.hint="for adjusting coloring variable into 0.0 - 1.0 range for #index; index=slope*fn(var))+offset";


        colorby.caption="color by";
        colorby.default=2;
        colorby.enum="first\naverage\nlast\nfraction in range";


        scalefunc.caption="scaling function";
        scalefunc.default="ident" ;
        scalefunc.hint="for adjusting coloring variable into 0.0 - 1.0 range for #index; index=slope*fn(var))+offset";


        rangefunc.caption="range function";
        rangefunc.default="sqr" ;
        rangefunc.hint="for creating new range variables";


        colorfunc.caption="coloring function";
        colorfunc.default="sqr" ;
        colorfunc.hint="for creating new coloring variables";

    }
}


passthru { // Kerry Mitchell 02jul99
//
// Special purpose coloring, designed to be used with the
// "Quadrilateral Color" transform.
//
// Reads the x and y pixel flags, and determines an #index
// value from them.
//
int xflag;
int yflag;
parameter real xytype;
parameter real n;
parameter real m;

    void init(void)
    {
        xflag=round(real(pixel)/20);
        yflag=round(imag(pixel)/20);
    }
    void loop(void)
    {
    }
    void final(void)
    {
        if ((xytype==1))
        {
            index=yflag/n;
        }
        else if ((xytype==2))
        {
            index=(m*xflag+yflag)/n;
        }
        else if ((xytype==3))
        {
            index=(m*yflag+xflag)/n;
        }
        else
        {

            index=xflag/n;
        }
    }
    void description(void)
    {
        this.title="Passthru";

        xytype.caption="flag type";
        xytype.enum="x\ny\nm*x+y\nm*y+x";
        xytype.default=0;


        m.caption="m: flag factor";
        m.default=1.0;
        m.hint="use with 'm*x+y' and 'm*y+x' flag types";


        n.caption="flag scale";
        n.default=1.0;
        n.hint="#index = flag/scale";

    }
}


distance-point { // Kerry Mitchell 06feb00
//
// Colors by the distance from the orbit to a specified
// point.  Offers 7 different ways to determine the distance,
// and colors by the closest approach, the furthest approach,
// or combinations thereof.
//
// See the helpfile for more information.
//
real Var_x;
real Var_y;
real d;
real dfirst;
real dlast;
real dmin;
real dmax;
real dsum;
real dproduct;
real temp;
real temp2;
real twooverpi;
real pio2;
complex tempc;
int iter;
real x0e;
parameter complex point;
parameter real r;
real cosy0e;
real siny0e;
real x0h;
real y0h;
real r0h;
parameter int metric;
parameter real power;
parameter int colorby;

    void init(void)
    {
        Var_x=0.0;
        Var_y=0.0;
        d=0.0;
        dfirst=0.0;
        dlast=0.0;
        dmin=1e20;
        dmax=0.0;
        dsum=0.0;
        dproduct=1.0;
        temp=0.0;
        temp2=0.0;
        twooverpi=2.0/pi;
        pio2=0.5*pi;
        tempc=(0,0);
        iter=0;
//
// reference point for elliptical geometry:
//   x becomes theta (longitude)
//   y becomes phi (latitude)
//   (spherical coordinates)
//
        x0e=real(point)/r*pi;
        temp=imag(point)/r*pio2;
        cosy0e=cos(temp);
        siny0e=sin(temp);
//
// reference point for hyperbolic geometry
//
        x0h=real(point)/r;
        y0h=imag(point)/r;
        r0h=sqrt(1-sqr(x0h)-sqr(y0h));

    }
    void loop(void)
    {
        iter=iter+1;
//
// determine distance based on metric choice
//
        if ((metric==1))
        {// elliptic geometry
//
// spherical coordinates:
//   x becomes longitude, y becomes latitude
//
            Var_x=real(z)/r*pi;
            Var_y=imag(z)/r*pio2;
            temp=cos(Var_y)*cosy0e*cos(Var_x-x0e)+sin(Var_y)*siny0e;
//
// use atan to find angle between reference and field points
//
            temp2=sqrt(1.0-temp*temp);
            if ((temp>0.0))
            {
                d=atan(temp2/temp);
            }
            else if ((temp<0.0))
            {
                d=atan(temp2/temp)+pi;
            }
            else
            {

                d=pio2;
            }
//
// distance = radius * angle
//
            d=r*d;
        }
        else if ((metric==2))
        {// hyperbolic geometry
            Var_x=real(z)/r;
            Var_y=imag(z)/r;
            temp=(1.0-Var_x*x0h-Var_y*y0h);
            temp=temp/(r0h*sqrt(1-sqr(Var_x)-sqr(Var_y)));
//
//  distance = radius * angle
//
            d=r*acosh(temp);
        }
        else if ((metric==3))
        {// minimum
//
// uses minimum of delta-x and delta-y
//   not a true distance, but included anyway
//
            Var_x=abs(real(z-point));
            Var_y=abs(imag(z-point));
            if ((Var_x<Var_y))
            {
                d=Var_x;
            }
            else
            {

                d=Var_y;
            }
        }
        else if ((metric==4))
        {// maximum
//
// uses maximum of delta-x and delta-y
//
            Var_x=abs(real(z-point));
            Var_y=abs(imag(z-point));
            if ((Var_x>Var_y))
            {
                d=Var_x;
            }
            else
            {

                d=Var_y;
            }
        }
        else if ((metric==5))
        {// sum
//
// uses sum of delta-x and delta-y
//
            Var_x=abs(real(z-point));
            Var_y=abs(imag(z-point));
            d=Var_x+Var_y;
        }
        else if ((metric==6))
        {// product
//
// uses product of delta-x and delta-y
//   not a true distance, but included anyway
//
            Var_x=abs(real(z-point));
            Var_y=abs(imag(z-point));
            d=sqrt(Var_x*Var_y);
        }
        else
        {
                       // euclidean
//
// generalized pythagorean theorem
//
            Var_x=abs(real(z-point));
            Var_y=abs(imag(z-point));
            d=(Var_x^power+Var_y^power)^(1/power);
        }
//
// fill variables for different coloring types
//
        if ((iter==1))
        {
            dfirst=d;
        }
        if ((d<dmin))
        {
            dmin=d;
        }
        if ((d>dmax))
        {
            dmax=d;
        }
        dsum=dsum+d;
        dproduct=dproduct*d;

    }
    void final(void)
    {
        dlast=d;
        if ((colorby==0))
        {// first
            index=dfirst;
        }
        else if ((colorby==1))
        {// last
            index=dlast;
        }
        else if ((colorby==2))
        {// first/last angle
            tempc=dfirst+flip(dlast);
            temp=atan2(tempc);
            index=temp*twooverpi;
        }
        else if ((colorby==3))
        {// minimum
            index=dmin;
        }
        else if ((colorby==4))
        {// maximum
            index=dmax;
        }
        else if ((colorby==5))
        {// min/max angle
            tempc=dmin+flip(dmax);
            temp=atan2(tempc);
            index=temp*twooverpi;
        }
        else if ((colorby==6))
        {// arithmetic mean
            index=dsum/iter;
        }
        else if ((colorby==7))
        {// geometric mean
            temp=log(dproduct)/iter;
            index=exp(temp);
        }
        else if ((colorby==8))
        {// amean/gmean angle
            temp=dsum/iter;
            temp2=log(dproduct)/iter;
            temp2=exp(temp2);
            tempc=temp+flip(temp2);
            temp=atan2(tempc);
            index=temp*twooverpi;
        }
    }
    void description(void)
    {
        this.title="Distance to a Point";
        this.helpfile="lkm-help\lkm-distance.html";

        point.caption="reference point";
        point.default=(0,0);
        point.hint="Point from which distance is measured.";


        metric.caption="distance metric";
        metric.default=0;
        metric.enum="Euclidean\nelliptic\nhyperbolic\nminimum\nmaximum\nsum\nproduct";
        metric.hint="What type of formula is used to figure the distance.";


        colorby.caption="color by:";
        colorby.default=3;
        colorby.enum="1st distance\nlast distance\n1st/last combo\nminimum\nmaximum\nmin/max combo\narithmetic mean\ngeometric mean\na/g combo";
        colorby.hint="Combinations plot the angle using the first variable as the x, and the 2nd as the y.";


        power.caption="power";
        power.default=2.0;
        power.min=1.0;
        power.hint="Exponent used with 'power' metric.  Must be at least 1.  Use 2 for standard Pythagorean distance.";


        r.caption="radius";
        r.default=4.0;
        r.hint="Used with 'elliptic' & 'hyperbolic' metrics.  Should be larger than the bailout.";

    }
}


pythagorean-triple { // Kerry Mitchell 06feb00
//
// Colors by the approach of the orbit to a Pythagorean
// triple:  when x, y, and r are all integers.  Uses
// several different ways to find the closest integer
// to determine the distance.  Colors according to
// the closest approach, the furthest approach, or
// combinations thereof.
//
// See the helpfile for more information.
//
real rx;
real ry;
real rr;
real r;
real rmin;
real rmax;
real rsum;
real rproduct;
real temp;
real temp2;
int iter;
int itermin;
int itermax;
complex zmin;
complex zmax;
parameter int inttype;
parameter int rtype;
parameter real power;
parameter int colorby;

    void init(void)
    {
        rx=0.0;
        ry=0.0;
        rr=0.0;
        r=0.0;
        rmin=1.0e12;
        rmax=0.0;
        rsum=0.0;
        rproduct=1.0;
        temp=0.0;
        temp2=0.0;
        iter=0;
        itermin=0;
        itermax=0;
        zmin=(0.0,0.0);
        zmax=(0.0,0.0);

    }
    void loop(void)
    {
        iter=iter+1;
//
// generate closest integer
//
        if ((inttype==0))
        {// round
            temp=cabs(z);
            rr=abs(temp-round(temp));
            temp=real(z);
            rx=abs(temp-round(temp));
            temp=imag(z);
            ry=abs(temp-round(temp));
        }
        else if ((inttype==1))
        {// trunc
            temp=cabs(z);
            rr=abs(temp-trunc(temp));
            temp=real(z);
            rx=abs(temp-trunc(temp));
            temp=imag(z);
            ry=abs(temp-trunc(temp));
        }
        else if ((inttype==2))
        {// ceiling
            temp=cabs(z);
            rr=abs(temp-ceil(temp));
            temp=real(z);
            rx=abs(temp-ceil(temp));
            temp=imag(z);
            ry=abs(temp-ceil(temp));
        }
        else if ((inttype==3))
        {// floor
            temp=cabs(z);
            rr=abs(temp-floor(temp));
            temp=real(z);
            rx=abs(temp-floor(temp));
            temp=imag(z);
            ry=abs(temp-floor(temp));
        }
//
// find distance to integer
//
        if ((rtype==1))
        {// minimum
            r=rr;
            if ((rx<r))
            {
                r=rx;
            }
            if ((ry<r))
            {
                r=ry;
            }
        }
        else if ((rtype==2))
        {// maximum
            r=rr;
            if ((rx>r))
            {
                r=rx;
            }
            if ((ry>r))
            {
                r=ry;
            }
        }
        else if ((rtype==3))
        {// sum
            r=rr+rx+ry;
        }
        else if ((rtype==4))
        {// product
            r=(rr*rx*ry)^(1/3);
        }
        else
        {
                       // euclidean
            r=(rx^power+ry^power+rr^power)^(1/power);
        }
//
// find min, max & build tallies for means
//
        rsum=rsum+r;
        rproduct=rproduct*r;
        if ((r<rmin))
        {
            rmin=r;
            zmin=z;
            itermin=iter;
        }
        if ((r>rmax))
        {
            rmax=r;
            zmax=z;
            itermax=iter;
        }

    }
    void final(void)
    {
        if ((colorby==0))
        {// minimum distance
            index=rmin;
        }
        else if ((colorby==1))
        {// iteration @ min
            index=0.01*itermin;
        }
        else if ((colorby==2))
        {// angle @ min
            temp=atan2(zmin)/pi;
            if ((temp<0.0))
            {
                temp=temp+2.0;
            }
            index=0.5*temp;
        }
        else if ((colorby==3))
        {// maximum distance
            index=rmax;
        }
        else if ((colorby==4))
        {// iteration @ max
            index=0.01*itermax;
        }
        else if ((colorby==5))
        {// angle @ max
            temp=atan2(zmax)/pi;
            if ((temp<0.0))
            {
                temp=temp+2.0;
            }
            index=0.5*temp;
        }
        else if ((colorby==6))
        {// min/max distance angle
            zmin=rmin+flip(rmax);
            temp=atan2(zmin)/pi;
            if ((temp<0.0))
            {
                temp=temp+2.0;
            }
            index=0.5*temp;
        }
        else if ((colorby==7))
        {// arithmetic mean
            index=rsum/iter;
        }
        else if ((colorby==8))
        {// geometric mean
            temp=log(rproduct)/iter;
            index=exp(temp);
        }
        else if ((colorby==9))
        {// amean/gmean angle
            temp=rsum/iter;
            temp2=log(rproduct)/iter;
            temp2=exp(temp2);
            zmax=temp+flip(temp2);
            temp=atan2(zmax)/pi;
            if ((temp<0.0))
            {
                temp=temp+2.0;
            }
            index=0.5*temp;
        }
    }
    void description(void)
    {
        this.title="Pythagorean Triple";
        this.helpfile="lkm-help\lkm-pythagorean.html";

        inttype.caption="integer type";
        inttype.default=0;
        inttype.enum="round\ntrunc\nceil\nfloor";


        rtype.caption="distance type";
        rtype.default=0;
        rtype.enum="Euclidean\nminimum\nmaximum\nsum\nproduct";
        rtype.hint="How the distance to the integer is calculated:  'Euclidean' is like the regular Pythagorean distance, 'minimum' takes the smallest of the 3 components, 'maximum' takes the largest, 'sum' adds them, and 'product' multiplies them.";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="minimum distance\niteration @ min\nangle @ min\nmaximum distance\niteration @ max\nangle @ max\nmin/max combo\narithmetic mean\ngeometric mean\namean/gmean combo";
        colorby.hint="How the color is set; the combos use the 1st quantity for the x and the 2nd for the y, and color by the angle.";


        power.caption="power";
        power.default=2.0;
        power.hint="Exponent for the 'Euclidean' type. Use 2 for standard Pythagorean distance.";

    }
}


bezier-curve { // Kerry Mitchell 08apr00
//
// Colors by the orbit's closest approach to a user-defined Bezier curve.
// The curve is determined by specifying beginning and ending anchor points,
// through which the curve passes, and 2 control points, which influence the
// shape of the curve.
//
real x0;
parameter complex z0;
real y0;
real x1;
parameter complex z1;
real y1;
real x2;
parameter complex z2;
real y2;
real x3;
parameter complex z3;
real y3;
real cx;
real bx;
real ax;
real cy;
real by;
real ay;
real t;
real r;
real Var_x;
real Var_y;
real u;
real v;
int iter;
real rmin;
int itermin;
complex zmin;
parameter real dt;
parameter int colorby;
parameter real nexp;

    void init(void)
    {
        x0=real(z0);
        y0=imag(z0);
        x1=real(z1);
        y1=imag(z1);
        x2=real(z2);
        y2=imag(z2);
        x3=real(z3);
        y3=imag(z3);
        cx=3*(x1-x0);
        bx=3*(x2-x1)-cx;
        ax=x3-x0-cx-bx;
        cy=3*(y1-y0);
        by=3*(y2-y1)-cy;
        ay=y3-y0-cy-by;
        t=0.0;
        r=0.0;
        Var_x=0.0;
        Var_y=0.0;
        u=0.0;
        v=0.0;
        iter=0;
        rmin=1.0e20;
        itermin=0;
        zmin=(0.0,0.0);
    }
    void loop(void)
    {
        iter=iter+1;
        u=real(z);
        v=imag(z);
//
// The curve is parameterized with x(t) and y(t).  Step through several t
// values to find the nearest approach of the orbit to the curve.
//
        t=0.0;
        while ((t<=1.0))
        {
            Var_x=((ax*t+bx)*t+cx)*t+x0;
            Var_y=((ay*t+by)*t+cy)*t+y0;
            r=(Var_x-u)*(Var_x-u)+(Var_y-v)*(Var_y-v);
            if ((r<rmin))
            {
                rmin=r;
                itermin=iter;
                zmin=z;
            }
            t=t+dt;
        }
    }
    void final(void)
    {
        if ((colorby==1))
        {// iteration @ min
            index=0.01*itermin;
        }
        else if ((colorby==2))
        {// angle @ min
            t=atan2(zmin);
            t=t/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
        else if ((colorby==3))
        {// draw section
            u=real(pixel);
            v=imag(pixel);
            t=0.0;
            rmin=1e20;
            while ((t<=1.0))
            {
                Var_x=((ax*t+bx)*t+cx)*t+x0;
                Var_y=((ay*t+by)*t+cy)*t+y0;
                r=(Var_x-u)*(Var_x-u)+(Var_y-v)*(Var_y-v);
                if ((r<rmin))
                {
                    rmin=r;
                }
                t=t+dt;
            }
            index=rmin^nexp;
        }
        else
        {
                       // minimum distance
            index=rmin^nexp;
        }
    }
    void description(void)
    {
        this.title="Bezier Curve";
        this.helpfile="lkm-help\lkm-bezier.html";

        z0.caption="1st anchor point";
        z0.default=(1.0,0.0);
        z0.hint="Curve starts at this point.";


        z1.caption="1st control point";
        z1.default=(1.0,1.0);
        z1.hint="Influences the shape of the curve.";


        z2.caption="2nd control point";
        z2.default=(0.0,0.0);
        z2.hint="Influences the shape of the curve.";


        z3.caption="2nd anchor point";
        z3.default=(0.0,1.0);
        z3.hint="Curve ends at this point.";


        dt.caption="step size";
        dt.default=0.1;
        dt.hint="Decrease for smoother line, increase to see dots.  Should be between 0 & 1.";
        dt.min=0.0;
        dt.max=1.0;


        nexp.caption="power";
        nexp.default=0.1;
        nexp.min=0.0;
        nexp.hint="Decrease to make thinner lines. Use with 'minimum distance' coloring.";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="minimum distance\niteration @ min\nangle @ min\nshow curve";

    }
}


rose-range-lite {
int rangeiter;
real rminfac;
parameter real scale;
parameter real rangewidth;
real rmaxfac;
real rz;
real tz;
real twopi;
real rmin;
real rmax;
real r1;
real rn;
real t1;
real tn;
real rotangle;
parameter real rot;
parameter real ac;
parameter real bc;
parameter real as;
parameter real bs;
parameter complex curvecenter;
parameter int colorby;

    void init(void)
    {
        rangeiter=0;
        rminfac=scale*(1.0-0.5*rangewidth);
        rmaxfac=scale*(1.0+0.5*rangewidth);
        rz=0.0;
        tz=0.0;
        twopi=2.0*pi;
        rmin=0.0;
        rmax=0.0;
        r1=0.0;
        rn=0.0;
        t1=0.0;
        tn=0.0;
        rotangle=-rot/180.0*pi;
    }
    void loop(void)
    {
        tz=atan2(z)+rotangle;
        rz=ac*cos(bc*tz)+as*sin(bs*tz);
        rmin=rz*rminfac;
        rmax=rz*rmaxfac;
        rz=cabs(z-curvecenter);
        if (((rz>=rmin)&&(rz<=rmax)))
        {
            rangeiter=rangeiter+1;
            rn=(rz-rmin)/(rmax-rmin);
            tn=tz;
            if ((rangeiter==1))
            {
                r1=rn;
                t1=tn;
            }
        }
    }
    void final(void)
    {
        if ((rangeiter==0))
        {
            solid=true;
        }
        else
        {

            if ((colorby==0))
            {// first mag
                index=r1;
            }
            else if ((colorby==1))
            {// last mag
                index=rn;
            }
            else if ((colorby==2))
            {// first angle
                if ((t1<0.0))
                {
                    t1=t1+twopi;
                }
                t1=t1/twopi;
                index=t1;
            }
            else if ((colorby==3))
            {// last angle
                if ((tn<0.0))
                {
                    tn=tn+twopi;
                }
                tn=tn/twopi;
                index=tn;
            }
        }
    }
    void description(void)
    {
        this.title="Rose Range Lite";
        this.helpfile="lkm-help\lkm-roserange.html";

        scale.caption="range scale";
        scale.default=1.0;


        rangewidth.caption="range width";
        rangewidth.default=0.1;


        colorby.caption="color by";
        colorby.default=3;
        colorby.enum="first magnitude\nlast magnitude\nfirst angle\nlast angle";


        ac.caption="cos amplitude";
        ac.default=1.0;


        bc.caption="cos frequency";
        bc.default=3.0;


        as.caption="sin amplitude";
        as.default=0.0;


        bs.caption="sin frequency";
        bs.default=0.0;


        curvecenter.caption="curve center";
        curvecenter.default=(0.0,0.0);


        rot.caption="rotation angle";
        rot.default=0.0;
        rot.hint="degrees";

    }
}


spiral { // Kerry Mitchell 13may00
//
// Colors by the orbit's closest approach to a user-defined spiral
//
// The spiral is given as theta = f(r), where r is the distance from the
// center of the spiral.
//
real r;
real rmin;
real tspiral;
complex zspiral;
real rotangle;
parameter real rot;
int iter;
int itermin;
complex zmin;
parameter complex spiralcenter;
parameter real Parm_rmin;
parameter real rmax;
parameter real rpower;
parameter real spiralfactor;
parameter int colorby;

    void init(void)
    {
        r=0.0;
        rmin=1e20;
        tspiral=0.0;
        zspiral=(0,0);
        rotangle=rot/180.0*pi;
        iter=0;
        itermin=0;
        zmin=(0,0);
    }
    void loop(void)
    {
        iter=iter+1;
//
// find r, the distance from the center of the spiral
//   then determine the angle as a function of r
//
        r=cabs(z-spiralcenter);
        if (((r>=Parm_rmin)&&(r<=rmax)))
        {
            tspiral=real(rfunc(r))^rpower;
            tspiral=spiralfactor*tspiral+rotangle;
//
// find the x- and y-coordinates of the point on the spiral
//   with the same r value as the current iterate
//
            zspiral=r*(cos(tspiral)+flip(sin(tspiral)))+spiralcenter;
//
// check the difference between the spiral and the current point
//   update the minimum values
//
            r=cabs(z-zspiral);
            if ((r<rmin))
            {
                rmin=r;
                itermin=iter;
                zmin=z;
            }
        }
    }
    void final(void)
    {
        if ((colorby==1))
        {// iteration number at minimum
            index=itermin/maxiter;
        }
        else if ((colorby==2))
        {// angle at minimum
            tspiral=atan2(zmin)/pi;
            if ((tspiral<0.0))
            {
                tspiral=tspiral+2.0;
            }
            index=0.5*tspiral;
        }
        else if ((colorby==3))
        {// magnitude at minimum
            index=cabs(zmin);
        }
        else if ((colorby==4))
        {// draw spiral
            r=cabs(pixel-spiralcenter);
            if (((r>=Parm_rmin)&&(r<=rmax)))
            {
                tspiral=real(rfunc(r))^rpower;
                tspiral=spiralfactor*tspiral+rotangle;
                zspiral=r*(cos(tspiral)+flip(sin(tspiral)))+spiralcenter;
                r=cabs(pixel-zspiral);
            }
            index=r;
        }
        else
        {
                       // minimum distance
            index=rmin;
        }
    }
    void description(void)
    {
        this.title="Spiral";
        this.helpfile="lkm-help\lkm-spiral.html";

        spiralfactor.caption="spiral factor";
        spiralfactor.default=8.0;
        spiralfactor.hint="Increase to make the spiral tighter; use negative values to reverse the direction.";


        rpower.caption="spiral power";
        rpower.default=1.0;
        rpower.hint="Increase in magnitude to make the spiral tighter; try negative values, too.";


        spiralcenter.caption="spiral center";
        spiralcenter.default=(0,0);


        rot.caption="rotation angle";
        rot.default=0.0;
        rot.hint="degrees";


        Parm_rmin.caption="minimum r";
        Parm_rmin.default=0.0;
        Parm_rmin.hint="To keep the spiral from going too far in.";


        rmax.caption="maximum r";
        rmax.default=1e20;
        rmax.hint="To keep the spiral from going too far out.";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="min. distance\niteration @ min\nangle @ min\nmagnitude @ min\ndraw spiral";


        rfunc.caption="r function";
        rfunc.default="ident" ;
        rfunc.hint="theta = size * [f(r)]^power";

    }
}


grid { // Kerry Mitchell 30aug2000
//
// Colors by the orbit's closest approach to a rectangular grid
//
real cerr;
real cerrmin;
real Var_x;
real Var_y;
real t;
int iter;
int itermin;
complex zmin;
parameter real tdeg;
complex rot;
complex temp;
parameter complex gridcenter;
parameter real Parm_width;
parameter real Parm_height;
parameter int colorby;

    void init(void)
    {
        cerr=0.0;
        cerrmin=1.0e12;
        Var_x=0.0;
        Var_y=0.0;
        t=0.0;
        iter=0;
        itermin=0;
        zmin=(0.0,0.0);
        t=tdeg*pi/180;
        rot=cos(t)-flip(sin(t));
    }
    void loop(void)
    {
        iter=iter+1;
        temp=(z-gridcenter)*rot;
        Var_x=(real(temp)%Parm_width)/Parm_width;
        Var_x=(Var_x+2)%1;
        Var_x=2*abs(0.5-Var_x);
        Var_y=(imag(temp)%Parm_height)/Parm_height;
        Var_y=(Var_y+2)%1;
        Var_y=2*abs(0.5-Var_y);
        if ((Var_x<Var_y))
        {
            cerr=Var_y;
        }
        else
        {

            cerr=Var_x;
        }
        if ((cerr<cerrmin))
        {
            cerrmin=cerr;
            zmin=z;
            itermin=iter;
        }
    }
    void final(void)
    {
        if ((colorby==0))
        {// minimum distance
            index=cerrmin;
        }
        else if ((colorby==1))
        {// iteration @ min
            index=0.01*itermin;
        }
        else if ((colorby==2))
        {// angle @ min
            t=atan2(zmin)/pi;
            if ((t<0.0))
            {
                t=t+2.0;
            }
            index=0.5*t;
        }
    }
    void description(void)
    {
        this.title="Grid";

        gridcenter.caption="center";
        gridcenter.default=(0.0,0.0);


        Parm_height.caption="height";
        Parm_height.default=1.0;
        Parm_height.hint="height of each cell";


        Parm_width.caption="width";
        Parm_width.default=1.0;
        Parm_width.hint="width of each cell";


        tdeg.caption="rotation, deg";
        tdeg.default=0.0;


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="minimum distance\niteration @ min\nangle @ min";

    }
}


single-truchet { // Kerry Mitchell 08nov00
//
// Colors the fractal randomly with truchet tiles
//
real Var_x;
real Var_y;
int intx;
int inty;
real r;
real r0;
real r1;
real rmax;
parameter real n;
parameter real Parm_width;
int tiletype;
int randint;
int modbase;
int multiplier;
int offset;
int niter;
int randiter;
parameter real tdeg;
complex scale;
parameter real scalefac;
complex temp;
parameter complex shift;
parameter real seed;

    void init(void)
    {
        Var_x=0.0;
        Var_y=0.0;
        intx=0;
        inty=0;
        r=0.0;
        r0=0.0;
        r1=0.0;
        rmax=0.5*(2^(1/n)-1)*Parm_width;
        tiletype=0;
        randint=0;
        modbase=2147483647;
        multiplier=1103515245;
        offset=12345;
        niter=0;
        randiter=0;
        r=-tdeg*pi/180;
        scale=(cos(r)+flip(sin(r)))/scalefac;
    }
    void loop(void)
    {
    }
    void final(void)
    {
//
// find fractional x and y from final z
//
        temp=z*scale-shift;
        Var_x=real(temp)+0.5;
        Var_y=imag(temp)+0.5;
        intx=round(Var_x);
        inty=round(Var_y);
        r=Var_x-intx;
        if ((r<0.0))
        {
            Var_x=1+r;
        }
        else
        {

            Var_x=r;
        }
        r=Var_y-inty;
        if ((r<0.0))
        {
            Var_y=1+r;
        }
        else
        {

            Var_y=r;
        }
//
// pick tile type
//
        if ((seed==0))
        {// force type 0
            tiletype=0;
        }
        else if ((seed==1))
        {// force type 1
            tiletype=1;
        }
        else
        {
                       // choose randomly
            intx=round(real(temp));
            inty=round(imag(temp));
            niter=abs(intx+inty+intx*inty);
            randint=seed+niter;
            randiter=0;
            while ((randiter<niter))
            {
                randiter=randiter+1;
                randint=(randint*multiplier+offset)%modbase;
            }
            tiletype=randint%2;
        }
//
// determine radial distances based on tile type
//
        if ((tiletype==1))
        {
            r0=(abs(Var_x)^n+abs(Var_y)^n)^(1/n);
            r1=(abs(Var_x-1)^n+abs(Var_y-1)^n)^(1/n);
        }
        else
        {

            r0=(abs(Var_x-1)^n+abs(Var_y)^n)^(1/n);
            r1=(abs(Var_x)^n+abs(Var_y-1)^n)^(1/n);
        }
//
// find radial distance and set index color
//
//   solid background
//
        solid=true;
//
//   gradation inside arc 0
//
        r=abs(r0-0.5)/rmax;
        if ((r<1))
        {
            solid=false;
            index=r;
        }
//
//   gradation inside arc 1
//
        r=abs(r1-0.5)/rmax;
        if ((r<1))
        {
            solid=false;
            index=r;
        }
    }
    void description(void)
    {
        this.title="Single Truchet";

        shift.caption="shift";
        shift.default=(0,0);
        shift.hint="Use to shift the center of the tiling pattern.";


        scalefac.caption="size";
        scalefac.default=1.0;
        scalefac.hint="Increase or decrease the size of the pattern.";


        tdeg.caption="rotation";
        tdeg.default=0.0;
        tdeg.hint="Rotate the pattern, +ccw, -cw, in degrees.";


        Parm_width.caption="arc width";
        Parm_width.default=0.5;
        Parm_width.min=0.0;
        Parm_width.max=1.0;
        Parm_width.hint="a fraction from 0 to 1";


        n.caption="exponent";
        n.default=2.0;
        n.min=0.0;
        n.hint="Use 2 for circular arcs, 1 for straight lines, and less than 1 for astroid arcs.";


        seed.caption="random seed";
        seed.default=12345;
        seed.min=0;
        seed.hint="Use 0 or 1 to force either tile; otherwise, use a large, odd integer.";

    }
}


double-truchet { // Kerry Mitchell 12nov2000
//
// Colors the fractal randomly with truchet tiles using 4 arcs
//
real Var_x;
real Var_y;
int intx;
int inty;
real r;
real r0;
real r1;
real r2;
real r3;
real ra;
real rb;
real rc;
real rd;
real rmax;
parameter real Parm_width;
real rcutoff;
int tiletype;
int randint;
int modbase;
int multiplier;
int offset;
int niter;
int randiter;
int seedy;
parameter real tdeg;
complex scale;
parameter real scalefac;
real recipn;
parameter real n;
real r2o2;
real omr2o2;
int overtype;
complex temp;
parameter complex shift;
parameter real seed;
parameter bool nozero;

    void init(void)
    {
        Var_x=0.0;
        Var_y=0.0;
        intx=0;
        inty=0;
        r=0.0;
        r0=0.0;
        r1=0.0;
        r2=0.0;
        r3=0.0;
        ra=0.0;
        rb=0.0;
        rc=0.0;
        rd=0.0;
        rmax=Parm_width;
        rcutoff=0.0;
        tiletype=0;
        randint=0;
        modbase=32768;
        multiplier=13821;
        offset=0;
        niter=0;
        randiter=0;
        seedy=0;
        r=-tdeg*pi/180;
        scale=(cos(r)+flip(sin(r)))/scalefac;
        recipn=1/n;
        r2o2=0.5*sqrt(2);
        omr2o2=1-0.5*sqrt(2);
        overtype=0;
    }
    void loop(void)
    {
    }
    void final(void)
    {
//
// find fractional x and y from final z
//
        temp=z*scale-shift;
        Var_x=real(temp)+0.5;
        Var_y=imag(temp)+0.5;
        intx=round(Var_x);
        inty=round(Var_y);
        r=Var_x-intx;
        if ((r<0.0))
        {
            Var_x=1+r;
        }
        else
        {

            Var_x=r;
        }
        r=Var_y-inty;
        if ((r<0.0))
        {
            Var_y=1+r;
        }
        else
        {

            Var_y=r;
        }
//
// pick tile type
//
        if ((seed==0))
        {// force type 0
            tiletype=0;
        }
        else if ((seed==1))
        {// force type 1
            tiletype=1;
        }
        else if ((seed==2))
        {// force type 2
            tiletype=2;
        }
        else if ((seed==3))
        {// force type 3
            tiletype=3;
        }
        else
        {
                       // choose randomly
            niter=0;
            seedy=2*seed+1;
            intx=round(real(temp))*2;
            if ((intx<=0))
            {
                intx=1-intx;
                niter=niter+1;
            }
            inty=round(imag(temp))*2;
            if ((inty<=0))
            {
                inty=1-inty;
                niter=niter+2;
            }
            niter=niter+intx+inty+intx*inty;
            randint=seedy+niter;
            randiter=0;
            while ((randiter<niter))
            {
                randiter=randiter+1;
                randint=(randint*multiplier+offset)%modbase;
            }
            tiletype=(randint%seedy)%4;
            if ((nozero==true))
            {
                tiletype=((randint%seedy)%3)+1;
            }
            randint=(randint*multiplier+offset)%modbase;
            overtype=(randint%seedy)%4;
        }
//
// determine radial distances based on tile type
//
        if ((tiletype==1))
        {
            ra=((abs(Var_x-0.5))^n+(abs(Var_y))^n)^recipn;
            rb=((abs(Var_x-1))^n+(abs(Var_y-0.5))^n)^recipn;
            rc=((abs(Var_x-0.5))^n+(abs(Var_y-1))^n)^recipn;
            rd=((abs(Var_x))^n+(abs(Var_y-0.5))^n)^recipn;
            rcutoff=r2o2-0.5;
        }
        else if ((tiletype==2))
        {
            ra=((abs(Var_x))^n+(abs(Var_y))^n)^recipn;
            rb=((abs(Var_x-1))^n+(abs(Var_y))^n)^recipn;
            rc=((abs(Var_x-1))^n+(abs(Var_y-1))^n)^recipn;
            rd=((abs(Var_x))^n+(abs(Var_y-1))^n)^recipn;
            rcutoff=omr2o2;
        }
        else if ((tiletype==3))
        {
            ra=((abs(Var_x))^n+(abs(Var_y))^n)^recipn;
            rb=((abs(Var_x-1))^n+(abs(Var_y-1))^n)^recipn;
            rc=((abs(Var_x-1))^n+(abs(Var_y))^n)^recipn;
            rd=((abs(Var_x))^n+(abs(Var_y-1))^n)^recipn;
            rcutoff=r2o2;
        }
        else
        {

            ra=abs(Var_x-omr2o2);
            rb=abs(Var_y-omr2o2);
            rc=abs(Var_x-r2o2);
            rd=abs(Var_y-r2o2);
            rcutoff=0;
        }
//
// determine which arc is over which
//
        if ((overtype==1))
        {
            r1=ra;
            r2=rb;
            r3=rc;
            r0=rd;
        }
        else if ((overtype==2))
        {
            r2=ra;
            r3=rb;
            r0=rc;
            r1=rd;
        }
        else if ((overtype==2))
        {
            r3=ra;
            r0=rb;
            r1=rc;
            r2=rd;
        }
        else
        {

            r0=ra;
            r1=rb;
            r2=rc;
            r3=rd;
        }
//
// find radial distance and set index color
//
//   solid background
//
        solid=true;
//
//   gradation inside arc 0
//
        r=abs(r0-rcutoff)/rmax;
        if ((r<=1))
        {
            solid=false;
            index=r;
        }
//
//   gradation inside arc 1
//
        r=abs(r1-rcutoff)/rmax;
        if ((r<=1))
        {
            solid=false;
            index=r;
        }
//
//   gradation inside arc 2
//
        r=abs(r2-rcutoff)/rmax;
        if ((r<=1))
        {
            solid=false;
            index=r;
        }
//
//   gradation inside arc 3
//
        r=abs(r3-rcutoff)/rmax;
        if ((r<=1))
        {
            solid=false;
            index=r;
        }
    }
    void description(void)
    {
        this.title="Double Truchet";

        shift.caption="shift";
        shift.default=(0,0);
        shift.hint="Use to shift the center of the tiling pattern.";


        scalefac.caption="size";
        scalefac.default=1.0;
        scalefac.hint="Increase or decrease the size of the pattern.";


        tdeg.caption="rotation";
        tdeg.default=0.0;
        tdeg.hint="Rotate the pattern, +ccw, -cw, in degrees.";


        Parm_width.caption="arc width";
        Parm_width.default=0.15;
        Parm_width.min=0.0;
        Parm_width.max=1.0;
        Parm_width.hint="From 0 to 1.";


        n.caption="exponent";
        n.default=2.0;
        n.min=0.0;
        n.hint="Use 2 for circular arcs, 1 for straight lines, and less than 1 for astroid arcs.";


        seed.caption="random seed";
        seed.default=123;
        seed.min=0;
        seed.hint="Use 0, 1, 2, or 3 to force either tile; otherwise, use a large integer.";


        nozero.caption="exclude tile 0?";
        nozero.default=false;
        nozero.hint="Check to exclude the straight arc tile.";

    }
}


string-art-rose { // Kerry Mitchell 21mar01
//
// "Sting Art" coloring, in which the "nails" are placed around
// a rose curve.
//
real Var_x;
real Var_y;
real theta;
real theta0;
parameter real t0deg;
real dtheta;
parameter real m;
parameter real n;
real r;
real rmin;
real rlast;
real rfirst;
complex t;
complex z1;
complex z2;
complex zfirst;
complex zlast;
parameter real trdeg;
complex rot;
complex zrot;
int itheta;
int nin;
parameter real ac;
parameter real bc;
parameter real as;
parameter real bs;
parameter complex curvecenter;
parameter real rangewidth;
parameter int colorby;
parameter real power;

    void init(void)
    {
        Var_x=0.0;
        Var_y=0.0;
        theta=0.0;
        theta0=t0deg*pi/180;
        dtheta=m/n*2*pi;
        r=0.0;
        rmin=1e20;
        rlast=0.0;
        rfirst=0.0;
        t=(0,0);
        z1=(0,0);
        z2=(0,0);
        zfirst=(0,0);
        zlast=(0,0);
        theta=trdeg*pi/180;
        rot=cos(theta)-flip(sin(theta));
        zrot=(0,0);
        itheta=0;
        nin=0;
    }
    void loop(void)
    {
        zrot=z*rot;// rotate z instead of the curve
//
// set up 1st line endpoint
//
        theta=theta0;
        r=ac*cos(bc*theta)+as*sin(bs*theta);
        z1=r*(cos(theta)+flip(sin(theta)))+curvecenter;
//
// loop through the rest of the points on the curve:
// next point in the loop becomes the new 2nd line endpoint
//
        itheta=0;
        while ((itheta<n))
        {
            itheta=itheta+1;
            theta=itheta*dtheta+theta0;
            r=ac*cos(bc*theta)+as*sin(bs*theta);
            z2=r*(cos(theta)+flip(sin(theta)))+curvecenter;
//
// find distance from z to the line
//
            t=(zrot-z1)/(z2-z1);
            Var_x=real(t);
            Var_y=sqr(imag(t));
            if ((Var_x<0.0))
            {
                r=sqr(Var_x)+Var_y;
            }
            else if ((Var_x>1.0))
            {
                r=sqr(Var_x-1.0)+Var_y;
            }
            else
            {

                r=Var_y;
            }
            r=sqrt(r)*cabs(z2-z1);// scale lines to have same width
//
// check to see if this point is closer to line, or falls in range
//
            if ((r<rmin))
            {
                rmin=r;
            }
            if ((r<=rangewidth))
            {
                nin=nin+1;
                rlast=r;
                zlast=z;
                if ((nin==1))
                {
                    rfirst=r;
                    zfirst=z;
                }
            }
            z1=z2;// update 1st endpoint
        }
    }
    void final(void)
    {
        if ((colorby==1))
        {// first mag
            if ((nin==0))
            {
                solid=true;
            }
            else
            {

                index=rfirst/rangewidth;
            }
        }
        else if ((colorby==2))
        {// first angle
            if ((nin==0))
            {
                solid=true;
            }
            else
            {

                theta=atan2(zfirst)/pi;
                if ((theta<0.0))
                {
                    theta=theta+2;
                }
                index=theta/2;
            }
        }
        else if ((colorby==3))
        {// last mag
            if ((nin==0))
            {
                solid=true;
            }
            else
            {

                index=rlast/rangewidth;
            }
        }
        else if ((colorby==4))
        {// last angle
            if ((nin==0))
            {
                solid=true;
            }
            else
            {

                theta=atan2(zlast)/pi;
                if ((theta<0.0))
                {
                    theta=theta+2;
                }
                index=theta/2;
            }
        }
        else
        {
                       // minimum dist.
            index=rmin^power;
        }
    }
    void description(void)
    {
        this.title="Rose Curve String Art";

        n.caption="# of points";
        n.default=7;
        n.hint="Both #'s should be coprime, like 7 & 3 or 16 & 9.";


        m.caption="# skipped";
        m.default=3;
        m.hint="Both #'s should be coprime, like 7 & 3 or 16 & 9.";


        curvecenter.caption="curve center";
        curvecenter.default=(0,0);


        t0deg.caption="line rotation";
        t0deg.default=0.0;
        t0deg.hint="Rotates the start of the lines, in degrees, +ccw.";


        trdeg.caption="curve rotation";
        trdeg.default=0.0;
        trdeg.hint="Rotates the curve, in degrees, +ccw.";


        ac.caption="cosine amplitude";
        ac.default=1.0;


        bc.caption="cosine frequency";
        bc.default=2.0;


        as.caption="sine amplitude";
        as.default=0.0;


        bs.caption="sine frequency";
        bs.default=0.0;


        colorby.caption="color by:";
        colorby.default=0;
        colorby.enum="minimum dist.\nfirst mag\nfirst angle\nlast mag\nlast angle";


        rangewidth.caption="range width";
        rangewidth.default=0.1;
        rangewidth.hint="Use with all but 'minimum dist.' 'color by:' choices.";


        power.caption="power";
        power.default=1.0;
        power.hint="Use with color density to control line thickness.  Use with 'minimum dist.' 'color by:' choice.";

    }
}


string-art-astroid { // Kerry Mitchell 21mar2001
//
// "Sting Art" coloring, in which the "nails" are placed around
// an astroid curve.
//
real Var_x;
real Var_y;
real theta;
real theta0;
parameter real t0deg;
real dtheta;
parameter real m;
parameter real n;
real r;
real rmin;
real rlast;
real rfirst;
real fac;
real recipn;
parameter real astroidn;
complex t;
complex z1;
complex z2;
complex zfirst;
complex zlast;
parameter real trdeg;
complex rot;
complex zrot;
int itheta;
int nin;
parameter real astroida;
parameter complex curvecenter;
parameter real rangewidth;
parameter int colorby;
parameter real power;

    void init(void)
    {
        Var_x=0.0;
        Var_y=0.0;
        theta=0.0;
        theta0=t0deg*pi/180;
        dtheta=m/n*2*pi;
        r=0.0;
        rmin=1e20;
        rlast=0.0;
        rfirst=0.0;
        fac=0.0;
        recipn=-1/astroidn;
        t=(0,0);
        z1=(0,0);
        z2=(0,0);
        zfirst=(0,0);
        zlast=(0,0);
        theta=trdeg*pi/180;
        rot=cos(theta)-flip(sin(theta));
        zrot=(0,0);
        itheta=0;
        nin=0;
    }
    void loop(void)
    {
        zrot=z*rot;// rotate z instead of the curve
//
// set up 1st line endpoint
//
        theta=theta0;
        Var_x=abs(cos(theta))^astroidn;
        Var_y=abs(sin(theta))^astroidn;
        fac=(Var_x+Var_y)^recipn;
        r=fac*astroida;
        z1=r*(cos(theta)+flip(sin(theta)))+curvecenter;
//
// loop through the rest of the points on the curve:
// next point in the loop becomes the new 2nd line endpoint
//
        itheta=0;
        while ((itheta<n))
        {
            itheta=itheta+1;
            theta=itheta*dtheta+theta0;
            Var_x=abs(cos(theta))^astroidn;
            Var_y=abs(sin(theta))^astroidn;
            fac=(Var_x+Var_y)^recipn;
            r=fac*astroida;
            z2=r*(cos(theta)+flip(sin(theta)))+curvecenter;
//
// find distance from z to the line
//
            t=(zrot-z1)/(z2-z1);
            Var_x=real(t);
            Var_y=sqr(imag(t));
            if ((Var_x<0.0))
            {
                r=sqr(Var_x)+Var_y;
            }
            else if ((Var_x>1.0))
            {
                r=sqr(Var_x-1.0)+Var_y;
            }
            else
            {

                r=Var_y;
            }
            r=sqrt(r)*cabs(z2-z1);// scale lines to have same width
//
// check to see if this point is closer to line, or falls in range
//
            if ((r<rmin))
            {
                rmin=r;
            }
            if ((r<=rangewidth))
            {
                nin=nin+1;
                rlast=r;
                zlast=z;
                if ((nin==1))
                {
                    rfirst=r;
                    zfirst=z;
                }
            }
            z1=z2;// update 1st endpoint
        }
    }
    void final(void)
    {
        if ((colorby==1))
        {// first mag
            if ((nin==0))
            {
                solid=true;
            }
            else
            {

                index=rfirst/rangewidth;
            }
        }
        else if ((colorby==2))
        {// first angle
            if ((nin==0))
            {
                solid=true;
            }
            else
            {

                theta=atan2(zfirst)/pi;
                if ((theta<0.0))
                {
                    theta=theta+2;
                }
                index=theta/2;
            }
        }
        else if ((colorby==3))
        {// last mag
            if ((nin==0))
            {
                solid=true;
            }
            else
            {

                index=rlast/rangewidth;
            }
        }
        else if ((colorby==4))
        {// last angle
            if ((nin==0))
            {
                solid=true;
            }
            else
            {

                theta=atan2(zlast)/pi;
                if ((theta<0.0))
                {
                    theta=theta+2;
                }
                index=theta/2;
            }
        }
        else
        {
                       // nearest
            index=rmin^power;
        }
    }
    void description(void)
    {
        this.title="Astroid String Art";

        n.caption="# of points";
        n.default=8;
        n.hint="Both #'s should be coprime, like 8 & 3 or 16 & 9.";


        m.caption="# skipped";
        m.default=3;
        m.hint="Both #'s should be coprime, like 8 & 3 or 16 & 9.";


        curvecenter.caption="curve center";
        curvecenter.default=(0,0);


        t0deg.caption="line rotation";
        t0deg.default=0.0;
        t0deg.hint="Rotates the start of the lines, in degrees, +ccw.";


        trdeg.caption="curve rotation";
        trdeg.default=0.0;
        trdeg.hint="Rotates the curve, in degrees, +ccw.";


        astroidn.caption="astroid power";
        astroidn.default=0.5;
        astroidn.min=0.0;
        astroidn.max=256.0;
        astroidn.hint="< 1 for a true astroid, 1 for a diamond, 2 for a circle, 256 for a near square.";


        astroida.caption="astroid scale";
        astroida.default=1.0;
        astroida.min=0.0;
        astroida.hint="Increase to make the curve larger.";


        colorby.caption="color by:";
        colorby.default=0;
        colorby.enum="minimum dist.\nfirst mag\nfirst angle\nlast mag\nlast angle";


        rangewidth.caption="range width";
        rangewidth.default=0.1;
        rangewidth.hint="Use with all but 'minimum dist.' 'color by:' choices.";


        power.caption="power";
        power.default=1.0;
        power.hint="Use with color density to control line thickness.  Use with 'minimum dist.' 'color by:' choice.";

    }
}


crosshatch { // Kerry Mitchell 23may2001
//
// "Colors" the fractal by drawing lines, grids, squares or triangles
// with the variation in spacing relaying the "color" information.
//
real space;
real size;
parameter real sizefac;
real Var_x;
real Var_y;
real t;
real r1;
real r2;
real r3;
real cost;
real sint;
real p;
parameter real randseed;
real colorz;
real count;
int iter;
complex zvalue;
complex zfirst;
complex zmean;
parameter int ztype;
parameter real tdeg;
parameter int shape;
parameter int colorby;
parameter real numsect;
parameter real density;
parameter real offset;
parameter int transfertype;

    void init(void)
    {
        space=0.0;
        size=0.01*sizefac/magn;
        Var_x=real(pixel);
        Var_y=imag(pixel);
        t=0.0;
        r1=0.0;
        r2=0.0;
        r3=0.0;
        cost=0.0;
        sint=0.0;
        p=randseed;
        colorz=0.0;
        count=0.0;
        iter=0;
        zvalue=(0,0);
        zfirst=(0,0);
        zmean=(0,0);
    }
    void loop(void)
    {
        iter=iter+1;
        if ((iter==1))
        {
            zfirst=z;
        }
        zmean=zmean+z;
        p=4*p*(1-p);
    }
    void final(void)
    {
//
// determine count
//
        if ((ztype==1))
        {// first z
            zvalue=zfirst;
        }
        else if ((ztype==2))
        {// last z
            zvalue=z;
        }
        else if ((ztype==3))
        {// average
            zvalue=zmean/iter;
        }
        else
        {
                       // pixel
            zvalue=pixel;
        }
//
// set up lines
//
        Var_x=real(zvalue);
        Var_y=imag(zvalue);
        t=tdeg*pi/180;
        cost=cos(t);
        sint=sin(t);
        r1=cost*Var_y-sint*Var_x;
        r2=sint*Var_y+cost*Var_x;
        if (((shape==2)||(shape==4)))
        {// triangular grid
            t=(tdeg+120)*pi/180;
            cost=cos(t);
            sint=sin(t);
            r2=cost*Var_y-sint*Var_x;
            t=(tdeg+240)*pi/180;
            cost=cos(t);
            sint=sin(t);
            r3=cost*Var_y-sint*Var_x;
        }
        if ((colorby==1))
        {// log(last magnitude)
            count=log(cabs(z));
            if ((numsect>0))
            {
                count=trunc(count*numsect)/numsect;
            }
        }
        else if ((colorby==2))
        {// last angle
            count=atan2(z);
            if ((numsect>0))
            {
                count=trunc(count*numsect)/numsect;
            }
        }
        else
        {
                       // iterations
            count=iter;
        }
//
// determine index (colorz)
//
        t=count*density+offset;
        if ((transfertype==1))
        {// sawtooth
            if ((t<0.0))
            {
                t=t+round(2-t);
            }
            colorz=t%1.0;
        }
        else
        {
                       // sinusoidal
            colorz=(1-cos(t))/2;
        }
//
// build lines
//
        if ((colorz>0.0))
        {
            if ((shape==1))
            {// square grid
                space=size/colorz;
                space=(space-size)*2+size;
                if ((r1<0.0))
                {
                    t=round(2-r1/space)+p;
                }
                else
                {

                    t=p;
                }
                r1=(r1+t*space)%space;
                if ((r2<0.0))
                {
                    t=round(2-r2/space)+p;
                }
                else
                {

                    t=p;
                }
                r2=(r2+t*space)%space;
                if (((r1<=size)||(r2<=size)))
                {
                    r1=r1/size;
                    r2=r2/size;
                    if ((r2>r1))
                    {
                        index=r2;
                    }
                    else
                    {

                        index=r1;
                    }
                }
                else
                {

                    solid=true;
                }
            }
            else if ((shape==2))
            {// triangular grid
                space=size/colorz;
                space=(space-size)*3+size;
                if ((r1<0.0))
                {
                    t=round(2-r1/space)+p;
                }
                else
                {

                    t=p;
                }
                r1=(r1+t*space)%space;
                if ((r2<0.0))
                {
                    t=round(2-r2/space)+p;
                }
                else
                {

                    t=p;
                }
                r2=(r2+t*space)%space;
                if ((r3<0.0))
                {
                    t=round(2-r3/space)+p;
                }
                else
                {

                    t=p;
                }
                r3=(r3+t*space)%space;
                if (((r1<=size)||(r2<=size)||(r3<=size)))
                {
                    r1=r1/size;
                    r2=r2/size;
                    r3=r3/size;
                    if (((r1>=r2)&&(r1>=r3)))
                    {
                        index=r1;
                    }
                    if (((r2>=r1)&&(r2>=r3)))
                    {
                        index=r2;
                    }
                    if (((r3>=r1)&&(r3>=r2)))
                    {
                        index=r3;
                    }
                }
                else
                {

                    solid=true;
                }
            }
            else if ((shape==3))
            {// squares
                space=size/colorz;
                space=(space-size)/3+size;
                if ((r1<0.0))
                {
                    t=round(2-r1/space)+p;
                }
                else
                {

                    t=p;
                }
                r1=(r1+t*space)%space;
                if ((r2<0.0))
                {
                    t=round(2-r2/space)+p;
                }
                else
                {

                    t=p;
                }
                r2=(r2+t*space)%space;
                if (((r1<=size)&&(r2<=size)))
                {
                    r1=r1/size;
                    r2=r2/size;
                    if ((r2>r1))
                    {
                        index=r2;
                    }
                    else
                    {

                        index=r1;
                    }
                }
                else
                {

                    solid=true;
                }
            }
            else if ((shape==4))
            {// triangles
                space=size/colorz;
                space=(space-size)/5+size;
                if ((r1<0.0))
                {
                    t=round(2-r1/space)+p;
                }
                else
                {

                    t=p;
                }
                r1=(r1+t*space)%space;
                if ((r2<0.0))
                {
                    t=round(2-r2/space)+p;
                }
                else
                {

                    t=p;
                }
                r2=(r2+t*space)%space;
                if ((r3<0.0))
                {
                    t=round(2-r3/space)+p;
                }
                else
                {

                    t=p;
                }
                r3=(r3+t*space)%space;
                if (((r1<=size)&&(r2<=size)&&(r3<=size)))
                {
                    r1=r1/size;
                    r2=r2/size;
                    r3=r3/size;
                    if (((r1>=r2)&&(r1>=r3)))
                    {
                        index=r1;
                    }
                    if (((r2>=r1)&&(r2>=r3)))
                    {
                        index=r2;
                    }
                    if (((r3>=r1)&&(r3>=r2)))
                    {
                        index=r3;
                    }
                }
                else
                {

                    solid=true;
                }
            }
            else
            {
                     // line
                space=size/colorz;
                if ((r1<0.0))
                {
                    t=round(2-r1/space)+p;
                }
                else
                {

                    t=p;
                }
                r1=(r1+t*space)%space;
                if ((r1<=size))
                {
                    r1=r1/size;
                    index=r1;
                }
                else
                {

                    solid=true;
                }
            }
        }
        else
        {

            solid=true;
        }
    }
    void description(void)
    {
        this.title="Crosshatch";
        this.helpfile="lkm-help\lkm-crosshatch.html";

        ztype.caption="variable";
        ztype.default=0;
        ztype.enum="pixel\nfirst z\nlast z\naverage";
        ztype.hint="The plane in which the crosshatching is done.";


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="iterations\nlog(last mag)\nlast angle";
        colorby.hint="What to use to determine color.";


        numsect.caption="# sections";
        numsect.default=0;
        numsect.min=0;
        numsect.hint="Use to discretize magnitude or angle; use '0' to turn off.";


        shape.caption="shape";
        shape.default=0;
        shape.enum="line\nsquare grid\ntriangular grid\nsquare\ntriangle";
        shape.hint="What shape to use to color the regions.";


        sizefac.caption="size";
        sizefac.default=1.0;
        sizefac.min=0.0;
        sizefac.hint="The relative size of the lines.";


        tdeg.caption="rotation, deg";
        tdeg.default=45.0;
        tdeg.hint="Rotation of the shape, + ccw.";


        density.caption="density";
        density.default=0.25;
        density.hint="How quickly the 'colors' change; like the 'color density' parameter.";


        offset.caption="offset";
        offset.default=0.0;
        offset.hint="Use to rotate the 'color' distribution; like the 'gradient offset' parameter.";


        transfertype.caption="transfer function";
        transfertype.default=0;
        transfertype.enum="sinusoidal\nsawtooth";
        transfertype.hint="How the count gets reduced to a value between 0 & 1.";


        randseed.caption="random seed";
        randseed.default=0.1;
        randseed.min=0.0;
        randseed.hint="Use to offset the lines randomly in each band; use '0' to turn off.";

    }
}


hilbert-curve { // Kerry Mitchell 15jul2001
//
// Draws a Hilbert curve in the #z plane and colors by the
// orbit's relationship to the curve.  Only uses the last
// iterate.  Allows you to change the locations of the 4
// sub-square centers, and the point where the curve enters
// and exits the block of 4 squares.
//
real Var_x;
real Var_y;
real r;
real rmin;
real u;
real v;
real msb;
complex t;
int iter;
int nlast;
int n;
int fac3;
int fac4;
int power4;
complex zh;
complex zh0;
real h;
parameter real fourcorners;
real ooh;
real ooomh;
real hoomh;
complex z0a;
parameter real enterexit;
complex z0b;
parameter complex centerll;
complex z1;
parameter complex centerul;
complex z2;
parameter complex centerur;
complex z3;
parameter complex centerlr;
complex z4;
complex z5;
parameter real niter;
parameter int colorby;
complex z0;

    void init(void)
    {
        Var_x=0.0;
        Var_y=0.0;
        r=0.0;
        rmin=1e20;
        u=0.0;
        v=0.0;
        msb=0.0;
        t=(0,0);
        iter=0;
        nlast=0;
        n=0;
        fac3=0;
        fac4=0;
        power4=0;
        zh=(0,0);
        zh0=(0,0);
//
// set up endpoints for the lines
//
        h=fourcorners;
        ooh=1/h;
        ooomh=1/(1-h);
        hoomh=h/(1-h);
        z0a=enterexit*h;
        z0b=flip(enterexit*h);
        Var_x=real(centerll)*h;
        Var_y=imag(centerll)*h;
        z1=Var_x+flip(Var_y);
        Var_x=real(centerul)*h;
        Var_y=imag(centerul)*(1-h)+h;
        z2=Var_x+flip(Var_y);
        Var_x=real(centerur)*(1-h)+h;
        Var_y=imag(centerur)*(1-h)+h;
        z3=Var_x+flip(Var_y);
        Var_x=real(centerlr)*(1-h)+h;
        Var_y=imag(centerlr)*h;
        z4=Var_x+flip(Var_y);
        z5=1+flip(enterexit*h);
    }
    void loop(void)
    {
    }
    void final(void)
    {
        iter=0;
        zh0=z*0.25+(0.5,0.5);
        zh=zh0;
        while ((iter<niter))
        {
            iter=iter+1;
            Var_x=real(zh);
            Var_y=imag(zh);
//
// lower left sub-square:  shrink, flip horizontally, rotate by -90 degrees
//
            if (((Var_x<h)&&(Var_y<h)))
            {
                nlast=3;
                u=ooh*Var_y;
                v=ooh*Var_x;
//
// upper left sub-square:  just shrink
//
            }
            else if (((Var_x<h)&&(Var_y>=h)))
            {
                nlast=2;
                u=ooh*Var_x;
                v=ooomh*Var_y-hoomh;
//
// upper right sub-square:  shrink & flip horizontally
//
            }
            else if (((Var_x>=h)&&(Var_y>=h)))
            {
                nlast=1;
                u=ooomh-ooomh*Var_x;
                v=ooomh*Var_y-hoomh;
//
// lower right sub-square:  shrink, rotate by 90 degrees
//
            }
            else if (((Var_x>=h)&&(Var_y<h)))
            {
                nlast=4;
                u=ooh*Var_y;
                v=ooomh-ooomh*Var_x;
            }
            else
            {

                nlast=-1;
                u=Var_x;
                v=Var_y;
            }
            zh=u+flip(v);
            msb=(msb+nlast-1)/4;
            n=n*4+nlast-1;
        }
        if ((n<0))
        {
            solid=true;
        }
        else
        {

            if ((colorby==1))
            {// least significant bit
                index=n/(4^niter);
            }
            else if ((colorby==2))
            {// most significat bit
                index=msb;
            }
            else if ((colorby==3))
            {// distance from initial z
                index=cabs(zh-zh0);
            }
            else if ((colorby==4))
            {// angle from initial z
                r=atan2(zh-zh0)/pi;
                if ((r<0.0))
                {
                    r=r+2;
                }
                index=r/2;
            }
            else
            {
                     // distance to curve
//
// determine how to start line 0-1
//
                z0=z0a;
                iter=0;
                fac3=2;
                fac4=1;
                power4=4;
                while ((iter<niter))
                {
                    iter=iter+1;
                    if ((iter%2==1))
                    {
                        if ((n%power4==fac3))
                        {
                            z0=z0b;
                        }
                    }
                    else
                    {

                        if ((n%power4==fac3))
                        {
                            z0=z0a;
                        }
                    }
                    power4=4*power4;
                    fac4=4*fac4;
                    fac3=fac3+2*fac4;
                }
//
// line from enter to lower left sub-square
//
                t=(zh-z0)/(z1-z0);
                Var_x=real(t);
                Var_y=imag(t);
                if ((Var_x<0.0))
                {
                    r=sqr(Var_x)+sqr(Var_y);
                }
                else if ((Var_x>1.0))
                {
                    r=sqr(Var_x-1.0)+sqr(Var_y);
                }
                else
                {

                    r=sqr(Var_y);
                }
                r=sqrt(r)*cabs(z1-z0);
                if ((r<rmin))
                {
                    rmin=r;
                }
//
// line from lower left to upper left sub-squares
//
                t=(zh-z1)/(z2-z1);
                Var_x=real(t);
                Var_y=imag(t);
                if ((Var_x<0.0))
                {
                    r=sqr(Var_x)+sqr(Var_y);
                }
                else if ((Var_x>1.0))
                {
                    r=sqr(Var_x-1.0)+sqr(Var_y);
                }
                else
                {

                    r=sqr(Var_y);
                }
                r=sqrt(r)*cabs(z2-z1);
                if ((r<rmin))
                {
                    rmin=r;
                }
//
// line from upper left to upper right sub-squares
//
                t=(zh-z2)/(z3-z2);
                Var_x=real(t);
                Var_y=imag(t);
                if ((Var_x<0.0))
                {
                    r=sqr(Var_x)+sqr(Var_y);
                }
                else if ((Var_x>1.0))
                {
                    r=sqr(Var_x-1.0)+sqr(Var_y);
                }
                else
                {

                    r=sqr(Var_y);
                }
                r=sqrt(r)*cabs(z3-z2);
                if ((r<rmin))
                {
                    rmin=r;
                }
//
// line from upper right to lower right sub-squares
//
                t=(zh-z3)/(z4-z3);
                Var_x=real(t);
                Var_y=imag(t);
                if ((Var_x<0.0))
                {
                    r=sqr(Var_x)+sqr(Var_y);
                }
                else if ((Var_x>1.0))
                {
                    r=sqr(Var_x-1.0)+sqr(Var_y);
                }
                else
                {

                    r=sqr(Var_y);
                }
                r=sqrt(r)*cabs(z4-z3);
                if ((r<rmin))
                {
                    rmin=r;
                }
//
// line from lower right sub-square to exit
//
                t=(zh-z4)/(z5-z4);
                Var_x=real(t);
                Var_y=imag(t);
                if ((Var_x<0.0))
                {
                    r=sqr(Var_x)+sqr(Var_y);
                }
                else if ((Var_x>1.0))
                {
                    r=sqr(Var_x-1.0)+sqr(Var_y);
                }
                else
                {

                    r=sqr(Var_y);
                }
                r=sqrt(r)*cabs(z5-z4);
                if ((r<rmin))
                {
                    rmin=r;
                }
//
                index=rmin;
            }
        }
    }
    void description(void)
    {
        this.title="Hilbert curve";

        centerll.caption="lower left center";
        centerll.default=(0.5,0.5);
        centerll.hint="Center of the lower left sub-square. Make both coordinates between 0 & 1; use (0.5,0.5) for standard Hilbert curve.";


        centerul.caption="upper left center";
        centerul.default=(0.5,0.5);
        centerul.hint="Center of the upper left sub-square. Make both coordinates between 0 & 1; use (0.5,0.5) for standard Hilbert curve.";


        centerur.caption="upper right center";
        centerur.default=(0.5,0.5);
        centerur.hint="Center of the upper right sub-square. Make both coordinates between 0 & 1; use (0.5,0.5) for standard Hilbert curve.";


        centerlr.caption="lower right center";
        centerlr.default=(0.5,0.5);
        centerlr.hint="Center of the lower right sub-square. Make both coordinates between 0 & 1; use (0.5,0.5) for standard Hilbert curve.";


        enterexit.caption="enter/exit";
        enterexit.default=0.5;
        enterexit.min=0.0;
        enterexit.max=1.0;
        enterexit.hint="Where the curve enters and exits the block of 4 sub-squares. Between 0 & 1; use 0.5 for standard Hilbert curve.";


        fourcorners.caption="4 corners";
        fourcorners.default=0.5;
        fourcorners.hint="Where the 4 corners meet.  Between 0 & 1; use 0.5 for standard Hilbert curve.";


        niter.caption="iterations";
        niter.default=0;
        niter.min=0;


        colorby.caption="color by";
        colorby.default=0;
        colorby.enum="distance\nlast = lsb\nlast = msb\ndistance from z\nangle from z";

    }
}


line { // Kerry Mitchell 15jul2001
//
// Just draws a line between 2 points
//
real Var_x;
real Var_y;
real r;
complex t;
parameter complex z1;
parameter complex z2;

    void init(void)
    {
        Var_x=0.0;
        Var_y=0.0;
        r=0.0;
        t=(0,0);
    }
    void loop(void)
    {
    }
    void final(void)
    {
        t=(z-z1)/(z2-z1);
        Var_x=real(t);
        Var_y=imag(t);
        if ((Var_x<0.0))
        {
            r=sqr(Var_x)+sqr(Var_y);
        }
        else if ((Var_x>1.0))
        {
            r=sqr(Var_x-1.0)+sqr(Var_y);
        }
        else
        {

            r=sqr(Var_y);
        }
        r=sqrt(r)*cabs(z2-z1);
        index=r;
    }
    void description(void)
    {
        this.title="Line";

        z1.caption="1st endpoint";
        z1.default=(0,0);


        z2.caption="2nd endpoint";
        z2.default=(1,1);

    }
}