' GALAXY - Ported to METAL by Allan Crossman
' Major additions by Lisa Lovelace
'
' A galaxy collision program adapted from the listing
' in Astronomy magazine, Dec 1988, by M. C. Schroeder
' Neil F. Comins. This program allows two galaxies
' to interact via their gravitational fields. The 
' nuclei of the TARGET galaxy and the INTRUDER galaxy
' are point masses, and they move according to their
' gravitational attraction. The stars respond to the
' gravity of the nuclei but do not influence their
' motions.
'
' With the help of Lisa, this version allows you to:
'   Set the galaxies as slanting at the start.
'   Rotate around the view at runtime.
'
' By the way, the speed is dependent on how many colours
' your monitor is set to - try using only 256 to make
' this program fly...

'' You can change these starting values and options:

NRR = 5              ' Rings of stars in each galaxy.
NRS = 25             ' Stars per ring in each galaxy.

M2 = 1               ' Mass of intruder (as proportion of target's mass)

X2 = 60              ' X position of intruder (relative to target)
Y2 = 160             ' Y position of intruder (relative to target)
Z2 = 20              ' Z position of intruder (relative to target)

VX2 = -0.2           ' X velocity of intruder (relative to target)
VY2 = -0.15          ' Y velocity of intruder (relative to target)
VZ2 = -0.12          ' Z velocity of intruder (relative to target)

theta1 = 0           ' Slant angle of target plane along Y axis
omega1 = 0           ' Rotate (after slant) along Z axis for target

theta2 = 60          ' Slant angle of intruder plane along Y axis
omega2 = 20          ' Rotate (after slant) along Z axis for intruder

innerR = 6           ' Radius of innermost ring of stars
outerR = 35          ' Radius of outermost ring of stars
                     ' (note: if only one ring, then outerR is used)

same_rotation = 0    ' Set to 1 to make both galaxies rotate clockwise
                     ' Set to 0 to make one clockwise, one anti-clockwise

intruder_stars = 1   ' Does the intruder object get stars as well?

zoom = 1             ' The initial zoom (can be adjusted at runtime)
zoomvalue = 0.02     ' How much a single zoom in/out moves you

' Set METAL console...

width = 800          ' You can also safely change these.
depth = 400

left = screen width / 2 - (width / 2)
top = screen height / 2 - (depth / 2)
right = screen width / 2 + (width / 2)
bottom = screen height / 2 + (depth / 2)

resize console left, top, right, bottom
set console title to "    Q : Quit     P : Pause     I / T : Centre view     + / - : Zoom     <- / -> : Rotate    "
cls

if intruder_stars then
   NS = NRR * NRS * 2
else
   NS = NRR * NRS
end if

if NRR > 1 then DR = (outerR - innerR) / (NRR - 1)

dim X(NS)
dim Y(NS)
dim Z(NS)
dim VX(NS)
dim VY(NS)
dim VZ(NS)

' The position of the TARGET galaxy
' is assigned by the computer.

' Initialize TARGET galaxy mass,
' position, and velocity.

M1 = 5
X1 = 0
Y1 = 0
Z1 = 0
VX1 = 0
VY1 = 0
VZ1 = 0
SF2 = 2     ' See note below.

' Scale the INTRUDER galaxy mass
' position and velocities.

M2 = M2 * M1
X2 = X2 + X1
Y2 = Y2 + Y1
Z2 = Z2 + Z1

' Initialize Lisa's rotation stuff...

dim TM(3,3)                   ' transfer matrix for view rotation
dim TM2(6,3,3)                ' 6 series of matrixes for additional rotation. (X-plus, X-minus, Y-.. Z...)
dim wkM(3,3)                  ' work area for matrix calc.
cosVphi = cos(1/180*3.14159)  ' unit angle for 1 step rotation.  (used for all (x,y,z) axis)  
sinVphi = sin(1/180*3.14159)
gosub initTMs                 ' Initialize matrixes

' Set up initial load,
' place stars in concentric rings.

I = 0

FOR IR = 1 TO NRR
   if NRR > 1 then R = innerR + (IR - 1) * DR else R = outerR
   V = sqr(M1 / R)
   TH = (.5 * V / R) * (180 / 3.14159)
   IF R = 10 THEN V = .9 * V
   FOR IT = 1 TO NRS
      T = (IT - 1) * 360 / NRS
      T1 = 3.14159 * (T - TH) / 180
      I = I + 1
 
' Initialize star positions (relative coordinates to Galaxy center at first).
 
      X(I) = R * cos(T / 57.2958)
      Y(I) = R * sin(T / 57.2958)
      Z(I) = 0
      
' Slant / rotate along Y axis (theta deg).
 
      workZ = cos(theta1 / 57.2958) * Z(I) - sin(theta1 / 57.2958) * X(I)
      workX = sin(theta1 / 57.2958) * Z(I) + cos(theta1 / 57.2958) * X(I)
      Z(I) = workZ
      X(I) = workX
 
' Then rotate around Z axis.

      workX = cos(omega1 / 57.2958) * X(I) - sin(omega1 / 57.2958) * Y(I)
      workY = sin(omega1 / 57.2958) * X(I) + cos(omega1 / 57.2958) * Y(I)
      X(I) = workX
      Y(I) = workY
 
' Then add center coordinate of galaxy.
      
      X(I) = X(I) + X1
      Y(I) = Y(I) + Y1
      Z(I) = Z(I) + Z1
 
' Initialize star velocities
' to place them in circular
' orbits about the target galaxy.
 
      VX(I) = -V * sin(T1)
      VY(I) = V * cos(T1)
      VZ(I) = 0
 
' Slant / rotate along Y axis (theta deg)
 
      workZ = cos(theta1 / 57.2958) * VZ(I) - sin(theta1 / 57.2958) * VX(I)
      workX = sin(theta1 / 57.2958) * VZ(I) + cos(theta1 / 57.2958) * VX(I)
      VZ(I) = workZ
      VX(I) = workX
 
' Then rotate around Z axis
 
      workX = cos(omega1 / 57.2958) * VX(I) - sin(omega1 / 57.2958) * VY(I)
      workY = sin(omega1 / 57.2958) * VX(I) + cos(omega1 / 57.2958) * VY(I)
      VX(I) = workX
      VY(I) = workY
 
' Then add velocity of galaxy core.
 
      VX(I) = VX(I) + VX1
      VY(I) = VY(I) + VY1
      VZ(I) = VZ(I) + VZ1
      
   NEXT IT
NEXT IR
 
' And do the same thing for the intruder galaxy,
' if the flag intruder_stars is set (above).
 
if intruder_stars then

   FOR IR = 1 TO NRR
      if NRR > 1 then R = innerR + (IR - 1) * DR else R = outerR
      V = sqr(M2 / R)
      TH = (.5 * V / R) * (180 / 3.14159)
      IF R = 10 THEN V = .9 * V
      FOR IT = 1 TO NRS
         T = (IT - 1) * 360 / NRS
         T1 = 3.14159 * (T - TH) / 180
         I = I + 1
 
' Initialize star positions (relative coordinates to Galaxy center at first).
 
         X(I) = R * cos(T / 57.2958)
         Y(I) = R * sin(T / 57.2958)
         Z(I) = 0
      
' Slant / rotate along Y axis (theta deg).
 
         workZ = cos(theta2 / 57.2958) * Z(I) - sin(theta2 / 57.2958) * X(I)
         workX = sin(theta2 / 57.2958) * Z(I) + cos(theta2 / 57.2958) * X(I)
         Z(I) = workZ
         X(I) = workX
 
' Then rotate around Z axis.
 
         workX = cos(omega2 / 57.2958) * X(I) - sin(omega2 / 57.2958) * Y(I)
         workY = sin(omega2 / 57.2958) * X(I) + cos(omega2 / 57.2958) * Y(I)
         X(I) = workX
         Y(I) = workY
 
' Then add center coordinate of galaxy.
      
         X(I) = X(I) + X2
         Y(I) = Y(I) + Y2
         Z(I) = Z(I) + Z2
 
' Initialize star velocities
' to place them in circular
' orbits about the intruder galaxy.
 
         if same_rotation then
            VX(I) = -V * sin(T1)
            VY(I) = V * cos(T1)
         else
            VX(I) = V * sin(T1)
            VY(I) = -V * cos(T1)
         end if
         VZ(I) = 0
 
' Slant / rotate along Y axis (theta deg).
 
         workZ = cos(theta2 / 57.2958) * VZ(I) - sin(theta2 / 57.2958) * VX(I)
         workX = sin(theta2 / 57.2958) * VZ(I) + cos(theta2 / 57.2958) * VX(I)
         VZ(I) = workZ
         VX(I) = workX
 
' Then rotate around Z axis.
 
         workX = cos(omega2 / 57.2958) * VX(I) - sin(omega2 / 57.2958) * VY(I)
         workY = sin(omega2 / 57.2958) * VX(I) + cos(omega2 / 57.2958) * VY(I)
         VX(I) = workX
         VY(I) = workY
 
' Then add velocity of galaxy core.
 
         VX(I) = VX(I) + VX2
         VY(I) = VY(I) + VY2
         VZ(I) = VZ(I) + VZ2
         
      NEXT IT
   NEXT IR
 
end if

polar = init screen(0, 0, width / 2, depth)
set screen to polar
cls

edgeon = init screen(0, 0, width / 2, depth)
set screen to edgeon
cls

GOSUB UpdateScreen



Main:

' Press P to pause...

keymap scan
if keymap key ("p") then goto Main

FOR I = 1 TO NS

' Determine the force on a star
' due to the galactic centers.
' SF2 below, called the softening factor,
' is used to prevent overflows
' during force calculations.
' SF2 is assigned a value above.
' You may experiment with its value
' but it should be kept as small as possible
' to better approximate
' a true 1/r**2 force.

   R1 = M1 / ((X(I) - X1) ^ 2 + (Y(I) - Y1) ^ 2 + (Z(I) - Z1) ^ 2 + SF2) ^ 1.5
   R2 = M2 / ((X(I) - X2) ^ 2 + (Y(I) - Y2) ^ 2 + (Z(I) - Z2) ^ 2 + SF2) ^ 1.5

' Calculate star's x,y,z accelerations.

   AX = R1 * (X1 - X(I)) + R2 * (X2 - X(I))
   AY = R1 * (Y1 - Y(I)) + R2 * (Y2 - Y(I))
   AZ = R1 * (Z1 - Z(I)) + R2 * (Z2 - Z(I))
			
' Update star velocities and positions
' using a time-centered
' leap-frog algorithm.

   VX(I) = VX(I) + AX
   VY(I) = VY(I) + AY
   VZ(I) = VZ(I) + AZ

   X(I) = X(I) + VX(I)
   Y(I) = Y(I) + VY(I)
   Z(I) = Z(I) + VZ(I)

NEXT I

' Update positions
' and velocities of the galactic centers.

S = ((X1 - X2) ^ 2 + (Y1 - Y2) ^ 2 + (Z1 - Z2) ^ 2 + SF2) ^ 1.5
AX = (X2 - X1) / S
AY = (Y2 - Y1) / S
AZ = (Z2 - Z1) / S
VX1 = VX1 + M2 * AX
VY1 = VY1 + M2 * AY
VZ1 = VZ1 + M2 * AZ
VX2 = VX2 - M1 * AX
VY2 = VY2 - M1 * AY
VZ2 = VZ2 - M1 * AZ
X1 = X1 + VX1
Y1 = Y1 + VY1
Z1 = Z1 + VZ1
X2 = X2 + VX2
Y2 = Y2 + VY2
Z2 = Z2 + VZ2

GOSUB UpdateScreen
IF keymap key ("q") or keymap key (chr$(27)) then
   kill screen virtue
   disable done
   end
end if
goto Main



UpdateScreen:

set screen to polar
cls
set screen to edgeon
cls

' Calculate center of mass of system
' for use as center of output display.

' User can press T or I to centre
' on target or intruder.

if keymap key ("t") and keymap key ("i") = 0 then
   XC = X1
   YC = Y1
   ZC = Z1
else
   if keymap key ("i") and keymap key ("t") = 0 then
      XC = X2
      YC = Y2
      ZC = Z2
   else
      XC = (M1 * X1 + M2 * X2) / (M1 + M2)
      YC = (M1 * Y1 + M2 * Y2) / (M1 + M2)
      ZC = (M1 * Z1 + M2 * Z2) / (M1 + M2)
   end if
end if


' Rotate with numeric keypad or arrows...

if keymap key (chr$(29)) or keymap key ("6") then RotSel = 6 : gosub doAdditionalRotation  ' Right arrow/ Z-
if keymap key (chr$(28)) or keymap key ("4") then RotSel = 5 : gosub doAdditionalRotation  ' Left arrow / Z+
if keymap key (chr$(30)) or keymap key ("8") then RotSel = 2 : gosub doAdditionalRotation  ' Up arrow   / X-
if keymap key (chr$(31)) or keymap key ("5") then RotSel = 1 : gosub doAdditionalRotation  ' Down arrow / X+

if keymap key ("7") then RotSel = 3 : gosub doAdditionalRotation ' Y+
if keymap key ("9") then RotSel = 4 : gosub doAdditionalRotation ' Y-

if keymap key (" ") then gosub initTM  ' Reset view angle


' Calculate position of galactic centers
' and display on screen.

if keymap key ("-") and zoom > zoomvalue then zoom = zoom - zoomvalue
if keymap key ("+") then zoom = zoom + zoomvalue

XX1 = (X1 - XC) * zoom
YY1 = (Y1 - YC) * zoom
ZZ1 = (Z1 - ZC) * zoom
XX2 = (X2 - XC) * zoom
YY2 = (Y2 - YC) * zoom
ZZ2 = (Z2 - ZC) * zoom


' View Rotation

workX = TM(1,1) * XX1 + TM(1,2) * YY1 + TM(1,3) * ZZ1
workY = TM(2,1) * XX1 + TM(2,2) * YY1 + TM(2,3) * ZZ1
ZZ1 = TM(3,1) * XX1 + TM(3,2) * YY1 + TM(3,3) * ZZ1
XX1 = workX
YY1 = workY

workX = TM(1,1) * XX2 + TM(1,2) * YY2 + TM(1,3) * ZZ2
workY = TM(2,1) * XX2 + TM(2,2) * YY2 + TM(2,3) * ZZ2
ZZ2 = TM(3,1) * XX2 + TM(3,2) * YY2 + TM(3,3) * ZZ2
XX2 = workX
YY2 = workY
 
set screen to polar
forecolor 0, 65535, 0
FCIRCLE (XX1 + (width / 4)), (YY1 + (depth / 2)), 2
forecolor 65535, 65535, 0
FCIRCLE (XX2 + (width / 4)), (YY2 + (depth / 2)), 2
forecolor 65535, 65535, 65535
 
set screen to edgeon
forecolor 0, 65535, 0
FCIRCLE (XX1 + (width / 4)), (ZZ1 + (depth / 2)), 2
forecolor 65535, 65535, 0
FCIRCLE (XX2 + (width / 4)), (ZZ2 + (depth / 2)), 2
forecolor 65535, 65535, 65535
 
 
' Place stars on screen.
 

FOR I = 1 TO NS

   XP = (X(I) - XC) * zoom
   YP = (Y(I) - YC) * zoom
   ZP = (Z(I) - ZC) * zoom
   
   workX = TM(1,1) * XP + TM(1,2) * YP + TM(1,3) * ZP
   workY = TM(2,1) * XP + TM(2,2) * YP + TM(2,3) * ZP
   ZP = TM(3,1) * XP + TM(3,2) * YP + TM(3,3) * ZP
   XP = workX
   YP = workY

   set screen to polar
   plot (XP + (width / 4)), (YP + (depth / 2))
   set screen to edgeon
   plot (XP + (width / 4)), (ZP + (depth / 2))

NEXT I
 
' Print out display headings.
 
set screen to polar
forecolor 65535, 65535, 0
line (width / 2) - 1, 0, (width / 2) - 1, depth
forecolor 65535, 30000, 0
text 10, 20, "Polar view"
forecolor 0, 65535, 0
text 10, depth - 13, "Zoom: " + str$(zoom)
set screen to edgeon
forecolor 65535, 30000, 0
text 10, 20, "Edge-on view"
forecolor 0, 65535, 0
 
' Copy virtual screen to console...
 
copyrect 0, 0, width / 2, depth, 0, 0, width / 2, depth, 0, polar, 0
copyrect 0, 0, width / 2, depth, width / 2, 0, width, depth, 0, edgeon, 0
 
RETURN



''''''''''''''''''''''''''''''
' Lisa's funky matrix stuff...

initTM:  ' initialize transfer matrix. This is also called when user press space bar. 
   for ti=1 to 3 : for tj=1 to 3
      if ti=tj then TM(ti,tj) = 1 else TM(ti,tj) = 0
   next : next
return


initTMs: ' Provide rotation matrixes ( + call 'initTM')   This is called once.

   gosub initTM

   ' Make unit matrix at first.
   for ti=1 to 3 : for tj=1 to 3
      for k = 1 to 6
         if ti=tj then TM2(k,ti,tj) = 1 else TM2(k,ti,tj) = 0
      next
   next : next

   ' Then arrange them for each rotation.
   ' X plus / TM2(1, ...
      TM2(1, 2,2) =  cosVphi : TM2(1, 2,3) = -sinVphi
      TM2(1, 3,2) =  sinVphi : TM2(1, 3,3) =  cosVphi
   ' X minus / TM2(2, ...
      TM2(2, 2,2) =  cosVphi : TM2(2, 2,3) =  sinVphi
      TM2(2, 3,2) = -sinVphi : TM2(2, 3,3) =  cosVphi
   ' Y plus
      TM2(3, 1,1) =  cosVphi : TM2(3, 1,3) =  sinVphi
      TM2(3, 3,1) = -sinVphi : TM2(3, 3,3) =  cosVphi
   ' Y minus
      TM2(4, 1,1) =  cosVphi : TM2(4, 1,3) = -sinVphi
      TM2(4, 3,1) =  sinVphi : TM2(4, 3,3) =  cosVphi
   ' Z plus
      TM2(5, 1,1) =  cosVphi : TM2(5, 1,2) = -sinVphi
      TM2(5, 2,1) =  sinVphi : TM2(5, 2,2) =  cosVphi
   ' Z minus
      TM2(6, 1,1) =  cosVphi : TM2(6, 1,2) =  sinVphi
      TM2(6, 2,1) = -sinVphi : TM2(6, 2,2) =  cosVphi

return


doAdditionalRotation: 
 ' Multiply TM by TM2.
 ' TM2 is placed at left. (TM2 * TM)
 ' loop counter is not used.
 ' needed matrix is selected by 'RotSel'. (1..X plus, 2..X minus ...)

   wkM(1,1) = TM2(RotSel, 1,1) * TM(1,1) + TM2(RotSel, 1,2) * TM(2,1) + TM2(RotSel, 1,3) * TM(3,1)
   wkM(1,2) = TM2(RotSel, 1,1) * TM(1,2) + TM2(RotSel, 1,2) * TM(2,2) + TM2(RotSel, 1,3) * TM(3,2)
   wkM(1,3) = TM2(RotSel, 1,1) * TM(1,3) + TM2(RotSel, 1,2) * TM(2,3) + TM2(RotSel, 1,3) * TM(3,3)

   wkM(2,1) = TM2(RotSel, 2,1) * TM(1,1) + TM2(RotSel, 2,2) * TM(2,1) + TM2(RotSel, 2,3) * TM(3,1)
   wkM(2,2) = TM2(RotSel, 2,1) * TM(1,2) + TM2(RotSel, 2,2) * TM(2,2) + TM2(RotSel, 2,3) * TM(3,2)
   wkM(2,3) = TM2(RotSel, 2,1) * TM(1,3) + TM2(RotSel, 2,2) * TM(2,3) + TM2(RotSel, 2,3) * TM(3,3)

   wkM(3,1) = TM2(RotSel, 3,1) * TM(1,1) + TM2(RotSel, 3,2) * TM(2,1) + TM2(RotSel, 3,3) * TM(3,1)
   wkM(3,2) = TM2(RotSel, 3,1) * TM(1,2) + TM2(RotSel, 3,2) * TM(2,2) + TM2(RotSel, 3,3) * TM(3,2)
   wkM(3,3) = TM2(RotSel, 3,1) * TM(1,3) + TM2(RotSel, 3,2) * TM(2,3) + TM2(RotSel, 3,3) * TM(3,3)

   TM(1,1) = wkM(1,1) : TM(1,2) = wkM(1,2) : TM(1,3) = wkM(1,3)
   TM(2,1) = wkM(2,1) : TM(2,2) = wkM(2,2) : TM(2,3) = wkM(2,3)
   TM(3,1) = wkM(3,1) : TM(3,2) = wkM(3,2) : TM(3,3) = wkM(3,3)

return