pgd_HRing_Mandelbrot_c {
; Herman Ring Mandelbrot set sliced in c-plane.
; Default contains the default Julia set's location.
init:
  c = #pixel
  ; This section selects the critical point differently than the sqrt function would. The sqrt function causes seams;
  ; this method does not.
  bool flag = (@cp == 1)
  z = (3+sqr(c)+sqrt(9-10*sqr(c)+sqr(sqr(c))))/(4*c)
  IF (real(c)*real(c)-imag(c)*imag(c)-5 >= 0)
    flag = !flag
  ENDIF
  IF (imag(c) < 0)
    flag = !flag
  ENDIF
  IF (flag)
    z = (3+sqr(c)-sqrt(9-10*sqr(c)+sqr(sqr(c))))/(4*c)
  ENDIF
  complex z0 = z
  complex a
  float zm
  IF (@aexp)
    a = exp(2*#pi*(0,1)*@a)
  ELSE
    a = @a
  ENDIF
loop:
  z = a*sqr(z)*(z-c)/(1-z*c)
  zm = |z|
  IF (@what == 1)
    IF (zm >= @bailout)
      z = z0 ; Fudge to prevent divergent points mysteriously showing as "outside"
      zm = |z|
    ENDIF
  ENDIF
bailout:
  ((zm < @bailout) || (@what == 1)) && ((zm > 1/@bailout) || (@what == 0))
default:
  magn = 0.3
  method = onepass
  maxiter = 10000000
  title = "Herman Ring Mandelbrot Set (c-plane)"
  param cp
    caption = "Critical point"
    hint = "There are two critical points with parameter-dependent fates that give slightly different M-sets when iterated."
    enum = "Critical point #1" "Critical point #2"
    default = 0
  endparam
  param what
    caption = "What to draw"
    hint = "Points may converge to zero, infinity, or one or more parameter-dependent finite attractors. Select whether points going to zero, infinity, or both count as 'outside'."
    enum = "Points diverging" "Points convergent to 0" "Points doing both"
    default = 2
  endparam
  param aexp
    caption = "Use a as an angle"
    hint = "The fractal will use exp(2*pi*i*alpha) in place of alpha. Real numbers between 0 and 1 will then give Herman rings."
    default = true
  endparam
  param a
    caption = "Parameter a"
    hint = "Herman rings occur for a = exp(2*pi*i*alpha) for real numbers alpha between 0 and 1."
    default = (0.61803398,0)
  endparam  
  param bailout
    caption = "Bailout"
    hint = "Escape radius (reciprocal is used for trap radius about 0)"
    default = 100000.0
    min = 0
  endparam
switch:
  type = "pgd_HRing_Julia"
  what = what
  aexp = aexp
  a = a
  c = #pixel
  bailout = bailout
;  magn = 0.3
}

pgd_HRing_Mandelbrot_c_2 {
; Herman Ring Mandelbrot set sliced in c-plane.
; Default contains the default Julia set's location.
; Slightly different treatment of critical points.
init:
  c = #pixel
  ; This section selects the critical point differently than the sqrt function would. The sqrt function causes seams;
  ; this method does not.
  bool flag = (@cp == 1)
  z = (3+sqr(c)+sqrt(9-10*sqr(c)+sqr(sqr(c))))/(4*c)
  IF (real(c)*real(c)-imag(c)*imag(c)-5 >= 0)
    flag = !flag
  ENDIF
  IF (flag)
    z = (3+sqr(c)-sqrt(9-10*sqr(c)+sqr(sqr(c))))/(4*c)
  ENDIF
  complex z0 = z
  complex a
  float zm
  IF (@aexp)
    a = exp(2*#pi*(0,1)*@a)
  ELSE
    a = @a
  ENDIF
loop:
  z = a*sqr(z)*(z-c)/(1-z*c)
  zm = |z|
  IF (@what == 1)
    IF (zm >= @bailout)
      z = z0 ; Fudge to prevent divergent points mysteriously showing as "outside"
      zm = |z|
    ENDIF
  ENDIF
bailout:
  ((zm < @bailout) || (@what == 1)) && ((zm > 1/@bailout) || (@what == 0))
default:
  magn = 0.3
  method = onepass
  maxiter = 10000000
  title = "Herman Ring Mandelbrot Set (c-plane) altcp"
  param cp
    caption = "Critical point"
    hint = "There are two critical points with parameter-dependent fates that give slightly different M-sets when iterated."
    enum = "Critical point #1" "Critical point #2"
    default = 0
  endparam
  param what
    caption = "What to draw"
    hint = "Points may converge to zero, infinity, or one or more parameter-dependent finite attractors. Select whether points going to zero, infinity, or both count as 'outside'."
    enum = "Points diverging" "Points convergent to 0" "Points doing both"
    default = 2
  endparam
  param aexp
    caption = "Use a as an angle"
    hint = "The fractal will use exp(2*pi*i*alpha) in place of alpha. Real numbers between 0 and 1 will then give Herman rings."
    default = true
  endparam
  param a
    caption = "Parameter a"
    hint = "Herman rings occur for a = exp(2*pi*i*alpha) for real numbers alpha between 0 and 1."
    default = (0.61803398,0)
  endparam  
  param bailout
    caption = "Bailout"
    hint = "Escape radius (reciprocal is used for trap radius about 0)"
    default = 100000.0
    min = 0
  endparam
switch:
  type = "pgd_HRing_Julia"
  what = what
  aexp = aexp
  a = a
  c = #pixel
  bailout = bailout
;  magn = 0.3
}

pgd_HRing_Mandelbrot_a {
; Herman Ring Mandelbrot set sliced in a-plane.
; Default contains the default Julia set's location.
init:
  complex a = #pixel
  bool flag = (@cp == 1)
  z = (3+sqr(@c)+sqrt(9-10*sqr(@c)+sqr(sqr(@c))))/(4*@c)
  IF (real(@c)*real(@c)-imag(@c)*imag(@c)-5 >= 0)
    flag = !flag
  ENDIF
  IF (imag(@c) < 0)
    flag = !flag
  ENDIF
  IF (flag)
    z = (3+sqr(@c)-sqrt(9-10*sqr(@c)+sqr(sqr(@c))))/(4*@c)
  ENDIF
  complex z0 = z
  float zm
loop:
  z = a*sqr(z)*(z-@c)/(1-z*@c)
  zm = |z|
  IF (@what == 1)
    IF (zm >= @bailout)
      z = z0 ; Fudge to prevent divergent points mysteriously showing as "outside"
      zm = |z|
    ENDIF
  ENDIF
bailout:
  ((zm < @bailout) || (@what == 1)) && ((zm > 1/@bailout) || (@what == 0))
default:
  center = (1,0)
  magn = 0.6
  method = onepass
  maxiter = 10000000
  title = "Herman Ring Mandelbrot Set (a-plane)"
  param cp
    caption = "Critical point"
    hint = "There are two critical points with parameter-dependent fates that give slightly different M-sets when iterated."
    enum = "Critical point #1" "Critical point #2"
    default = 0
  endparam
  param what
    caption = "What to draw"
    hint = "Points may converge to zero, infinity, or one or more parameter-dependent finite attractors. Select whether points going to zero, infinity, or both count as 'outside'."
    enum = "Points diverging" "Points convergent to 0" "Points doing both"
    default = 2
  endparam
  param c
    caption = "Parameter c"
    hint = "Herman rings occur for c imaginary or real, bigger in size than about 3.5, and some other points too."
    default = (3.7,0)
  endparam  
  param bailout
    caption = "Bailout"
    hint = "Escape radius (reciprocal is used for trap radius about 0)"
    default = 100000.0
    min = 0
  endparam
switch:
  type = "pgd_HRing_Julia"
  what = what
  aexp = aexp
  a = #pixel
  c = c
  bailout = bailout
;  magn = 0.3
}

pgd_HRing_Julia {
; Herman Ring Julia sets.
; Can come out rather shoddy with solid guessing depending.
; Default is, in fact, a Herman ring.
init:
  z = #pixel
  complex a
  float zm
  IF (@aexp)
    a = exp(2*#pi*(0,1)*@a)
  ELSE
    a = @a
  ENDIF
loop:
  z = a*sqr(z)*(z-@c)/(1-z*@c)
  zm = |z|
  IF (@what == 1)
    IF (|z| >= @bailout)
      z = #pixel ; Fudge to prevent divergent points mysteriously showing as "outside"
      zm = |z|
    ENDIF
  ENDIF
bailout:
  ((zm < @bailout) || (@what == 1)) && ((zm > 1/@bailout) || (@what == 0))
default:
  center = (1.5,0)
  magn = 0.3
  method = onepass
  maxiter = 1000
  title = "Herman Ring Julia sets"
  param what
    caption = "What to draw"
    hint = "Points may converge to zero, infinity, or one or more parameter-dependent finite attractors. Select whether points going to zero, infinity, or both count as 'outside'."
    enum = "Points diverging" "Points convergent to 0" "Points doing both"
    default = 2
  endparam
  param aexp
    caption = "Use a as an angle"
    hint = "The fractal will use exp(2*pi*i*alpha) in place of alpha. Real numbers between 0 and 1 will then give Herman rings."
    default = true
  endparam
  param a
    caption = "Parameter a"
    hint = "Herman rings occur for a = exp(2*pi*i*alpha) for real numbers alpha between 0 and 1."
    default = (0.61803398,0)
  endparam  
  param c
    caption = "Parameter c"
    hint = "Another complex parameter. It should be either real or imaginary and larger than about 3.5 for Herman rings to appear."
    default = (3.7,0.0)
  endparam
  param bailout
    caption = "Bailout"
    hint = "Escape radius (reciprocal is used for trap radius about 0)"
    default = 100000.0
    min = 0
  endparam
}

pgd_GenericCubicJulia {
global:
  complex ta = -3*@a*@a
  bool foundatt1 = false
  bool foundatt2 = false
  complex att1 = -@a
  complex att2 = @a
  IF (@what != 0)
    int i = 0
    foundatt1 = true
    foundatt2 = true
    WHILE ((i < #maxit) && (foundatt1 || foundatt2))
      IF (foundatt1)
        att1 = (att1*att1 + ta)*att1 + @c
        IF (|att1| > @bailout)
          foundatt1 = false
        ENDIF
      ENDIF
      IF (foundatt2)
        att2 = (att2*att2 + ta)*att2 + @c
        IF (|att2| > @bailout)
          foundatt2 = false
        ENDIF
      ENDIF
      i = i + 1
    ENDWHILE
    IF (foundatt1 && foundatt2)
      i = 0
      WHILE ((i < #maxit) && foundatt2)
        att2 = (att2*att2 + ta)*att2 + @c
        IF (|att2 - att1| < (1/@bailout))
          foundatt2 = false
        ENDIF
        i = i + 1
      ENDWHILE
    ENDIF
  ENDIF
init:
  z = #pixel
  bool escaped = false
  bool bail = false
  IF ((@what == 1) && (!foundatt1))
    bail = true
  ENDIF
  IF ((@what == 2) && (!foundatt2))
    bail = true
  ENDIF
loop:
  IF (!escaped)
    z = (z*z + ta)*z + @c
    IF (|z| > @bailout)
      bail = true
    ENDIF
    IF (foundatt1)
      IF (|z - att1| < (1/@bailout))
        escaped = true
        bail = (@what == 2)
      ENDIF
    ENDIF
    IF (foundatt2)
      IF (|z - att2| < (1/@bailout))
        escaped = true
        bail = (@what == 1)
      ENDIF
    ENDIF    
  ENDIF    
bailout:
  !bail
default:
  title = "Generic Cubic Julia"
  param a
    caption = "Julia seed a"
    default = (-0.5,-0.2)
  endparam
  param c
    caption = "Julia seed c"
    default = (0.13,-0.341)
  endparam
  param what
    caption = "Draw what"
    enum = "Escaping points" "Finite attractor 1" "Finite attractor 2"
    default = 0
  endparam
  param bailout
    caption = "Bailout"
    hint = "Escape radius (inverse is used for finite attractors)"
    default = 128.0
    min = 0
  endparam
}

pgd_GenericCubicMandelbrot_a {
init:
  complex a = #pixel
  complex ta = -3*a*a
  IF ((@what == 1) || (@what == 3) || (@what == 5))
    z = -a
  ELSE
    z = a
  ENDIF
  complex w = -z ; We iterate also the critical point not being studied.
  bool bail = false
  bool zesc = false
  bool wesc = false
  bool same = true
  bool div = false
  int i = 0
  int j = 0
  bool odd = false
  complex z_values[@zmax]
  complex w_values[@zmax]
loop:
  IF (!zesc)
    z = (z*z + ta)*z + @c
    IF (|z| > @bailout)
      zesc = true
      same = false
      div = true
      bail = true
      IF (@allout && (@what != 0) && (@what != 1))
        z = (501,10)
      ENDIF
    ENDIF
  ENDIF
  IF (@what > 1)
    IF (i < @zmax)
      z_values[i] = z
    ENDIF
    IF (!wesc)
      w = (w*w + ta)*w + @c
      IF (i < @zmax)
        w_values[i] = w
      ENDIF
      IF (|w| > @bailout)
        IF ((@what == 4) || (@what == 5))
          bail = true
          IF (@allout)
            z = (501,10)
          ENDIF
        ENDIF
        wesc = true
        same = false
        div = true
      ELSEIF (odd)
        IF (|w - w_values[j]| < (1/@bailout))
          wesc = true
        ENDIF
      ENDIF
    ENDIF
    IF (odd)
      IF (!zesc)
        IF (|z - z_values[j]| < (1/@bailout))
          zesc = true
          IF (@allout && div && ((@what == 2) || (@what == 3)))
            bail = true
          ENDIF
        ENDIF
      ENDIF
      j = j + 1
    ENDIF
    i = i + 1
    IF (i < @zmax)
      odd = !odd
    ELSE
      odd = false
    ENDIF
    IF (((zesc && wesc) || (i == #maxit)) && (!bail))
      ; We're done
      IF (same)
        ; Iterate z some more and see if it becomes equal to itself again or to w.
        zesc = false
        complex z0 = z
        int i0 = i
        odd = false
        WHILE (!zesc)
          z = (z*z + ta)*z + @c
          IF (|z - w| < (1/@bailout))
            zesc = true
            odd = true
          ELSEIF (|z| > @bailout)
            same = false
            zesc = true
          ELSE
            IF (|z - z0| < (1/@bailout))
              same = false
              zesc = true
            ENDIF
          ENDIF
          i = i + 1
          IF (i >= #maxit)
            same = odd
            zesc = true
          ENDIF
        ENDWHILE
        i = i - i0
      ENDIF
      IF (@allout)
        bail = true
        IF ((@what == 2) || (@what == 3))
          IF (same)
            z = (501,10)
          ENDIF
        ELSE
          IF (!same)
            z = (501,10)
          ENDIF
        ENDIF
      ELSEIF ((@what == 2) || (@what == 3))
        bail = same
      ELSE
        bail = !same && !div
      ENDIF
    ENDIF
  ENDIF
bailout:
  !bail
default:
  title = "Generic Cubic Mandelbrot (a plane)"
  param what
    caption = "Draw what"
    enum = "CP1 escapes" "CP2 escapes" "CP1 own attractor" "CP2 own attractor" "CP1 shared attractor" "CP2 shared attractor"
    default = 0
  endparam
  param allout
    caption = "All outside"
    default = false
  endparam
  param c
    caption = "Slice location"
    hint = "Enter a complex number here"
    default = (0.13,-0.341)
  endparam
  param zmax
    caption = "Maximum iterations to record"
    hint = "Maximum iterations is best, but values over 1,000,000 may cause low memory errors."
    default = 1000
    visible = ((@what != 0) && (@what != 1))
  endparam
  param bailout
    caption = "Bailout"
    hint = "Escape radius (inverse is used for finite attractors)"
    default = 128.0
    min = 0
  endparam
switch:
  type = "pgd_GenericCubicJulia"
  a = #pixel
  c = c
  bailout = bailout  
}

pgd_GenericCubicMandelbrot_c {
init:
  complex c = #pixel
  complex ta = -3*@a*@a
  IF ((@what == 1) || (@what == 3) || (@what == 5))
    z = -@a
  ELSE
    z = @a
  ENDIF
  complex w = -z ; We iterate also the critical point not being studied.
  bool bail = false
  bool zesc = false
  bool wesc = false
  bool same = true
  bool div = false
  int i = 0
  int j = 0
  bool odd = false
  complex z_values[@zmax]
  complex w_values[@zmax]
loop:
  IF (!zesc)
    z = (z*z + ta)*z + c
    IF (|z| > @bailout)
      zesc = true
      same = false
      div = true
      bail = true
      IF (@allout && (@what != 0) && (@what != 1))
        z = (501,10)
      ENDIF
    ENDIF
  ENDIF
  IF (@what > 1)
    IF (i < @zmax)
      z_values[i] = z
    ENDIF
    IF (!wesc)
      w = (w*w + ta)*w + c
      IF (i < @zmax)
        w_values[i] = w
      ENDIF
      IF (|w| > @bailout)
        IF ((@what == 4) || (@what == 5))
          bail = true
          IF (@allout)
            z = (501,10)
          ENDIF
        ENDIF
        wesc = true
        same = false
        div = true
      ELSEIF (odd)
        IF (|w - w_values[j]| < (1/@bailout))
          wesc = true
        ENDIF
      ENDIF
    ENDIF
    IF (odd)
      IF (!zesc)
        IF (|z - z_values[j]| < (1/@bailout))
          zesc = true
          IF (@allout && div && ((@what == 2) || (@what == 3)))
            bail = true
          ENDIF
        ENDIF
      ENDIF
      j = j + 1
    ENDIF
    i = i + 1
    IF (i < @zmax)
      odd = !odd
    ELSE
      odd = false
    ENDIF
    IF (((zesc && wesc) || (i == #maxit)) && (!bail))
      ; We're done
      IF (same)
        ; Iterate z some more and see if it becomes equal to itself again or to w.
        zesc = false
        complex z0 = z
        int i0 = i
        odd = false
        WHILE (!zesc)
          z = (z*z + ta)*z + c
          IF (|z - w| < (1/@bailout))
            zesc = true
            odd = true
          ELSEIF (|z| > @bailout)
            same = false
            zesc = true
          ELSE
            IF (|z - z0| < (1/@bailout))
              same = false
              zesc = true
            ENDIF
          ENDIF
          i = i + 1
          IF (i >= #maxit)
            same = odd
            zesc = true
          ENDIF
        ENDWHILE
        i = i - i0
      ENDIF
      IF (@allout)
        bail = true
        IF ((@what == 2) || (@what == 3))
          IF (same)
            z = (501,10)
          ENDIF
        ELSE
          IF (!same)
            z = (501,10)
          ENDIF
        ENDIF
      ELSEIF ((@what == 2) || (@what == 3))
        bail = same
      ELSE
        bail = !same && !div
      ENDIF
    ENDIF
  ENDIF
bailout:
  !bail
default:
  title = "Generic Cubic Mandelbrot (c plane)"
  param what
    caption = "Draw what"
    enum = "CP1 escapes" "CP2 escapes" "CP1 own attractor" "CP2 own attractor" "CP1 shared attractor" "CP2 shared attractor"
    default = 0
  endparam
  param allout
    caption = "All outside"
    default = false
  endparam
  param a
    caption = "Slice location"
    hint = "Enter a complex number here"
    default = (-0.5,-0.2)
  endparam
  param zmax
    caption = "Maximum iterations to record"
    hint = "Maximum iterations is best, but values over 1,000,000 may cause low memory errors."
    default = 1000
    visible = ((@what != 0) && (@what != 1))
  endparam
  param bailout
    caption = "Bailout"
    hint = "Escape radius (inverse is used for finite attractors)"
    default = 128.0
    min = 0
  endparam
switch:
  type = "pgd_GenericCubicJulia"
  a = #pixel
  c = c
  bailout = bailout  
}

; Generic rotated slices for 4d parameter spaces formula:
; a = (xcos t + iycos u + zsin t + iwsin u)
; c = (-xsin t - iysin u + zcos t + iwcos u)
; (x + iy) is pixel, (z + iw) is slice offset
; t and u together determine slice orientation.
; When t = u we get the special case
; a = (cos t)pixel + (sin t)offset
; c = -(sin t)pixel + (cos t)offset
; This is a C-linear combination, so conformal features of
; 2-D Mandelbrots such as circular buds will guaranteeably
; be present. When t != u wackier things might occur.
; For example, to get the a-plane t = u = 0, to get the
; c-plane t = u = pi/2, and to get real(a), real(c) mapped,
; t = 0, u = -pi/2 resulting in a = x - iw, c = z + iy.

pgd_GenericCubicMandelbrot_anyslice {
init:
  complex a
  complex c
  float v = real(#pixel)
  float vv = imag(#pixel)
  float p = real(@offset)
  float q = imag(@offset)
  float r = @t*#pi/2
  float s = @u*#pi/2
  a = v*cos(r) + (0,1)*vv*cos(s) + p*sin(r) + (0,1)*q*sin(s)
  c = p*cos(r) + (0,1)*q*cos(s) - v*sin(r) - (0,1)*vv*sin(s)
  complex ta = -3*a*a
  IF ((@what == 1) || (@what == 3) || (@what == 5))
    z = -a
  ELSE
    z = a
  ENDIF
  complex w = -z ; We iterate also the critical point not being studied.
  bool bail = false
  bool zesc = false
  bool wesc = false
  bool same = true
  bool div = false
  int i = 0
  int j = 0
  bool odd = false
  complex z_values[@zmax]
  complex w_values[@zmax]
loop:
  IF (!zesc)
    z = (z*z + ta)*z + c
    IF (|z| > @bailout)
      zesc = true
      same = false
      div = true
      bail = true
      IF (@allout && (@what != 0) && (@what != 1))
        z = (501,10)
      ENDIF
    ENDIF
  ENDIF
  IF (@what > 1)
    IF (i < @zmax)
      z_values[i] = z
    ENDIF
    IF (!wesc)
      w = (w*w + ta)*w + c
      IF (i < @zmax)
        w_values[i] = w
      ENDIF
      IF (|w| > @bailout)
        IF ((@what == 4) || (@what == 5))
          bail = true
          IF (@allout)
            z = (501,10)
          ENDIF
        ENDIF
        wesc = true
        same = false
        div = true
      ELSEIF (odd)
        IF (|w - w_values[j]| < (1/@bailout))
          wesc = true
        ENDIF
      ENDIF
    ENDIF
    IF (odd)
      IF (!zesc)
        IF (|z - z_values[j]| < (1/@bailout))
          zesc = true
          IF (@allout && div && ((@what == 2) || (@what == 3)))
            bail = true
          ENDIF
        ENDIF
      ENDIF
      j = j + 1
    ENDIF
    i = i + 1
    IF (i < @zmax)
      odd = !odd
    ELSE
      odd = false
    ENDIF
    IF (((zesc && wesc) || (i == #maxit)) && (!bail))
      ; We're done
      IF (same)
        ; Iterate z some more and see if it becomes equal to itself again or to w.
        zesc = false
        complex z0 = z
        int i0 = i
        odd = false
        WHILE (!zesc)
          z = (z*z + ta)*z + c
          IF (|z - w| < (1/@bailout))
            zesc = true
            odd = true
          ELSEIF (|z| > @bailout)
            same = false
            zesc = true
          ELSE
            IF (|z - z0| < (1/@bailout))
              same = false
              zesc = true
            ENDIF
          ENDIF
          i = i + 1
          IF (i >= #maxit)
            same = odd
            zesc = true
          ENDIF
        ENDWHILE
        i = i - i0
      ENDIF
      IF (@allout)
        bail = true
        IF ((@what == 2) || (@what == 3))
          IF (same)
            z = (501,10)
          ENDIF
        ELSE
          IF (!same)
            z = (501,10)
          ENDIF
        ENDIF
      ELSEIF ((@what == 2) || (@what == 3))
        bail = same
      ELSE
        bail = !same && !div
      ENDIF
    ENDIF
  ENDIF
bailout:
  !bail
default:
  title = "Generic Cubic Mandelbrot (any way you want to slice it)"
  param what
    caption = "Draw what"
    enum = "CP1 escapes" "CP2 escapes" "CP1 own attractor" "CP2 own attractor" "CP1 shared attractor" "CP2 shared attractor"
    default = 0
  endparam
  param allout
    caption = "All outside"
    default = false
  endparam
  param offset
    caption = "Slice offset"
    hint = "Enter a complex number here"
    default = (0.1,0)
  endparam
  param t
    caption = "Slice angle 1"
    hint = "0-1 is a quarter rotation."
    default = 0.5
  endparam
  param u
    caption = "Slice angle 2"
    hint = "0-1 is a quarter rotation."
    default = 0.5
  endparam
  param zmax
    caption = "Maximum iterations to record"
    hint = "Maximum iterations is best, but values over 1,000,000 may cause low memory errors."
    default = 1000
    visible = ((@what != 0) && (@what != 1))
  endparam
  param bailout
    caption = "Bailout"
    hint = "Escape radius (inverse is used for finite attractors)"
    default = 128.0
    min = 0
  endparam
switch:
  type = "pgd_GenericCubicJulia"
  a = #pixel
  c = c
  bailout = bailout  
}

pgd_Mandelnewton {
; Iterates Newton's method for roots of (((c^2+c)^2+c)^2+c)^2+c... = 0.
; The number of nested Mandelbrot iterations is determined via parameter.
; The higher the number the "busier" the fractal, and the more the
; outlines of the M-set are revealed. The fractals themselves are
; Julia sets with 2^n basins of attraction.
; All points bail out. Use decomp, iterations, or similar color methods.
; Orbit traps might be interesting too.
init:
  z = #pixel
  float b = 1/@bailout
  float b2 = b*b
  complex zz = z
loop:
  complex oz = z
  complex zd = 1
  int i = 1
  WHILE (i < @numi)
    i = i + 1
    zd = 2*zd*z + 1
    z = sqr(z) + oz
    IF (|zd| > @bailout)
      i = @numi
    ENDIF
  ENDWHILE
  ; z is now m(oz) and zd is m'(oz), with m(x) = (x^2+x)^2+x...
  z = oz - (z/zd) ; Apply Newton's method once
  IF (@zeroinside)
    IF (|z| < b2)
      z = zz ; Trap z
    ENDIF
  ENDIF
bailout:
  |z - oz| > b
default:
  title = "Mandelnewton Julia sets"
  periodicity = 0
  maxiter = 1000
  center = (-0.5,0)
  param numi
    caption = "Order"
    default = 6
    min = 1
  endparam
  param zeroinside
    caption = "A(0) inside"
    default = false
  endparam
  param bailout
    caption = "Bailout"
    default = 100000
    min = 0
  endparam
}

pgd_TwofoldJulia {
  ; Two-fold rotational symmetry and complementary critical
  ; point fates make for some very intertwining (and
  ; entertaining) Julia sets.
global:
  bool zeroattracts = (|@seed| > 1)
  bool foundatt1 = false
  bool foundatt2 = false
  complex att1 = 1
  complex att2 = -1
  float b = 1/@bailout
  IF (@what != 0)
    int j = 0
    foundatt1 = true
    ; Locate attractor 1 goes to, if any.
    WHILE ((j < #maxit) && foundatt1)
      att1 = att1/(@seed * (sqr(att1) + 1))
      IF (|att1| > @bailout)
        att1 = 0
      ENDIF
      IF (zeroattracts)
        IF (|att1| < b)
          foundatt1 = false   ; it went to zero.
        ENDIF
      ENDIF
      j = j + 1
    ENDWHILE
    ; Is it really an attractor? Or did it just wander off
    ; and die somewhere? Is the attractor self-complementary
    ; or is there a distinct mirror attractor?    
    IF (foundatt1)
      j = 0
      foundatt1 = false
      foundatt2 = true
      bool nd = true
      complex a = att1
      WHILE ((j < #maxit) && nd)
        att1 = att1/(@seed * (sqr(att1) + 1))
        IF (|att1| > @bailout)
          att1 = 0
        ENDIF
        IF (zeroattracts)
          IF (|att1| < b)
            foundatt1 = false   ; it went to zero.
            foundatt2 = false   ; so will the other one.
            nd = false     
          ENDIF
        ENDIF
        IF (|att1 - a| < b)
          foundatt1 = true      ; it repeated, yay!
        ENDIF        
        att2 = att2/(@seed * (sqr(att2) + 1))
        IF (|att2| > @bailout)
          att2 = 0
        ENDIF
        IF (|att2 - a| < b)
          foundatt2 = false      ; other cp went to same.
        ENDIF
        j = j + 1
      ENDWHILE
    ENDIF
  ENDIF
init:
  z = #pixel
  bool notdone = true
  bool landedright = false
  IF ((@what == 0) && (!zeroattracts))
    notdone = false
  ENDIF
  IF ((@what == 1) && (!foundatt1))
    notdone = false
  ENDIF
  IF ((@what == 2) && (!foundatt2))
    notdone = false
  ENDIF
  i = 0
loop:
  z = z/(@seed * (sqr(z) + 1))
  IF (|z| > @bailout)
    z = 0
  ENDIF
  IF (zeroattracts)
    IF (|z| < b)
      notdone = false
      landedright = (@what == 0)
    ENDIF
  ENDIF
  IF (foundatt1)
    IF (|z - att1| < b)
      notdone = false
      landedright = (@what == 1)
    ENDIF
  ENDIF
  IF (foundatt2)
    IF (|z - att2| < b)
      notdone = false
      landedright = (@what == 2)
    ENDIF
  ENDIF
  IF (i == #maxit)
    IF (@what != 3)
      notdone = false
      landedright = false
    ENDIF
  ENDIF
  i = i + 1
  IF (!notdone)
    IF (!landedright)
      z = (501,10)  ; Tell Smooth (Generalized) 2 to color this thing solid.
    ENDIF
  ENDIF
bailout:
  notdone
default:
  title = "Twofold Julia"
  periodicity = 0
  maxiter = 1000
  method = onepass
  magn = 0.5
  param seed
    caption = "Julia seed"
    default = (0.09583333334,0.3041666667)
  endparam
  param what
    caption = "Draw what"
    enum = "Basin of 0" "Other attractor #1" "Other attractor #2" "Trapped points"
    default = 1
  endparam
  param bailout
    caption = "Bailout"
    hint = "Escape radius (inverse is used for finite attractors)"
    default = 100000.0
    min = 0
  endparam
switch:
  type = "pgd_TwofoldMandel"
  bailout = bailout
}

pgd_TwofoldMandel {
  ; Two-fold rotational symmetry and complementary critical
  ; point fates make for some very intertwining (and
  ; entertaining) Julia sets. This Mandelbrot catalogues
  ; them and is fascinating in its own right. Buds come in
  ; five types -- "zero", "merge", "split", "shuffle", and
  ; "mix" -- at least in principle. This formula attempts to
  ; detect and render them separately. "Zero" definitely
  ; exists and is the region surrounding the set -- outside
  ; the big circle zero is an attractor, and captures the
  ; critical points +/-i, as well as the point at infinity.
  ; Interior buds and cardioids come in the other flavors,
  ; defined by symmetries. Each is associated with an even
  ; number of attracting points, which may form a single
  ; orbit or divide into two in a symmetrical manner. Thus,
  ; odd length cycles can occur. Buds with two odd period
  ; attractors, one grabbing 1 and the other -1, are
  ; "merge" buds. Rotate the fractal 180 degrees and you get
  ; mirror image "split" buds, in which the same points form
  ; a single orbit that grabs both 1 and -1. "Shuffle" buds
  ; appear opposite other "shuffle" buds and have a single
  ; attracting cycle of length divisible by four. But the
  ; order the points are visited differs between the mirrored
  ; buds. "Mix" buds appear opposite other "mix" buds and
  ; have two even, symmetric attracting cycles. Half the points
  ; in each cycle swap with half the points in the other if
  ; you turn the fractal 180 degrees though, hence the name.
  ; It is provable that mix buds can only touch one another
  ; and shuffle buds, while the others can touch one another
  ; and the rim of the "zero" circle. The two large central
  ; buds are a "merge" bud with a symmetric pair of attracting
  ; fixed points and a corresponding "split" bud with an
  ; attracting 2-cycle.
  ; The associated Julia sets have 2-fold rotational
  ; symmetry which either preserves a single attracting basin
  ; (of zero or an even-length attracting cycle) or swaps two
  ; basins (of symmetrical, equal length cycles). They are
  ; also symmetric under inversion, since infinity is a pre-
  ; image of zero. Inside the Mandelbrot circle, there can
  ; only be one or two attractors; i, -i, 0, and infinity all
  ; belong to the Julia set in this region of the parameter
  ; space! 1 and -1 may go to a single, even-length cycle,
  ; two cycles that may have odd length, or wander endlessly.
  ; In the latter case the Julia set grows to consume the
  ; entire Riemann sphere. Otherwise, the Julia set is very
  ; complex and intricate, separating as it does the even and
  ; odd points in a single cycle, or the two separate basins
  ; of separate cycles. If there are two cycles each of period
  ; greater tha one, each cycle's basin is intruded into by
  ; fingers of Julia set, as in all such fractals; but the
  ; only other basin is also divided into cyclic basins, which
  ; makes the boundary wildly convoluted with intricate spirals.
  ; There's almost always a spiral centered at zero.
  ; Outside the M-set circle, zero seems to suck up everything,
  ; and the Julia set falls apart into a dust that contains
  ; none of the points 1, -1, i, -i, 0, and infinity.
  ; Explore away!
init:
  z = 1
  bool notdone = true
  bool landedright = false
  float b = 1/@bailout
  complex z_values[@zmax]
  int i = 0
  int j = 0
  bool odd = false
  bool zeroattracts = (|#pixel| > 1)
  int zm = @zmax - 1
loop:
  z = z/(#pixel*(sqr(z) + 1))
  IF (!@allinside)
    IF (zeroattracts)
      IF (|z| < b)
        notdone = false
        landedright = (@what == 0)
      ENDIF
    ENDIF
    IF (i < zm)
      z_values[i] = z
      IF (odd)
        IF (|z - z_values[j]| < b)
          ; caught a periodic cycle
          bool same = false
          int k = 0
          complex w = -1
          WHILE ((k < #maxit) && (!same))
            w = w/(#pixel*(sqr(w) + 1))
            IF (|w - z| < b)
              same = true
            ENDIF
            k = k + 1
          ENDWHILE
          IF (same)
            ; Splitter or shuffler
            ; See if iterating using -#pixel produces a cycle
            ; containing its negatives.
            k = 0
            complex u = w
            same = false
            WHILE ((k < #maxit) && (!same))
              w = -w/(#pixel*(sqr(w) + 1))
              IF (|w + u| < b)
                same = true
              ENDIF
              k = k + 1
            ENDWHILE
            IF (same)
              ; Shuffler
              landedright = (@what == 4)
            ELSE
              ; Splitter
              landedright = (@what == 2)
            ENDIF
          ELSE
            ; Merger or mixer
            ; See if iterating using -#pixel produces a cycle
            ; containing its negatives.
            k = 0
            complex u = w
            WHILE ((k < #maxit) && (!same))
              w = -w/(#pixel*(sqr(w) + 1))
              IF (|w + u| < b)
                same = true
              ENDIF
              k = k + 1
            ENDWHILE
            IF (same)
              ; Merger
              landedright = (@what == 1)
            ELSE
              ; Mixer
              landedright = (@what == 5)
            ENDIF
          ENDIF
          notdone = false
        ENDIF
        j = j + 1
      ENDIF
      odd = !odd
    ELSE
      IF (@what != 3)
        notdone = false
        landedright = false
      ENDIF
    ENDIF
    IF (@ztrap && (@what == 3))
      IF (|z| > @ztbail)
        notdone = false
        landedright = true
      ENDIF
    ENDIF
    i = i + 1
    IF (!notdone)
      IF (!landedright)
        z = (501,10)  ; Tell Smooth (Generalized) 2 to color this thing solid.
      ENDIF
    ENDIF
  ENDIF
bailout:
  notdone
default:
  title = "Twofold Mandelbrot"
  periodicity = 0
  maxiter = 1000
  method = onepass
  magn = 1.5
  param allinside
    caption = "All inside"
    default = false
  endparam
  param what
    caption = "Draw what"
    enum = "Zero attracts" "Odd cycle" "Splitter cycle" "Trapped points" "Shuffler cycle" "Mixer cycle"
    default = 1
    visible = !@allinside
  endparam
  param ztrap
    caption = "Trap points"
    default = false
    visible = ((@what == 3) && !@allinside)
  endparam
  param ztbail
    caption = "Trap distance"
    default = 4.0
    visible = ((@what == 3) && !@allinside)
  endparam
  param zmax
    caption = "Max iters recorded"
    default = 1000
    visible = !@allinside
  endparam
  param bailout
    caption = "Bailout"
    hint = "Escape radius (inverse is used for finite attractors)"
    default = 100000.0
    min = 0
    visible = !@allinside
  endparam
switch:
  type = "pgd_TwofoldJulia"
  seed = #pixel
  bailout = bailout
}

pgd_Twofold2_Julia {
  ; Two-fold rotational symmetry and complementary critical
  ; point fates make for some very intertwining (and
  ; entertaining) Julia sets.
  ; f(z) = (z^2 - 1/3)/cz^3
global:
  bool foundatt1 = false
  bool foundatt2 = true
  complex att1 = 1
  complex att2 = -1
  float b = 1/@bailout
  int j = 0
  ; Locate attractor 1 goes to, if any.
  WHILE (j < #maxit)
    att1 = (sqr(att1) - 1/3)/(@seed * sqr(att1) * att1)
    j = j + 1
  ENDWHILE
  ; Is it really an attractor? Or did it just wander off
  ; and die somewhere? Is the attractor self-complementary
  ; or is there a distinct mirror attractor?    
  j = 0
  complex a = att1
  WHILE (j < #maxit)
    att1 = (sqr(att1) - 1/3)/(@seed * sqr(att1) * att1)
    IF (|att1 - a| < b)
      foundatt1 = true      ; it repeated, yay!
    ENDIF        
    att2 = (sqr(att2) - 1/3)/(@seed * sqr(att2) * att2)
    IF (|att2 - a| < b)
      foundatt2 = false      ; other cp went to same.
    ENDIF
    j = j + 1
  ENDWHILE
init:
  z = #pixel
  bool notdone = true
  bool landedright = false
  IF ((@what == 1) && (!foundatt1))
    notdone = false
  ENDIF
  IF ((@what == 2) && (!foundatt2))
    notdone = false
  ENDIF
  i = 0
loop:
  z = (sqr(z) - 1/3)/(@seed * sqr(z) * z)
  IF (|z| > @bailout)
    notdone = false
    landedright = (@what == 0)
  ENDIF
  IF (foundatt1)
    IF (|z - att1| < b)
      notdone = false
      landedright = (@what == 1)
    ENDIF
  ENDIF
  IF (foundatt2)
    IF (|z - att2| < b)
      notdone = false
      landedright = (@what == 2)
    ENDIF
  ENDIF
  IF (i == #maxit)
    IF (@what != 3)
      notdone = false
      landedright = false
    ENDIF
  ENDIF
  i = i + 1
  IF (!notdone)
    IF (!landedright)
      z = (501,10)  ; Tell Smooth (Generalized) 2 to color this thing solid.
    ENDIF
  ENDIF
bailout:
  notdone
default:
  title = "Twofold2 Julia"
  periodicity = 0
  maxiter = 1000
  method = onepass
  magn = 0.5
  param seed
    caption = "Julia seed"
    default = (0.75, 0)
  endparam
  param what
    caption = "Draw what"
    enum = "Main basin" "Other attractor #1" "Other attractor #2" "Trapped points"
    default = 0
  endparam
  param bailout
    caption = "Bailout"
    hint = "Escape radius (inverse is used for finite attractors)"
    default = 100000.0
    min = 0
  endparam
switch:
  type = "pgd_Twofold2_Mandel"
  bailout = bailout
}

pgd_Twofold2_Mandel {
  ; Two-fold rotational symmetry and complementary critical
  ; point fates make for some very intertwining (and
  ; entertaining) Julia sets.
  ; f(z) = (z^2 - 1/3)/cz^3
init:
  z = 1
  bool notdone = true
  bool landedright = false
  int i = 0
  int j = 0
  bool odd = false
  complex z_values[@zmax]
  float b = 1/@bailout
loop:
  z = (sqr(z) - 1/3)/(#pixel * sqr(z) * z)
  IF (|z| > @bailout)
    notdone = false
    landedright = (@what == 0)
  ENDIF
  z_values[i] = z
  IF (odd)
    IF (|z - z_values[j]| < b)
      notdone = false
      complex w = -1
      j = 0
      bool sep = true
      WHILE ((j < #maxit) && sep)
        w = (sqr(w) - 1/3)/(#pixel * sqr(w) * w)
        IF (|w - z| < b)
          sep = false
        ENDIF
        j = j + 1
      ENDWHILE
      IF (sep)
        landedright = (@what == 2)
      ELSE
        landedright = (@what == 1)
      ENDIF
    ENDIF
    j = j + 1
  ENDIF
  odd = !odd
  i = i + 1
  IF (i == @zmax)
    IF (@what != 3)
      notdone = false
      landedright = false
    ENDIF
  ENDIF
  IF (!notdone)
    IF (!landedright)
      z = (501,10)  ; Tell Smooth (Generalized) 2 to color this thing solid.
    ENDIF
  ENDIF
bailout:
  notdone
default:
  title = "Twofold2 Mandelbrot"
  periodicity = 0
  maxiter = 1000
  method = onepass
  param what
    caption = "Draw what"
    enum = "Main basin" "One finite attractor" "Two attractors" "Trapped points"
    default = 0
  endparam
  param zmax
    caption = "Max iters to record"
    default = 1000
  endparam
  param bailout
    caption = "Bailout"
    hint = "Escape radius (inverse is used for finite attractors)"
    default = 100000.0
    min = 0
  endparam
switch:
  type = "pgd_Twofold2_Julia"
  seed = #pixel
  bailout = bailout
}

pgd_ApollonianTwist {
  ; Each iteration one of a collection of fractional linear
  ; transformations is applied, selected randomly. Points that
  ; stray near infinity escape.
  ; Fractional linear transformation is (az + b)/(cz + d);
  ; we force ad - bc = 1 and a = 1, do d = 1 + bc.
  ; Fixed points are (+/-sqrt(b(bc + 4)) - bc^2)/2c^2
  ; b = 0, bc = -4, or c = 0 give a single fixed point (at
  ; infinity in the latter case); otherwise there are two.
  ;
  ; Use Closest Approach of Attractor to Point of Origin
  ; coloring with this for good results. Try other inside
  ; options too. No points bail out!!
  ;
  ; To find an interesting fractal, use Orbit Density
  ; coloring and set the "replace" parameter to the name
  ; of one of the other parameters. It determines a parameter
  ; to vary with the pixel, producing a "mandelbrot" catalogue
  ; dimension. With Orbit Density coloring and the Lighting
  ; gradient, darker greys correspond to parameters with
  ; more curvelike than dustlike or space-filling attractors.
  ; Use say 0.01 tolerance and power 10 in the Orbit Density
  ; settings. Find a dark spot and right click the parameter
  ; you picked to replace, pick eyedropper, and click the
  ; dark spot. The image will regenerate unchanged. Move the
  ; "replace" option to another parameter. Do this a few times
  ; cycling through all the parameters and then switch to
  ; "replace" = "None" and go to Closest Approach to Point of
  ; Origin coloring, power = 10, etc. -- and locate the fractal.
  ; Outzoom or use Riemann Sphere Mapping Projections in
  ; "Orthographic" or "Mollweide" mode to show the whole
  ; Riemann sphere (set location to center 0, 0 magnification
  ; 1 to get the full sphere shown) and then zoom in where
  ; there are dark areas. These are the fractal. If it's
  ; generating too slowly lower maximum iterations; if the
  ; shape of the attractor seems pleasing, use Switch to go
  ; to the faster formula that precomputes the attractor
  ; and set inside to Basic from lkm, option real. Pick
  ; a more colorful gradient.
init:
  complex b[@numt]
  complex c[@numt]
  complex d[@numt]
  int j = 0
  complex bb = @bstart
  IF (@rep == 1)
    bb = #pixel
  ENDIF
  complex cc = @cstart
  IF (@rep == 3)
    cc = #pixel
  ENDIF
  WHILE (j < @numt)
    b[j] = bb
    c[j] = cc
    d[j] = 1 + bb*cc
    IF (@rep == 2)
      bb = (bb * #pixel)
    ELSE
      bb = (bb * @bmul)
    ENDIF
    IF (@rep == 4)
      cc = (cc * #pixel)
    ELSE
      cc = (cc * @cmul)
    ENDIF
    j = j + 1
  ENDWHILE
  IF (@rep == 0)
    z = #pixel
  ELSE
    z = 0
  ENDIF
  float w = 0.61803398
loop:
  w = 3.99999*w*(1 - w)
  int tt = trunc(w*@numt)
  IF (tt == @numt)
    tt = @numt - 1
  ENDIF
  z = (z + b[tt])/(c[tt]*z + d[tt])
bailout:
  true
default:
  title = "Apollonian Twist"
  periodicity = 0
  param numt
    caption = "Number of transformations"
    default = 4
    min = 3
  endparam
  param bstart
    caption = "Start value 1"
    default = (2,0)
  endparam
  param bmul
    caption = "Value 1 growth factor"
    default = (0,0.2)
  endparam
  param cstart
    caption = "Start value 2"
    default = (0,1)
  endparam
  param cmul
    caption = "Value 2 growth factor"
    default = (0.5,0)
  endparam
  param rep
    caption = "Which to replace"
    hint = "Setting this to other than 'none' gives Mandelbrot-like slice info."
    enum = "None" "Start value 1" "Value 1 growth factor" "Start value 2" "Value 2 growth factor"
    default = 0
  endparam
switch:
  type = "pgd_ApollonianTwist2"
  numt = numt
  bstart = bstart
  bmul = bmul
  cstart = cstart
  cmul = cmul
}

pgd_ApollonianTwist2 {
  ; Fast final-rendering version that precomputes attractor.
  ; Ends with real(z) related to distance from pixel to
  ; attractor. So, use color option "real".
global:
  complex b[@numt]
  complex c[@numt]
  complex d[@numt]
  int j = 0
  complex bb = @bstart
  complex cc = @cstart
  WHILE (j < @numt)
    b[j] = bb
    c[j] = cc
    d[j] = 1 + bb*cc
    bb = (bb * @bmul)
    cc = (cc * @cmul)
    j = j + 1
  ENDWHILE
  complex q = 0
  float w = 0.61803398
  j = 0
  complex att[#maxit]
  int ttt = 0
  WHILE (j < 50)
    w = 3.99999*w*(1 - w)
    int tt = trunc(w*@numt)
    IF (tt == @numt)
      tt = @numt - 1
    ENDIF
    IF (j == 0)
      ttt = tt
    ENDIF
    q = (q + b[tt])/(c[tt]*q + d[tt])
    j = j + 1
  ENDWHILE
  j = 0
  WHILE (j < #maxit)
    w = 3.99999*w*(1 - w)
    int tt = trunc(w*@numt)
    IF (tt == @numt)
      tt = @numt - 1
    ENDIF
    q = (q + b[tt])/(c[tt]*q + d[tt])
    att[j] = q
    j = j + 1
  ENDWHILE
init:
  int i = 0
  float p = 0
  complex r = (#pixel + b[ttt])/(c[ttt]*#pixel + d[ttt])
loop:
  IF (i == 0)
    p = |r - att[0]|
    z = p^(1/@power)
  ELSE
    float s = |r - att[i]|
    IF (s < p)
      p = s
      z = p^(1/@power)
    ENDIF
  ENDIF
  i = i + 1
bailout:
  true
default:
  title = "Apollonian Twist Fast"
  periodicity = 0
  param numt
    caption = "Number of transformations"
    default = 4
    min = 3
  endparam
  param bstart
    caption = "Start value 1"
    default = (2,0)
  endparam
  param bmul
    caption = "Value 1 growth factor"
    default = (0,0.2)
  endparam
  param cstart
    caption = "Start value 2"
    default = (0,1)
  endparam
  param cmul
    caption = "Value 2 growth factor"
    default = (0.5,0)
  endparam
  param power
    caption = "Power (scales coloring)"
    default = 10.0
  endparam
switch:
  type = "pgd_ApollonianTwist"
  numt = numt
  bstart = bstart
  bmul = bmul
  cstart = cstart
  cmul = cmul
}

pgd_ParabolicWorkshop {
init:
  z = #pixel
  complex lambda = @mag*exp(2*#pi*(0,1)*@num/@denom)
  complex crit = -lambda/(@power*@alpha)
  IF ((|z^(@power - 1) - crit|) < 0.002)
    z = @bailout + 1
  ENDIF
loop:
  z = lambda*z + @alpha*z^@power
bailout:
  |z| < @bailout
default:
  title = "Parabolic Workshop"
  param num
    caption = "Multiplier Numerator"
    default = 2
  endparam
  param denom
    caption = "Multiplier Denominator"
    default = 3
  endparam
  param power
    caption = "Power"
    default = 2
    min = 2
  endparam
  param alpha
    caption = "Coefficient"
    default = (1,0)
  endparam
  param mag
    caption = "Magnitude (1 for parabolic)"
    default = 1.0
  endparam
  param bailout
    caption = "Bailout"
    default = 100000.0
  endparam
}

pgd_FunkyFixedCubeM {
init:
  IF (real(#pixel) >= 0)
    z = -#pixel - sqrt(sqr(#pixel) - 3)
  ELSE
    z = -#pixel + sqrt(sqr(#pixel) - 3)
  ENDIF
  z = z/3
loop:
  z = z + (#pixel + z)*sqr(z)
bailout:
  |z| < 100
default:
  title = "Funky Fixed Points Cubic Mset"
}

pgd_HenonM {
  ; There's apparently no analogue in systems of two real
  ; variables of critical points as dynamical determinants
  ; (no, points where the Jacobian is singular don't cut
  ; it -- I tested that, and besides this particular
  ; system has none!) so you'll have to make several layers
  ; with different start values to get something approaching
  ; a complete picture. Use Lyapunov coloring inside,
  ; "brightest" merge modes, black or dark outside coloring,
  ; and you might get exceptionally nice images however.
  ;
  ; Check "Lyapunov" to pass Lyapunov information to the
  ; PGD Lyapunov coloring algorithm in z.
init:
  z = @start
  int i = 0
  complex dxx = 1 ; Components of Jacobian of nth iterate
  complex dxy = 0 ; wind up in here.
  complex dyx = 0 ; We start with the identity matrix.
  complex dyy = 1
  bool notdone = true
loop:
  z = 1 - real(#pixel)*sqr(real(z)) + imag(z) + (0,1)*imag(#pixel)*real(z)
  IF (|z| >= @bailout)
    notdone = false
  ENDIF
  IF (@lya)
    IF ((i == #maxit - 1) || !notdone)
      ; Pass Lyapunov information in z
      IF (@lw == 2)
        ; We want the larger eigenvalue of the Jacobian.
        complex alpha = dxx + dyy
        complex beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
        z = log(real(alpha + beta)/2)/i
      ELSEIF (@lw == 3)
        ; We want the smaller eigenvalue.
        alpha = dxx + dyy
        beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
        z = log(real(alpha - beta)/2)/i
      ELSEIF(@lw == 0)
        ; We want the determinant.
        z = log(real(dxx*dyy - dxy*dyx))/i
      ELSEIF(@lw == 1)
        ; We want the ratio of the eigenvalues.
        alpha = dxx + dyy
        beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
        z = log(real(alpha + beta)/real(alpha - beta))/i
      ELSEIF(@lw == 4)
        ; We want the imaginary part.
        beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
        z = log(imag(beta))/i
      ENDIF
    ELSE
      i = i + 1
      ; Constructing Lyapunov information
      ; We multiply our Jacobian matrix on the left by the
      ; current iteration's Jacobian.
      float d = -2*real(#pixel)*real(z)
      complex nxx = dxx*d + dyx
      complex nxy = dxy*d + dyy
      dyx = dxx*imag(#pixel)
      dyy = dxy*imag(#pixel)
      dxx = nxx
      dxy = nxy
    ENDIF
  ENDIF
bailout:
  notdone
default:
  title = "Henon Mandelbrotlike Set"
  center = (1,0)
  param start
    caption = "Starting (x,y)"
    default = (0,0)
  endparam
  param lya
    caption = "Lyapunov"
    default = false
  endparam
  param lw
    caption = "Show what"
    enum = "Expansivity" "Stretch" "Larger" "Smaller" "Twist"
    default = 0
  endparam
  param bailout
    caption = "Bailout"
    default = 100000.0
    min = 0.0
  endparam
switch:
  type = "pgd_HenonJ"
  seed = #pixel
  lya = lya
  lw = lw
  bailout = bailout
}

pgd_HenonJ {
  ; There's potentially no limit to the number of distinct
  ; attracting cycles, so they can't be selected by critical
  ; point. Instead, they are selected by periodicity.
  ; Set one layer to "outside" and if any areas are solid,
  ; set a new layer to "periodic" and try various integers.
  ; The periodic layer can be made exact, or made to show
  ; any multiple of a period. Any remaining dark areas lie
  ; in another cycle, or they may be ergodic. An "ergodic"
  ; layer catches any apparently aperiodic points. Periodic
  ; and ergodic make all points bail out. Outside makes
  ; only divergent points bail out.
  ;
  ; Check "Lyapunov" to pass Lyapunov information to the
  ; PGD Lyapunov coloring algorithm in z.
init:
  z = #pixel
  bool notdone = true
  bool landedright = false
  int per = @period
  IF (@what == 2)
    per = 1
  ENDIF
  complex z_values[@period]
  int i = 0
  complex dxx = 1 ; Components of Jacobian of nth iterate
  complex dxy = 0 ; wind up in here.
  complex dyx = 0 ; We start with the identity matrix.
  complex dyy = 1
  int stopwi = trunc(#maxit/10)
  IF (stopwi < per)
    stopwi = per
  ENDIF
loop:
  z = 1 - real(@seed)*sqr(real(z)) + imag(z) + (0,1)*imag(@seed)*real(z)
  IF (@lya)
    ; Constructing Lyapunov information
    ; We multiply our Jacobian matrix on the left by the
    ; current iteration's Jacobian.
    float d = -2*real(#pixel)*real(z)
    complex nxx = dxx*d + dyx
    complex nxy = dxy*d + dyy
    dyx = dxx*imag(#pixel)
    dyy = dxy*imag(#pixel)
    dxx = nxx
    dxy = nxy
  ENDIF
  IF (|z| > @bailout)
    notdone = false
    landedright = (@what == 0)
  ELSEIF (@what != 0)
    IF (i < stopwi)
      IF (i >= per)
        IF (|z - z_values[i % per]| < 1/@bailout)
          notdone = false
          IF (@what == 1)
            landedright = true
            int j = 1
            WHILE (j < per)
              IF(|z - z_values[(i + j)%per]| < @tol/@bailout)
                landedright = false
              ENDIF
              j = j + 1
            ENDWHILE
          ENDIF
        ENDIF
      ENDIF
      z_values[i % per] = z
    ELSE
      IF (|z - z_values[i % per]| < 1/@bailout)
        notdone = false
        IF ((@what == 1) && !@exact)
          landedright = true
          j = 1
          WHILE (j < per)
            IF(|z - z_values[(i + j)%per]| < @tol/@bailout)
              landedright = false
            ENDIF
            j = j + 1
          ENDWHILE
        ENDIF
      ENDIF
    ENDIF
    i = i + 1
    IF (i == (#maxit - 1))
      IF ((@what == 0) && @lya)
        ; Pass Lyapunov information in z
        IF (@lw == 2)
          ; We want the larger eigenvalue of the Jacobian.
          complex alpha = dxx + dyy
          complex beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(real(alpha + beta)/2)/i
        ELSEIF (@lw == 3)
          ; We want the smaller eigenvalue.
          alpha = dxx + dyy
          beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(real(alpha - beta)/2)/i
        ELSEIF(@lw == 0)
          ; We want the determinant.
          z = log(real(dxx*dyy - dxy*dyx))/i
        ELSEIF(@lw == 1)
          ; We want the ratio of the eigenvalues.
          alpha = dxx + dyy
          beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(real(alpha + beta)/real(alpha - beta))/i
        ELSEIF(@lw == 4)
          ; We want the imaginary part.
          beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(imag(beta))/i
        ENDIF
      ELSE
        notdone = false
        landedright = (@what == 2)
      ENDIF
    ENDIF
  ENDIF
  IF (!notdone && @lya)
    ; Pass Lyapunov information in z
    IF (@lw == 2)
      ; We want the larger eigenvalue of the Jacobian.
      complex alpha = dxx + dyy
      complex beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(real(alpha + beta)/2)/i
    ELSEIF (@lw == 3)
      ; We want the smaller eigenvalue.
      alpha = dxx + dyy
      beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(real(alpha - beta)/2)/i
    ELSEIF(@lw == 0)
      ; We want the determinant.
      z = log(real(dxx*dyy - dxy*dyx))/i
    ELSEIF(@lw == 1)
      ; We want the ratio of the eigenvalues.
      alpha = dxx + dyy
      beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(real(alpha + beta)/real(alpha - beta))/i
    ELSEIF(@lw == 4)
      ; We want the imaginary part.
      beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(imag(beta))/i
    ENDIF
  ENDIF  
  IF (!notdone && !landedright)
    z = (501,10)
  ENDIF
bailout:
  notdone
default:
  title = "Henon Julia Sets"
  param seed
    caption = "Julia seed"
    default = (1.4,0.3)
  endparam
  param what
    caption = "Draw what"
    enum = "Outside" "Periodic" "Ergodic"
    default = 0
  endparam
  param period
    caption = "Period"
    default = 2
    visible = (@what == 1)
  endparam
  param tol
    caption = "Tolerance"
    default = 100.0
    hint = "Tolerance when checking not really a shorter \
            period."
    visible = (@what == 1)
  endparam
  param exact
    caption = "Exact"
    default = false
    visible = (@what == 1)
  endparam
  param lya
    caption = "Lyapunov"
    default = false
  endparam
  param lw
    caption = "Show what"
    enum = "Expansivity" "Stretch" "Larger" "Smaller" "Twist"
    default = 0
    visible = @lya
  endparam
  param bailout
    caption = "Bailout"
    default = 100000.0
  endparam
switch:
  type = "pgd_HenonM"
  lya = lya
  lw = lw
  bailout = bailout
}

pgd_HenonMS {
; Henon Mandelbrot that samples multiple points. Passes
; Lyapunov information if desired, and can be restricted
; to particular periods.
init:
  float xx[@samples, @samples]
  float yy[@samples, @samples]
  float ls[@samples, @samples]
  float tx[@samples, @samples]
  float ty[@samples, @samples]
  int i = 0
  int j = 0
  int k = 0
  complex dxx[@samples, @samples]
  complex dxy[@samples, @samples]
  complex dyx[@samples, @samples]
  complex dyy[@samples, @samples]
  bool done[@samples, @samples]
  bool notdone = true
  float samplespacing = 0.0
  IF (@samples > 1)
    samplespacing = 2.0/(@samples - 1)
  ENDIF
  float fpx = 0.0
  float fpy = 0.0
  float rt = sqr(imag(#pixel)) - 2*imag(#pixel) + 1 + 4*real(#pixel)
  float rt2 = 4*real(#pixel) - 3*sqr(1 - imag(#pixel))
  IF (rt2 >= 0)
    fpx = (1 - imag(#pixel) + sqrt(rt2))/(2*real(#pixel))
    fpy = fpx*imag(#pixel)
  ELSEIF (rt >= 0)
    fpx = (imag(#pixel) - 1 + sqrt(rt))/(2*real(#pixel))
    fpy = fpx*imag(#pixel)
  ENDIF
  float scl = 1
  IF ((rt > 0) && (rt2 < 0))
    scl = 1/((sqr(real(#pixel)) + 1))
    IF (rt/2 < scl)
      scl = rt/2
    ENDIF
  ELSEIF (rt2 > 0)
     scl = scl/((10*sqr(imag(#pixel)) + 1))
;     scl = scl * (0.001 + (2*real(#pixel) + 3*imag(#pixel) - 2*sqr(imag(#pixel)) - 3))
  ENDIF
  WHILE (j < @samples)
    WHILE (k < @samples)
      xx[j, k] = -1.0 + j*samplespacing
      yy[j, k] = -1.0 + k*samplespacing
      IF (xx[j, k] < 0)
        xx[j, k] = -sqr(xx[j, k])
      ELSE
        xx[j, k] = sqr(xx[j, k])
      ENDIF
      IF (yy[j, k] < 0)
        yy[j, k] = -sqr(yy[j, k])
      ELSE
        yy[j, k] = sqr(yy[j, k])
      ENDIF
      xx[j, k] = scl*xx[j, k] + fpx + 1/@bailout
      yy[j, k] = scl*yy[j, k] + fpy + 1/@bailout
      ls[j, k] = sqr(xx[j, k])
      dxx[j, k] = 1
      dxy[j, k] = 0
      dyx[j, k] = 0
      dyy[j, k] = 1
      done[j, k] = false
      k = k + 1
    ENDWHILE
    j = j + 1
    k = 0
  ENDWHILE
  bool foundit = false
  int numtrapped = @samples*@samples
  float avgiters = 0.0
loop:
  j = 0
  k = 0
  notdone = false
  WHILE (j < @samples)
    WHILE (k < @samples)
      IF (!done[j, k])
        float t = 1 - real(#pixel)*ls[j, k] + yy[j, k]
        yy[j, k] = imag(#pixel)*xx[j, k]
        xx[j, k] = t
        ls[j, k] = t*t
        IF ((ls[j, k] + sqr(yy[j, k])) > @bailout)
          done[j, k] = true
          numtrapped = numtrapped - 1
          avgiters = avgiters + i
        ELSEIF (i == @preiters)
          tx[j, k] = xx[j, k]
          ty[j, k] = yy[j, k]
        ELSEIF (i > @preiters)
          IF ((sqr(tx[j, k] - xx[j, k]) + sqr(ty[j, k] - yy[j, k])) < 1/@bailout)
            done[j, k] = true
            IF (@what == 1)
              IF ((@per == i - @preiters) || ((@per < i - @preiters) && @lrgr))
                foundit = true
                notdone = false
                IF (@lya)
                  z = log(real(dxx[j, k]*dyy[j, k] - dxy[j, k]*dyx[j, k]))/i
                  IF (real(z) > 0)
                    foundit = false ; Looked periodic but the Lyapunov exponent tells a different tale.
                  ENDIF
                ENDIF
              ENDIF
            ENDIF
          ENDIF
        ENDIF
        IF (@lya)
          IF (!done[j, k])
            float d = -2*real(#pixel)*xx[j, k]
            complex nxx = dxx[j, k]*d + dyx[j, k]
            complex nxy = dxy[j, k]*d + dyy[j, k]
            dyx[j, k] = dxx[j, k]*imag(#pixel)
            dyy[j, k] = dxy[j, k]*imag(#pixel)
            dxx[j, k] = nxx
            dxy[j, k] = nxy
          ENDIF
        ENDIF
        IF (!done[j, k] && !foundit)
          notdone = true
        ENDIF
      ENDIF
      k = k + 1
    ENDWHILE
    j = j + 1
    k = 0
  ENDWHILE
  i = i + 1
  IF (!notdone)
    IF (@what == 0)
      IF (numtrapped == 0)
        foundit = true
        z = 1.5*log(avgiters/(@samples*@samples))/log(@repiters)
      ENDIF
    ENDIF
  ELSEIF (i == (#maxit - 1))
    notdone = false
    IF (@what == 2)
      j = 0
      k = 0
      float biggest = 0.0
      WHILE (j < @samples)
        WHILE (k < @samples)
          IF ((j == 0) && (k == 0))
            biggest = log(real(dxx[j, k]*dyy[j, k] - dxy[j, k]*dyx[j, k]))/i
          ELSE
            float q = log(real(dxx[j, k]*dyy[j, k] - dxy[j, k]*dyx[j, k]))/i
            IF (q > biggest)
              biggest = q
            ENDIF
          ENDIF
          k = k + 1
        ENDWHILE
        j = j + 1
        k = 0
      ENDWHILE
      foundit = (biggest > 0)
      z = biggest
    ENDIF
  ENDIF
  IF (!notdone && !foundit)
    z = (501,10)
  ENDIF
bailout:
  notdone
default:
  title = "Henon Mandelbrotlike Set (Multisampled)"
  center = (1,0)
  method = onepass
  periodicity = 0
  maxiter = 10000
  param samples
    caption = "Samples"
    hint = "More samples makes for a more accurate but slower picture."
    default = 4
  endparam
  param preiters
    caption = "Preiterations"
    hint = "Start looking for periodic attractor after this many iterations. \
            Higher is more accurate. Too high or too low slows down calculation."
    default = 100
  endparam
  param what
    caption = "Draw what"
    enum = "All diverge" "Periodic" "Chaotic"
    default = 0
  endparam
  param repiters
    caption = "Repeat iterations"
    hint = "Gradient repeats after this many iterations."
    default = 1000
    visible = (@what == 0)
  endparam
  param per
    caption = "Period"
    hint = "Period to draw"
    default = 1
    visible = (@what == 1)
  endparam
  param lrgr
    caption = "Draw all higher periods too"
    default = false
    visible = (@what == 1)
  endparam
  param lya
    caption = "Lyapunov"
    default = true
  endparam  
  param bailout
    caption = "Bailout"
    default = 100000.0
    min = 0.0
  endparam
switch:
  type = "pgd_HenonJ"
  seed = #pixel
  lya = lya
  bailout = bailout
}

pgd_ForcedLogistic {
  ; Forced Logistic mapping alternates in a pattern between
  ; two parameter values a and b for r in rx(1 - x).
  ; The pattern is specified in binary -- if you want to
  ; use abaab, for example, you put "10110" in the sequence
  ; field. In the sequence length field you put 5.
  ;
  ; Use PGD Lyapunov coloring for nice yummy smooth
  ; Lyapunov striping. 
init:
  float der = 1.0
  bool notdone = true
  int spow = 1
  float parm[@slength]
  int i = 0
  WHILE (i < @slength)
    int sseg = trunc(@sequence/spow)
    int sdig = sseg - 10*trunc(sseg/10)
    IF (sdig > 0)
      parm[i] = real(#pixel)
    ELSE
      parm[i] = imag(#pixel)
    ENDIF
    spow = spow*10
    i = i + 1
  ENDWHILE
  i = 0
  z = @start
  bool landedright = false
  float t = 0.0
loop:
  float r = parm[i % @slength]
  z = r*z*(1 - z)
  IF (@lya)
    der = der * r * (1 - 2*real(z))
  ENDIF
  IF (abs(real(z)) > @bailout)
    notdone = false
    landedright = (@what == 2)
    IF (@lya)
      z = log(der)/i - 0.01
      IF (der == 0)
        z = -0.75
      ENDIF
    ENDIF
  ELSEIF (i == @preiters)
    t = real(z)
  ELSEIF (i > @preiters)
    IF (abs(real(z) - t) < 1/@bailout)
      notdone = false
      IF (@what == 0)
        IF ((@per == i - @preiters) || ((@per < i - @preiters) && @lrgr))
          landedright = true
          IF (@lya)
            z = log(der)/i - 0.01
            IF (real(z) > 0)
              landedright = false
            ENDIF
            IF (der == 0)
              z = -0.75
            ENDIF
          ENDIF
        ENDIF
      ELSE
        landedright = false
      ENDIF
    ENDIF
  ENDIF
  i = i + 1
  IF ((i == (#maxit - 1)) && notdone)
    notdone = false
    landedright = (@what == 1)
    IF (@lya)
      z = log(der)/i - 0.01
      IF (real(z) <= 0)
        landedright = (@what == 0)
      ENDIF
      IF (der == 0)
        z = -0.75
      ENDIF
    ENDIF
  ENDIF
  IF (!notdone && !landedright)
    z = (501,10)
  ENDIF
bailout:
  notdone
default:
  title = "Forced Logistic Mapping"
  center = (2,2.5)
  method = onepass
  periodicity = 0
  maxiter = 10000
  param sequence
    caption = "Sequence"
    hint = "Enter it in binary -- e.g., if you want 'abaab' \
            you use 10110."
    default = 10
  endparam
  param slength
    caption = "Sequence Length"
    default = 2
  endparam
  param what
    caption = "Draw what"
    enum = "Periodic" "Chaotic" "All diverge"
    default = 0
  endparam
  param per
    caption = "Period"
    hint = "Period to draw"
    default = 1
    visible = (@what == 0)
  endparam
  param lrgr
    caption = "Draw all higher periods too"
    default = true
    visible = (@what == 0)
  endparam
  param start
    caption = "Starting point"
    default = 0.5
  endparam
  param lya
    caption = "Lyapunov"
    default = true
  endparam  
  param preiters
    caption = "Preiterations"
    hint = "Start looking for periodic attractor after this many iterations. \
            Higher is more accurate. Too high or too low slows down calculation."
    default = 100
  endparam
  param bailout
    caption = "Bailout"
    default = 100000.0
    min = 0.0
  endparam
switch:
  type = "pgd_ForcedLogisticJ"
  sequence = sequence
  slength = slength
  bailout = bailout
  seed = #pixel
}

pgd_ForcedLogisticJ {
  ; Forced Logistic mapping alternates in a pattern between
  ; two parameter values a and b for r in rx(1 - x).
  ; The pattern is specified in binary -- if you want to
  ; use abaab, for example, you put "10110" in the sequence
  ; field. In the sequence length field you put 5.
  ;
  ; The Julia sets are rendered by letting x be a complex
  ; number.
init:
  int spow = 1
  float parm[@slength]
  int i = 0
  WHILE (i < @slength)
    int sseg = trunc(@sequence/spow)
    int sdig = sseg - 10*trunc(sseg/10)
    IF (sdig > 0)
      parm[i] = real(@seed)
    ELSE
      parm[i] = imag(@seed)
    ENDIF
    spow = spow*10
    i = i + 1
  ENDWHILE
  z = #pixel
  i = 0
loop:
  z = parm[i%@slength]*z*(1 - z)
  i = i + 1
bailout:
  |z| < @bailout
default:
  title = "Forced Logistic Mapping (Julia sets)"
  center = (0.5,0)
  param sequence
    caption = "Sequence"
    hint = "Enter it in binary -- e.g., if you want 'abaab' \
            you use 10110."
    default = 10
  endparam
  param slength
    caption = "Sequence Length"
    default = 2
  endparam
  param seed
    caption = "Julia seed"
    default = (3.514, 3.213)
  endparam
  param bailout
    caption = "Bailout"
    default = 100000.0
    min = 0.0
  endparam
switch:
  type = "pgd_ForcedLogistic"
  sequence = sequence
  slength = slength
  bailout = bailout
}

pgd_VLotkaJ {
  ; Another system in two real variables. Like Henon, there's
  ; all sorts of periodicities possible, as well as strange
  ; attractors; unlike Henon, non-chaotic aperiodic orbits
  ; (which tend to fill in a simple curve, and in which
  ; pairs of close-by points neither get far apart as on a
  ; strange attractor nor draw closer) are possible too.
  ;
  ; Check "Lyapunov" to pass Lyapunov information to the
  ; PGD Lyapunov coloring algorithm in z.
  ;
  ; Possible sections to select and color include escaping
  ; points, points that go to periodic points, chaotic
  ; points, and non-chaotic aperiodic points. To distinguish
  ; the latter the formula will perform lyapunov calculations
  ; even if "Lyapunov" is unchecked if you select either of
  ; those; Lyapunov values close to zero will be taken to
  ; indicate that an aperiodic orbit is not chaotic, while
  ; values significantly above zero will be taken to
  ; indicate chaos.
  ;
  ; Outside coloring suggested: smoothed iterations.
  ; Inside coloring sugegsted: "Closest Approach of Attractor
  ; to Point of Origin", especially for aperiodic attractor
  ; basins. This option, with suitable settings, makes the
  ; attractor itself visible.
  ;
  ; Volterra-Lotka differential equations:
  ; dx/dt = f(x,y) = x - xy
  ; dy/dt = g(x,y) = -y + xy
  ; Our discretization (also used by Fractint):
  ; x(new) = x + h/2 * [ f(x,y) + f[x + pf(x,y), y + pg(x,y)] ]
  ; y(new) = y + h/2 * [ g(x,y) + g[x + pf(x,y), y + pg(x,y)] ]  
  ; f(x,y) and g(x,y) are as above; h and p are the parameters
  ; of the discretization. For this type, they constitute
  ; the Julia seed.
init:
  z = #pixel
  bool notdone = true
  bool landedright = false
  int per = @period
  IF ((@what == 2) || (@what == 3))
    per = 1
  ENDIF
  complex z_values[@period]
  int i = 0
  float dxx = 1 ; Components of Jacobian of nth iterate
  float dxy = 0 ; wind up in here.
  float dyx = 0 ; We start with the identity matrix.
  float dyy = 1
  int stopwi = trunc(#maxit/10)
  IF (stopwi < per)
    stopwi = per
  ENDIF
  IF (@dtype)
    stopwi = @preit
  ENDIF
  float h = real(@seed)/2
  float p = imag(@seed)
  float pp2 = p + 2
  float pm2 = p - 2
  float pp1 = p + 1
  float pm1 = p - 1
  float pm2pp1 = pm2*pp1
  float pp21mp = -(pp2*pm1)
  float p1mp = -p*pm1
  float p1pp = p*pp1
  float tp1mp = 2*p1mp
  float tp1pp = 2*p1pp
  float tp2 = 2*sqr(p)
  bool lya2 = (@lya || (@what == 2) || (@what == 3))
  complex tar = 0
loop:
  float qx = real(z)
  float qy = imag(z)
  float xy = qx*qy
  float f = qx - xy
  float g = xy - qy
  float xx = qx + p*f
  float yy = qy + p*g
  float xyzzy = xx*yy
  float ff = xx - xyzzy
  float gg = xyzzy - yy
  z = z + h*(f + ff + (0,1)*(g + gg))
  IF (lya2 && (!@dtype || (i > stopwi)))
    ; Constructing Lyapunov information
    ; We multiply our Jacobian matrix on the left by the
    ; current iteration's Jacobian.
    ;
    ; xx = 1 + h/2 * [ (p + 2) + (p - 2)(p + 1)y - 2p(1 + p)xy + p(1 - p)y^2 + 2p^2xy^2 ]
    ; xy =     h/2 * [           (p - 2)(p + 1)x - p(1 + p)x^2 + 2p(1 - p)xy + 2p^2x^2y ]
    ; yx =     h/2 * [ (p + 2)(1 - p)y           + 2p(1 + p)xy + p(p - 1)y^2 - 2p^2xy^2 ]
    ; yy = 1 + h/2 * [ (p + 2)(1 - p)x + (p - 2) + p(1 + p)x^2 + 2p(p - 1)xy - 2p^2x^2y ]
    ;
    float x2 = sqr(qx)
    float y2 = sqr(qy)
    float x2y = qy*x2
    float xy2 = qx*y2
    float cxx = 1 + h*(pp2 + pm2pp1*qy - tp1pp*xy + p1mp*y2 + tp2*xy2)
    float cxy = h*(pm2pp1*qx - p1pp*x2 + tp1mp*xy + tp2*x2y)
    float cyx = h*(pp21mp*qy + tp1pp*xy - p1mp*y2 - tp2*xy2)
    float cyy = 1 + h*(pp21mp*qx + pm2 + p1pp*x2 - tp1mp*xy - tp2*x2y)
    float nxx = cxx*dxx + cxy*dyx
    float nxy = cxx*dxy + cxy*dyy
    float nyx = cyx*dxx + cyy*dyx
    dyy = cyx*dxy + cyy*dyy
    dxx = nxx
    dxy = nxy
    dyx = nyx
  ENDIF
  IF (|z| > @bailout)
    notdone = false
    landedright = (@what == 0)
  ELSEIF (@what != 0)
    IF (@dtype)
      IF (i == stopwi)
        tar = z
      ELSEIF (i > stopwi)
        IF (|z - tar| < 1/@bailout)
          int j = i - stopwi
          IF (j < @maxp)
            notdone = false
            landedright = false
            IF (@what == 1)
              IF (j == per)
                landedright = true
              ELSEIF ((j > per) && !@exact)
                landedright = true
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF
    ELSE
      IF (i < stopwi)
        IF (i >= per)
          IF (|z - z_values[i % per]| < 1/@bailout)
            notdone = false
            IF (@what == 1)
              landedright = true
              int j = 1
              WHILE (j < per)
                IF(|z - z_values[(i + j)%per]| < @tol/@bailout)
                  landedright = false
                ENDIF
                j = j + 1
              ENDWHILE
            ENDIF
          ENDIF
        ENDIF
        z_values[i % per] = z
      ELSE
        IF ((|z - z_values[i % per]| < 1/@bailout) && ((i - stopwi) < @maxp))
          notdone = false
          IF ((@what == 1) && !@exact)
            landedright = true
            j = 1
            WHILE (j < per)
              IF(|z - z_values[(i + j)%per]| < @tol/@bailout)
                landedright = false
              ENDIF
              j = j + 1
            ENDWHILE
          ENDIF
        ENDIF
      ENDIF
    ENDIF
    i = i + 1
    IF (i == (#maxit - 1))
      IF (@dtype)
        IF (i < stopwi)
          i = 0
        ELSE
          i = i - stopwi
        ENDIF
      ENDIF
      IF ((@what == 0) && @lya)
        ; Pass Lyapunov information in z
        IF (@lw == 2)
          ; We want the larger eigenvalue of the Jacobian.
          complex alpha = dxx + dyy
          complex beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(real(alpha + beta)/2)/i
        ELSEIF (@lw == 3)
          ; We want the smaller eigenvalue.
          alpha = dxx + dyy
          beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(real(alpha - beta)/2)/i
        ELSEIF(@lw == 0)
          ; We want the determinant.
          z = log(real(dxx*dyy - dxy*dyx))/i
        ELSEIF(@lw == 1)
          ; We want the ratio of the eigenvalues.
          alpha = dxx + dyy
          beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(real(alpha + beta)/real(alpha - beta))/i
        ELSEIF(@lw == 4)
          ; We want the imaginary part.
          beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(imag(beta))/i
        ENDIF
      ELSE
        notdone = false
        IF ((@what == 0) || (@what == 1))
          landedright = false
        ELSE
          float ly = log(real(dxx*dyy - dxy*dyx))/i
          IF (ly > @chaostol)
            landedright = (@what == 3)
          ELSE
            landedright = (@what == 2)
          ENDIF
        ENDIF
      ENDIF
    ENDIF
  ENDIF
  IF (!notdone && @lya)
    ; Pass Lyapunov information in z
    IF (@dtype)
      IF (i < stopwi)
        i = 0
      ELSE
        i = i - stopwi
      ENDIF
    ENDIF
    ly = log(real(dxx*dyy - dxy*dyx))/i
    IF (@what == 1)
      IF (ly > @chaostol)
        landedright = false ; Chaos alert!
      ENDIF
    ENDIF
    IF (@lw == 2)
      ; We want the larger eigenvalue of the Jacobian.
      complex alpha = dxx + dyy
      complex beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(real(alpha + beta)/2)/i
    ELSEIF (@lw == 3)
      ; We want the smaller eigenvalue.
      alpha = dxx + dyy
      beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(real(alpha - beta)/2)/i
    ELSEIF(@lw == 0)
      ; We want the determinant.
      z = ly
    ELSEIF(@lw == 1)
      ; We want the ratio of the eigenvalues.
      alpha = dxx + dyy
      beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(real(alpha + beta)/real(alpha - beta))/i
    ELSEIF(@lw == 4)
      ; We want the imaginary part.
      beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(imag(beta))/i
    ENDIF
  ENDIF  
  IF (!notdone && !landedright)
    z = (501,10)
  ENDIF
bailout:
  notdone
default:
  title = "Discrete Volterra-Lotka Julia Sets"
  center = (3,2.25)
  magn = 0.66666666
  method = onepass
  maxiter = 1000
  periodicity = 0
  param seed
    caption = "Julia seed"
    default = (1.4,0.3)
  endparam
  param what
    caption = "Draw what"
    enum = "Outside" "Periodic" "Limit Circle" "Ergodic"
    default = 0
  endparam
  param period
    caption = "Period"
    default = 2
    visible = (@what == 1)
  endparam
  param exact
    caption = "Exact"
    default = false
    hint = "Turned on, only the exact period is colored. \
            Turned off, multiples of the period are colored \
            (regular detection type) or all periods equal or \
            greater (alternate detection type)."
    visible = (@what == 1)
  endparam
  param dtype
    caption = "Use alternate periodicity detection"
    default = false
    visible = (@what != 0)
  endparam
  param preit
    caption = "Preiterations"
    default = 100
    min = 1
    visible = ((@what != 0) && @dtype)
  endparam
  param tol
    caption = "Tolerance"
    default = 100.0
    hint = "Tolerance when checking not really a shorter \
            period."
    visible = ((@what == 1) && !@dtype)
  endparam
  param maxp
    caption = "Max Period"
    default = 50
    min = 1
    hint = "If a period is detected that is longer than this, \
            it is assumed the orbit is really aperiodic but \
            bounded, and has settled to its attractor."
    visible = (@what != 0)
  endparam
  param chaostol
    caption = "Chaos tolerance"
    default = 0.01
    min = 0.0
    hint = "Aperiodic orbits with a Lyapunov exponent estimated \
            to exceed this value are classified as chaotic."
    visible = ((@what == 2) || (@what == 3))
  endparam
  param lya
    caption = "Lyapunov"
    default = false
  endparam
  param lw
    caption = "Show what"
    enum = "Expansivity" "Stretch" "Larger" "Smaller" "Twist"
    default = 0
    visible = @lya
  endparam
  param bailout
    caption = "Bailout"
    default = 100000.0
  endparam
switch:
  type = "pgd_VLotkaM"
  lya = lya
  lw = lw
  bailout = bailout
}

pgd_VLotkaM {
  ; Another system in two real variables. Like Henon, there's
  ; all sorts of periodicities possible, as well as strange
  ; attractors; unlike Henon, non-chaotic aperiodic orbits
  ; (which tend to fill in a simple curve, and in which
  ; pairs of close-by points neither get far apart as on a
  ; strange attractor nor draw closer) are possible too.
  ;
  ; Check "Lyapunov" to pass Lyapunov information to the
  ; PGD Lyapunov coloring algorithm in z.
  ;
  ; Possible sections to select and color include escaping
  ; points, points that go to periodic points, chaotic
  ; points, and non-chaotic aperiodic points. To distinguish
  ; the latter the formula will perform lyapunov calculations
  ; even if "Lyapunov" is unchecked if you select either of
  ; those; Lyapunov values close to zero will be taken to
  ; indicate that an aperiodic orbit is not chaotic, while
  ; values significantly above zero will be taken to
  ; indicate chaos.
  ;
  ; Suggested use: one outside layer (escaping points),
  ; another escaping points layer with inside colored by
  ; closest approach of attractor to point of origin,
  ; and tinting layers for periodic points (one, or several
  ; for different periods), aperiodic points, and chaotic
  ; points. Inside coloring by PGD Lyapunov (check the
  ; "Lyapunov" box) is also good.
  ;
  ; Use an exponent of 4 with outside smoothed iterations
  ; coloring.
  ;
  ; Volterra-Lotka differential equations:
  ; dx/dt = f(x,y) = x - xy
  ; dy/dt = g(x,y) = -y + xy
  ; Our discretization (also used by Fractint):
  ; x(new) = x + h/2 * [ f(x,y) + f[x + pf(x,y), y + pg(x,y)] ]
  ; y(new) = y + h/2 * [ g(x,y) + g[x + pf(x,y), y + pg(x,y)] ]  
  ; f(x,y) and g(x,y) are as above; h and p are the parameters
  ; of the discretization. For this type, they are
  ; determined by the pixel coordinates.
init:
  z = @start
  bool notdone = true
  bool landedright = false
  int per = @period
  IF ((@what == 2) || (@what == 3))
    per = 1
  ENDIF
  complex z_values[@period]
  int i = 0
  float dxx = 1 ; Components of Jacobian of nth iterate
  float dxy = 0 ; wind up in here.
  float dyx = 0 ; We start with the identity matrix.
  float dyy = 1
  int stopwi = trunc(#maxit/10)
  IF (stopwi < per)
    stopwi = per
  ENDIF
  IF (@dtype)
    stopwi = @preit
  ENDIF
  float h = real(#pixel)/2
  float p = imag(#pixel)
  float pp2 = p + 2
  float pm2 = p - 2
  float pp1 = p + 1
  float pm1 = p - 1
  float pm2pp1 = pm2*pp1
  float pp21mp = -(pp2*pm1)
  float p1mp = -p*pm1
  float p1pp = p*pp1
  float tp1mp = 2*p1mp
  float tp1pp = 2*p1pp
  float tp2 = 2*sqr(p)
  bool lya2 = (@lya || (@what == 2) || (@what == 3))
  complex tar = 0
loop:
  float qx = real(z)
  float qy = imag(z)
  float xy = qx*qy
  float f = qx - xy
  float g = xy - qy
  float xx = qx + p*f
  float yy = qy + p*g
  float xyzzy = xx*yy
  float ff = xx - xyzzy
  float gg = xyzzy - yy
  z = z + h*(f + ff + (0,1)*(g + gg))
  IF (lya2 && (!@dtype || (i > stopwi)))
    ; Constructing Lyapunov information
    ; We multiply our Jacobian matrix on the left by the
    ; current iteration's Jacobian.
    ;
    ; xx = 1 + h/2 * [ (p + 2) + (p - 2)(p + 1)y - 2p(1 + p)xy + p(1 - p)y^2 + 2p^2xy^2 ]
    ; xy =     h/2 * [           (p - 2)(p + 1)x - p(1 + p)x^2 + 2p(1 - p)xy + 2p^2x^2y ]
    ; yx =     h/2 * [ (p + 2)(1 - p)y           + 2p(1 + p)xy + p(p - 1)y^2 - 2p^2xy^2 ]
    ; yy = 1 + h/2 * [ (p + 2)(1 - p)x + (p - 2) + p(1 + p)x^2 + 2p(p - 1)xy - 2p^2x^2y ]
    ;
    float x2 = sqr(qx)
    float y2 = sqr(qy)
    float x2y = qy*x2
    float xy2 = qx*y2
    float cxx = 1 + h*(pp2 + pm2pp1*qy - tp1pp*xy + p1mp*y2 + tp2*xy2)
    float cxy = h*(pm2pp1*qx - p1pp*x2 + tp1mp*xy + tp2*x2y)
    float cyx = h*(pp21mp*qy + tp1pp*xy - p1mp*y2 - tp2*xy2)
    float cyy = 1 + h*(pp21mp*qx + pm2 + p1pp*x2 - tp1mp*xy - tp2*x2y)
    float nxx = cxx*dxx + cxy*dyx
    float nxy = cxx*dxy + cxy*dyy
    float nyx = cyx*dxx + cyy*dyx
    dyy = cyx*dxy + cyy*dyy
    dxx = nxx
    dxy = nxy
    dyx = nyx
  ENDIF
  IF (|z| > @bailout)
    notdone = false
    landedright = (@what == 0)
  ELSEIF (@what != 0)
    IF (@dtype)
      IF (i == stopwi)
        tar = z
      ELSEIF (i > stopwi)
        IF (|z - tar| < 1/@bailout)
          int j = i - stopwi
          IF (j < @maxp)
            notdone = false
            landedright = false
            IF (@what == 1)
              IF (j == per)
                landedright = true
              ELSEIF ((j > per) && !@exact)
                landedright = true
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF
    ELSE
      IF (i < stopwi)
        IF (i >= per)
          IF (|z - z_values[i % per]| < 1/@bailout)
            notdone = false
            IF (@what == 1)
              landedright = true
              int j = 1
              WHILE (j < per)
                IF(|z - z_values[(i + j)%per]| < @tol/@bailout)
                  landedright = false
                ENDIF
                j = j + 1
              ENDWHILE
            ENDIF
          ENDIF
        ENDIF
        z_values[i % per] = z
      ELSE
        IF ((|z - z_values[i % per]| < 1/@bailout) && ((i - stopwi) < @maxp))
          notdone = false
          IF ((@what == 1) && !@exact)
            landedright = true
            j = 1
            WHILE (j < per)
              IF(|z - z_values[(i + j)%per]| < @tol/@bailout)
                landedright = false
              ENDIF
              j = j + 1
            ENDWHILE
          ENDIF
        ENDIF
      ENDIF
    ENDIF
    i = i + 1
    IF (i == (#maxit - 1))
      IF (@dtype)
        IF (i < stopwi)
          i = 0
        ELSE
          i = i - stopwi
        ENDIF
      ENDIF
      IF ((@what == 0) && @lya)
        ; Pass Lyapunov information in z
        IF (@lw == 2)
          ; We want the larger eigenvalue of the Jacobian.
          complex alpha = dxx + dyy
          complex beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(real(alpha + beta)/2)/i
        ELSEIF (@lw == 3)
          ; We want the smaller eigenvalue.
          alpha = dxx + dyy
          beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(real(alpha - beta)/2)/i
        ELSEIF(@lw == 0)
          ; We want the determinant.
          z = log(real(dxx*dyy - dxy*dyx))/i
        ELSEIF(@lw == 1)
          ; We want the ratio of the eigenvalues.
          alpha = dxx + dyy
          beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(real(alpha + beta)/real(alpha - beta))/i
        ELSEIF(@lw == 4)
          ; We want the imaginary part.
          beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          z = log(imag(beta))/i
        ENDIF
      ELSE
        notdone = false
        IF ((@what == 0) || (@what == 1))
          landedright = false
        ELSE
          float ly = log(real(dxx*dyy - dxy*dyx))/i
          $define DEBUG
          alpha = dxx + dyy
          beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
          print (real(log(real(alpha + beta)/2)/i), "   ", real(log(real(alpha - beta)/2)/i), "   ", ly)
          IF (ly > @chaostol)
            landedright = (@what == 3)
          ELSE
            landedright = (@what == 2)
          ENDIF
        ENDIF
      ENDIF
    ENDIF
  ENDIF
  IF (!notdone && @lya)
    ; Pass Lyapunov information in z
    IF (@dtype)
      IF (i < stopwi)
        i = 0
      ELSE
        i = i - stopwi
      ENDIF
    ENDIF
    ly = log(real(dxx*dyy - dxy*dyx))/i
    IF (@what == 1)
      IF (ly > @chaostol)
        landedright = false ; Chaos alert!
      ENDIF
    ENDIF
    IF (@lw == 2)
      ; We want the larger eigenvalue of the Jacobian.
      complex alpha = dxx + dyy
      complex beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(real(alpha + beta)/2)/i
    ELSEIF (@lw == 3)
      ; We want the smaller eigenvalue.
      alpha = dxx + dyy
      beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(real(alpha - beta)/2)/i
    ELSEIF(@lw == 0)
      ; We want the determinant.
      z = ly
    ELSEIF(@lw == 1)
      ; We want the ratio of the eigenvalues.
      alpha = dxx + dyy
      beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(real(alpha + beta)/real(alpha - beta))/i
    ELSEIF(@lw == 4)
      ; We want the imaginary part.
      beta = sqrt(dyy*dyy - 2*dxx*dyy + dxx*dxx + 4*dxy*dyx)
      z = log(imag(beta))/i
    ENDIF
  ENDIF  
  IF (!notdone && !landedright)
    z = (501,10)
  ENDIF
bailout:
  notdone
default:
  title = "Discrete Volterra-Lotka MSet"
  center = (2,1.5)
  method = onepass
  maxiter = 1000
  periodicity = 0
  param start
    caption = "Starting (x,y)"
    hint = "Points close to the fixed point (1,1) are good."
    default = (1.1,1.1)
  endparam
  param what
    caption = "Draw what"
    enum = "Outside" "Periodic" "Limit Circle" "Ergodic"
    default = 0
  endparam
  param period
    caption = "Period"
    default = 2
    visible = (@what == 1)
  endparam
  param exact
    caption = "Exact"
    default = false
    hint = "Turned on, only the exact period is colored. \
            Turned off, multiples of the period are colored \
            (regular detection type) or all periods equal or \
            greater (alternate detection type)."
    visible = (@what == 1)
  endparam
  param dtype
    caption = "Use alternate periodicity detection"
    default = false
    visible = (@what != 0)
  endparam
  param preit
    caption = "Preiterations"
    default = 100
    min = 1
    visible = ((@what != 0) && @dtype)
  endparam
  param tol
    caption = "Tolerance"
    default = 100.0
    hint = "Tolerance when checking not really a shorter \
            period."
    visible = ((@what == 1) && !@dtype)
  endparam
  param maxp
    caption = "Max Period"
    default = 50
    min = 1
    hint = "If a period is detected that is longer than this, \
            it is assumed the orbit is really aperiodic but \
            bounded, and has settled to its attractor."
    visible = (@what != 0)
  endparam
  param chaostol
    caption = "Chaos tolerance"
    default = 0.01
    min = 0.0
    hint = "Aperiodic orbits with a Lyapunov exponent estimated \
            to exceed this value are classified as chaotic."
    visible = ((@what == 2) || (@what == 3))
  endparam
  param lya
    caption = "Lyapunov"
    default = false
  endparam
  param lw
    caption = "Show what"
    enum = "Expansivity" "Stretch" "Larger" "Smaller" "Twist"
    default = 0
    visible = @lya
  endparam
  param bailout
    caption = "Bailout"
    default = 100000.0
  endparam
switch:
  type = "pgd_VLotkaJ"
  lya = lya
  lw = lw
  bailout = bailout
  seed = #pixel
}

pgd_VLotkaMS {
  ; Another system in two real variables. Like Henon, there's
  ; all sorts of periodicities possible, as well as strange
  ; attractors; unlike Henon, non-chaotic aperiodic orbits
  ; (which tend to fill in a simple curve, and in which
  ; pairs of close-by points neither get far apart as on a
  ; strange attractor nor draw closer) are possible too.
  ;
  ; Check "Lyapunov" to pass Lyapunov information to the
  ; PGD Lyapunov coloring algorithm in z.
  ;
  ; This version samples multiple points near (1,1) and
  ; combines information from them.
  ;
  ; Possible sections to select and color include parameter
  ; points where all samples escape, parameter points where
  ; at least one sample is periodic (and you can select
  ; specific periods, or a specific period and all larger
  ; periods), parameter points where at least one sample
  ; has a limit circle attractor, and parameter points where
  ; at least one sample has a strange attractor. The latter
  ; are distinguished from one another via Lyapunov
  ; exponent calculations (whether or not "Lyapunov" is
  ; checked) and from periodic attractors by not having a
  ; short period before returning to some neighborhood of
  ; an earlier point (smaller with higher bailout).
  ;
  ; Lyapunov coloring is recommended. For periodic, chaotic,
  ; and limit circle points if "Lyapunov" is checked the
  ; exponent handed to the PGD Lyapunov coloring method is
  ; from a sample that met the criterion -- the right
  ; periodicity, or a chaotic attractor, or a limit circle
  ; attractor caught that sample point. Parameter points
  ; where all samples escape have the averaged number of
  ; iterations stuffed in z for PGD Lyapunov to use -- the
  ; result isn't "really" Lyapunov coloring but it gives
  ; a nice "soft iterations" effect outside the Mandelbrot.
  ;
  ; Volterra-Lotka differential equations:
  ; dx/dt = f(x,y) = x - xy
  ; dy/dt = g(x,y) = -y + xy
  ; Our discretization (also used by Fractint):
  ; x(new) = x + h/2 * [ f(x,y) + f[x + pf(x,y), y + pg(x,y)] ]
  ; y(new) = y + h/2 * [ g(x,y) + g[x + pf(x,y), y + pg(x,y)] ]  
  ; f(x,y) and g(x,y) are as above; h and p are the parameters
  ; of the discretization. For this type, they are
  ; determined by the pixel coordinates.
init:
  float samplespacing = 0.0
  IF (@samples > 1)
    samplespacing = 2.0/(@samples - 1)
  ENDIF
  complex w[@samples, @samples]
  complex tar[@samples, @samples]
  bool done[@samples, @samples]
  float dxx[@samples, @samples]
  float dxy[@samples, @samples]
  float dyx[@samples, @samples]
  float dyy[@samples, @samples]
  int j = 0
  int k
  WHILE (j < @samples)
    k = 0
    WHILE (k < @samples)
      w[j, k] = 0.1*(samplespacing*j + (0,1)*samplespacing*k) + (0.901,0.901)
      done[j, k] = false
      dxx[j, k] = 1
      dxy[j, k] = 0
      dyx[j, k] = 0
      dyy[j, k] = 1
      k = k + 1
    ENDWHILE
    j = j + 1
  ENDWHILE
  bool notdone = true
  bool landedright = false
  int per = @period
  IF ((@what == 2) || (@what == 3))
    per = 1
  ENDIF
  int i = 0
  float h = real(#pixel)/2
  float p = imag(#pixel)
  float pp2 = p + 2
  float pm2 = p - 2
  float pp1 = p + 1
  float pm1 = p - 1
  float pm2pp1 = pm2*pp1
  float pp21mp = -(pp2*pm1)
  float p1mp = -p*pm1
  float p1pp = p*pp1
  float tp1mp = 2*p1mp
  float tp1pp = 2*p1pp
  float tp2 = 2*sqr(p)
  bool lya2 = (@lya || (@what == 2) || (@what == 3))
  float avgiters = 0.0
  int numtrapped = @samples*@samples
loop:
  notdone = false
  landedright = false
  j = 0
  WHILE (j < @samples)
    k = 0
    WHILE (k < @samples)
      IF (!done[j, k])
        float qx = real(w[j, k])
        float qy = imag(w[j, k])
        float xy = qx*qy
        float f = qx - xy
        float g = xy - qy
        float xx = qx + p*f
        float yy = qy + p*g
        float xyzzy = xx*yy
        float ff = xx - xyzzy
        float gg = xyzzy - yy
        w[j, k] = w[j, k] + h*(f + ff + (0,1)*(g + gg))
        IF (lya2 && (i > @preiters))
          ; Constructing Lyapunov information
          ; We multiply our Jacobian matrix on the left by the
          ; current iteration's Jacobian.
          float x2 = sqr(qx)
          float y2 = sqr(qy)
          float x2y = qy*x2
          float xy2 = qx*y2
          float cxx = 1 + h*(pp2 + pm2pp1*qy - tp1pp*xy + p1mp*y2 + tp2*xy2)
          float cxy = h*(pm2pp1*qx - p1pp*x2 + tp1mp*xy + tp2*x2y)
          float cyx = h*(pp21mp*qy + tp1pp*xy - p1mp*y2 - tp2*xy2)
          float cyy = 1 + h*(pp21mp*qx + pm2 + p1pp*x2 - tp1mp*xy - tp2*x2y)
          float nxx = cxx*dxx[j, k] + cxy*dyx[j, k]
          float nxy = cxx*dxy[j, k] + cxy*dyy[j, k]
          float nyx = cyx*dxx[j, k] + cyy*dyx[j, k]
          dyy[j, k] = cyx*dxy[j, k] + cyy*dyy[j, k]
          dxx[j, k] = nxx
          dxy[j, k] = nxy
          dyx[j, k] = nyx
        ENDIF
        IF (|w[j, k]| > @bailout)
          done[j, k] = true
          avgiters = avgiters + i
          numtrapped = numtrapped - 1
        ELSEIF (i == @preiters)
          tar[j, k] = w[j, k]
        ELSEIF (i > @preiters)
          IF (|w[j, k] - tar[j, k]| < 1/@bailout)
            int l = i - @preiters
            IF (l < @maxp)
              done[j, k] = true
              IF (@what == 1)
                IF (l == per)
                  landedright = true
                ELSEIF ((l > per) && !@exact)
                  landedright = true
                ENDIF
                IF (landedright && @lya)
                  ; Pass Lyapunov information in z
                  IF (i < @preiters)
                    i = 0
                  ELSE
                    i = i - @preiters
                  ENDIF
                  z = log(real(dxx[j, k]*dyy[j, k] - dxy[j, k]*dyx[j, k]))/i
                ENDIF
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      IF (!done[j, k] && !landedright)
        notdone = true
      ENDIF
      k = k + 1
    ENDWHILE
    j = j + 1
  ENDWHILE
  i = i + 1
  IF ((@what == 0) && !notdone)
    IF (numtrapped > 0)
      landedright = false
    ELSE
      landedright = true
      z = 1.5*log(avgiters/(@samples*@samples))/log(@repiters)
    ENDIF
  ELSEIF ((i == (#maxit - 1)) && notdone)
    notdone = false
    IF (i < @preiters)
      i = 0
    ELSE
      i = i - @preiters
    ENDIF
    IF ((@what == 0) || (@what == 1))
      landedright = false
    ELSEIF ((@what == 2) || (@what == 3))
      landedright = false
      float biggest = 0
      j = 0
      WHILE (j < @samples)
        k = 0
        WHILE (k < @samples)
          float ly = log(real(dxx[j, k]*dyy[j, k] - dxy[j, k]*dyx[j, k]))/i
          IF (ly > biggest)
            biggest = ly
          ENDIF
          IF ((ly > @chaostol) && (@what == 3))
            landedright = true
          ELSEIF ((ly <= @chaostol) && (@what == 2) && (!done[j, k]))
            landedright = true
          ENDIF
          k = k + 1
        ENDWHILE
        j = j + 1
      ENDWHILE
      z = biggest
    ENDIF
  ENDIF
  IF (!notdone && !landedright)
    z = (501,10)
  ENDIF
bailout:
  notdone
default:
  title = "Discrete Volterra-Lotka MSet (Multisampled)"
  center = (2,1.5)
  method = onepass
  maxiter = 1000
  periodicity = 0
  param samples
    caption = "Samples"
    default = 4
    min = 1
  endparam
  param what
    caption = "Draw what"
    enum = "All escape" "Periodic" "Limit Circle" "Ergodic"
    default = 0
  endparam
  param repiters
    caption = "Repeat iterations"
    hint = "Gradient repeats after this many iterations."
    default = 1000
    visible = (@what == 0)
  endparam
  param period
    caption = "Period"
    default = 2
    visible = (@what == 1)
  endparam
  param exact
    caption = "Exact"
    default = false
    hint = "Turned on, only the exact period is colored. \
            Turned off, multiples of the period are colored \
            (regular detection type) or all periods equal or \
            greater (alternate detection type)."
    visible = (@what == 1)
  endparam
  param preiters
    caption = "Preiterations"
    default = 100
    min = 1
    visible = (@what != 0)
  endparam
  param maxp
    caption = "Max Period"
    default = 50
    min = 1
    hint = "If a period is detected that is longer than this, \
            it is assumed the orbit is really aperiodic but \
            bounded, and has settled to its attractor."
    visible = (@what != 0)
  endparam
  param chaostol
    caption = "Chaos tolerance"
    default = 0.01
    min = 0.0
    hint = "Aperiodic orbits with a Lyapunov exponent estimated \
            to exceed this value are classified as chaotic."
    visible = ((@what == 2) || (@what == 3))
  endparam
  param lya
    caption = "Lyapunov"
    default = false
  endparam
  param bailout
    caption = "Bailout"
    default = 100000.0
  endparam
switch:
  type = "pgd_VLotkaJ"
  lya = lya
  lw = lw
  bailout = bailout
  seed = #pixel
}

pgd_ComplexHenonJulia {
; The Henon mapping with both variables complex. Parameters
; allow varying the slice of the 4d dynamic space shown,
; and whether escaping/nonescaping points or points of
; certain periodicities are shown. Closest approach of
; attractor to point of origin coloring may show a
; projection of the finite attractor onto the plane, and
; there are many other ways to get interesting coloring
; too.
init:
  complex w
  float v = real(#pixel)
  float vv = imag(#pixel)
  float pp = real(@offset)
  float q = imag(@offset)
  float r = @t*#pi/2
  float s = @u*#pi/2
  float ss = @uu*#pi/2
  z = v*cos(r) + (0,1)*vv*cos(s) + pp*sin(r) + (0,1)*q*sin(s)
  w = pp*cos(r) + (0,1)*q*cos(s) - v*sin(r) - (0,1)*vv*sin(s)
  w = w*exp((0,1)*ss)
  IF (abs(imag(z)) < 1/@bailout)
    z = real(z)
  ENDIF
  IF (abs(imag(w)) < 1/@bailout)
    w = real(w)
  ENDIF
  int i = 0
  bool notdone = true
  bool landedright = false
  complex t1 = 0
  complex t2 = 0
  complex iz = z
  complex iw = w
  float d = 0
loop:
  t = @seed2*z
  z = 1 - @seed1*sqr(z) + w
  w = t
  IF ((|z| + |w|) > @bailout)
    notdone = false
    landedright = (@what == 0)
  ELSEIF (@what != 0)
    IF (i == @preit)
      t1 = z
      t2 = w
      IF (@catt)
        d = |z - iz| + |w - iw|
      ENDIF
    ELSEIF (i > @preit)
      IF (@catt)
        float dd = |z - iz| + |w - iw|
        IF (dd < d)
          d = dd
        ENDIF
      ENDIF
      IF ((|z - t1| + |w - t2|) < (1/@bailout))
        notdone = false
        landedright = (@what == 1)
        int p = i - @preit
        IF (p < @per)
          landedright = false
        ELSEIF ((@maxp > 0) && (p > @maxp))
          landedright = (@what == 2)
        ELSEIF ((!@lrgr) && (p > @per))
          landedright = false
        ENDIF
      ENDIF
    ENDIF
    i = i + 1
    IF (i == (#maxit - 1))
      notdone = false
      landedright = (@what == 2)
    ENDIF
  ENDIF
  IF (!notdone)
    IF (@catt)    
      z = d^(1/@power)
    ENDIF
    IF (!landedright)
      z = (501,10)
    ENDIF
  ENDIF
bailout:
  notdone
default:
  title = "Complex Henon Julia"
  maxiter = 1000
  method = onepass
  heading
    caption = "Julia seeds"
  endheading
  param seed1
    caption = "Julia Seed 1"
    default = (1.4,0)
  endparam
  param seed2
    caption = "Julia Seed 2"
    default = (0.3,0)
  endparam
  heading
    caption = "Image settings"
  endheading
  param what
    caption = "Draw what"
    enum = "Outside/Inside" "Periodic" "Aperiodic"
    default = 0
  endparam
  param preit
    caption = "Preiterations"
    default = 100
    min = 0
    hint = "Number of iterations before looking for periodicity. \
            As with maxiter, larger numbers give slower, \
            more accurate calculations."
    visible = (@what != 0)
  endparam
  param per
    caption = "Period"
    default = 1
    min = 1
    visible = (@what == 1)
  endparam
  param lrgr
    caption = "and higher periods"
    default = true
    visible = (@what == 1)
  endparam
  param maxp
    caption = "Max Period"
    default = 100
    min = 0
    hint = "If this is not zero, orbits that don't come close \
            to themselves after this many iterations after \
            the preiterations are treated as aperiodic."
    visible = (@what != 0)
  endparam
  param catt
    caption = "Pass distance to attractor"
    default = false
    hint = "If set, the real component of z will be left \
            indicating the distance from the initial point \
            to the attractor in 4-space. Use with real(z) \
            coloring. The attractor is assumed to be all \
            the points hit after the preiterations."
    visible = (@what != 0)
  endparam
  param power
    caption = "Power"
    default = 10.0
    hint = "Affects distance to attractor coloring. Higher \
            brings the attractor points in focus more."
    visible = ((@what != 0) && @catt)
  endparam
  param bailout
    caption = "Bailout"
    hint = "Escape radius (inverse is used for finite attractors)"
    default = 100000.0
    min = 0
  endparam
  heading
    caption = "Image slice"
  endheading
  param offset
    caption = "Slice offset"
    hint = "Shifts the slice in directions orthogonal to it."
    default = (0.01,0)
  endparam
  param t
    caption = "Slice angle 1"
    hint = "0-1 is a quarter rotation."
    default = 0.0
  endparam
  param u
    caption = "Slice angle 2"
    hint = "0-1 is a quarter rotation."
    default = 1.0
  endparam
  param uu
    caption = "Rotation"
    hint = "0-1 is a quarter rotation."
    default = 1.0
  endparam
switch:
  type = "pgd_ComplexHenonMandelbrot"
  bailout = bailout  
}

pgd_ComplexHenonMandelbrot {
; The Henon mapping with both variables complex. Parameters
; allow varying the slice of the 4d parameter space shown,
; and whether escaping/nonescaping points or points of
; certain periodicities are shown. Closest approach of
; attractor to point of origin coloring may be interesting,
; as might Lyapunov coloring or the built-in alternate
; closest approach coloring. (For Lyapunov coloring, check
; "Lyapunov" and use PGD Lyapunov coloring.
;
; The default starting point for iteration gives fairly
; good results, but cannot show the whole Mandelbrot object.
; To see it "all" requires checking out numerous starting
; points and slices.
;
; Unfortunately, the Switch function will not work properly
; unless the slice angles are all zero. (It will work properly
; even if the offset is non-zero however.)
;
init:
  complex w
  float v = real(#pixel)
  float vv = imag(#pixel)
  float pp = real(@offset)
  float q = imag(@offset)
  float r = @t*#pi/2
  float s = @u*#pi/2
  float ss = @uu*#pi/2
  complex s1 = v*cos(r) + (0,1)*vv*cos(s) + pp*sin(r) + (0,1)*q*sin(s)
  complex s2 = pp*cos(r) + (0,1)*q*cos(s) - v*sin(r) - (0,1)*vv*sin(s)
  s2 = s2*exp((0,1)*ss)
  IF (abs(imag(s1)) < 1/@bailout)
    s1 = real(s1)
  ENDIF
  IF (abs(imag(s2)) < 1/@bailout)
    s2 = real(s2)
  ENDIF
  int i = 0
  bool notdone = true
  bool landedright = false
  complex t1 = 0
  complex t2 = 0
  z = @start1
  complex w = @start2
  complex iz = z
  complex iw = w
  float d = 0
  complex dxx = 1
  complex dxy = 0
  complex dyx = 0
  complex dyy = 1
loop:
  IF (@lya)
    complex ddd = -2*s1*z
    complex nxx = dxx*ddd + dyx
    complex nxy = dxy*ddd + dyy
    dyx = dxx*s2
    dyy = dxy*s2
    dxx = nxx
    dxy = nxy
  ENDIF  
  t = s2*z
  z = 1 - s1*sqr(z) + w
  w = t
  IF ((|z| + |w|) > @bailout)
    notdone = false
    landedright = (@what == 0)
  ELSEIF (@what != 0)
    IF (i == @preit)
      t1 = z
      t2 = w
      IF (@catt)
        d = |z - iz| + |w - iw|
      ENDIF
    ELSEIF (i > @preit)
      IF (@catt)
        float dd = |z - iz| + |w - iw|
        IF (dd < d)
          d = dd
        ENDIF
      ENDIF
      IF ((|z - t1| + |w - t2|) < (1/@bailout))
        notdone = false
        $define DEBUG
        print ("Here!")
        landedright = (@what == 1)
        int p = i - @preit
        IF (p < @per)
          landedright = false
        ELSEIF ((@maxp > 0) && (p > @maxp))
          landedright = (@what == 2)
        ELSEIF ((!@lrgr) && (p > @per))
          landedright = false
        ENDIF
      ENDIF
    ENDIF
    i = i + 1
    IF (i == (#maxit - 1))
      notdone = false
      landedright = (@what == 2)
    ENDIF
  ELSE
    i = i + 1
  ENDIF
  IF (!notdone)
    IF (@catt)    
      z = d^(1/@power)
    ELSEIF (@lya)
      z = log(real(dxx*dyy - dxy*dyx))/i
    ENDIF
    IF (!landedright)
      z = (501,10)
    ENDIF
  ENDIF
bailout:
  notdone
default:
  title = "Complex Henon Mandelbrot"
  maxiter = 1000
  method = onepass
  heading
    caption = "Starting point"
  endheading
  param start1
    caption = "Starting x"
    default = (0.5,0)
  endparam
  param start2
    caption = "Starting y"
    default = (0.5,0)
  endparam
  heading
    caption = "Image settings"
  endheading
  param what
    caption = "Draw what"
    enum = "Outside/Inside" "Periodic" "Aperiodic"
    default = 0
  endparam
  param preit
    caption = "Preiterations"
    default = 100
    min = 0
    hint = "Number of iterations before looking for periodicity. \
            As with maxiter, larger numbers give slower, \
            more accurate calculations."
    visible = (@what != 0)
  endparam
  param per
    caption = "Period"
    default = 1
    min = 1
    visible = (@what == 1)
  endparam
  param lrgr
    caption = "and higher periods"
    default = true
    visible = (@what == 1)
  endparam
  param maxp
    caption = "Max Period"
    default = 100
    min = 0
    hint = "If this is not zero, orbits that don't come close \
            to themselves after this many iterations after \
            the preiterations are treated as aperiodic."
    visible = (@what != 0)
  endparam
  param lya
    caption = "Lyapunov"
    hint = "Lyapunov information will be passed in z if this \
            is set. Use with PGD Lyapunov coloring."
    default = false
    visible = (!@catt)
  endparam
  param catt
    caption = "Pass distance to attractor"
    default = false
    hint = "If set, the real component of z will be left \
            indicating the distance from the initial point \
            to the attractor in 4-space. Use with real(z) \
            coloring. The attractor is assumed to be all \
            the points hit after the preiterations."
    visible = ((@what != 0) && !@lya)
  endparam
  param power
    caption = "Power"
    default = 10.0
    hint = "Affects distance to attractor coloring. Higher \
            brings the attractor points in focus more."
    visible = ((@what != 0) && @catt && !@lya)
  endparam
  param bailout
    caption = "Bailout"
    hint = "Escape radius (inverse is used for finite attractors)"
    default = 100000.0
    min = 0
  endparam
  heading
    caption = "Image slice"
  endheading
  param offset
    caption = "Slice offset"
    hint = "Shifts the slice in directions orthogonal to it."
    default = (0.01,0)
  endparam
  param t
    caption = "Slice angle 1"
    hint = "0-1 is a quarter rotation."
    default = 0.0
  endparam
  param u
    caption = "Slice angle 2"
    hint = "0-1 is a quarter rotation."
    default = 1.0
  endparam
  param uu
    caption = "Rotation"
    hint = "0-1 is a quarter rotation."
    default = 1.0
  endparam
switch:
  type = "pgd_ComplexHenonJulia"
  seed1 = #pixel
  seed2 = offset
  bailout = bailout
}

pgd_Freaky01Md {
; Mandelbrot from a two-parameter family of rational functions.
init:
     complex d = pixel
     complex a = @a
     complex b = (1 - 3*a - d)/2
     complex c = (1 + a - d)/2
     z = -1
loop:
     z = a*sqr(z) + b/z + c*z + d
bailout:
        |z| < @bail
default:
        title = "Freaky 01 Mandelbrot (d-plane)"
        maxiter = 1000
        param a
              caption = "Parameter a"
              default = (1.0,0.0)
        endparam
        param bail
              caption = "Bailout"
              default = 1000.0
        endparam
}

pgd_Freaky02Md {
; Mandelbrot from a two-parameter family of rational functions.
global:
     complex r = @r
     complex r2 = sqr(r)
init:
     complex a = pixel
     z = r
     IF (@what == 1)
          z = -z
     ENDIF
loop:
     z = 1 + a*(r2 - z)*(1/z - 1)
bailout:
        |z| < @bail
default:
        title = "Freaky 02 Mandelbrot (a-plane)"
        maxiter = 1000
        param r
              caption = "Parameter r"
              default = (1.4,0.0)
        endparam
        param bail
              caption = "Bailout"
              default = 1000.0
        endparam
        param what
              caption = "What critical point"
              enum = "r" "-r"
              default = 0
        endparam
}

pgd_Freaky03Md {
; Mandelbrot from a two-parameter family of rational functions.
global:
     complex r = @r
     complex dr = 2*r
     complex r2 = sqr(r)
init:
     complex s = pixel
     z = r
     IF (@what != 0)
          z = sqrt(9 - 8*s)
          IF (@what != 1)
            z = -z
          ENDIF
          z = r/4 * (3 + z)
     ENDIF
loop:
     z = ((z - dr)*sqr(z) + r2*z)/(r - s*z)
bailout:
        |z| < @bail
default:
        title = "Freaky 03 Mandelbrot (s-plane)"
        maxiter = 1000
        param r
              caption = "Parameter r"
              default = (0.4,0.0)
        endparam
        param bail
              caption = "Bailout"
              default = 1000.0
        endparam
        param what
              caption = "What critical point"
              enum = "r" "conj" "-conj"
              default = 0
        endparam
}

pgd_Freaky03MsMod {
; Mandelbrot from a two-parameter family of rational functions.
global:
     complex r = @r
     complex dr = 2*r
     complex r2 = sqr(r)
init:
     complex s = pixel
     z1 = r
     z = sqrt(9 - 8*s)
     IF (@which == 1)
        z = -z
     ENDIF
     z = r/4 * (3 + z)
     bool b
     float b2 = @bail
     IF (@what == 1)
        b2 = 1/b2
     ENDIF
loop:
     z1 = ((z1 - dr)*sqr(z1) + r2*z1)/(r - s*z)
     z = ((z - dr)*sqr(z) + r2*z)/(r - s*z)
     IF (@what == 1)
        b = |z - z1| > b2
     ELSE
        b = |z| < b2
     ENDIF
bailout:
        b
default:
        title = "Freaky 03 Mandelbrot (s-plane) 2"
        maxiter = 1000
        param r
              caption = "Parameter r"
              default = (0.4,0.0)
        endparam
        param bail
              caption = "Bailout"
              default = 1000.0
        endparam
        param which
              caption = "Which critical point"
              enum = "conj" "-conj"
              default = 0
        endparam
        param what
              caption = "What to draw as outside"
              enum = "z -> inf" "z -> f.att"
              default = 0
        endparam
}

pgd_Freaky03MrMod {
; Mandelbrot from a two-parameter family of rational functions.
init:
     complex r = pixel
     complex dr = 2*r
     complex r2 = sqr(r)
     complex s = @s
     z1 = r
     z = sqrt(9 - 8*s)
     IF (@which == 1)
        z = -z
     ENDIF
     z = r/4 * (3 + z)
     IF (@which == 2)
       z = 0
     ENDIF
     bool b
     b2 = 1/@bail
loop:
     IF(@which != 2)
               z1 = ((z1 - dr)*sqr(z1) + r2*z1)/(r - s*z1)
     ENDIF
     z = ((z - dr)*sqr(z) + r2*z)/(r - s*z)
     IF (@what == 1 && @which != 2)
        b = |z - z1| > b2 && |z| < @bail
     ELSE
        b = |z| < @bail
     ENDIF
bailout:
        b
default:
        title = "Freaky 03 Mandelbrot (r-plane) 2"
        maxiter = 1000
        param s
              caption = "Parameter s"
              default = (2.0,0.0)
        endparam
        param bail
              caption = "Bailout"
              default = 1000.0
        endparam
        param which
              caption = "Which critical point"
              enum = "conj" "-conj" "r"
              default = 0
        endparam
        param what
              caption = "What to draw as outside"
              enum = "z -> inf" "z -> f.att"
              default = 0
        endparam
}

pgd_Freaky03MsModB {
; Mandelbrot from a two-parameter family of rational functions.
global:
     complex r = @r
     complex dr = 2*r
     complex r2 = sqr(r)
init:
     complex s = pixel
     z1 = r
     z = sqrt(9 - 8*s)
     IF (@which == 1)
        z = -z
     ENDIF
     z = r * (3 + z) / (4 * s)
     bool b
     float b2 = @bail
     IF (@what == 1)
        b2 = 1/b2
     ENDIF
loop:
     z1 = ((z1 - dr)*sqr(z1) + r2*z1)/(r - s*z)
     z = ((z - dr)*sqr(z) + r2*z)/(r - s*z)
     IF (@what == 1)
        b = |z - z1| > b2
     ELSE
        b = |z| < b2
     ENDIF
bailout:
        b
default:
        title = "Freaky 03 Mandelbrot (s-plane) 2b"
        maxiter = 1000
        param r
              caption = "Parameter r"
              default = (0.4,0.0)
        endparam
        param bail
              caption = "Bailout"
              default = 1000.0
        endparam
        param which
              caption = "Which critical point"
              enum = "conj" "-conj"
              default = 0
        endparam
        param what
              caption = "What to draw as outside"
              enum = "z -> inf" "z -> f.att"
              default = 0
        endparam
}

pgd_Freaky03MrModB {
; Mandelbrot from a two-parameter family of rational functions.
init:
     complex r = pixel
     complex dr = 2*r
     complex r2 = sqr(r)
     complex s = @s
     z1 = r
     z = sqrt(9 - 8*s)
     IF (@which == 1)
        z = -z
     ENDIF
     z = r * (3 + z) / (4 * s)
     IF (@which == 2)
       z = 0
     ENDIF
     bool b
     b2 = 1/@bail
loop:
     IF(@which != 2)
               z1 = ((z1 - dr)*sqr(z1) + r2*z1)/(r - s*z1)
     ENDIF
     z = ((z - dr)*sqr(z) + r2*z)/(r - s*z)
     IF (@what == 1 && @which != 2)
        b = |z - z1| > b2 && |z| < @bail
     ELSE
        b = |z| < @bail
     ENDIF
bailout:
        b
default:
        title = "Freaky 03 Mandelbrot (r-plane) 2b"
        maxiter = 1000
        param s
              caption = "Parameter s"
              default = (2.0,0.0)
        endparam
        param bail
              caption = "Bailout"
              default = 1000.0
        endparam
        param which
              caption = "Which critical point"
              enum = "conj" "-conj" "r"
              default = 0
        endparam
        param what
              caption = "What to draw as outside"
              enum = "z -> inf" "z -> f.att"
              default = 0
        endparam
}

pgd_Freaky03JMod {
; Mandelbrot from a two-parameter family of rational functions.
init:
     complex r = @r
     complex dr = 2*r
     complex r2 = sqr(r)
     complex s = @s
     z = pixel
     z1 = r
     zr = sqrt(9 - 8*s)
     za = r/4 * (3 + zr)
     zb = r/4 * (3 - zr)
     bool b
     b2 = 1/@bail
     bool d = |z - z1| > 0.1 && |z - za| > 0.1 && |z - zb| > 0.1
     d = |real(z) - 3*r/4| > 0.1
loop:
     IF (d)
          z1 = ((z1 - dr)*sqr(z1) + r2*z1)/(r - s*z1)
          z = ((z - dr)*sqr(z) + r2*z)/(r - s*z)
          za = ((za - dr)*sqr(za) + r2*za)/(r - s*za)
          zb = ((zb - dr)*sqr(zb) + r2*zb)/(r - s*zb)
          b = |z - z1| > b2 && |z - za| > b2 && |z - zb| > b2 && |z| < @bail
     ELSE
         z = sqr(z)
     ENDIF
bailout:
        b
default:
        title = "Freaky 03 Julia 2"
        maxiter = 1000
        param r
              caption = "Parameter r"
              default = (1.5,0.1)
        endparam
        param s
              caption = "Parameter s"
              default = (2.0,0.0)
        endparam
        param bail
              caption = "Bailout"
              default = 1000.0
        endparam
}