!READ_ME_FIRST{; press F2 to read this!
; 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 one/two radii is 
;  required. If these are zero, the division by zero error will give
;  you a blank screen.
;
; The formulas named Mcol$$ do not have mathematical explanations -
; I was guesssing to see if I could come up with anything interesting.
; It will probably be easiest to start with one of my par files for 
; these types.  
;
; 
; Here is a brief key to the formula names
;
; Type 
; ~~~~
; M    - mandelbrot
; J    - julia
; test - just what it says - only included for you to see my
;	   train of thought.  I used these to develop the Fn types
;
; Trap
; ~~~~
; Cir  - indicates circular trap
; Ell  - ellipse
; Prb  - parabola
; Hyp  - hyperbola
;
; 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:
; ~~~~~~~~
; L    - only last value of z that meets criterion is used
; A    - uses all values of z that meet criterion
; d    - uses 'distance' from center fro traps
; e    - uses 'distance' from Edge - the best kind usually
; de   - uses 'distance from end'
; 1    - adds constant if it hits
; z    - uses value of z for colouring, not just distance
; m    - subtracts normal colouring (no iteration bands)
;
; NB/ for the options, e-m are generally the best - colouring from the
; edge means that points in the center are most strongly coloured, 
; points at edge are almost the same colour as the 'background'.
; By choosing 'm' you are also eliminating the normal iteration bands.  
; You can still see them if you select outside = iter instead of real
;
; Also, the colour factor can be positive or negative and this will 
; give different results - colouring goes opposite way
;
; Have a look at this file in NotePad  to read the other
; comments and  finally - I hope you enjoy them
;
}


 
MCol1(XAXIS) {
c=pixel, done=0:

z= done
z=z*z+c
done = z
z= fn1(z) + imag(p1)*|fn2(c)| ;
|done| < real(p1)

}


JCol1(ORIGIN) {
c=p1, done=pixel:

z= done
z=z*z+c
done = z
z= fn1(z) + imag(p2)*|fn2(c)| ;
|done| < real(p2)

}


MCol2(XAXIS) {
; try 4,0.5,flip,conj
c=pixel, done=0:

z= done
z=z*z+c
done = z
a= fn1(atan(imag(z)/real(z))) + fn2(z)
z=  imag(p1)* cos(a) + flip (imag(p1) * sin(a));

|done| < real(p1)

}



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
}


MCol5(XAXIS) { ; re(p2)=bailout, try fn1=tan
                        ; p1=1,0  p2 = 8,1
c=pixel, w=0, s=0:

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


MCol5-m {
c=pixel, z=0, s=0, fin=0, count=0:

z=z*z+c
count=count+1
s=s+p1*1/z
If (|z| > real(p2))
	s=fn1(|s|)*imag(p12)
	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)
}



MCol7(XAXIS) { 
; try 2, 0, 8, 0.1, flip, sqrt
c=pixel, temp=0, q=0:

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

|temp| <= real(p2)
}

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

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

|temp| <= real(p2)
}

MCol9(XAXIS) { ; rP3-bailout, iP3-strength of colouring
               ; rP2-how close bands are, iP2-pallette shift
c=pixel, temp=0, q=0:

z=temp*temp+c
temp=z
q=q+imag(p3)/(|z+p1|-real(p2))*z

z=q+256*imag(p2)
|temp| <= real(p3)

}

;---------------------------------------------------------------------------
; trap 1 - circle
; uses only the last z which lies in the circle

M-Cir-L-z(XAXIS) { ; trap: p1 = centre of circle
; rP2 = sqr(radius), rP3 = bailout
; iP3 = colour factor
; try 1, 1, 0.5, 0, 16, 20

c=pixel, z=0, q=0, fin=0:
z=z*z+c
IF (|z-p1|<real(p2))
	q=real(z)*imag(p3)
END IF
IF (|z| > real(p3))		; ie going to bailout/quit
	z=q-255*floor((q-1)/real(255))
	fin=-1
END IF

fin == 0
}

; trap 2 - circle
; uses the last z which lands in circle, multiplied by the distance of that z from 
; the centre of the circle

M-Cir-L-z*d(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
IF (|z-p1|<real(p2))
	q=real(z)*imag(p3)*|z-p1|
END IF

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

fin ==0
}



:-------------------------------------------------------
; trap 3, circle (last (z)n: distance from edge of circle* z) 

M-Cir-L-z*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= z * test*imag(p3)
END IF
IF (|z| > real(p3))
	z=q-255*floor((q-1)/real(255))
	fin=-1
END IF

fin ==0

}


M-Cir-L-z*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= z * test*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 4- circle, (last z(n), distance from edge of circle)

M-Cir-L-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=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-L-em(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=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 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-em(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 6- rectangular hyperbola (adds 1 if it hits at all)

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

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

fin == 0
}


:-------------------------------------------------------
; trap 7- rectangular hyperbola, (add 1 each time it hits)

M-Hyp-A-1(XAXIS) { ; trap- p1 = shift, rP2=distance
		; rP3 = bailout
c=pixel, z=0, q=0, fin=0:
z=z*z+c
IF (|(real(z)-real(p1))*(imag(z)-imag(p1))| < real(p2))
	q=q+imag(p3)
END IF

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

fin == 0
}


:-------------------------------------------------------
; trap 8- rectangular hyperbola, (z(n) of last time it hit)

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

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

fin == 0
}

}

:-------------------------------------------------------
; trap 9- rectangular hyperbola, (adds all z(n) when it hits)

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

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

fin == 0
}


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

IF (|(real(z)-real(p1))*(imag(z)-imag(p1))| < real(p2))
	q=q+z
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 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-em(XAXIS) { ; trap- p1 = shift, rP2=distance
		; iP2=palette shift, rP3 = bailout
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 11 - scrapped

:-------------------------------------------------------
; trap 12- ellipse (all z(n): adds distance from centre of ellipse)

M-Ell-A-d(XAXIS) { ; trap- p1 = centre of ellipse, 
; rP2 = sqr(x radius), iP2= sqr(y radius)
; rP3 = bailout, iP3 = colour factor
c=pixel, z=0, q=0, fin=0:
z=z*z+c

test=sqr(real(z-p1))/real(p2)+sqr(imag(z-p1))/imag(p2)

IF (test<1)
	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-Ell-A-dm(XAXIS) { ; trap- 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+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 13- ellipse -(all z(n): adds distance from end of ellipse)

M-Ell-A-de(XAXIS) { ; trap- p1 = centre of ellipse
; rP2 = sqr(x radius), iP2= sqr(y radius)
; rP3 = bailout, iP3 = colour factor
; sign of radii determines which 'end' to colour from

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

test1=real(z-p1)
test2=imag(z-p1)
IF (((sqr(test1)/|real(p2)|) + (sqr(test2)/|imag(p2)|)) < 1 )
	q=q + test1+real(p2)+test2+imag(p2)
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-Ell-A-dem(XAXIS) { ; trap- p1 = centre of ellipse
; rP2 = sqr(x radius), iP2= sqr(y radius)
; rP3 = bailout, iP3 = colour factor
; sign of radii determines which 'end' to colour from

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

test1=real(z-p1)
test2=imag(z-p1)
IF (((sqr(test1)/|real(p2)|) + (sqr(test2)/|imag(p2)|)) < 1 )
	q=q + test1+real(p2)+test2+imag(p2)
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 14- ellipse -(last z(n): adds distance from end of ellipse)

M-Ell-L-de(XAXIS) { ; trap- p1 = centre of ellipse
; rP2 = sqr(x radius), iP2= sqr(y radius)
; rP3 = bailout, iP3 = colour factor
; sign of radii determines which 'end' to colour from

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

test1=real(z-p1)
test2=imag(z-p1)
IF (((sqr(test1)/|real(p2)|) + (sqr(test2)/|imag(p2)|)) < 1 )
	q= test1+real(p2)+test2+imag(p2)
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-Ell-L-dem(XAXIS) { ; trap- p1 = centre of ellipse
; rP2 = sqr(x radius), iP2= sqr(y radius)
; rP3 = bailout, iP3 = colour factor
; sign of radii determines which 'end' to colour from
c=pixel, z=0, q=0, count=0, fin=0:
z=z*z+c
count=count+1

test1=real(z-p1)
test2=imag(z-p1)
IF (((sqr(test1)/|real(p2)|) + (sqr(test2)/|imag(p2)|)) < 1 )
	q= test1+real(p2)+test2+imag(p2)
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 15- ellipse -(all z(n): adds distance from end of ellipse, scaled)


M-Ell-A-de2(XAXIS) { ; trap- p1 = centre of ellipse
; rP2 = sqr(x radius), iP2= sqr(y radius)
; rP3 = bailout, iP3 = colour factor
; sign of radii determines which 'end' to colour from
c=pixel, z=0, q=0, fin=0:
z=z*z+c

test1=real(z-p1)
test2=imag(z-p1)
IF (((sqr(test1)/|real(p2)|) + (sqr(test2)/|imag(p2)|)) < 1 )
	q=q + (test1+real(p2)+test2+imag(p2))/(real(p2)+imag(p2))
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-Ell-A-de2m(XAXIS) { ; trap- p1 = centre of ellipse
; rP2 = sqr(x radius), iP2= sqr(y radius)
; rP3 = bailout, iP3 = colour factor
; sign of radii determines which 'end' to colour from

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

test1=real(z-p1)
test2=imag(z-p1)
IF (((sqr(test1)/|real(p2)|) + (sqr(test2)/|imag(p2)|)) < 1 )
	q=q + (test1+real(p2)+test2+imag(p2))/(real(p2)+imag(p2))
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:
;
; p1 is focus, P2 =gives eqn of directrix:
; y=mx + k (m=rP2, k=iP2)
; it can be shown that if PN is perp to directrix, where
; P=(xp,yp) (general point) 
; N=(xn,yn), point on directrix
; then
; (yp-yn)/(xp-xn)=-1/m
; and
; yn=m*xn+k
; giving:
;
; xn = m*(yp-k)+xp
;      ----------
;         1+m^2
; yn = m*xn+k
;
; therefore P is found, PN is known therefore parabola defined (PF=PN)
; note: below, N= (xn,yn), (instead of using separate variables)

M-Prb-A-e(XAXIS) { ; p1 = focus, rP2 = directrix:
                ; y=rP2 * x + iP2
		; rP3 = bailout, iP3 = colour factor
c=pixel, z=0, q=0, fin=0, m=real(p2), k=imag(p2):
z=z*z+c

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))
	fin=-1
END IF

fin == 0

}

; instead of using y=mx+c, use another point, G(xg,yg) such that FG is perp to 
; directrix, G is on directrix.  Then it can be shown that
; m= (xg-xf)/(yf-yg) and
; k= -m*xg+yg
; provided that yf <> yg  (otherwise division by zero gives blank screen)


M-Prb2-A-e(XAXIS) { ; p1 = focus, P2 = point 'opposite' focus
                ; on directrix
		; rP3 = bailout, iP3 = colour factor
c=pixel, z=0, q=0, fin=0
m=(real(p2)-real(p1))/(imag(p1)-imag(p2))
k=-m*real(p2)+imag(p2):
z=z*z+c
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))
	fin=-1
END IF

fin == 0
}


M-Prb2-A-em(XAXIS) { ; p1 = focus, P2 = point 'opposite' focus
                ; on directrix
		; rP3 = bailout, iP3 = colour factor
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(ORIGIN) { ; focus=0 , rP2 = point 'opposite' focus
                ; on directrix
		; rP3 = bailout, iP3 = colour factor
		; P1 = c
c=P1, z=pixel, q=0, fin=0
m=real(p2)/(-imag(p2))
k=-m*real(p2)+imag(p2):
z=z*z+c
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))
	fin=-1
END IF
fin ==0

}

J-Prb2-A-em(ORIGIN) { ; focus=0 , rP2 = point 'opposite' focus
                ; on directrix
		; rP3 = bailout, iP3 = colour factor
		; P1 = c
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(XAXIS) { ; trap- p1 = centre of ellipse
; rP2 = sqr(x radius), iP2= sqr(y radius)
; rP3 = bailout, iP3 = colour factor
; ellipse (all z(n): adds distance from edge of ellipse)

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))
	fin=-1
END IF

fin == 0
}

M-Ell-A-em(XAXIS) { ; trap- p1 = centre of ellipse
; rP2 = sqr(x radius), iP2= sqr(y radius)
; rP3 = bailout, iP3 = colour factor
; ellipse (all z(n): adds distance from edge of ellipse)

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 (iP2,0) but can be tilted
; tilted ellipse:

; ex= sin(a)*(y-yc)+cos(a)*(x-xc)
; ey= (y-yc-exsin(a))/cos(a)
; where a is angle of tilt, (xc,yc) is centre of ellipse
; ex and ey are x and y coordinates relative to the elipses (tilted) axes
; the eqn of ellipse then becomes 
; sqr(ex/r1)+sqr(ey/r2)=1
; where r1 and r2 are the respective radii
; In the first formula, yc is 0 always.  in the second,
; yc can be set ,but the bailout cannot



M-Ell2-A-e(XAXIS) { ;  p1 = radii of ellipse, rP2 = tilt
                ; iP2=centre of ellipse(on real axis)
                ; rP3 = bailout, iP3 = colour factor
; ellipse (all z(n): adds distance from edge of ellipse)

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

fin == 0
}

M-Ell2-A-em(XAXIS) { ;  p1 = radii of ellipse, rP2 = tilt
                ; iP2=centre of ellipse(on real axis)
                ; rP3 = bailout, iP3 = colour factor
; ellipse (all z(n): adds distance from edge of ellipse)

c=pixel, z=0, q=0, count=0, fin=0
cost=cos(real(p2)), sint=sin(real(p2)):
z=z*z+c
count=count+1
ex=cost*(real(z)-imag(p2))+sint*(imag(z))
ey=(imag(z)-ex*sint)/cost
test=sqr(ex/real(p1))+sqr(ey/imag(p1))
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(XAXIS) { ;  p1 = radii of ellipse, 
           ; P2=centre of ellipse
           ; rP3 = tilt, iP3 = colour factor
; ellipse (all z(n): adds distance from edge of ellipse)
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))
	fin=-1
END IF

fin == 0
}

M-Ell3-A-em(XAXIS) { ;  p1 = radii of ellipse, 
           ; P2=centre of ellipse
           ; rP3 = tilt, iP3 = colour factor
; ellipse (all z(n): adds distance from edge of ellipse)
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
}


;-----------------------------------------------------------------------------
; the following produces a 'pond ripple' effect around the desired point.  
; d is the distance from desired point. The curve should
; exponentionally decreases while oscillating (!)
; For best results it should be a maximum at the bgeinning for large central maximum:
;
; z(d)=ae^(-bd)*cos(cd+E)                    ; first attempt
; z'(d)=ae^(-bd)[-b*cos(cd+E)-c*sin(cd+E)]
; we require z'(0)=0 for a maximum at d=0
; ae^(-bd)<> 0 for d real
; -b*cos(E) = c*sin(E)
; tan(E) = -b/c
;
;In fact, abs(cos(cd+E)) is used instead, but E still applies
;
;

test1  {;p1 = point
a=real(p2)  				; colour factor
b=real(p3)  				; decay factor
c=imag(p2)  				; oscillator factor
ep=atan(-b/c) 				; see above for E
d=cabs(pixel-p1)			; distance from point

z=a*(exp(-d*b))*abs(cos(d*c+ep))

1>2					; alway bailout immediately
}


; this is a by product of the above!
Diffraction1{; try 0,0,100,500,0.1,0 for brilliant pattern - it
; is just *concentric circles* !
d=|pixel-p1|

z=100+real(p2)*(exp(-d*real(p3))) *abs(cos(d*imag(p2)))

1>2
}


Diffraction2{; try 0,0,100,500,0.1,0 for brilliant pattern - it
; is just *concentric circles* !
d=sqrt(|pixel-p1|)

z=100+real(p2)*(exp(-d*real(p3))) *abs(cos(d*imag(p2)))

1>2
}


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

M-Fn1-A-m {; this has a bug - use M-Fn2-A-M instead
; p1 is center of circles, iP3 = bailout
a=real(p2)  				; colour factor
b=real(p3)  				; decay factor
o=imag(p2)  				; oscillator factor
ep=atan(-b/c) 				; 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
}

;First try has a mistake  - ep=atan(-b/c) should be atan(-b/o).  Old version
; kept because I already made a fractal with it!

M-Fn2-A-m {; 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
}


;----------------------------------------------------------------------
; a different way of doing Fn2 - not an exponential decay though
M-Fn3-A-m{;  p1 is center of circles
; try 1,1,10,10,16,0, fn1=ident, fn2=ident
; fn1 is the oscillating function, fn2 is the decay function
; rP3 = bailout
a=real(p2)     		; colour factor
b=imag(p2)		; oscillate factor
z=0, c=pixel, q=0, count=0, fin=0:
z=z*z+c
count=count+1
d=cabs(z-p1)
q=q+abs(sin(b*fn1(d)))/fn2(d)
IF (|z| > real(p3))
	q=q*a
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1
END IF
fin==0
}		

; the test formula for developing the above
test3 {; try 1,1,10,10,16,0, fn1=ident, fn2=ident
; fn1 modifies the oscillating function, fn2 modifies the decay function
a=real(p2)     		; colour factor
b=imag(p2)		; oscillate factor
d=cabs(pixel-p1)
z=a*abs(sin(b*fn1(d)))/fn2(d)
1>2
}

;----------------------------------------------------------------------
; modified test1 to give elliptical rings

test4  {;p1 = point, iP3=aspect ratio (height/width) for rings
; try 0,0,10,10,1,2
a=real(p2)  				; colour factor
c=imag(p2)  				; oscillator factor
b=real(p3)  				; decay factor
ep=atan(-b/c) 				; see above for E
d=pixel-p1				; distance from point
d=cabs(imag(p3)*real(d)+flip(imag(d)))	; but distorted for ellipses


z=a*(exp(-d*b))*abs(cos(d*c+ep)):

1>2					; alway bailout immediately for test
}

; implementing test4
M-Fn4-A-m {; p1 is point which is center of colouring
;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
}
;----------------------------------------------------------------------
; attempt a simple exponential decay, so the formula can be taken apart
test5.1 {; p1 is point
a= real(p2) 			; colour factor
b= imag(p2)			; decay factor
d= cabs(pixel-p1)		; distance from point
z= a*exp(-d*b):
1>2				;bailout
}

test5.2 {; p1 is point
a= real(p2) 			; colour factor
b= imag(p2)			; decay factor
d= |pixel-p1| 			; sq of distance from point
z= a*exp(-d*b):
1>2				;bailout
}

; 5.2 gives a better effect than  5.1 - the peak for 5.1 is 'sharp'
; but for 5.2 is smooth.  Compare graphs of y=e^(-abs(x)) (5.1)
; with y=e^(-x*x) (5.2)
;  Also, 5.2 can be separated:
; d=x*x + y*y
; e^(x*x + y*y)=e^(x*x)*e^(y*y)
; 
; In test5.3, x and y decay with separate decay factors - giving ellipse
; as you can deduce

test5.3 {; p1 is point
; bx=by gives test5.2
a= real(p2) 			; colour factor
bx= real(p3)			; decay factor in x direction
by= imag(p3)			; decay factor in y direction
x=real(pixel-p1)
y=imag(pixel-p1)
z= a*(exp(-x*x*bx)*exp(-y*y*by)):
1>2				;bailout
}

; Adding functions for decay in different directions gives more
; possibilties of shapes.  Good combinations include (either way round):
; sqr/sqrt, sqr/ident, sqrt/ident

test5.4 {; p1 is point
; fn1=fn1=ident gives test5.3
; 
a= real(p3) 			; colour factor
bx= real(p2)			; decay factor in x direction
by= imag(p2)			; decay factor in y direction
x=real(pixel-p1)
y=imag(pixel-p1)
z= a*(exp(-fn1(x*x)*bx)*exp(-fn2(y*y)*by)):
1>2				;bailout
}

;Try test 5.4 in a fractal

M-Fn5-A-m {; use M-Fn6-A-m for corrected version!
; 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
a= imag(p3) 			; colour factor
bx= real(p2)			; decay factor in x direction
by= imag(p2)			; decay factor in y direction
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
}

;----------------------------------------------------------------------
; unlike with earlier functions, this one has the possibility to be negative - 
; producing colour 0 only.  To wrap around to end:
; 0 --> 255
; -1 --> 254
; -255 --> 0, but 0 must go to 255
; can be done by using floor instead of trunc at end of formula
; e.g.(-5/255)=-0.01960
; floor (-5/255)= -1
; trunc (-5/255)= 0
; therefore floor will wrap all colours correctly, not just positive ones
; and will always be used in future


M-Fn6-A-m {; 
; 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
a= imag(p3) 			; colour factor
bx= real(p2)			; decay factor in x direction
by= imag(p2)			; decay factor in y direction
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 {; use Fn7b - simpler, better, faster
; 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.
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

}

; oops!  the above variable base for exponentiation is unnecessary:
; a^bx  = exp(ln(a^bx)) = exp(bx*ln(a))
; = exp(cx) where c = b ln(a)
; ie variables a and b can be replaced with c, one variable.  Therefore use the below:

M-Fn7b-A-m {; use Fn7b - simpler, better, faster
; rP1 is decay factor, iP1 is colour factor, rP2 is bailout
; try 0.1,10,128,0  for super smooth colouring.
a= imag(p1) 			; colour factor
b= real(p1)				; decay factor
bail=real(p2)			; bailout
z=0, c=pixel, q=0, count=0, fin=0:
z=z*z+c
count=count+1
q=q+exp(-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

}

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

; some 'funny' behaviour was noticed using M-Fn6 and fn1=fn2=flip:
; exp(-flip(x*x)*bx)*exp(-flip(y*y)*by)
; = exp(-bx*i*x*x)*exp(-by*i*y*y)
; = exp(i*(-bx*x*x -by*y*y))
; = cos(-bx*x*x - by*y*y)+ i sin(-bx*x*x - by*y*y)
; since only real part is used for colouring, the following should
; produce the same result (and does, slightly faster)
;

M-Fn8-A-m {; 
; P1 is point, center of colouring, rP3 is bailout
; try 5,5,0.005,0.001,128,10 for some nice diffraction patterns
a= imag(p3) 			; colour factor
bx= real(p2)			; oscillate factor in x direction
by= imag(p2)			; oscillate factor in y direction
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+cos(-x*x*bx-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
}

;----------------------------------------------------------------------
;Extending what was tried in M-Fn8 - use any function of x and y inside
; the cos() bit.  Fn1=Fn2=sqr will give Fn8
M-Fn9-A-m {; 
; P1 is point, center of colouring, rP3 is bailout
; try 5,5,0.05,0.01,128,10 for some nice diffraction patterns
; use fn1=fn2=sqr for M-Fn8
bx= real(p2)			; oscillate factor in x direction
by= imag(p2)			; oscillate 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+cos(-fn1(x)*bx-fn2(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
}


;--------------------------------------------------------------------
; lets mess with a newton!  First try ripples (M-Fn1)

Newton1 { ; P1 is center or colouring, iP3 bailout
;-----------Initialise----------------------------------------
a=real(p2)  				; colour factor
o=imag(p2)  				; oscillator factor
b=real(p3)  				; decay factor
ep=atan(-b/o) 				; epsilon to give maximum
z=pixel, q=0, count=0, fin=0:
; -----------Newton iteration---------------------------------
z2=z*z
z3=z2*z
fz=z3-1
z=z-fz/(3*z2)
; -----------REM ripple function bit!-------------------------
count=count+1
d=cabs(z-p1)				; distance from point
q=q+(exp(-d*b))*abs(cos(d*o+ep))
IF (|fz| < imag(p3))
	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 - 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)

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
	
}

;--------------------------------------------------------------------
; previous value of z to give angle of centre of circle
; - achieved just by putting angle calculations before iteration

M-EllM2-A-e-m (XAXIS){;
; ellipse - rP1, iP1 = radii of ellipse
; rP2=distance of circle  from origin
; iP2=start angle 
; 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:

;----------- GET ANGLE FROM LAST Z -------------------------
IF (real(z) != 0)			; avoid DBZ error
	angle=atan(imag(z)/real(z))	; use last value of z to get angle
	IF (real(z)<0)		; reflex angle
		angle=angle+pi
	END IF
ELSE
	angle=0
END IF
angle=angle+imag(p2)

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

z=sqr(z)+c				; Mandelbrot iteration
count=count+1

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

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 {;
; 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 {;
; 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:
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
	
}





;------------------------------------------------------------------
; Now try variable radius factor, depending on fn1(|z|)
; The position of centre can be fixed
; (other possibilities include just varying one radius, so
; changing the aspect ratio of ellipse)

M-EllM4-A-e-m {; variable 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)*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
	
}

;----------------------------------------------------------------
;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 

M-P001-m {;
; real(p1) = colour factor 1
; imag(p1) = colour factor 2
; real(p2) = bailout
; imag(p2) = offset

; try 10,10,20,0,log
z=pixel, c=pixel, count=0, small=1000000, fin=0:
z=z*z+c
count=count+1
IF (abs(real(z)-imag(p2))< small)
	small = abs(real(z))
END IF

IF (|z| > real(p2))
	small=fn1(small*real(p1))*imag(p1)
	z= small - 255*floor((small-1)/real(255))
	z=z-count
	fin=1
END IF

fin==0

}

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


M-P002-m {;

z=pixel, c=pixel, count=0, small=1000000, fin=0:
lastz=z
z=z*z+c
count=count+1
IF ( abs(real(z)/real(lastz)) < imag(p2))
	small = abs(real(z)/real(lastz))
END IF

IF (|z| > real(p2))
	small=fn1(small*real(p1))*imag(p1)
	z= small - 255*floor((small-1)/real(255))
	z=z-count
	fin=1
END IF

fin==0

}



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


M-P003-m { ;
; P1 is a point - smallest dis from P1 is used
z=0, c=pixel, count=0, fin=0, sm=100000:

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

d=|z-P1|
IF (d<sm)
	sm=d
END IF

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

fin==0

}



M-P004-m { ;
; P1 is a point - real(z) of smallest dis from P1 is used
z=0, c=pixel, count=0, fin=0, small=100000, q=0:

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

d=|z-P1|
IF (d<small)
	small=d
	q=real(z)
END IF

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

fin==0

}


M-P005-m { ;
; P1 is a point - imag(z) of smallest dis from P1 is used
z=0, c=pixel, count=0, fin=0, small=100000, q=0:

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

d=|z-P1|
IF (d<small)
	small=d
	q=imag(z)
END IF

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

fin==0

}

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

M-P006-m { ;
; smallest |z| is used
z=0, c=pixel, count=0, fin=0, small=100000, q=0, mz=0:

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

mz=|z|
IF (mz<small)
	small=mz
END IF

IF (mz > real(p2))
	small=small*imag(p2)
	z= small - 255*floor((small-1)/real(255))
	z=z-count
	fin=1
END IF

fin==0

}


M-P007-m { ;
; real(z) of smallest |z| is used
z=0, c=pixel, count=0, fin=0, small=100000, q=0, mz=0:

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

mz=|z|
IF (mz<small)
	small=mz
	rz=real(z)
END IF

IF (mz > real(p2))
	rz=rz*imag(p2)
	z= rz - 255*floor((rz-1)/real(255))
	z=z-count
	fin=1
END IF

fin==0

}

M-P008-m { ;
; imag(z) of smallest |z| is used
z=0, c=pixel, count=0, fin=0, small=100000, mz=0:

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

mz=|z|
IF (mz<small)
	small=mz
	iz=imag(z)
END IF

IF (mz > real(p2))
	iz=iz*imag(p2)
	z= iz - 255*floor((iz-1)/real(255))
	z=z-count
	fin=1
END IF

fin==0

}

;-------------------------------------------------------------
; 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


Test-Ast { ;
c=pixel, z=0:
d= (abs(real(c-P1)/real(P2)))^(2/3) + (abs(imag(c-P1)/imag(P2)))^(2/3)
IF (d < 1)
	z= (1-d)*imag(p3)
END IF

1==0
}

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

test-fan { ;
z=pixel:
temp=z-P1
d=|temp|
th=atan(imag(temp)/real(temp)) 
z=imag(p3)*exp(-d*real(p2))*cos(imag(p2)* th)
1==0
}

test-fan2 { ;
z=pixel:
zm=z-P1

z=imag(p3)*exp( (-1)*|zm|*real(p2)) * cos(imag(p2)* atan(imag(zm)/real(zm)) )
1==0
}


M-FN-fan-m { ;
; 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

}

test-fan3 { ; using a quadratic for the deacy
z=pixel:
zm=z-P1

h =1-real(p2)*|zm|
IF (h>0)
	z=imag(p3)*h * cos(imag(p2)* atan(imag(zm)/real(zm)) )
ELSE 
	z=0
END IF

1==0
}

test-fan4 { ; using a quadratic for the deacy, dips at centre
; the quadratic has roots at x=0 and x=real(p2).
; the maximum value should also be fixed ie independent of
; real(p2)
; width=real(p2)
z=pixel:
zm=z-P1

d=cabs(zm)
h =d*(real(p2)-d)
IF (h>0)
	z=imag(p3)/sqr(real(p2))*h * cos(imag(p2)* atan(imag(zm)/real(zm)) )+100
ELSE 
	z=100
END IF

1==0
}


test-fan5 { ; as per fan4 but use |zm| instead of cabs(zm)
z=pixel:
zm=z-P1

d=|zm|
h =d*(real(p2)-d)
IF (h>0)
	z=imag(p3)/sqr(real(p2))*h * cos(imag(p2)* atan(imag(zm)/real(zm)) )+100
ELSE 
	z=100
END IF

1==0
}


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

;use test-fan4 in a fractal

M-FN-fan2-m { ;
; P1 = center of fan.
; real(p2) = width
; imag (p2) = number of arms
; real(p3) = bailout
; imag(p3) = colour factor
; try 0.5, 0.5, 1, 6, 16, 80
z=0, c=pixel, count=0, fin=0, q=0:

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

zm=z-P1
d=cabs(zm)
h=d*(real(p2)-d)
IF (h>0)
	q=q+ imag(p3)/sqr(real(p2))*h * cos(imag(p2)* atan(imag(zm)/real(zm)) )
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

}


;----------------------------------------------------------------------
; lemniscate(?) shape
test-lem1 {
z=pixel:
z2=z-P1
a2=sqr(real(p2))

theta=atan(imag(z2)/real(z2))
test=a2*cos(2*theta)-|z2|
IF (test>0)
	col=test
ELSE
	col=0
END IF


z=imag(p3)*col+(100,0)

1==0
}

; try to get a figure-of-8 shape, based on lem1

test-lem2 {
z=pixel:
z2=z-P1
a2=sqr(real(p2))
b2=sqr(imag(p2))

theta=atan(imag(z2)/real(z2))
test=a2*cos(2*theta)-|z2|
IF (test>0)
	col=test
ELSE
	col=0
END IF

; after cut off limit, start decreasing
; sqr(b)=b2 used, since all other dimensions are squares
; tests show that if a2:b2 is fixed, the realative 'width' of band
; produced is constant

IF (test > b2) 			;start decreasing colouring again
	col=col-2*(test- b2)	;need to get colour function to decrease
				; ie 1*(test-b) would get it to be level
				; 2*(test-b) gets it to go down again, making sure
				; the colours match at 'join'
END IF

IF (col<0)
	col=0
END IF

z=imag(p3)*col+(100,0)

1==0
}

; using test-lem2, if imag(p2)~0.6*real(p2), then
; a nice 8 shape is produced.  Use this fact to change
; imag(p2) automatically ie do not need to put in imag(p2)
; Also, the max colour obtained varies wit the size of lem.
; try to fix this by dividing by a2 (sqr(a)) or a

; imag(p2) = 0.6*real(p2)
; sqr(imag(p2)) = 0.36 * sqr(real(p2))
; b2 = 0.36 * a2

test-lem3 {
z=pixel:
z2=z-P1
a2=real(p2)
theta=atan(imag(z2)/real(z2))
test=a2*cos(2*theta)-|z2|
IF (test>0)
	col=test
ELSE
	col=0
END IF

b2=a2*0.36
IF (test > b2) 			;start decreasing colouring again
	col=col-2*(test- b2)
END IF

IF (col<0)
	col=0
END IF

z=imag(p3)/a2*col+(100,0)

1==0
}

; this can be optimized into lem4, which is actually easier to understand


test-lem4 {
z=pixel:
z2=z-P1
a2=real(p2)
b2=a2*0.36

theta=atan(imag(z2)/real(z2))
test=a2*cos(2*theta)-|z2|
IF (test>0)
	IF (test <b2)
		col=test
	ELSE
		col=2*b2-test
		IF (col<0)
			col=0
		END IF
	END IF
ELSE
	col=0
END IF

z=imag(p3)/a2*col+(100,0)

1==0
}




; try to use lem4 in a fractal.  imag(p2) is free. use this
; for number of arms

M-Lem1-A-e-m { ; figure of 8 shape, based
; on the leminsate shape
; P1 is center
; real(p2)=width
; imag(p2)=number of arms (=2 for 8 shape) (even numbers only)
; real(p2)=bailout
; imag(p3)=colour factor

; try 1, 1, 0.5, 2, 16, 10

z=0, c=pixel, q=0, fin=0, count=0
a2=real(p2)
b2=a2*0.36:

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

z2=z-P1
theta=atan(imag(z2)/real(z2))
test=a2*cos(imag(p2)*theta)-|z2|
IF (test>0)
	IF (test <b2)
		col=test
	ELSE
		col=2*b2-test
		IF (col<0)
			col=0
		END IF
	END IF
ELSE
	col=0
END IF

q=q+col

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

}


;-----------------------------------------------------------------------------
; generalize lem4 - use other functions instead of cos
; Also, fix the theta variable to give angle from -pi to pi
; since atan can't distingush between (-y/x) and (y/-x)



test-lem5 {
; try 0,0,1,2,16,80,cos
z=pixel:
z2=z-P1
a2=real(p2)
b2=a2*0.36

theta=atan(imag(z2)/real(z2))

IF (real(z2)<0)
	theta=theta+pi
END IF


test=a2*fn1(imag(p2)*theta)-|z2|
IF (test>0)
	IF (test <b2)
		col=test
	ELSE
		col=2*b2-test
		IF (col<0)
			col=0
		END IF
	END IF
ELSE
	col=0
END IF

z=imag(p3)/a2*col+(100,0)

1==0
}


test-theta { ; checking my theta function
z=pixel:
theta=atan(imag(z)/real(z))
IF (real(z)<0)
	theta=theta+pi
END IF
IF (theta > pi)
	theta=theta-(2*pi)
END IF
z=real(p1)*theta
1==0
}


;----------------------------------------------------
; try some new shapes

; sine based curve


test-sin1 { ;
z=pixel:
z2=pixel-P1

test1=imag(p2)*sin(real(p2)*real(z2)) *exp(-real(p2)/5*abs(real(z2)))
IF ( abs(imag(z2))<abs(test1) )
	col=imag(p3)
ELSE
	col=0
END IF

z=col
1==0
}


; a good 'aspect ratio' for sine shape was found with real(p2)=31, imag(p2)=1
; and using the (real(p2)/5) in the deacy bit - this will be used in next test

test-sin2 { ;
try 0,0,15,0.5,0,30
z=pixel:
z2=pixel-P1

test1=imag(p2)*sin(real(p2)*real(z2)) *exp(-real(p2)/5*abs(real(z2)))
IF ( abs(imag(z2))<abs(test1) )
	col=imag(p3)*(abs(test1)-abs(imag(z2)))
ELSE
	col=0
END IF

z=col
1==0
}

;try test-sin2 in a fractal
;-------------------------------------------------------------


M-Sin1-A-e-m {; a sine shape superimposed onto MSet
; P1 = center of shape
; real(p2)= inverse of 'width'- larger gives narrower shape
; imag(p2)= 'height' of shape
; real(p3)= bailout
; imag(p3)= colour factor
; try 0.5,0.5,10,0.3,16,10

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

z=z*z+c
z2=z-P1
count=count+1

test1=abs(imag(p2)*sin(real(p2)*real(z2)) *exp(-real(p2)/5*abs(real(z2))))
IF ( abs(imag(z2))< test1 )
	q=q+(test1-abs(imag(z2)))
END IF

q=q+col

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

}

;------------------------------------------------------------
; generalize the above - not just sine, any function
M-shape-A-e-m {; a 'shape' superimposed onto MSet
; ie finite trap, not a function
; P1 = center of shape
; real(p2)= inverse of 'width'- larger gives narrower shape
; imag(p2)= 'height' of shape
; real(p3)= bailout
; imag(p3)= colour factor
; try 0.5,0.5,10,0.3,16,10, cos

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

z=z*z+c
z2=z-P1
count=count+1

test1=abs(imag(p2)*fn1(real(p2)*real(z2)) *exp(-real(p2)/5*abs(real(z2))))
IF ( abs(imag(z2))< test1 )
	q=q+(test1-abs(imag(z2)))
END IF

q=q+col

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

}

;--------------------------------------------------------------------------
; using Shape formula above along with a smooth algorithm in a PHC formula
; to create a HC image by shrinking.  NB, specifc fractal in mind, so paramters
; are altered accordingly

M-shape2-A-e-m {; a 'shape' superimposed onto MSet
; ie finite trap, not a function
; P1 = center of shape
; real(p2)= inverse of 'width'- larger gives narrower shape
; imag(p2)= 'height' of shape
; real(p3)= bailout
; imag(p3)= colour factor
; try 0.5,0.5,10,0.3,16,10, cos

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

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


IF (whitesq==0)
	smooth=exp(-2*|z|)
	q=q+smooth*15
ELSE
	z2=z-P1
	test1=abs(imag(p2)*fn1(real(p2)*real(z2)) *exp(-real(p2)/5*abs(real(z2))))
	IF ( abs(imag(z2))< test1 )
		q=q+(test1-abs(imag(z2)))*imag(p3)
	END IF
END IF


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

fin==0

}


;-----------------------------
; using a julibrot, my formula 1 therefore JB1 is the type
; attempt using smoothed colours

JB1-Smooth(ORIGIN){;real p1 = real of c, imag p1=imag of z0
; rP2 is decay factor, iP2 is colour factor, rP3 is bailout
; try -,-,0.1,10,128,0  for super smooth colouring.
a= imag(p2) 			; colour factor
b= real(p2)				; decay factor
bail=real(p3)			; bailout
z = real(pixel)+ flip(imag(p1))
c = real(p1)+flip(imag(pixel))
q=0, count=0, fin=0:

z=z*z+c
count=count+1
q=q+exp(-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


}

JB1-Ripples(ORIGIN){;real p1 = real of c, imag p1=imag of z0
; P2 is center of colouring
o=10		  				; oscillator factor
b=real(p3)  				; decay factor
a=imag(p3)  				; colour factor
ep=atan(-b/o) 				; 

z = real(pixel)+ flip(imag(p1))
c = real(p1)+flip(imag(pixel))
q=0, count=0, fin=0:

z=z*z+c
count=count+1
d=cabs(z-p2)				; distance from point
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

ENDIF

fin==0
}


; specialised version of above for Ailen_Skull to shift colours about

JB1-Ripples-b(ORIGIN){;real p1 = real of c, imag p1=imag of z0
; P2 is center of colouring
o=10		  				; oscillator factor
b=real(p3)  				; decay factor
a=imag(p3)  				; colour factor
ep=atan(-b/o) 				; 

z = real(pixel)+ flip(imag(p1))
c = real(p1)+flip(imag(pixel))
q=0, count=0, fin=0:

z=z*z+c
count=count+1
d=cabs(z-p2)				; distance from point
q=q+(exp(-d*b))*abs(cos(d*o+ep))
IF ( |z| > 16)
	q=q*a+60
	z=q-255*floor((q-1)/real(255))
	z=z-count
	fin=-1

ENDIF

fin==0
}