PROGRAM Newtons_Basin   !© W.van Duyn 29th April 2004
!Plot the cube roots of 1
!F(z)=z^3-1
!Use Newton's method
!F(z_i+1)=z_i-f(z_i)/f'(z_i)
!z1=1
!z2=-0.5+i*SQR(3)/2
!z3=-0.5-i*SQR(3)/2

SET MODE "COLOR"
SET window 0,1199,0,903
SET background color "white"
INPUT PROMPT "side try 1... ": side
INPUT PROMPT "size try 400... ": f
CLEAR
SET COLOR 16
FLOOD 1,1
LET n=50
LET lim=0.001
LET k=20
LET m=f/(2*side)

FOR x=-f/2 TO f/2
FOR y=-f/2 TO f/2
LET a=x/m   !real part
LET b=y/m   !virtual part
FOR i=1 TO n

LET p=a^2-b^2   !square of a+bi separated into real and virtual parts
LET q=2*a*b   !trivial for cubic but useful for higher powers.
IF p=0 AND q=0 THEN EXIT FOR
    !For higher powers add 1 to the 2* and 3*
LET at=(2*a*p^2+2*a*q^2+p)/(3*(p^2+q^2))
LET b=(2*b*p^2+2*b*q^2-q)/(3*(p^2+q^2))
LET a=at

IF abs(a-1) < lim AND abs(b) < lim THEN
SET COLOR MIX(i+100) 1,1-i/k,i/k
LET col=13
CALL create2(i,x,y,col,k)
EXIT FOR
END IF

IF abs(a+0.5) < lim AND abs(b-SQR(3)/2) < lim THEN
SET COLOR MIX(i+100) 1-i/k,i/k,1
LET col=11
CALL create2(i,x,y,col,k)
EXIT FOR
END IF

IF abs(a+0.5) < lim AND abs(b+SQR(3)/2) < lim THEN
SET COLOR MIX(i+100) i/k,1,1-i/k
LET col=14
CALL create2(i,x,y,col,k)
EXIT FOR
END IF

NEXT i
NEXT y
NEXT x
END

SUB create2(i,x,y,col,k)
IF i<=k THEN SET COLOR i+100
IF i>k THEN SET COLOR col
PLOT x+600,y+450
END SUB