'============================================================================
'  G.H. George                                                 1993 SEP 02
'                                ATTRACT.BAS         Modified  1998 APR 30
'  Dynamic illustration of mathematical attractors.
'  From the Mathematical Gazette, v. 77, # 479, 1993 July, Note 77.10,
'  pp. 247-250, "Extraordinary pictures of some simple sums", N. MacKinnon.
'============================================================================

'   Procedures:
'   GetValue:    prints welcome message; gets the value of the exponent  d .
'   DrawPicture: plots the points.
DECLARE SUB GetValue (d!)
DECLARE SUB DrawPicture (d!)
DECLARE SUB DrawAxes (zoom!)

'   Global constants:
'    False:  logical false (= 0)
'     True:  logical true (= -1)
'       Pi:  pi
CONST False = 0, True = NOT False, Pi = 3.141593

'   Variables:
'        d:  the exponent in the summations.

CALL GetValue(d)
CALL DrawPicture(d)

END

SUB DrawAxes (zoom)
'---------------------------------------------------------------------------
'   Ask the user for the zoom factor.  Draw axes and guide circle.
'---------------------------------------------------------------------------

'   Procedures:
'       <none>

'   Parameter:
'     zoom:  the zoom factor for the plot.

'   Local variables:
'        i:  loop counter = step of five degrees in angle
'    theta:  angle in radians on the guide circle, inside [0, 2Pi).
'        x:  x coordinate of the current point on the guide circle.
'        y:  y coordinate of the current point on the guide circle.

LOCATE 5, 45
INPUT "Zoom factor: ", zoom             ' Ask the user for the zoom factor.
                                        '   (any accepted!)
LOCATE 2, 53
PRINT USING "(Guide circle radius = ####)"; 50 / zoom

LINE (320, 48)-(320, 199), 3            ' Draw y axis (cyan colour)
LINE (0, 120)-(639, 120), 3             ' Draw x axis

PSET (320 + 2.4 * 50, 120), 3           ' Plot first point on guide circle.

FOR i = 1 TO 72                         ' Draw guide circle.
    LET theta = i * Pi / 36
    LET x = 50 * COS(theta)             ' radius of guide circle = 50 pixels
    LET y = 50 * SIN(theta)
    LINE -(320 + 2.4 * x, 120 - y), 3      ' Aspect ratio = 12/5 in screen 8
NEXT i

END SUB

SUB DrawPicture (d)
'---------------------------------------------------------------------------
'   Plots the points  (x_n, y_n).   The user controls the zoom factor.
'---------------------------------------------------------------------------

'   Procedures:
'       DrawAxes:   ask the user for a zoom factor; draw axes & guide circle.

'   Local variables:
'        x:  current value of the partial sum of the  x_n .
'        y:  current value of the partial sum of the  y_n .
'     zoom:  zoom factor (user's choice)
'        i:  loop counter = current point to plot
'    theta:  current value of i^d = argument for sin and cos functions
'   Colour:  the current colour in which to plot the curve.

'   Parameter:
'        d:  the exponent in the summations.

DIM i AS LONG, theta AS DOUBLE, x AS DOUBLE, y AS DOUBLE, Colour AS INTEGER

CALL DrawAxes(zoom)

' Plot all points.

LET Colour = 1

FOR i = 1 TO 25000
    LET theta = i ^ d
    LET theta = theta - 2 * Pi * INT(theta / (2 * Pi))
                                                ' Keep angle inside [0, 2Pi)
    LET x = x + COS(theta)
    LET y = y + SIN(theta)
    PSET (2.4 * zoom * x + 320, 120 - zoom * y), Colour
                        '^-- aspect ratio = 12/5 

    IF i MOD 1000 = 0 THEN                      ' Update printed value of  i
        LOCATE 6, 65                            '  every 1000th time
        PRINT USING "i = #####"; i
        LET Colour = (i / 1000 MOD 15) + 1      ' Change colour every 1000
    END IF                                      '  points.

NEXT i

END SUB

SUB GetValue (d)
'---------------------------------------------------------------------------
'   Prints welcome message; gets the value of the exponent  d .
'---------------------------------------------------------------------------

'   No procedures or local variables.

'   Parameter:
'        d:  the exponent in the summations.

SCREEN 8
CLS
PRINT "This program plots the points  (x_n, y_n)  where"
PRINT "             n                                n"
PRINT "    x_n  =  Sum  cos(i^d)     and    y_n  =  Sum  sin(i^d)"
PRINT "            i=1                              i=1"
PRINT
INPUT "Enter a value for the parameter  d: ", d

' Data validation:  [deactivated here; *any* value of  d  is accepted]

'DO WHILE d < 0
'    PRINT "Error:  d  must be positive.   Try again."
'    INPUT "Enter a value for the parameter  d: ", d
'LOOP

END SUB