Skip to content

Instantly share code, notes, and snippets.

@jgamble
Last active October 23, 2018 22:44
Show Gist options
  • Select an option

  • Save jgamble/b5d3d504e085b8f6e1b6a879a894e0a0 to your computer and use it in GitHub Desktop.

Select an option

Save jgamble/b5d3d504e085b8f6e1b6a879a894e0a0 to your computer and use it in GitHub Desktop.
Spirograph simulation -- rediscovered the code from an August 1983 printout. There are two routines, SPIRO and XPIRO
C
C Subroutine spiro simulates the designs of the toy spirograph.
C Where:
C R1 -- The radius of the stationary wheel.
C R2 -- The radius of the moving 'pen' wheel.
C P -- The pen's distance from the center of the 'pen' wheel.
C ANGLE -- The angle by which the drawn figure is rotated, in degrees.
C STEP -- The increment of rotation of the pen wheel around the
C stationary wheel, in degrees.
C X, Y -- Arrays to hold the X and Y coordinates of the design.
C N -- The dimension of the arrays X and Y. it should equal
C 360 * (the number of rotations of the 'pen' wheel) / STEP.
C
C If R1 and R2 are both positive, the 'pen' wheel is on the outer edge
C of the stationary wheel.
C If R1 is positive and R2 negative, the 'pen' wheel is on the inner
C edge of the stationary wheel.
C If R1 is negative and R2 positive, the 'pen' wheel encircles the
C stationary wheel.
C
SUBROUTINE SPIRO(R1, R2, P, ANGLE, STEP, X, Y, N)
INTEGER I, N
REAL X(N), Y(N), R1, R2, R, P, D, LEG
REAL ROTATE, ANGLE, THETA1, THETA2, ALPHA, GAMMA, TOOTH, STEP, DR
C
THETA1 = 0.0
THETA2 = 0.0
ALPHA = 0.0
DR = ATAN(1.0)/45.0
TOOTH = STEP * DR
ROTATE = ANGLE * DR
R = R1 + R2
C
IF (ABS(P) .GT. ABS(R)) THEN
LEG = SQRT(P*P - R*R)
ELSE
LEG = -1.0
ENDIF
C
DO 10 I = 1, N
D = SQRT(R*R + P*P - 2*P*R*COS(THETA1))
C
IF (D .NE. 0.0) THEN
GAMMA = ASIN(SIN(THETA2)*P/D)
IF (D .LT. LEG) GAMMA = 180.0 * DR - GAMMA
ALPHA = THETA1 - GAMMA * ROTATE
ENDIF
C
X(I) = D * COS(ALPHA)
Y(I) = D * SIN(ALPHA)
THETA1 = FLOAT(I) * TOOTH
THETA2 = THETA1 * R1/R2
10 CONTINUE
RETURN
END
C
C Subroutine xpiro is an extended verion of subroutine spiro.
C Where:
C F1 -- The radius of the stationary wheel.
C F2 -- The radius of the moving 'pen' wheel.
C SPIN -- The function that determines the rate of spin of the
C function F2. It is found by comparing the line
C length functions of F1 and F2. The line length
C function of a polar equation is:
C
C LENGTH(ANGLE) = INTEGRAL( SQRT(F(ANGLE)**2 + F'(ANGLE)**2)
C
C Then if ALEN is the inverse line length function,
C SPIN = ALEN2(LENGTH1(ANGLE))
C
C Obviously, if F1 = F2, then SPIN(ANGLE) = ANGLE.
C
C P -- The pen's distance from the center of function F2.
C ANGLE -- The angle by which the drawn figure is rotated, in degrees.
C STEP -- The increment of rotation of F2 around the stationary
C wheel, in degrees.
C X, Y -- Arrays to hold the X and Y coordinates of the design.
C N -- The dimension of the arrays X and Y. it should equal
C 360 * (the number of rotations of the 'pen' wheel) / STEP.
C SIDE -- A REAL that determines whether F2 is moving on the
C inside (SIDE = -1.0) or the outside (side = 1.0) of F1.
C
C
SUBROUTINE XPIRO(F1, F2, SPIN, P, ANGLE, STEP, X, Y, N)
INTEGER I, N
REAL X(N), Y(N), R, P, D, LEG, SIDE
REAL ROTATE, ANGLE, THETA1, THETA2, ALPHA, GAMMA, TOOTH, STEP, DR
REAL F1, F2, SPIN
C
THETA1 = 0.0
THETA2 = 0.0
ALPHA = 0.0
DR = ATAN(1.0)/45.0
TOOTH = STEP * DR
ROTATE = ANGLE * DR
C
DO 10 I = 1, N
R = F1(THETA1) + F2(THETA2) * SIDE
D = SQRT(R*R + P*P - 2*P*R*COS(THETA1))
C
IF (ABS(P) .GT. ABS(R)) THEN
LEG = SQRT(P*P - R*R)
ELSE
LEG = 0.0
ENDIF
IF (D .NE. 0.0) THEN
GAMMA = ASIN(SIN(THETA2)*P/D)
IF (D .LT. LEG) GAMMA = 180.0 * DR - GAMMA
ALPHA = THETA1 - GAMMA * ROTATE
ENDIF
C
X(I) = D * COS(ALPHA)
Y(I) = D * SIN(ALPHA)
THETA1 = FLOAT(I) * TOOTH
THETA2 = SPIN(THETA1) * SIDE
10 CONTINUE
RETURN
END
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment