# julia set for z -> z^2+c
par am=.25,q=0
cx=am*cos(q)
cy=am*sin(q)
r=sqrt(sqrt((x-cx)^2+(y-cy)^2))
th=atan2(y-cy,x-cx)/2
s=sign(ran(1)-.5)
init x=.1
x'=s*r*cos(th)
y'=s*r*sin(th)
@ total=2000, meth=disc,trans=100
@ xp=x,yp=y,xlo=-2,xhi=2,ylo=-2,yhi=2
done