dim as double theta, phi, xa1, xa2, xa3, ya1, ya2, ya3
dim as double x1, y1, z1, x2, y2, z2, x3, y3, z3
dim as double px1, px2, px3, px4, px5, px6, px7, px8
dim as double py1, py2, py3, py4, py5, py6, py7, py8
dim as double pz1, pz2, pz3, pz4, pz5, pz6, pz7, pz8

dim as double ax1, ax2, ax3, ax4, ax5, ax6, ax7, ax8
dim as double ay1, ay2, ay3, ay4, ay5, ay6, ay7, ay8


screen 18,, 2

theta = 0

phi = 0

d = 4
z = 1

n = 0




do


IF n = 0 THEN
SCREENset 1, 2
ELSEIF n = 1 THEN
SCREENset 2, 1

END IF



WINDOW (-z, z)-(z, -z)


CLS

REM First Rotations (in x-y plane)

REM Unrotated, point1 = (1,0,0)

xa1 = COS(theta)
ya1 = SIN(theta)

REM Unrotated, point3 = (001)

xa3 = 0
ya3 = 0

REM Unrotated, point2 = (010)

xa2 = -SIN(theta)
ya2 = COS(theta)

REM Second rotation, in x'-z plane

x1 = xa1 * COS(phi)
y1 = ya1
z1 = xa1 * SIN(phi)

x3 = -SIN(phi)
y3 = ya3
z3 = COS(phi)

x2 = xa2 * COS(phi)
y2 = ya2
z2 = xa2 * SIN(phi)

REM Now that the unit vectors have been calculated, we calculate
REM the position of the cube's corners


px1 = -x1 - x2 + x3 'This point was originally (-1, -1, +1)
py1 = -y1 - y2 + y3
pz1 = -z1 - z2 + z3

px2 = -x1 + x2 + x3 'This point was originally (-1, +1, +1)
py2 = -y1 + y2 + y3
pz2 = -z1 + z2 + z3

px3 = x1 + x2 + x3
py3 = y1 + y2 + y3
pz3 = z1 + z2 + z3

px4 = -x1 - x2 - x3
py4 = -y1 - y2 - y3
pz4 = -z1 - z2 - z3

px5 = -x1 + x2 - x3
py5 = -y1 + y2 - y3
pz5 = -z1 + z2 - z3

px6 = x1 + x2 - x3
py6 = y1 + y2 - y3
pz6 = z1 + z2 - z3

px7 = x1 - x2 - x3
py7 = y1 - y2 - y3
pz7 = z1 - z2 - z3

px8 = x1 - x2 + x3
py8 = y1 - y2 + y3
pz8 = z1 - z2 + z3


REM Now, using perspective, we reduce the above 3D points to points in
REM the plane

ax1 = py1 / (-px1 + d)
ay1 = pz1 / (-px1 + d)

ax2 = py2 / (-px2 + d)
ay2 = pz2 / (-px2 + d)


ax3 = py3 / (-px3 + d)
ay3 = pz3 / (-px3 + d)

ax4 = py4 / (-px4 + d)
ay4 = pz4 / (-px4 + d)

ax5 = py5 / (-px5 + d)
ay5 = pz5 / (-px5 + d)

ax6 = py6 / (-px6 + d)
ay6 = pz6 / (-px6 + d)

ax7 = py7 / (-px7 + d)
ay7 = pz7 / (-px7 + d)

ax8 = py8 / (-px8 + d)
ay8 = pz8 / (-px8 + d)


REM Now, we connect the points

COLOR 10

LINE (ax1, ay1)-(ax2, ay2) 'Going clockwise around original
LINE (ax2, ay2)-(ax5, ay5) 'back face
LINE (ax5, ay5)-(ax4, ay4)
LINE (ax4, ay4)-(ax1, ay1)

COLOR 11

LINE (ax1, ay1)-(ax8, ay8) 'Connecting front to back face
LINE (ax2, ay2)-(ax3, ay3)
LINE (ax5, ay5)-(ax6, ay6)
LINE (ax4, ay4)-(ax7, ay7)

COLOR 12

LINE (ax8, ay8)-(ax3, ay3) 'Going clockwise around original
LINE (ax3, ay3)-(ax6, ay6) 'front face
LINE (ax6, ay6)-(ax7, ay7)
LINE (ax7, ay7)-(ax8, ay8)


theta = theta + .02
phi = phi - .02

IF n = 1 THEN
n = 0
ELSEIF n = 0 THEN
n = 1
END IF

sleep 12




loop while inkey$=""




REM Following is derivation of rotation coordinate transformation

'X = r*cos(t1 + t2)
'Y = r*sin(t1 + t2)

'x = r*cos(t1)
'y = r*sin(t1)

'X = r*(cos(t1)*cos(t2) - sin(t1)*sin(t2))
'Y = r*(sin(t1)*cos(t2) + cos(t1)*sin(t2))

'X = x*cos(t2) - y*sin(t2)
'Y = y*cos(t2) + x*sin(t2)