!READ_ME_FIRST{; press F2 to read this!
; This file: colour-c.frm
; Cleaned up version of colour.frm, available from:
; http://members.xoom.com/lukeplant/
;
; This formula file is (c) Luke Plant 1998, lukeplant@aol.com
; To understand and use the formulae:
; please refer to the actual formula file itself.  Some comments
; are written in the formulae themselves, but most are written BETWEEN
; the formulae to explain the maths behind them and how to use them.  I
; have not written any foolproofing into the formulae - it is very easy
; to end up with a blank screen! A few notes on how to avoid this -  
;  1) ALMOST ALL THESE FORMULAE WORK WITH OUTSIDE=REAL ONLY
;  2) For the type M-Cir-xxx and M-Ell-xxx, at least one/two radii are 
;  required. If these are zero, the division by zero error will give
;  you a blank screen.
;
;
; Here is a brief key to the formula names
;
; Type 
; ~~~~
; M    - mandelbrot
; J    - julia
;
; Trap
; ~~~~
; Cir  - indicates circular trap
; Ell  - ellipse
; Prb  - parabola
; Hyp  - hyperbola
; Ast  - astroid
; MEll - moving ellipse. the position/size of the trap is
;        designed to vary
;
; Numbers after the type indicate improved/different possibilities
; in the formulae, such as rotated ellipses.
;
; Fn$  - not strictly a trap, but uses some function of distance
;        from point, giving ripples or even very smooth colouring
;        This type is usually inifinite - all points are affected 
;	   to some degree.
;
; options:
; ~~~~~~~~
; A    - uses (A)ll values of z that meet criterion
; e    - uses 'distance' from Edge - the best kind usually
; m    - subtracts normal colouring (no iteration bands)
;
; NB/ other options are available in the full formula file
; The ones shown above generally give the best results - only these have
; been kept for this version
;
; Have a look at this file in NotePad  to read the other
; comments and  finally - I hope you enjoy them
;
}


 


MCol3(XAXIS) {
; try 4,1,tan
c=pixel, done=0, s=0:

z=done
z=z*z+c
s=s+fn1(z)
done=z
z=s * imag(p1)
|done| <= real(p1)
}

MCol3b(XAXIS) { ; re(p2)=bail,
; try 1,2,4,2,recip
c=pixel, done=0, s=0:

z=done
z=z*z+c
s=s+fn1(z)
done=z
z=s * p1
|done| <= real(p2)
}

MCol3-shift(XAXIS) { ; re(p2)=bailout
; try 4,1,tan
c=pixel, w=0, s=0:
w=w*w+c
s=s+fn1(w)
z=s*p1+(2048,2048)
|w| <= real(p2)
}


MCol3b-shift(XAXIS) { ; re(p2)=bail
; try 1,2,4,2,recip
c=pixel, done=0, s=0:

z=done
z=z*z+c
s=s+fn1(z)
done=z
z=s * p1 +imag(p2)*(256,256)
|done| <= real(p2)
}


MCol4(XAXIS) {
; try 16,1,tanh,flip
c=pixel, w=0, s=0:

w=w*w+c
s=fn2(s+fn1(w))
z=s * imag(p1)
|w|<= real(p1)
}

MCol4-m {
; try 16,1,tanh,flip
c=pixel, z=0, s=0, fin=0, count=0:

z=z*z+c
count=count+1
s=fn2(s+fn1(z))
If (|z| > real(p1))
	s=s*imag(p1)
	z=s-255*floor((s-1)/real(255))
	z=z-count-flip(count)
	fin=1
End If

fin==0
}


MCol6(XAXIS) { 
;try 2, 0, 8, 0.1, asinh
c=pixel, temp=0, q=0:

z=temp*temp+c
q=q+p1*fn1(temp/(z+imag(p2)))
temp=z
z=q+(256,256)

|temp| <= real(p2)
}

MCol6-m(XAXIS) { 
;try 2, 0, 8, 0.1, asinh

c=pixel, temp=0, q=0, count=0, fin=0:

z=temp*temp+c
count=count+1
q=q+p1*fn1(temp/(z+imag(p2)))
temp=z

If (|z| > real(p2))
	z=q-255*floor((q-1)/real(255))
	z=z-count-flip(count)
	fin=1
End If

fin==0
}


JCol6(ORIGIN) { ; P1 is normal Julia parameter
;try -0.1, 0.75, 2, 0, 8, 0.1, asinh

temp=pixel, c=p1, q=0:
z=temp*temp+c
q=q+p2*fn1(temp/(z+imag(p3)))
temp=z
z=q

|temp| <= real(p3)
}


JCol6-m(ORIGIN) { ; P1 is normal Julia parameter
;try -0.1, 0.75, 2, 0, 8, 0.1, asinh 
temp=z=pixel, c=p1, q=0, fin=0, count=0:

z=z*z+c
count=count+1
q=q+p2*fn1(temp/(z+imag(p3)))
temp=z
If (|z| > real(p3))
	z=q-255*floor((q-1)/real(255))
	z=z-count-flip(count)
	fin=1
End If

fin == 0
}



MCol6b(XAXIS) { 
c=pixel, temp=0, q=0:

z=temp*temp+c
q=q+p1*fn1(temp/z)
temp=z
z=q+(256,256)

|temp| <= real(p2)
}

JCol6b(ORIGIN) { 
temp=pixel, c=p1, q=0:

z=temp*temp+c
q=q+p2*fn1(temp/z)
temp=z
z=q

|temp| <= real(p3)
}



:-------------------------------------------------------
; trap 5- circle - (all z(n): adds distance from edge of circle)

M-Cir-A-e(XAXIS) { ; trap- p1 = centre of circle
; rP2 = sqr(radius), rP3 = bailout
; iP3 = colour factor

c=pixel, z=0, q=0, fin=0:
z=z*z+c
test=(|z-p1| - real(p2))
If (test<0)
	q=q+ test/real(p2)*imag(p3)
END IF
IF (|z| > real(p3))
	z=q-255*floor((q-1)/real(255))
	fin=-1
END IF

fin == 0
}

M-Cir-A-e-m(XAXIS) { ; trap- p1 = centre of circle
; rP2 = sqr(radius), rP3 = bailout
; iP3 = colour factor

c=pixel, z=0, q=0, count=0, fin=0:
z=z*z+c
count=count+1
test=(|z-p1| - real(p2))
If (test<0)
	q=q+ test/real(p2)*imag(p3)
END IF
IF (|z| > real(p3))
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF

fin == 0
}


;-------------------------------------------------------
; trap 10- rectangular hyperbola, (all z(n): 'distance' from edge)

M-Hyp-A-e(XAXIS) { ; trap- p1 = shift, rP2=distance
		; iP2=palette shift, rP3 = bailout
c=pixel, z=0, q=0, fin=0:
z=z*z+c

test=(|(real(z)-real(p1))*(imag(z)-imag(p1))| - real(p2))
IF (test<0) 
	q=q+test
END IF

IF (|z| > real(p3))
	q=q*imag(p3)
	z=q-255*floor((q-1)/real(255))
	fin=-1
END IF

fin == 0
}


M-Hyp-A-e-m(XAXIS) { ; hyperbola:
; P1 = center, rP2='distance'
; rP3 = bailout, iP3= colour factor
; try-   -0.75,0.5,1,0,16,50

c=pixel, z=0, q=0, count=0, fin=0:
z=z*z+c
count=count+1

test=(|(real(z)-real(p1))*(imag(z)-imag(p1))| - real(p2))
IF (test<0) 
	q=q+test
END IF

IF (|z| > real(p3))
	q=q*imag(p3)
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF

fin == 0
}

;--------------------------------------------------------
; trap - parabola:  parabola is defined using the focus
; and the point 'opposite' the focus on the directrix


M-Prb2-A-e-m(XAXIS) { ; parabola:
; p1 = focus, P2 = point 'opposite' focus on directrix
; rP3 = bailout, iP3 = colour factor
; try-   1,1,1.05,1.05,16,20

c=pixel, z=0, q=0, count=0, fin=0
m=(real(p2)-real(p1))/(imag(p1)-imag(p2))
k=-m*real(p2)+imag(p2):
z=z*z+c
count=count+1
N=(m*(imag(z)-k) + real(z))/(1+sqr(m))
N=N+flip(m*N+k)
test1 = |z-p1|-|z-N|
IF (test1<0)
	q=q +  test1
END IF

IF (|z| > real(p3))
	q=q*imag(p3)
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF

fin == 0
}

;--- Julia versions have focus as origin always - not enough parameters

J-Prb2-A-e-m(ORIGIN) { ; parabola
; P1 = julia parameters i.e. = c
; focus=0 , P2 = point 'opposite' focus on directrix
; rP3 = bailout, iP3 = colour factor

c=P1, z=pixel, q=0, fin=0, count=0
m=real(p2)/(-imag(p2))
k=-m*real(p2)+imag(p2):
z=z*z+c
count=count+1
N=(m*(imag(z)-k) + real(z))/(1+sqr(m))
N=N+flip(m*N+k)
test1 = |z-p1|-|z-N|
IF (test1<0)
	q=q + test1
END IF
IF (|z| > real(p3))
	q=q*imag(p3)
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF
fin ==0

}



;-------------------------------------------------------
; ellipse (all z(n): adds distance from edge of ellipse)


M-Ell-A-e-m(XAXIS) { ; ellipse
; p1 = centre of ellipse
; rP2 = sqr(x radius), iP2= sqr(y radius)
; rP3 = bailout, iP3 = colour factor

c=pixel, z=0, q=0, count=0, fin=0:
z=z*z+c
count=count+1

test=sqr(real(z-p1))/real(p2)+sqr(imag(z-p1))/imag(p2)
IF (test<1)
	q=q+1-test
END IF

IF (|z| > real(p3))
	q=q*imag(p3)
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF

fin == 0
}


;-----------------------------------------------
;ellipse centered on (rP2,ip2) and can be tilted
; bailout set to 16

M-Ell3-A-e-m(XAXIS) { ; ellipse:
; p1 = radii of ellipse, 
; P2=centre of ellipse
; rP3 = tilt, iP3 = colour factor
; bailout=16

c=pixel, z=0, q=0, count=0, fin=0
sint=sin(real(p3)), cost=cos(real(p3)):
z=z*z+c
count=count+1
ex=cost*(real(z)-real(p2))+sint*(imag(z)-imag(p2))
ey=(imag(z)-imag(p2)-ex*sint)/cost
test=sqr(ex/real(p1))+sqr(ey/imag(p1))
IF (test<1)
	q=q+(1-test)
END IF

IF (|z| > 16)
	q=q*imag(p3)
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF

fin == 0
}


;-----------------------------------------------------------------------------
;-------------------------------------------------------------------
;Implementing test1 into MSet - NB 


M-Fn2-A-m { ; pond ripple effect
; P1 is center of circles, iP3 = bailout
; try 1,1,10,10,1,16
a=real(p2)  				; colour factor
o=imag(p2)  				; oscillator factor
b=real(p3)  				; decay factor
ep=atan(-b/o) 				; see above for E
z=0, c=pixel, q=0, count=0, fin=0:

z=z*z+c
count=count+1
d=cabs(z-p1)				; distance from point
q=q+(exp(-d*b))*abs(cos(d*o+ep))
IF ( |z| > imag(p3))
	q=q*a
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF
fin==0
}

;--------------------------------------------------------
; as above but with elliptical rings

M-Fn4-A-m {; pond ripple effect, with elliptical rings
; p1 is point which is center of rings
; imag(p3) is height/width for ellipses
; try 1,1,10,10,1,2 to start with

a=real(p2)  				; colour factor
o=imag(p2)  				; oscillator factor
b=real(p3)  				; decay factor
ep=atan(-b/o) 				; see above for E
z=0, c=pixel, q=0, count=0, fin=0:

z=z*z+c
count=count+1
d=z-p1						; distance from point
d=cabs(imag(p3)*real(d)+flip(imag(d)))	; but distorted for ellipses
q=q+(exp(-d*b))*abs(cos(d*o+ep))
IF ( |z| > 16)
	q=q*a
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF
fin==0
}

;----------------------------------------------------------------------
M-Fn6-A-m {; distorted pond ripple effect - uses
; functions to give variations
; P1 is point, center of colouring, rP3 is bailout
; Fn1 is decay function for x direction
; Fn2 is decay function for y direction
; NB - large decay factor gives small width/height
; try-   2, 2, 0.1, 0.1, 20, 20, ident, flip

bx= real(p2)			; decay factor in x direction
by= imag(p2)			; decay factor in y direction
a= imag(p3) 			; colour factor
z=0, c=pixel, q=0, count=0, fin=0:
z=z*z+c
count=count+1
x=real(z-p1)
y=imag(z-p1)
q=q+(exp(-fn1(x*x)*bx)*exp(-fn2(y*y)*by))
IF (|z| > real(p3))
	q=q*a
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF
fin==0
}

;----------------------------------------------------------------------
; special case of the above - gives very smooth colouring 
; the functions have been restricted but the base can be changed

M-Fn7-A-m {;
; rP1 is decay factor, iP1 is colour factor, rP2 is bailout
; iP2 = base for exponent
; try 0.1,30,100,100  for super smooth colouring.
; increase iP1 to put bands closer together

a= imag(p1) 			; colour factor
b= real(p1)				; decay factor
bail=real(p2)			; bailout
p=imag(p2)				; base for exponent
z=0, c=pixel, q=0, count=0, fin=0:
z=z*z+c
count=count+1
q=q+(p^(-b*|z|))
IF (|z| > bail)
	q=q*a
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF
fin==0

}


;--------------------------------------------------------------------
; Back to MSet traps, not functions.  More unusual patterns wanted,
; try some 'moving target' types.  Possibilities include moving
; position/size of target, depending on z,|z|, count or other

;--------------------------------------------------------------------
; using count to uniformly increase angle of centre of circle
M-EllM1-A-e-m {; ellipse - position variable
; ellipse - rP1, iP1 = radii of ellipse
; rP2=distance of circle  from origin 
; iP2 = step for changing angle, 
; (angle determined by count and iP2)
; rP3 = bailout, iP3 = colour factor
; (NB colour factor can be negative or positive - try both)
; try-    4, 2, 2, -0.1, 16, 15 

z=0, c=pixel, count=0, fin=0, q=0, angle=0:
z=sqr(z)+c
count=count+1

angle=count*imag(p2)
x=cos(angle)				; centre of circle
x=real(p2)*(x+flip(sin(angle)))	; = x (a,b)

test=z-x
test=sqr(real(test)/real(p1))+sqr(imag(test)/imag(p1))
					; test for points in ellipse
IF (test<1)
	q=q+1-test			; colour according to distance 
END IF				; from edge


IF (|z| > real(p3))		; bailout
	q=q*imag(p3)
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF

fin==0
	
}


;------------------------------------------------------------
; Try fixed angle of centre of circle. Position of circle is fixed
; by default distance, angle plus function of |z|
;
M-EllM3-A-e-m {; moving ellipse - the distance
; of the centre of the ellipse varies according to
; fn1(|z))
;
; try .1,2,2,0,16,10, fn1=sqr
; ellipse - rP1, iP1 = radii of ellipse
; rP2=default distance of centre of ellipse from origin
; iP2=angle of centre of ellipse
; rP3 = bailout, iP3 = colour factor
; (NB colour factor can be negative or positive - try both)

z=0, c=pixel, count=0, fin=0, q=0, angle=0:
z=sqr(z)+c
count=count+1

sqmodz=fn1(|z|)
x=cos(imag(p2))				; centre of circle=x
x=real(p2)*sqmodz*(x+flip(sin(imag(p2))))

test=z-x
test=sqr(real(test)/real(p1))+sqr(imag(test)/imag(p1))
					; test for points in ellipse
IF (test<1)
	q=q+1-test			; colour according to distance 
END IF				; from edge


IF (|z| > real(p3))		; bailout
	q=q*imag(p3)
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF

fin==0
	
}


;------------------------------------------------------------------
; identical to M-EllM3b-A-e-m, but use previous value of |z|

M-EllM3b-A-e-m {; as per EllM3b, but use previous value of z 
; try .1,2,2,0,16,10, fn1=sqr
; moving ellipse - rP1, iP1 = radii of ellipse
; rP2=default distance of centre of ellipse from origin
; iP2=angle of centre of ellipse
; rP3 = bailout, iP3 = colour factor
; (NB colour factor can be negative or positive - try both)

z=0, c=pixel, count=0, fin=0, q=0, angle=0:
sqmodz=fn1(|z|)

z=sqr(z)+c
count=count+1

x=cos(imag(p2))				; centre of circle=x
x=real(p2)*sqmodz*(x+flip(sin(imag(p2))))

test=z-x
test=sqr(real(test)/real(p1))+sqr(imag(test)/imag(p1))
					; test for points in ellipse
IF (test<1)
	q=q+1-test			; colour according to distance 
END IF				; from edge


IF (|z| > real(p3))		; bailout
	q=q*imag(p3)
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF

fin==0
	
}





;----------------------------------------------------------------
;Try varying one radius only

M-EllM5-A-e-m {; variable horizontal radius of ellipse
; try .4,1,-0.5,0.5,16,10, tan , recip
; ellipse - rP1, iP1 = radii of ellipse
; rP2,iP2 = centre of ellipse
; rP3 = bailout, iP3 = colour factor
; (NB colour factor can be negative or positive - try both)

z=0, c=pixel, count=0, fin=0, q=0, angle=0:
z=sqr(z)+c
count=count+1

rf=fn2(fn1(|z|))			; radius factor

test=z-p2
test=sqr(real(test)/(real(p1)*rf))+sqr(imag(test)/(imag(p1)))
					; test for points in ellipse
IF (test<1)
	q=q+1-test			; colour according to distance 
END IF				; from edge


IF (|z| > real(p3))		; bailout
	q=q*imag(p3)
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF

fin==0
	
}

M-EllM6-A-e-m {; variable vertical radius of ellipse
; try .4,1,-0.5,0.5,16,10, tan , recip
; ellipse - rP1, iP1 = radii of ellipse
; rP2,iP2 = centre of ellipse
; rP3 = bailout, iP3 = colour factor
; (NB colour factor can be negative or positive - try both)

z=0, c=pixel, count=0, fin=0, q=0, angle=0:
z=sqr(z)+c
count=count+1

rf=fn2(fn1(|z|))			; radius factor

test=z-p2
test=sqr(real(test)/(real(p1)))+sqr(imag(test)/(imag(p1)*rf))
					; test for points in ellipse
IF (test<1)
	q=q+1-test			; colour according to distance 
END IF				; from edge


IF (|z| > real(p3))		; bailout
	q=q*imag(p3)
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF

fin==0
	
}

;-----------------------------------------------------------

; smallest |real(z)| used.  Series of formula - M-Pxx,
; ie single point used
; start from z=c rather than z=0 to eliminate low values of
; real(z) for first iteration
; imag(p2)= offset ie smallest value of real(z) relative to 

;-------------------------------------------------------------
; new shape - astroid.  formula : (y/a)^(2/3) + (x/a)^(2/3) =1
; use (y/a)^(2/3) + (x/b)^(2/3) = 1 for variable aspect ratio


M-Ast-A-e-m { ;
; P1 = center of astroid.
; real(p2) = width
; imag (p2) = height
; real(p3) = bailout
; imag(p3) = colour factor
; try -1,1,1,2,16,30
z=0, c=pixel, count=0, fin=0, q=0:

z=z*z+c
count=count+1

d= (abs(real(z-P1)/real(P2)))^(.666666666) + (abs(imag(z-P1)/imag(P2)))^(.666666666)
IF (d <= 1)
	q=q+1-d
END IF

IF (|z| > real(p3))
	q=q*imag(p3)
	z= q - 255*floor((q-1)/real(255))
	z=z-count
	fin=1
END IF

fin==0

}

;--------------------------------------------------------------------------

; new FN type - doesn't give a circular shape, but a fan shape from center

M-FN-fan-m { ; fan/propellor shape
; P1 = center of fan.
; real(p2) = decay factor
; imag (p2) = number of arms
; real(p3) = bailout
; imag(p3) = colour factor
; try 1,1,1,4,16,30

z=0, c=pixel, count=0, fin=0, q=0:

z=z*z+c
count=count+1

zm=z-P1
q=q+ exp( (-1)*|zm|*real(p2) ) * cos( imag(p2)*atan(imag(zm)/real(zm)) )

IF (|z| > real(p3))
	q=q*imag(p3)
	z= q - 255*floor((q-1)/real(255))
	z=z-count
	fin=1
END IF

fin==0

}