comment {

This file contains the standard fractal types shipped with
Gnofract4D. They come from a variety of sources, which I've noted in
the comments for each fractal.

Each fractal includes its own Julia "twin" implicitly, so those aren't
listed as separate types.

If you modify these formulas, import of .fct files generated by previous
Gnofract4D versions may not work, so I suggest copying the formulas to
another file and renaming them before changing them.
 
In some cases, these have been adapted from Fractint formulas. 
These have had -G4 appended to avoid any confusion.


}

Mandelbrot {
; The classic Mandelbrot set
init:
	z = #zwpixel
loop:
	z = z * z + #pixel
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
}

HyperMandel {
; Hypercomplex Mandelbrot set
init:
	hyper h = @h0
	hyper c = (real(#pixel),imag(#pixel),real(#zwpixel), imag(#zwpixel)) 
loop:
	h = @hfunc(h) + c
	z = hyper_ri(h) ; for coloring algorithms to use
bailout:
	|h| < @bailout
default:
float param bailout
	default = 64.0
endparam
hyper func hfunc
	default = sqr
endfunc
hyper param h0
	default = (0.0,0.0,0.0,0.0)
endparam
}

HyperJulia {
; Hypercomplex Julia set
init:
	hyper h = (real(#pixel),imag(#pixel),real(#zwpixel), imag(#zwpixel)) 
loop:
	h = @hfunc(h) + @c
	z = hyper_ri(h) ; for coloring algorithms to use
bailout:
	|h| < @bailout
default:
float param bailout
	default = 64.0
endparam
hyper func hfunc
	default = sqr
endfunc
hyper param c
	default = (0.16,0.1,0.1,0.1)
endparam
}

Mandelbar {
; I first came across this at http://mathworld.wolfram.com/MandelbarSet.html
init:
	z = #zwpixel
loop:
	z = conj(z)^@a + #pixel
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
complex param a
	title = "Power"
	default = (2.0,0.0)
endparam
zcenter=1.0e-10
}

Quadratic {
; Included for backwards compatibility with earlier versions. There are more
; parameters here than are strictly warranted, since any combination of A,B, 
; and C is actually equivalent to moving the initial point around.

init:
	z = #zwpixel
loop:
	z = (@a * z + @b) * z + @c * #pixel
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
complex param a
	default = (1.0, 0.0)
endparam
complex param b
	default = (1.0, 0.0)
endparam
complex param c
	default = (1.0, 0.0)
endparam
}

Cubic Mandelbrot {
; z <- z^3 + c
; The cubic set actually has two critical values, but this formula just uses  
; zero - to be fixed later.
init:
	z = #zwpixel
	; nothing to do here
loop:
	z = z * z * (z - 3.0 * @a) + #pixel
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
complex param a
	default = (0.0,0.0)
endparam
}

ManZPower {
; The Mandelbrot set raised to the N'th power
init:
	z = #zwpixel
loop:
	z = z^@a + #pixel
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
complex param a
	title = "Power"
	default = (4.0,0.0)
endparam
float func bailfunc
	default = cmag
endfunc
zcenter=1.0e-10
}

Barnsley Type 1 {
; From Michael Barnsley's book Fractals Everywhere, via Fractint
init:
	z = #zwpixel
loop:
	float x_cy = real(z) * imag(#pixel)
	float x_cx = real(z) * real(#pixel)
	float y_cy = imag(z) * imag(#pixel)
	float y_cx = imag(z) * real(#pixel)

	if(real(z) >= 0)
		z = (x_cx - real(#pixel) - y_cy, y_cx - imag(#pixel) + x_cy)
	else 
		z = (x_cx + real(#pixel) - y_cy, y_cx + imag(#pixel) + x_cy)
	endif
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
}

Barnsley Type 2 {
; From Michael Barnsley's book Fractals Everywhere, via Fractint
init:
	z = #zwpixel
loop:
	float x_cy = real(z) * imag(#pixel)
	float x_cx = real(z) * real(#pixel)
	float y_cy = imag(z) * imag(#pixel)
	float y_cx = imag(z) * real(#pixel)

	if(real(z) * imag(#pixel) + imag(z) * real(#pixel) >= 0)
		z = (x_cx - real(#pixel) - y_cy, y_cx - imag(#pixel) + x_cy)
	else 
		z = (x_cx + real(#pixel) - y_cy, y_cx + imag(#pixel) + x_cy)
	endif
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
}

Barnsley Type 3 {
; From Michael Barnsley's book Fractals Everywhere, via Fractint
init:
	z = #zwpixel
loop:
	float x2 = real(z) * real(z)
	float y2 = imag(z) * imag(z)
	float xy = real(z) * imag(z)

	if(real(z) > 0)
		z = (x2 - y2 - 1.0, xy * 2.0)
	else
		z = (x2 - y2 - 1.0 + real(#pixel) * real(z), \
		     xy * 2.0 + imag(#pixel) * real(z))
	endif
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
}
 
XBuffalo {
init:
	z = #zwpixel
loop:
	z = (abs(real(z)), imag(z))
	z = (z - 1.0) * z + #pixel
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
magnitude=6.0
}     

YBuffalo {
init:
	z = #zwpixel
loop:
	z = (real(z), abs(imag(z)))
	z = (z - 1.0) * z + #pixel
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
magnitude=6.0
}     

Buffalo {
; From the web page http://www.theory.org/fracdyn/ 

; The Burning Ship is essentially a Mandelbrot variant where the real
; and imaginary parts of the current point are set to their absolute
; values after each iteration, ie z <- (|x| + i |y|)^2 + c. The
; Buffalo fractal uses the same method with the function z <- z^2 - z
; + c, making it equivalent to the Quadratic type with the "absolute
; value" modification.  

init:
	z = #zwpixel
loop:
	z = abs(z)
	z = (z - 1.0) * z + #pixel
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
magnitude=6.0
}     
		
Burning Ship {
; From the web page http://www.theory.org/fracdyn/ 

; The Burning Ship is essentially a Mandelbrot variant where the real
; and imaginary parts of the current point are set to their absolute
; values after each iteration, ie z <- (|x| + i |y|)^2 + c. The
; Buffalo fractal uses the same method with the function z <- z^2 - z
; + c, making it equivalent to the Quadratic type with the "absolute
; value" modification.  


init:
	z = #zwpixel
loop:
	z = abs(z)
	z = z*z + #pixel
bailout: 
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
xycenter = (-0.5,-0.5)
}

Cubic Burning Ship {
init:
	z = #zwpixel
loop:
	z = (abs(real(z)),abs(imag(z)))
	z = z*z*z + #pixel
bailout: 
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
}

Lambda {

; The Lambda function calculates lambda * z * ( 1 - z). The complex
; parameter lambda is set by the z and w parameters, so if lambda is
; zero, all you'll see is a blank screen.

init:
	z = #zwpixel
loop:
	t = z * (1.0 - z)
	z = t * #pixel
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
xcenter=1.0
zcenter=0.5
magnitude=8.0
}

Magnet {

; Magnet and Magnet 2 are from Fractint, but images generated by
; Gnofract 4D look a bit different because I don't look for a finite
; attractor.

init:
	z = #zwpixel
loop:
	z = (z * z + #pixel - 1.0)/(2.0 * z + #pixel - 2.0)
	z = z *z
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
xcenter=2.0
magnitude=8.0
}

Magnet 2 {

; Magnet and Magnet 2 are from Fractint, but images generated by
; Gnofract 4D look a bit different because I don't look for a finite
; attractor.

init:
	z = #zwpixel
loop:
	cm1 = #pixel - 1.0, cm2 = #pixel - 2.0
	z = (z * z * z + 3.0 * cm1 * z + cm1 * cm2)/ \
	     (3.0 * z * z + 3.0 * cm2 * cm2 + cm1 * cm2 + 1.0)
	z = z*z

bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
xcenter=2.0
magnitude=3.0
}

Newton {
; Newton is the Newton-Raphson method applied to z^a = root
init:
	z = #zwpixel
	nm1 = @a - 1.0
loop:
	last = z
	z = z - (z ^ @a - @root)/ (@a * z ^ nm1)
bailout:
	|z - last| > #tolerance
default:
xzangle=1.5707963267948966
ywangle=1.5707963267948966
xcenter=1.0
param a
	default = (3.0, 0.0)
endparam
param root
	default = (1.0, 0.0)
endparam

}


NewtBasin {
; Newton-Raphson method which determines which 'basin of attraction'
; each point ends up in.
init:
	z = #zwpixel
	nm1 = @power - 1.0
loop:
	last = z
	z = z - (z ^ @power - 1.0)/ (@power * z ^ nm1)
bailout:
	|z - last| > #tolerance
final:
	int i = 0
	int whichroot = 0
	float dist = 1.0e20
	while i < @power
		root = (cos(i*#pi * 2.0 / @power), sin(i*#pi * 2.0 / @power))
		float thisdist = |root - z|
		if(thisdist < dist)
			whichroot = i
			dist = thisdist
		endif
		i = i + 1
	endwhile
	#fate = whichroot

default:
xzangle=1.5707963267948966
ywangle=1.5707963267948966
xcenter=1.0
int param power
	default = 3
endparam
}

;NewtBasin {
;; Newton is the Newton-Raphson method applied to z^a = -1
;init:
;	z = #zwpixel
;	ca = (@a,0)
;	nm1 = @a - 1.0
;loop:
;	last = z
;	z = z - (z ^ ca)/ (@a * z ^ nm1)
;bailout:
;	|z - last| > #tolerance
;final:
;;	int i = 0
;;	int whichroot = 0
;;	float dist = 1.0e20
;;	while i < @a
;;		root = (cos(i*#pi * 2.0 / @a), sin(i*#pi * 2.0 / @a))
;;		float thisdist = |root - z|
;;		if(thisdist < dist)
;;			whichroot = i
;;		endif
;;	endwhile
;;	if whichroot > 0
;;		#fate = whichroot + 1
;;	endif		
;default:
;xzangle=1.5707963267948966
;ywangle=1.5707963267948966
;xcenter=1.0
;int param a
;	default = 3
;endparam
;}

Nova {
; Nova is Paul Derbyshire's Nova fractal.
init:
	z = #zwpixel
loop:
	last = z
	z = z - (@a * z * z * z - @b)/(@c * z * z) + #pixel
bailout:
	|z - last| > @epsilon
default:
param a
	default = (1.0, 0.0)
endparam
param b
	default = (1.0, 0.0)
endparam
param c
	default = (3.0, 0.0)
endparam
float param epsilon
	default = 0.01
endparam
zcenter=1.0
magnitude=3.0
}


Tetrate {
; Tetrate computes z <- z^c
init:
	z = #zwpixel
loop:
	z = #pixel^z
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default = cmag
endfunc
}


T02-01-G4 {; Modified for Gf4d by EY
	; V.1.1 - earlier versions may be discarded
        ; Copyright (c)1998,1999 Morgan L. Owens
        ; Chebyshev Types:
        ; Inspired by Clifford A. Pickover:
        ; Dynamic (Euler method)
        ;
        ; T(n+1) = 2xT(n)-T(n-1)
        ; T(0)  = 1
        ; T(1)  = x
        ;
        ; = 2zT01-T00
  t=#zwpixel, z=pixel:
  x=real(z), y=imag(z)
  Tx=(x+x)*x-1
  Ty=(y+y)*y-1
  x=x-t*Ty, y=y+t*Tx
  z=x+flip(y)
  |z| < @bailout
default:
float param bailout
	default = 4.0
endparam
zwcenter = (0.3,0.2)
}

T03-01-G4 {; based on T03-01 in CHBY1.FRM by Morgan L. Owens
        ; Modified for Gf4D by EY
        ; = 2zT02-T01
  t=#zwpixel, z=pixel:
  float x=real(z), float y=imag(z)
  float Tx=x*(4*x*x-3)
  float Ty=y*(4*y*y-3)
  cx=x-t*Ty, cy=y+t*Tx
  z=cx+flip(cy)
  |z| < @bailout
default:
float param bailout
	default = 4.0
endparam
zwcenter = (0.02,0.01)
}

Cubic Connectedness Locus {
; iterates both critical points of the cubic set,
; bails out only if both escape
; NB: this formula is incorrect, but is kept for backwards compatibility

init:
 k = #zwpixel
 z1 = k
 z2 = -k

loop:
 z1 = z1*z1*z1 - 3*k + #pixel
 z2 = z2*z2*z2 - 3*k + #pixel
 z = z1 + z2 ; so coloring algorithms have something to work ond
 
bailout:
 |z1| < @bailout || |z2| < @bailout

default:
 float param bailout
 	default = 4.0
 endparam
}

Cubic Connectedness Locus 2 {
; iterates both critical points of the cubic set,
; bails out if either escape
; NB: this formula is incorrect, but is kept for backwards compatibility
init:
k = #zwpixel
z1 = k
z2 = -k
loop:
z1 = z1*z1*z1 - 3*k + #pixel
z2 = z2*z2*z2 - 3*k + #pixel
z = z1 + z2 ; so coloring algorithms have something to work on

bailout:
|z1| < @bailout && |z2| < @bailout

final:
if |z1| < @bailout
   #fate = 1
elseif |z2| < @bailout
   #fate = 2
endif

default:
float param bailout
	default = 4.0
endparam
}

Cubic Connectedness Locus 3 {
; iterates both critical points of the cubic set,
; bails out only if both escape

init:
k = #zwpixel
z1 = k
z2 = -k
k2 = 3.0 * k * k
int z1_iter = -1
int z2_iter = -1

loop:
z1 = z1*z1*z1 - k2 * z1 + #pixel
z2 = z2*z2*z2 - k2 * z2 + #pixel
z = z1 + z2 ; so coloring algorithms have something to work ond

if |z1| > @bailout && z1_iter == -1
   z1_iter = #numiter
elseif |z2| > @bailout && z2_iter == -1
   z2_iter = #numiter
endif

bailout:
|z1| < @bailout || |z2| < @bailout

final:
if |z1| < @bailout
   if |z2| >= @bailout
      #fate = 1
      #inside = false
      #numiter = z2_iter
   endif
elseif |z2| < @bailout
   #fate = 2
   #inside = false
   #numiter = z1_iter
endif

default:
float param bailout
	default = 4.0
endparam
}

FnParts {
; A generalization of Burning Ship - apply separate functions to the 
; X and Y parts of Z, then another to Z itself
init:
	z = #zwpixel
loop:
	z = @fnComplex(@fnReal(real(z)), @fnImag(imag(z))) + #pixel
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4.0
endparam
float func bailfunc
	default =cmag
endfunc
float func fnReal
	argtype = float
	default = abs
endfunc
float func fnImag
	argtype = float
	default = abs
endfunc
func fnComplex
	default = sqr
endfunc
}

Chebyshev {
; a generalization of a number of Fractint Chebyshev functions
init:
  float s=sqrt(2)
  float a=0, float b=0, float c=0
  if @functype == "04-01"
    a = 8*s+7
  elseif @functype == "05-01"
    a = 8*s+7
  elseif @functype == "06-01"
    a = 33*s+12
  elseif @functype == "07-01"
    a=87*s+18, b=8*s-57
  elseif @functype == "08-01"
    a=185*s+25, b=41*s-141
  elseif @functype == "09-01"
    a=345*s+33, b=136*s-285
  elseif @functype == "10-01"
    a=588*s+42, b=321*s-510, c=41*s-1830
  endif
  float Tx, float Ty, float xx, float yy
  complex t=#zwpixel
  z=#pixel
loop:
  float x=real(z), float y=imag(z)
if @functype == "02-01"
  Tx=s*x*x-1
  Ty=s*y*y-1
elseif @functype == "03-01"
  Tx=x*(s*(x*x-2)-1)
  Ty=y*(s*(y*y-2)-1)
elseif @functype == "04-01"
  xx=x*x, yy=y*y
  Tx=xx*(s*(xx-5)-1)+3
  Ty=yy*(s*(yy-5)-1)+3
elseif @functype == "05-01"
  xx=x*x, yy=y*y
  Tx=x*(xx*(s*(xx-9)-1)+a)
  Ty=y*(yy*(s*(yy-9)-1)+a)
elseif @functype == "06-01"
  xx=x*x, yy=y*y
  Tx=xx*(xx*(s*(xx-14)-1)+a)-15
  Ty=yy*(yy*(s*(yy-14)-1)+a)-15
elseif @functype == "07-01"
  xx=x*x, yy=y*y
  Tx=x*(xx*(xx*(s*(xx-20)-1)+a)+b)
  Ty=y*(yy*(yy*(s*(yy-20)-1)+a)+b)
elseif @functype == "08-01"
  xx=x*x, yy=y*y
  Tx=xx*(xx*(xx*(s*(xx-27)-1)+a)+b)+105
  Ty=yy*(yy*(yy*(s*(yy-27)-1)+a)+b)+105
elseif @functype == "09-01"
  xx=x*x, yy=y*y
  Tx=x*(xx*(xx*(xx*(s*(xx-35)-1)+a)+b)+561)
  Ty=y*(yy*(yy*(yy*(s*(yy-35)-1)+a)+b)+561)
elseif @functype == "10-01"
  xx=x*x, yy=y*y
  Tx=xx*(xx*(xx*(xx*(s*(xx-44)-1)+a)+b)+c)-945
  Ty=yy*(yy*(yy*(yy*(s*(yy-44)-1)+a)+b)+c)-945
endif
  z= x-t*Ty + flip(y+t*Tx)
bailout:
	@bailfunc(z) < @bailout
default:
float param bailout
	default = 4000.0
endparam
float func bailfunc
	default =cmag
endfunc
int param functype
	enum = "02-01" "03-01" "04-01" "05-01" "06-01" "07-01" "08-01" "09-01" "10-01"
endparam
zcenter=1.0
wcenter=0.1
magnitude=8.0
maxiter=256
}