'                                IGMS22.BAS
'
'                INTERGALACTIC GRAVITATIONAL MOTION SIMULATOR
'
'                     v2.2    Last Updated 9 April 1999
'
'                          Programmed by D. Perley
'           http://www.people.cornell.edu/pages/dap29/programs.html
'
'    You are free to perform any modifications to this code that you retain
'    my name above and credit me where appropriate.  The code may be
'    redistributed free of restriction.
'
'    It is likely that you will have memory problems in trying to use
'    this software; the program size is simply too large and the memory
'    allocation requests too great for the QBasic software to handle.
'    If it will not compile when you try to run the program, try changing
'    the maxentry constant down to 50 or lower; this decreases the amount
'    of memory requested but does not impede program function unless you
'    are trying to save a file to disk (this reduces the maximum size of a
'    saved file.)

DECLARE SUB FindPixelStyle ()                 ' finds image for particle
DECLARE SUB AdjustColors ()                   ' changes doppler mimic stats
DECLARE SUB DopplerColor ()                   ' assigns colors to particles
DECLARE SUB Catalog ()                        ' prints all SIM/INT file names
DECLARE SUB ChgSpd ()                         ' alters timing loop
DECLARE SUB Core ()                           ' calculates forces, moves
DECLARE SUB DetermineTimestep ()              ' determines timestep
DECLARE SUB DisplayInfo ()                    ' displays simulation properties
DECLARE SUB DrawParticle ()                   ' draws given particle
DECLARE SUB GraphicCalc ()                    ' calculates graphic coords
DECLARE SUB IFGalaxies ()                     ' setup galaxy property page
DECLARE SUB IFModule ()                       ' setup module
DECLARE SUB IFParticles ()                    ' setup particle property page
DECLARE SUB IFProperties ()                   ' setup simulation property page
DECLARE SUB ImageRotate ()                    ' rotates view of image
DECLARE SUB ImageSlide ()                     ' moves view along a plane
DECLARE SUB ImageZoom ()                      ' changes image size
DECLARE SUB Initialize ()                     ' specifies defaults
DECLARE SUB Initialize2 ()                    ' initializes simulation
DECLARE SUB LimInput ()                       ' special string-input function
DECLARE SUB Pause ()                          ' <space> pause waiting area
DECLARE SUB RetroProject ()                   ' finds position in reverse time
DECLARE SUB Snapshot ()                       ' stores current graphic data
DECLARE SUB StoreInitial ()                   ' stores initial pos, vel
DECLARE SUB StoreSimulation ()                ' records positions on disk
DECLARE SUB TerminateSim ()                   ' empties pos/vel variables
DECLARE SUB TrackParticle ()                  ' moves screen to follow object
DECLARE SUB UserInput ()                      ' mid-program alteration center
                                            
CONST maxparticles = 302                      ' total number of particles
CONST maxentry = 60                           ' maximum entries handleable
CONST G = .1                                  ' gravitational constant
CONST rotateinc = 5 * .0174532                ' increment of rotation

DIM SHARED x(1 TO maxparticles) AS SINGLE     ' x-coordinate
DIM SHARED vx(1 TO maxparticles) AS SINGLE    ' x-axis velocity
DIM SHARED ax(1 TO maxparticles) AS SINGLE    ' x-axis acceleration
DIM SHARED y(1 TO maxparticles) AS SINGLE     ' y-coordinate
DIM SHARED vy(1 TO maxparticles) AS SINGLE    ' y-axis velocity
DIM SHARED ay(1 TO maxparticles) AS SINGLE    ' y-axis acceleration
DIM SHARED z(1 TO maxparticles) AS SINGLE     ' z-coordinate
DIM SHARED vz(1 TO maxparticles) AS SINGLE    ' z-axis velocity
DIM SHARED az(1 TO maxparticles) AS SINGLE    ' z-axis acceleration
DIM SHARED m(1 TO maxparticles) AS SINGLE     ' mass
DIM SHARED c(1 TO maxparticles) AS INTEGER    ' color of pixel
DIM SHARED c2(1 TO maxparticles) AS INTEGER   ' color of tracer
DIM SHARED xg(1 TO maxparticles) AS SINGLE    ' pixel y-coordinate
DIM SHARED yg(1 TO maxparticles) AS SINGLE    ' pixel x-coordinate
DIM SHARED rx(1 TO maxparticles) AS SINGLE    ' virtual x addition
DIM SHARED ry(1 TO maxparticles) AS SINGLE    ' virtual y addition
DIM SHARED rz(1 TO maxparticles) AS SINGLE    ' virtual z addition
DIM SHARED gx(1 TO 2) AS SINGLE               ' galaxy initial x-coord
DIM SHARED gvx(1 TO 2) AS SINGLE              ' galaxy initial x velocity
DIM SHARED gy(1 TO 2) AS SINGLE               ' galaxy initial y-coord
DIM SHARED gvy(1 TO 2) AS SINGLE              ' galaxy initial y velocity
DIM SHARED gz(1 TO 2) AS SINGLE               ' galaxy initial z-coord
DIM SHARED gvz(1 TO 2) AS SINGLE              ' galaxy initial z velocity
DIM SHARED gm(1 TO 2) AS SINGLE               ' galaxy mass
DIM SHARED gc(1 TO 2) AS INTEGER              ' galaxy color
DIM SHARED rings(1 TO 2) AS INTEGER           ' number of rings in halo
DIM SHARED rmin(1 TO 2) AS INTEGER            ' radius of innermost ring
DIM SHARED rmax(1 TO 2) AS INTEGER            ' radius of outermost ring
DIM SHARED particles(1 TO 2) AS INTEGER       ' particles per ring
DIM SHARED tilt(1 TO 2) AS INTEGER            ' tilt of halo (1 axis only)
DIM SHARED particle AS INTEGER                ' particle ID
DIM SHARED particle2 AS INTEGER               ' particle ID of interactor
DIM SHARED time AS SINGLE                     ' time
DIM SHARED timestep AS SINGLE                 ' minimum uniform time step
DIM SHARED exitflag AS INTEGER                ' end condition
DIM SHARED xd AS DOUBLE                       ' x distance between particles
DIM SHARED yd AS DOUBLE                       ' y distance between particles
DIM SHARED zd AS DOUBLE                       ' z distance between particles
DIM SHARED R AS DOUBLE                        ' distance between particles
DIM SHARED inkeyv AS STRING                   ' INKEY input
DIM SHARED speed AS LONG                      ' amount of timing-loop delay
DIM SHARED n AS LONG                          ' misc. counting variable
DIM SHARED inc AS INTEGER                     ' change to speed
DIM SHARED spdchg AS INTEGER                  ' direction to change speed
DIM SHARED np AS INTEGER                      ' new particle ID
DIM SHARED offx AS INTEGER                    ' graphic offset, x-axis
DIM SHARED offy AS INTEGER                    ' graphic offset, y-axis
DIM SHARED offa AS SINGLE                     ' x-y rotation
DIM SHARED offb AS SINGLE                     ' x-z rotation
DIM SHARED offc AS SINGLE                     ' y-z rotation
DIM SHARED tracers AS INTEGER                 ' tracers flag
DIM SHARED pauseflag AS INTEGER               ' pause flag
DIM SHARED pixelcolor AS INTEGER              ' pixel color in DrawParticle
DIM SHARED store AS INTEGER                   ' storage (none, RAM, disk)
DIM SHARED oldx AS SINGLE                     ' swap variable for x-coord
DIM SHARED oldy AS SINGLE                     ' swap variable for y-coord
DIM SHARED oldz AS SINGLE                     ' swap variable for z-coord
DIM SHARED display AS INTEGER                 ' display information flag
DIM SHARED winx1 AS INTEGER                   ' viewing window tlc x-coord
DIM SHARED winy1 AS INTEGER                   ' viewing window tlc y-coord
DIM SHARED winx2 AS INTEGER                   ' viewing window brc x-coord
DIM SHARED winy2 AS INTEGER                   ' viewing window brc y-coord
DIM SHARED dataentry AS LONG                  ' data storage entry #
DIM SHARED acceleration AS SINGLE             ' total acceleration on all axes
DIM SHARED subtime AS SINGLE                  ' fractional time component
DIM SHARED difftimestep AS INTEGER            ' differential timestep flag
DIM SHARED ia AS SINGLE                       ' initial particle angle
DIM SHARED ringparticles AS INTEGER           ' number of particles in ring
DIM SHARED rotation AS INTEGER                ' direction of rotation
DIM SHARED subtimestep AS SINGLE              ' non-differential subtimestep
DIM SHARED shotnum AS INTEGER                 ' snapshot ID
DIM SHARED simname AS STRING * 6              ' simulation name
DIM SHARED filename AS STRING * 12            ' file name
DIM SHARED filenum AS STRING * 2              ' snapshot ID as 2-char string
DIM SHARED autosnap AS INTEGER                ' frequency of auto-snapshot
DIM SHARED recfreq AS INTEGER                 ' frequency of data recording
DIM SHARED source AS INTEGER                  ' 1 if source is file
DIM SHARED fileindex AS INTEGER               ' file index
DIM SHARED rx2 AS SINGLE                      ' virtual x addition (particle2)
DIM SHARED ry2 AS SINGLE                      ' virtual y addition (particle2)
DIM SHARED rz2 AS SINGLE                      ' virtual z addition (particle2)
DIM SHARED xtrans0 AS SINGLE                  ' x axis transformation variable
DIM SHARED xtrans1 AS SINGLE                  ' x axis transformation variable
DIM SHARED xtrans2 AS SINGLE                  ' x axis transformation variable
DIM SHARED xtrans3 AS SINGLE                  ' x axis transformation variable
DIM SHARED ytrans0 AS SINGLE                  ' y axis transformation variable
DIM SHARED ytrans1 AS SINGLE                  ' y axis transformation variable
DIM SHARED ytrans2 AS SINGLE                  ' y axis transformation variable
DIM SHARED ytrans3 AS SINGLE                  ' y axis transformation variable
DIM SHARED ztrans0 AS SINGLE                  ' z axis transformation variable
DIM SHARED ztrans1 AS SINGLE                  ' z axis transformation variable
DIM SHARED ztrans2 AS SINGLE                  ' z axis transformation variable
DIM SHARED ztrans3 AS SINGLE                  ' z axis transformation variable
DIM SHARED optioncode AS INTEGER              ' allows diff. action in 1 sub
DIM SHARED numparticles AS INTEGER            ' total particles
DIM SHARED track AS INTEGER                   ' tracked particle ID
DIM SHARED scancode AS INTEGER                ' code for non-ASCII keypresses
DIM SHARED pnt AS INTEGER                     ' interface option ID
DIM SHARED sel AS INTEGER                     ' interface selection ID
DIM SHARED ylab AS INTEGER                    ' y-coord of option label
DIM SHARED xlab AS INTEGER                    ' x-coord of option label
DIM SHARED ysel AS INTEGER                    ' y-coord of selection label
DIM SHARED xsel AS INTEGER                    ' x-coord of selection label
DIM SHARED label AS STRING                    ' option label
DIM SHARED selection AS STRING                ' selection label
DIM SHARED inputstring AS STRING              ' LimInput argument
DIM SHARED lim AS INTEGER                     ' LimInput argument
DIM SHARED xpos AS INTEGER                    ' LimInput argument
DIM SHARED ypos AS INTEGER                    ' LimInput argument
DIM SHARED holdvar AS INTEGER                 ' miscellaneous variable holder
DIM SHARED motion AS INTEGER                  ' flags option change
DIM SHARED skipinput AS INTEGER               ' skips input for 1 turn
DIM SHARED monitor AS INTEGER                 ' view simulation flag
DIM SHARED gtrack AS INTEGER                  ' galaxy to track
DIM SHARED nextif AS INTEGER                  ' flags move to next IF page
DIM SHARED value AS SINGLE                    ' selection numerical label
DIM SHARED msy AS INTEGER                     ' LimInput argument
DIM SHARED check AS INTEGER                   ' flags variable conflict
DIM SHARED runsim AS INTEGER                  ' flags begin simulation
DIM SHARED oldvalue AS INTEGER                ' holds value
DIM SHARED errcode AS INTEGER                 ' identifies variable conflict
DIM SHARED viewcat AS INTEGER                 ' flags view catalog
DIM SHARED retro AS INTEGER                   ' flags retroproject
DIM SHARED perir AS DOUBLE                    ' R at perigalaction
DIM SHARED galaxy AS INTEGER                  ' galaxy ID
DIM SHARED ringspace AS SINGLE                ' distance between halo rings
DIM SHARED tiltrad AS SINGLE                  ' tilt in rad
DIM SHARED ulm AS INTEGER                     ' highest particle on if page
DIM SHARED page AS INTEGER                    ' page on if particle display
DIM SHARED centx AS SINGLE                    ' acts as effective 0 x-coord
DIM SHARED centy AS SINGLE                    ' acts as effective 0 y-coord
DIM SHARED centz AS SINGLE                    ' acts as effective 0 z-coord
DIM SHARED autoterminate AS INTEGER           ' flags automatic termination
DIM SHARED doppler AS INTEGER                 ' doppler shift flag
DIM SHARED vs AS SINGLE                       ' line of sight velocity
DIM SHARED pixelstyle AS INTEGER              ' type of display for particle
DIM SHARED pixelmode AS INTEGER               ' pixelstyles for all particles
DIM SHARED contrast AS INTEGER                ' doppler mimic contrast
DIM SHARED restvel AS INTEGER                 ' doppler mimic 0 point
DIM SHARED displaytime AS SINGLE              ' displayed time
DIM SHARED IGdistance AS SINGLE               ' perigalaction distance
DIM SHARED IGvelocity AS SINGLE               ' perigalaction velocity
DIM SHARED IGinclination AS SINGLE            ' kinematic inclination (y-axis)
DIM SHARED IGargument AS SINGLE               ' kinematic inclination (x-axis)
DIM SHARED IGt0 AS SINGLE                     ' initial time bf. perigalaction
DIM SHARED subtimestep(1 TO maxparticles) AS SINGLE     'timesteps / timestep
'$DYNAMIC
DIM SHARED sx(1 TO maxparticles, 0 TO maxentry)  AS SINGLE  'stores all x
DIM SHARED sy(1 TO maxparticles, 0 TO maxentry)  AS SINGLE  'stores all y
DIM SHARED sz(1 TO maxparticles, 0 TO maxentry)  AS SINGLE  'stores all z


CALL Initialize
DO
   CALL IFModule
   CALL Initialize2
   DO
      CALL Core
      CALL UserInput
   LOOP UNTIL exitflag = 1
   IF store = 2 THEN CALL StoreSimulation
   CALL TerminateSim
LOOP

Overflow:
   value = 0
RESUME NEXT

FileNotFound:
   check = 0
   errcode = 18
RESUME NEXT

MissingEntry:
RESUME NEXT

ExtNotFound:
RESUME NEXT

REM $STATIC
SUB AdjustColors
   'Changes contrast and rest velocity of doppler mimic colors

   IF (inkeyv = CHR$(78) OR inkeyv = CHR$(110)) AND contrast < 30 THEN contrast = contrast + 1
   IF (inkeyv = CHR$(77) OR inkeyv = CHR$(109)) AND contrast > 2 THEN contrast = contrast - 1
   IF (inkeyv = CHR$(74) OR inkeyv = CHR$(106)) AND restvel > -50 THEN restvel = restvel - 1
   IF (inkeyv = CHR$(75) OR inkeyv = CHR$(107)) AND restvel < 50 THEN restvel = restvel + 1

   IF monitor = 1 THEN
      IF tracers = 1 THEN tracers = 0
      IF tracers = 2 THEN CLS
      pixelcolor = 0
      FOR particle = 1 TO numparticles
         CALL FindPixelStyle
         CALL DrawParticle
      NEXT particle
   

      FOR particle = 1 TO numparticles
         CALL FindPixelStyle
         IF doppler = 1 THEN CALL DopplerColor
         pixelcolor = c(particle)
         CALL GraphicCalc
         CALL DrawParticle
      NEXT particle
   END IF
END SUB

SUB Catalog
   'Lists all loadable files
  
   CLS
   LOCATE 1, 4: PRINT "INTERGALACTIC GRAVITATIONAL MOTION SIMULATOR"
   LOCATE 1, 66: PRINT "Catalog"
   ON ERROR GOTO FileNotFound
   FILES "*.SET"
   FILES "*.SIM"
   DO
   LOOP WHILE INKEY$ = ""
   nextif = 1
   ON ERROR GOTO 0
END SUB

SUB ChgSpd
   'Alters the amount of time spent in the timing loop

   IF spdchg = -1 THEN
                            inc = 32000
      IF speed < 100000 THEN inc = 16000
      IF speed < 40000 THEN inc = 8000
      IF speed < 20000 THEN inc = 4000
      IF speed < 10000 THEN inc = 2000
      IF speed < 5000 THEN inc = 1000
      IF speed < 2500 THEN inc = 500
      IF speed < 1000 THEN inc = 200
      IF speed < 500 THEN inc = 100
      IF speed > 500000 THEN inc = 0
      speed = speed + inc
   END IF

   IF spdchg = 1 THEN
                            inc = 32000
      IF speed < 116000 THEN inc = 16000
      IF speed < 48000 THEN inc = 8000
      IF speed < 24000 THEN inc = 4000
      IF speed < 12000 THEN inc = 2000
      IF speed < 6000 THEN inc = 1000
      IF speed < 3000 THEN inc = 500
      IF speed < 1200 THEN inc = 200
      IF speed < 600 THEN inc = 100
      IF speed >= inc THEN speed = speed - inc
   END IF

   spdchg = 0
END SUB

SUB Core
   'Calculates new coordinates and velocities for all particles

   FOR particle = 1 TO numparticles
      IF monitor = 1 AND retro = 0 THEN
         CALL FindPixelStyle
         IF tracers = 1 THEN pixelcolor = c2(particle)
         IF tracers = 0 THEN pixelcolor = 0
         CALL DrawParticle
      END IF
   
      IF source = 0 OR retro = 1 THEN
         x(particle) = x(particle) + rx(particle)
         y(particle) = y(particle) + ry(particle)
         z(particle) = z(particle) + rz(particle)
         rx(particle) = 0
         ry(particle) = 0
         rz(particle) = 0
      ELSE
         x(particle) = sx(particle, dataentry)
         y(particle) = sy(particle, dataentry)
         z(particle) = sz(particle, dataentry)
         IF dataentry < maxentry THEN
            vx(particle) = (sx(particle, dataentry + 1) - sx(particle, dataentry)) / recfreq
            vy(particle) = (sy(particle, dataentry + 1) - sy(particle, dataentry)) / recfreq
            vz(particle) = (sz(particle, dataentry + 1) - sz(particle, dataentry)) / recfreq
         END IF
      END IF

      IF track <> 0 THEN CALL TrackParticle

      IF monitor = 1 AND retro = 0 THEN
         IF doppler = 1 THEN CALL DopplerColor
         pixelcolor = c(particle)
         CALL FindPixelStyle
         CALL GraphicCalc
         CALL DrawParticle
      END IF
   NEXT particle

   IF source = 0 OR retro = 1 THEN
      FOR particle = 1 TO numparticles
         FOR subtime = subtimestep(particle) TO timestep STEP subtimestep(particle)
            ax(particle) = 0
            ay(particle) = 0
            az(particle) = 0
            FOR particle2 = 1 TO numparticles
               IF particle <> particle2 AND m(particle2) <> 0 THEN
                  xd = x(particle2) - (x(particle) + rx(particle))
                  yd = y(particle2) - (y(particle) + ry(particle))
                  zd = z(particle2) - (z(particle) + rz(particle))
                  R = (xd ^ 2 + yd ^ 2 + zd ^ 2) ^ (1 / 2)
                  ax(particle) = ax(particle) + xd * G * m(particle2) / R ^ 3
                  ay(particle) = ay(particle) + yd * G * m(particle2) / R ^ 3
                  az(particle) = az(particle) + zd * G * m(particle2) / R ^ 3
               END IF
            NEXT particle2
            vx(particle) = vx(particle) + ax(particle) * subtimestep(particle)
            vy(particle) = vy(particle) + ay(particle) * subtimestep(particle)
            vz(particle) = vz(particle) + az(particle) * subtimestep(particle)
            rx(particle) = rx(particle) + vx(particle) * subtimestep(particle)
            ry(particle) = ry(particle) + vy(particle) * subtimestep(particle)
            rz(particle) = rz(particle) + vz(particle) * subtimestep(particle)
         NEXT subtime
        
         IF difftimestep = 1 AND retro = 0 THEN CALL DetermineTimestep
    
         IF store > 0 AND retro = 0 AND dataentry <= maxentry AND time / timestep MOD recfreq = 0 THEN
            sx(particle, dataentry) = x(particle)
            sy(particle, dataentry) = y(particle)
            sz(particle, dataentry) = z(particle)
         END IF
      NEXT particle
   END IF

   IF autosnap > 0 THEN
      IF (time - 1) MOD autosnap = 0 THEN CALL Snapshot
   END IF

   IF retro = 0 THEN
      FOR n = 1 TO speed * 4
      NEXT n
   END IF

   displaytime = time
   IF source = 0 THEN
      time = time + timestep
      IF time / timestep MOD recfreq = 0 THEN dataentry = dataentry + 1
   END IF
   IF source = 1 THEN
      time = time + timestep * recfreq
      dataentry = dataentry + 1
   END IF

   IF dataentry > maxentry AND source = 0 AND store = 2 AND retro = 0 THEN exitflag = 1: autoterminate = 1
   IF dataentry > maxentry AND source = 1 AND retro = 0 THEN exitflag = 1: autoterminate = 1
END SUB

SUB DetermineTimestep
   'Calculates particle's subtimestep based on acceleration

   acceleration = ((ax(particle) ^ 2 + ay(particle) ^ 2) ^ (1 / 2) * 8) / timestep
                               subtimestep(particle) = timestep
   IF acceleration > .0005 THEN subtimestep(particle) = 4
   IF acceleration > .005 THEN subtimestep(particle) = 2
   IF acceleration > .05 THEN subtimestep(particle) = 1
   IF acceleration > .1 THEN subtimestep(particle) = .5
   IF acceleration > .5 THEN subtimestep(particle) = .25
   IF acceleration > .8 THEN subtimestep(particle) = .125
   IF acceleration > 1.5 THEN subtimestep(particle) = .0625
   IF subtimestep(particle) > timestep THEN subtimestep(particle) = timestep
END SUB

SUB DisplayInfo
   'Gives information while paused during simulation
 
   LOCATE 1, 1: PRINT "    Name: "; simname
   LOCATE 2, 1: PRINT "    Time:"; displaytime
   LOCATE 3, 1: PRINT "         "; displaytime * 1.227 * 10 ^ 6; "yr"
   LOCATE 4, 1: PRINT "Timestep:"; timestep; : IF difftimestep = 1 THEN PRINT "(differential)"
   LOCATE 5, 1: PRINT "   Delay:"; speed - 1; "  "
   LOCATE 6, 1: PRINT "X offset:"; offx
   LOCATE 7, 1: PRINT "Y offset:"; offy
   LOCATE 8, 1: PRINT "    Zoom:"; INT((640 / (winx2 - winx1)) * 100); "%"
   LOCATE 9, 1: PRINT "C offset:"; CINT(offc * 10 / .174532); "  "
   LOCATE 10, 1: PRINT "B offset:"; CINT(offb * 10 / .174532); "  "
   LOCATE 11, 1: PRINT "A offset:"; CINT(offa * 10 / .174532); "  "
   LOCATE 12, 1: PRINT "Contrast:"; contrast
   LOCATE 13, 1: PRINT "Rest Vel:"; restvel
   LOCATE 14, 1: PRINT " Tracers: ";
                IF tracers = 0 THEN PRINT "off"
                IF tracers = 1 THEN PRINT "on "
   LOCATE 15, 1: PRINT " Storage: ";
                IF store = 0 THEN PRINT "none"
                IF store = 1 THEN PRINT "RAM "
                IF store = 2 THEN PRINT "disk"
   LOCATE 16, 1: PRINT "  Source: ";
                IF source = 0 THEN PRINT "program"
                IF source = 1 THEN PRINT "file   "
   LOCATE 17, 1: PRINT "Tracking:";
                IF track > 0 THEN PRINT track
                IF track = 0 THEN PRINT " manual"
                IF track = -1 THEN PRINT " center"
END SUB

SUB DopplerColor
   'Assigns each particle a color based on its line-of-sight velocity

   xtrans0 = vx(particle)
   ytrans0 = vy(particle)
   ztrans0 = vz(particle)

   xtrans1 = xtrans0 * COS(offa) + ytrans0 * SIN(offa)
   ytrans1 = ytrans0 * COS(offa) - xtrans0 * SIN(offa)
   ztrans1 = ztrans0

   xtrans2 = xtrans1 * COS(offb) + ztrans1 * SIN(offb)
   ytrans2 = ytrans1
   ztrans2 = ztrans1 * COS(offb) - xtrans1 * SIN(offb)

   xtrans3 = xtrans2
   ytrans3 = ytrans2 * COS(offc) + ztrans2 * SIN(offc)
   ztrans3 = ztrans2 * COS(offc) - ytrans2 * SIN(offc)

   vs = ztrans3 - (restvel / 10)

   IF c(particle) = 15 THEN EXIT SUB

                   c(particle) = 4
   IF vs > -.1 * contrast THEN c(particle) = 12
   IF vs > -.075 * contrast THEN c(particle) = 14
   IF vs > -.05 * contrast THEN c(particle) = 10
   IF vs > -.025 * contrast THEN c(particle) = 2
  
   IF vs > 0 THEN c(particle) = 3
   IF vs > .025 * contrast THEN c(particle) = 9
   IF vs > .05 * contrast THEN c(particle) = 1
   IF vs > .075 * contrast THEN c(particle) = 13
   IF vs > .1 * contrast THEN c(particle) = 5
END SUB

SUB DrawParticle
   'Draws particle in color pixelcolor and style pixelstyle
 
   IF xg(particle) >= winx1 AND xg(particle) <= winx2 AND yg(particle) >= winy1 AND yg(particle) <= winy2 THEN
      PSET (xg(particle), yg(particle)), pixelcolor
      IF pixelstyle = 1 THEN CIRCLE (xg(particle), yg(particle)), 1, pixelcolor
   END IF
END SUB

SUB FindPixelStyle
   'Finds image ID to use as particle representation
  
   IF pixelmode = 1 THEN
      IF c(particle) = 15 THEN pixelstyle = 1 ELSE pixelstyle = 0
   END IF
END SUB

SUB GraphicCalc
   'Calculates screen coordinates given actual coordinates and view variables

   xtrans0 = x(particle) - centx
   ytrans0 = y(particle) - centy
   ztrans0 = z(particle) - centz

   xtrans1 = xtrans0 * COS(offa) + ytrans0 * SIN(offa)
   ytrans1 = ytrans0 * COS(offa) - xtrans0 * SIN(offa)
   ztrans1 = ztrans0

   xtrans2 = xtrans1 * COS(offb) + ztrans1 * SIN(offb)
   ytrans2 = ytrans1
   ztrans2 = ztrans1 * COS(offb) - xtrans1 * SIN(offb)

   xtrans3 = xtrans2
   ytrans3 = ytrans2 * COS(offc) + ztrans2 * SIN(offc)
   ztrans3 = ztrans2 * COS(offc) - ytrans2 * SIN(offc)

   xg(particle) = xtrans3 + offx + 320 + centx
   yg(particle) = ytrans3 + offy + 240 + centy
END SUB

SUB IFGalaxies
   'Interface page 2: input for galactic and intergalactic properties

   COLOR 7, 0
   CLS
   LOCATE 1, 4: PRINT "INTERGALACTIC GRAVITATIONAL MOTION SIMULATOR"
   LOCATE 1, 66: PRINT "Galaxies"
   LOCATE 3, 3: PRINT "GALAXY 1"
   LOCATE 3, 28: PRINT "GALAXY 2"
   LOCATE 3, 46: PRINT "RETROTRAJECTORY PARAMETERS"

   FOR pnt = 1 TO 31
      GOSUB GetLabel2
      GOSUB GetSelection2
      COLOR 7, 0: LOCATE ylab, xlab: PRINT label
      COLOR 7, 0: LOCATE ysel, xsel - (value < 0): PRINT value
   NEXT pnt

   LOCATE 23, 1: COLOR 8: PRINT "F1 - Check"
   LOCATE 23, 14: COLOR 7: PRINT "F2 - Particles "
   LOCATE 23, 31: COLOR 7: PRINT "F3 - Exit"
   LOCATE 23, 63: COLOR 8: PRINT "F5 - Catalog"

   IF errcode = 0 THEN pnt = 1
   GOSUB ErrorMotion
   errcode = 0

   DO
      IF errcode = 19 THEN pnt = 28: errcode = 0
     
      IF IGdistance > 0 AND IGvelocity <> 0 AND gm(1) > 0 AND gm(2) > 0 THEN COLOR 7, 0 ELSE COLOR 8, 0
      LOCATE 23, 43: PRINT "F4 - Retroproject"

      IF scancode = 72 THEN pnt = pnt + -1
      IF scancode = 80 THEN pnt = pnt + 1

      IF pnt > 31 THEN pnt = 1
      IF pnt < 1 THEN pnt = 31

      GOSUB GetLabel2
      GOSUB GetSelection2

      LOCATE ylab, xlab:  COLOR 15, 0: PRINT label
      LOCATE ysel, xsel + 1: COLOR 15, 3: PRINT "             "
      ypos = ysel: xpos = xsel: lim = 13: inputstring = STR$(value): msy = 1: CALL LimInput

      GOSUB ChangeValues
                                           
      LOCATE ysel, xsel: COLOR 0, 0: PRINT "              "
      LOCATE ylab, xlab: COLOR 7, 0: PRINT label
      LOCATE ysel, xsel - (value < 0): COLOR 7, 0: PRINT value

      IF scancode = 59 THEN check = -1
      IF scancode = 60 THEN nextif = 1 ELSE nextif = 0
      IF scancode = 61 THEN exitflag = 1 ELSE exitflag = 0
      IF scancode = 62 AND IGdistance > 0 AND IGvelocity <> 0 AND gm(1) > 0 AND gm(2) > 0 THEN
         CALL RetroProject
         FOR pnt = 1 TO 31
            GOSUB GetLabel2
            GOSUB GetSelection2
            COLOR 7, 0: LOCATE ylab, xlab: PRINT label
            COLOR 7, 0: LOCATE ysel, xsel - (value < 0): PRINT "              "
            COLOR 7, 0: LOCATE ysel, xsel - (value < 0): PRINT value
         NEXT pnt
      END IF

      IF nextif = 1 THEN
         GOSUB HaloCheck
         GOSUB ErrorMotion
         IF errcode = 20 AND rings(2) * particles(2) > rings(1) * particles(1) THEN pnt = 24
         IF errcode = 20 AND rings(1) * particles(1) >= rings(2) * particles(2) THEN pnt = 11
         IF errcode > 0 THEN nextif = 0
      END IF
   LOOP UNTIL exitflag = 1 OR nextif = 1
EXIT SUB

'************************
GetLabel2:
   IF pnt = 1 THEN ylab = 4: xlab = 1: label = "        x:"
   IF pnt = 2 THEN ylab = 5: xlab = 1: label = "        y:"
   IF pnt = 3 THEN ylab = 6: xlab = 1: label = "        z:"
   IF pnt = 4 THEN ylab = 8: xlab = 1: label = "       Vx:"
   IF pnt = 5 THEN ylab = 9: xlab = 1: label = "       Vy:"
   IF pnt = 6 THEN ylab = 10: xlab = 1: label = "       Vz:"
   IF pnt = 7 THEN ylab = 12: xlab = 1: label = "     mass:"
   IF pnt = 8 THEN ylab = 14: xlab = 1: label = "    rings:"
   IF pnt = 9 THEN ylab = 15: xlab = 1: label = "     Rmin:"
   IF pnt = 10 THEN ylab = 16: xlab = 1: label = "     Rmax:"
   IF pnt = 11 THEN ylab = 17: xlab = 1: label = "particles:"
   IF pnt = 12 THEN ylab = 19: xlab = 1: label = "     tilt:"
   IF pnt = 13 THEN ylab = 21: xlab = 1: label = "    color:"

   IF pnt = 14 THEN ylab = 4: xlab = 26: label = "        x:"
   IF pnt = 15 THEN ylab = 5: xlab = 26: label = "        y:"
   IF pnt = 16 THEN ylab = 6: xlab = 26: label = "        z:"
   IF pnt = 17 THEN ylab = 8: xlab = 26: label = "       Vx:"
   IF pnt = 18 THEN ylab = 9: xlab = 26: label = "       Vy:"
   IF pnt = 19 THEN ylab = 10: xlab = 26: label = "       Vz:"
   IF pnt = 20 THEN ylab = 12: xlab = 26: label = "     mass:"
   IF pnt = 21 THEN ylab = 14: xlab = 26: label = "    rings:"
   IF pnt = 22 THEN ylab = 15: xlab = 26: label = "     Rmin:"
   IF pnt = 23 THEN ylab = 16: xlab = 26: label = "     Rmax:"
   IF pnt = 24 THEN ylab = 17: xlab = 26: label = "particles:"
   IF pnt = 25 THEN ylab = 19: xlab = 26: label = "     tilt:"
   IF pnt = 26 THEN ylab = 21: xlab = 26: label = "    color:"

   IF pnt = 27 THEN ylab = 4: xlab = 51: label = "     Distance:"
   IF pnt = 28 THEN ylab = 6: xlab = 51: label = "     Velocity:"
   IF pnt = 29 THEN ylab = 8: xlab = 51: label = "  Inclination:"
   IF pnt = 30 THEN ylab = 10: xlab = 51: label = "     Argument:"
   IF pnt = 31 THEN ylab = 12: xlab = 51: label = " Initial Time:"
RETURN

GetSelection2:
   value = 0
   IF pnt = 1 THEN ysel = 4: xsel = 11: value = gx(1)
   IF pnt = 2 THEN ysel = 5: xsel = 11: value = gy(1)
   IF pnt = 3 THEN ysel = 6: xsel = 11: value = gz(1)
   IF pnt = 4 THEN ysel = 8: xsel = 11: value = gvx(1)
   IF pnt = 5 THEN ysel = 9: xsel = 11: value = gvy(1)
   IF pnt = 6 THEN ysel = 10: xsel = 11: value = gvz(1)
   IF pnt = 7 THEN ysel = 12: xsel = 11: value = gm(1)
   IF pnt = 8 THEN ysel = 14: xsel = 11: value = rings(1)
   IF pnt = 9 THEN ysel = 15: xsel = 11: value = rmin(1)
   IF pnt = 10 THEN ysel = 16: xsel = 11: value = rmax(1)
   IF pnt = 11 THEN ysel = 17: xsel = 11: value = particles(1)
   IF pnt = 12 THEN ysel = 19: xsel = 11: value = tilt(1)
   IF pnt = 13 THEN ysel = 21: xsel = 11: value = gc(1)

   IF pnt = 14 THEN ysel = 4: xsel = 36: value = gx(2)
   IF pnt = 15 THEN ysel = 5: xsel = 36: value = gy(2)
   IF pnt = 16 THEN ysel = 6: xsel = 36: value = gz(2)
   IF pnt = 17 THEN ysel = 8: xsel = 36: value = gvx(2)
   IF pnt = 18 THEN ysel = 9: xsel = 36: value = gvy(2)
   IF pnt = 19 THEN ysel = 10: xsel = 36: value = gvz(2)
   IF pnt = 20 THEN ysel = 12: xsel = 36: value = gm(2)
   IF pnt = 21 THEN ysel = 14: xsel = 36: value = rings(2)
   IF pnt = 22 THEN ysel = 15: xsel = 36: value = rmin(2)
   IF pnt = 23 THEN ysel = 16: xsel = 36: value = rmax(2)
   IF pnt = 24 THEN ysel = 17: xsel = 36: value = particles(2)
   IF pnt = 25 THEN ysel = 19: xsel = 36: value = tilt(2)
   IF pnt = 26 THEN ysel = 21: xsel = 36: value = gc(2)

   IF pnt = 27 THEN ysel = 4: xsel = 65: value = IGdistance
   IF pnt = 28 THEN ysel = 6: xsel = 65: value = IGvelocity
   IF pnt = 29 THEN ysel = 8: xsel = 65: value = IGinclination
   IF pnt = 30 THEN ysel = 10: xsel = 65: value = IGargument
   IF pnt = 31 THEN ysel = 12: xsel = 65: value = IGt0
RETURN

ChangeValues:
   oldvalue = value
   ON ERROR GOTO Overflow
   value = VAL(inputstring)
   IF pnt = 1 THEN gx(1) = value
   IF pnt = 2 THEN gy(1) = value
   IF pnt = 3 THEN gz(1) = value
   IF pnt = 4 THEN gvx(1) = value
   IF pnt = 5 THEN gvy(1) = value
   IF pnt = 6 THEN gvz(1) = value
   IF pnt = 7 THEN value = ABS(value): gm(1) = value
   IF pnt = 8 THEN value = INT(ABS(value)) MOD 15: rings(1) = value
   IF pnt = 9 THEN value = INT(ABS(value)): rmin(1) = value
   IF pnt = 10 THEN value = INT(ABS(value)): rmax(1) = value
   IF pnt = 11 THEN value = INT(ABS(value)): particles(1) = value
   IF pnt = 12 THEN
      IF ABS(value) = 180 THEN value = 180 ELSE value = value MOD 180
      tilt(1) = value
   END IF
   IF pnt = 13 THEN value = INT(ABS(value)) MOD 16: gc(1) = value
                                       
   IF pnt = 14 THEN gx(2) = value
   IF pnt = 15 THEN gy(2) = value
   IF pnt = 16 THEN gz(2) = value
   IF pnt = 17 THEN gvx(2) = value
   IF pnt = 18 THEN gvy(2) = value
   IF pnt = 19 THEN gvz(2) = value
   IF pnt = 20 THEN value = ABS(value): gm(2) = value
   IF pnt = 21 THEN value = INT(ABS(value)) MOD 15: rings(2) = value
   IF pnt = 22 THEN value = INT(ABS(value)): rmin(2) = value
   IF pnt = 23 THEN value = INT(ABS(value)): rmax(2) = value
   IF pnt = 24 THEN value = INT(ABS(value)): particles(2) = value
   IF pnt = 25 THEN
      IF ABS(value) = 180 THEN value = 180 ELSE value = value MOD 180
      tilt(2) = value
   END IF
   IF pnt = 26 THEN value = INT(ABS(value)) MOD 16: gc(2) = value

   IF pnt = 27 THEN IGdistance = value
   IF pnt = 28 THEN IGvelocity = value
   IF pnt = 29 THEN
      IF ABS(value) = 180 THEN value = 180 ELSE value = value MOD 180
      IGinclination = value
   END IF
   IF pnt = 30 THEN
      IF value <> 180 THEN value = ABS(value) MOD 180
      IGargument = value
   END IF
   IF pnt = 31 THEN value = -ABS(value): IGt0 = value
RETURN

HaloCheck:
   errcode = 0
   IF -(gm(1) > 0) + -(gm(2) > 0) + rings(1) * particles(1) + rings(2) * particles(2) > maxparticles THEN errcode = 20
   IF rings(1) > 0 AND rmin(1) = 0 THEN errcode = 8
   IF rings(2) > 0 AND rmin(2) = 0 THEN errcode = 9
   IF rmax(1) < rmin(1) OR rings(1) > 1 AND rmax(1) = rmin(1) THEN errcode = 10
   IF rmax(2) < rmin(2) OR rings(2) > 1 AND rmax(2) = rmin(2) THEN errcode = 11
   IF rings(1) > 0 AND particles(1) = 0 THEN errcode = 12
   IF rings(2) > 0 AND particles(2) = 0 THEN errcode = 13
   IF gm(1) = 0 AND rings(1) > 0 THEN errcode = 6
   IF gm(2) = 0 AND rings(2) > 0 THEN errcode = 7
   IF rings(1) = 0 AND (tilt(1) <> 0 OR rmin(1) > 0 OR rmax(1) > 0 OR particles(1) > 0) THEN errcode = 14
   IF rings(2) = 0 AND (tilt(2) <> 0 OR rmin(2) > 0 OR rmax(2) > 0 OR particles(2) > 0) THEN errcode = 15
   IF gm(1) > 0 AND gm(2) > 0 AND gx(1) = gx(2) AND gy(1) = gy(2) AND gz(1) = gz(2) THEN errcode = 5
RETURN

ErrorMotion:
   IF errcode = 5 THEN pnt = 14
   IF errcode = 6 THEN pnt = 7
   IF errcode = 7 THEN pnt = 20
   IF errcode = 8 THEN pnt = 9
   IF errcode = 9 THEN pnt = 22
   IF errcode = 10 THEN pnt = 10
   IF errcode = 11 THEN pnt = 23
   IF errcode = 12 THEN pnt = 11
   IF errcode = 13 THEN pnt = 24
   IF errcode = 14 THEN pnt = 8
   IF errcode = 15 THEN pnt = 21
RETURN
END SUB

SUB IFModule
   'Center for interface operations

   numparticles = 0
   nextif = 1
   runsim = 0
   pnt = 1

   SCREEN 0
   CLS

   DO
      IF nextif = 1 THEN CALL IFProperties
      IF nextif = 1 THEN CALL IFGalaxies
      IF nextif = 1 THEN CALL IFParticles
      IF viewcat = 1 THEN CALL Catalog
   LOOP UNTIL exitflag = 1 OR runsim = 1

   IF exitflag = 1 THEN END
END SUB

SUB IFParticles
   'Interface page 3:  Calculates and prints initial particle conditions

   ON ERROR GOTO 0
   np = 1
   IF gm(1) > 0 THEN m(np) = gm(1): x(np) = gx(1): y(np) = gy(1): z(np) = gz(1): vx(np) = gvx(1): vy(np) = gvy(1): vz(np) = gvz(1): c(np) = 15: c2(np) = 8: np = np + 1
   IF gm(2) > 0 THEN m(np) = gm(2): x(np) = gx(2): y(np) = gy(2): z(np) = gz(2): vx(np) = gvx(2): vy(np) = gvy(2): vz(np) = gvz(2): c(np) = 15: c2(np) = 8: np = np + 1
  
   FOR galaxy = 1 TO 2
      IF rings(galaxy) > 0 AND gm(galaxy) > 0 THEN
         ringparticles = particles(galaxy)
         IF rings(galaxy) > 1 THEN ringspace = (rmax(galaxy) - rmin(galaxy)) / (rings(galaxy) - 1)
         IF rings(galaxy) = 1 THEN ringspace = 0
         R = rmin(galaxy)
         FOR n = 1 TO rings(galaxy)
            GOSUB CreateRing
            R = R + ringspace
         NEXT n
      END IF
   NEXT galaxy

   numparticles = np - 1
   CLS
   LOCATE 1, 4: PRINT "INTERGALACTIC GRAVITATIONAL MOTION SIMULATOR"
   LOCATE 1, 66: PRINT "Particles"
   LOCATE 3, 1
   PRINT TAB(1); " #"; TAB(5); " x, Vx"; TAB(22); " y, Vy"; TAB(39); " z, Vz"; TAB(56); " mass"; TAB(64); " color"
   LOCATE 23, 1: COLOR 8: PRINT "F1- Check"
   LOCATE 23, 14: COLOR 7: PRINT "F2 - Properties"
   LOCATE 23, 31: COLOR 7: PRINT "F3 - Exit"
   LOCATE 23, 43: COLOR 8: PRINT "F4 - Retroproject"
   LOCATE 23, 63: COLOR 8: PRINT "F5 - Catalog"

   COLOR 7
   page = 0
   DO
      LOCATE 4, 1
      IF numparticles > (page + 1) * 9 THEN ulm = (page + 1) * 9 ELSE ulm = numparticles
      FOR ypos = 4 TO 22
         PRINT "                                                                              "
      NEXT ypos
      LOCATE 4, 1
      FOR particle = page * 9 + 1 TO ulm
         PRINT TAB(1); particle; TAB(6); x(particle); TAB(23); y(particle); TAB(40); z(particle); TAB(57); m(particle); TAB(65); c(particle)
         PRINT TAB(6); vx(particle); TAB(23); vy(particle); TAB(40); vz(particle); TAB(57);
      NEXT particle

      DO
         inkeyv = INKEY$
      LOOP WHILE inkeyv = ""
      IF LEN(inkeyv) > 1 THEN scancode = ASC(MID$(inkeyv, 2, 1)) ELSE scancode = 0

      IF scancode = 72 AND page > 0 THEN page = page - 1
      IF scancode = 80 AND ulm < numparticles THEN page = page + 1
      IF scancode = 60 THEN nextif = 1 ELSE nextif = 0
      IF scancode = 61 THEN exitflag = 1

   LOOP UNTIL nextif = 1 OR exitflag = 1
   pnt = 1
EXIT SUB

CreateRing:
   ia = 0
   FOR np = np TO np + ringparticles - 1
      m(np) = 0
      c(np) = gc(galaxy)
      c2(np) = 0
     
      ia = ia + (2 * 3.14159) / ringparticles
      x(np) = gx(galaxy) + COS(ia) * R
      y(np) = gy(galaxy) + SIN(ia) * R
      z(np) = gz(galaxy)
      vx(np) = gvx(galaxy) - ABS(R) / R * (G * gm(galaxy) / ABS(R)) ^ (1 / 2) * SIN(ia)
      vy(np) = gvy(galaxy) + ABS(R) / R * (G * gm(galaxy) / ABS(R)) ^ (1 / 2) * COS(ia)
      vz(np) = gvz(galaxy)

      tiltrad = tilt(galaxy) * .174532 / 10

      xtrans0 = x(np) - gx(galaxy)
      ytrans0 = y(np) - gy(galaxy)
      ztrans0 = z(np) - gz(galaxy)
      GOSUB Transform2
      x(np) = xtrans1 + gx(galaxy)
      y(np) = ytrans1 + gy(galaxy)
      z(np) = ztrans1 + gz(galaxy)

      xtrans0 = vx(np) - gvx(galaxy)
      ytrans0 = vy(np) - gvy(galaxy)
      ztrans0 = vz(np) - gvz(galaxy)
      GOSUB Transform2
      vx(np) = xtrans1 + gvx(galaxy)
      vy(np) = ytrans1 + gvy(galaxy)
      vz(np) = ztrans1 + gvz(galaxy)
   NEXT np
RETURN

Transform2:
      xtrans1 = xtrans0 * COS(tiltrad) + ztrans0 * SIN(tiltrad)
      ytrans1 = ytrans0
      ztrans1 = ztrans0 * COS(tiltrad) - xtrans0 * SIN(tiltrad)
RETURN
END SUB

SUB IFProperties
   'Interface page 1:  Input for metaproperties and miscellaneous properties

   CLS
   LOCATE 1, 4: PRINT "INTERGALACTIC GRAVITATIONAL MOTION SIMULATOR"
   LOCATE 1, 66: PRINT "Properties"
   LOCATE 3, 1: PRINT "FILE"
   LOCATE 14, 1: PRINT "TIME"
   LOCATE 3, 40: PRINT "DISPLAY"
   LOCATE 23, 31: COLOR 7: PRINT "F3 - Exit"
   LOCATE 23, 43: COLOR 8: PRINT "F4 - Retroproject"
   LOCATE 23, 63: COLOR 7: PRINT "F5 - Catalog"

   viewcat = pnt
   COLOR 7, 0
   FOR pnt = 1 TO 13
      GOSUB GetLabel1
      GOSUB GetSel
      GOSUB GetSelection1
      LOCATE ylab, xlab: PRINT label
      LOCATE ysel, xsel: PRINT selection
   NEXT
   pnt = viewcat

   DO
      COLOR 7, 0
      LOCATE 23, 1: IF check = 1 THEN PRINT "F1 - Run  " ELSE PRINT "F1 - Check"
   
      IF selection <> "" THEN LOCATE ysel, xsel: COLOR 0, 0: PRINT "                "
      IF motion = 1 OR errcode > 0 THEN LOCATE ylab, xlab: COLOR 7, 0: PRINT label
      IF motion = 1 OR errcode > 0 THEN LOCATE ysel, xsel: COLOR 7, 0: PRINT selection

      IF errcode = 1 THEN pnt = 1: sel = 1
      IF errcode = 2 THEN pnt = 3: sel = store
      IF errcode = 3 THEN pnt = 4: sel = 0
      IF errcode = 4 THEN pnt = 9: sel = gtrack - (-1 = gtrack) * 4
      IF errcode = 18 THEN pnt = 1: sel = 0
      errcode = 0
    
      GOSUB GetLabel1
      GOSUB GetSelection1

      IF source = 0 THEN COLOR 7 ELSE COLOR 8
      LOCATE 23, 14: PRINT "F2 - Galaxies"

      COLOR 15, 3
      IF pnt = 1 AND sel = 1 THEN
         COLOR 15, 0: LOCATE ylab, xlab: PRINT label
         COLOR 15, 3: LOCATE ysel, xsel: PRINT "        "
         lim = 6: xpos = xsel: ypos = ysel: inputstring = "": msy = 0: CALL LimInput
         simname = inputstring
         sel = 0: skipinput = 1
         selection = " " + inputstring
      END IF
   
      IF (pnt = 12 OR pnt = 13) AND sel = 1 THEN
         COLOR 15, 0: LOCATE ysel, xsel - 4: PRINT CHR$(16)
         COLOR 15, 3: LOCATE ysel, xsel: PRINT "       " + ".SET "
         lim = 6: xpos = xsel: ypos = ysel: inputstring = "": msy = 0: CALL LimInput
         COLOR 0, 0: LOCATE ysel, xsel - 4: PRINT "                "
         IF pnt = 13 AND inputstring <> "" THEN
            filename = inputstring + ".SET"
            fileindex = FREEFILE
            OPEN "" + filename FOR OUTPUT AS #fileindex
            WRITE #fileindex, simname, source, store, recfreq, timestep, difftimestep, speed, monitor, gtrack, gx(1), gy(1), gz(1), gvx(1), gvy(1), gvz(1), gm(1), rings(1), rmin(1), rmax(1), particles(1), tilt(1), gc(1), gx(2), gy(2), gz(2), gvx(2), gvy(2), gvz(2), gm(2), rings(2), rmin(2), rmax(2), particles(2), tilt(2), gc(2), IGdistance, IGvelocity, IGinclination, IGargument, IGt0, doppler, pixelstyle, pixelmode
            CLOSE #fileindex
         END IF
         IF pnt = 12 AND inputstring <> "" THEN
            filename = inputstring + ".SET"
            fileindex = FREEFILE
            ON ERROR GOTO FileNotFound
            OPEN "" + filename FOR INPUT AS #fileindex
            INPUT #fileindex, simname, source, store, recfreq, timestep, difftimestep, speed, monitor, gtrack, gx(1), gy(1), gz(1), gvx(1), gvy(1), gvz(1), gm(1), rings(1), rmin(1), rmax(1), particles(1), tilt(1), gc(1), gx(2), gy(2), gz(2), gvx(2), gvy(2), gvz(2), gm(2), rings(2), rmin(2), rmax(2), particles(2), tilt(2), gc(2), IGdistance, IGvelocity, IGinclination, IGargument, IGt0
            ON ERROR GOTO MissingEntry
            INPUT #fileindex, doppler, pixelstyle, pixelmode
            CLOSE #fileindex
            COLOR 7, 0
            FOR pnt = 1 TO 11
               GOSUB GetSel
               GOSUB GetSelection1
               LOCATE ysel, xsel: PRINT "                "
               LOCATE ysel, xsel: PRINT selection
               errcode = 0
            NEXT pnt
            pnt = 12: selection = ""
         END IF
         sel = 0: skipinput = 1
      END IF
      LOCATE ylab, xlab:  COLOR 15, 0: PRINT label
      LOCATE ysel, xsel:  COLOR 15, 3: PRINT selection
      IF pnt = 1 THEN LOCATE ysel, xsel:  COLOR 15, 0: PRINT selection
  
      GOSUB IFUserInput
      
      IF motion <> 2 THEN GOSUB GetSel

      IF motion = 2 THEN check = 0

      IF pnt > 13 THEN pnt = 1
      IF pnt < 1 THEN pnt = 13

      IF (pnt = 2 OR pnt = 6 OR pnt = 8) AND sel = -1 THEN sel = 1
      IF (pnt = 2 OR pnt = 6 OR pnt = 8) AND sel = 2 THEN sel = 0
      IF pnt = 3 AND sel = -1 THEN sel = 2
      IF pnt = 3 AND sel = 3 THEN sel = 0
      IF pnt = 9 AND sel = -1 THEN sel = 3
      IF pnt = 9 AND sel = 4 THEN sel = 0
      IF pnt = 10 AND sel = 2 THEN sel = 0
      IF pnt = 10 AND sel = -1 THEN sel = 1
      IF pnt = 11 AND sel = 3 THEN sel = 0
      IF pnt = 11 AND sel = -1 THEN sel = 2
  
      IF check = -1 THEN GOSUB CheckSels
      COLOR 7, 0
   LOOP UNTIL exitflag = 1 OR nextif = 1 OR runsim = 1 OR viewcat = 1
   check = 0
EXIT SUB

'*******************
GetLabel1:
   IF pnt = 1 THEN ylab = 4: xlab = 4: label = "   Simulation:"
   IF pnt = 2 THEN ylab = 6: xlab = 4: label = "        Input:"
   IF pnt = 3 THEN ylab = 8: xlab = 4: label = "       Output:"
   IF pnt = 4 THEN ylab = 10: xlab = 4: label = "Rec. Freqency:"
   IF pnt = 5 THEN ylab = 15: xlab = 4: label = "     Timestep:"
   IF pnt = 6 THEN ylab = 17: xlab = 4: label = " Differential:"
   IF pnt = 7 THEN ylab = 19: xlab = 4: label = "        Delay:"
   IF pnt = 8 THEN ylab = 4: xlab = 43: label = "      Monitor:"
   IF pnt = 9 THEN ylab = 6: xlab = 43: label = "     Tracking:"
   IF pnt = 10 THEN ylab = 8: xlab = 43: label = "Doppler Mimic:"
   IF pnt = 11 THEN ylab = 10: xlab = 43: label = "Partcl. Style:"
   IF pnt = 12 THEN ylab = 15: xlab = 43: label = "Load Settings"
   IF pnt = 13 THEN ylab = 17: xlab = 43: label = "Save Settings"
RETURN

GetSelection1:
   selection = ""
   IF pnt = 1 AND simname <> "" THEN ysel = 4: xsel = 19: selection = " " + simname + " "
   IF pnt = 2 AND sel = 0 THEN ysel = 6: xsel = 19: selection = " program ": source = 0
   IF pnt = 2 AND sel = 1 THEN ysel = 6: xsel = 19: selection = " file ": source = 1
   IF pnt = 3 AND sel = 0 THEN ysel = 8: xsel = 19: selection = " none ": store = 0
   IF pnt = 3 AND sel = 1 THEN ysel = 8: xsel = 19: selection = " RAM ": store = 1
   IF pnt = 3 AND sel = 2 THEN ysel = 8: xsel = 19: selection = " disk ": store = 2
   IF pnt = 4 THEN
      holdvar = recfreq
      IF sel = -1 AND holdvar > 0 THEN recfreq = recfreq - 1
      IF sel = 1 AND holdvar >= 0 THEN recfreq = recfreq + 1
      IF sel = -1 AND holdvar > 30 THEN recfreq = recfreq - 1
      IF sel = 1 AND holdvar >= 30 THEN recfreq = recfreq + 1
      IF sel = -1 AND holdvar > 60 THEN recfreq = recfreq - 3
      IF sel = 1 AND holdvar >= 60 THEN recfreq = recfreq + 3
      IF recfreq > 100 THEN recfreq = 100
      ysel = 10: xsel = 19: selection = STR$(recfreq) + " "
      sel = 0
   END IF
   IF pnt = 5 THEN
      IF sel = -1 AND timestep > .03125 THEN timestep = timestep / 2
      IF sel = 1 AND timestep < 16 THEN timestep = timestep * 2
      ysel = 15: xsel = 19: selection = STR$(timestep) + " "
      sel = 0
   END IF
   IF pnt = 6 AND sel = 0 THEN ysel = 17: xsel = 19: selection = " disabled ": difftimestep = 0
   IF pnt = 6 AND sel = 1 THEN ysel = 17: xsel = 19: selection = " enabled ": difftimestep = 1
   IF pnt = 7 THEN
      spdchg = -sel: CALL ChgSpd
      ysel = 19: xsel = 19: selection = STR$(speed - 1) + " "
      sel = 0
   END IF
   IF pnt = 8 AND sel = 0 THEN ysel = 4: xsel = 58: selection = " disabled ": monitor = 0
   IF pnt = 8 AND sel = 1 THEN ysel = 4: xsel = 58: selection = " enabled ": monitor = 1
   IF pnt = 9 AND sel = 0 THEN ysel = 6: xsel = 58: selection = " manual ": gtrack = 0
   IF pnt = 9 AND sel = 1 THEN ysel = 6: xsel = 58: selection = " galaxy 1 ": gtrack = 1
   IF pnt = 9 AND sel = 2 THEN ysel = 6: xsel = 58: selection = " galaxy 2 ": gtrack = 2
   IF pnt = 9 AND sel = 3 THEN ysel = 6: xsel = 58: selection = " center of mass ": gtrack = -1
   IF pnt = 10 AND sel = 0 THEN ysel = 8: xsel = 58: selection = " disabled ": doppler = 0
   IF pnt = 10 AND sel = 1 THEN ysel = 8: xsel = 58: selection = " enabled ": doppler = 1
   IF pnt = 11 AND sel = 0 THEN ysel = 10: xsel = 58: selection = " small ": pixelstyle = 0: pixelmode = 0
   IF pnt = 11 AND sel = 1 THEN ysel = 10: xsel = 58: selection = " cross ": pixelstyle = 1: pixelmode = 0
   IF pnt = 11 AND sel = 2 THEN ysel = 10: xsel = 58: selection = " variable ": pixelstyle = 1: pixelmode = 1
   IF pnt = 12 AND sel = 1 THEN ysel = 15: xsel = 62: selection = ""
   IF pnt = 13 AND sel = 1 THEN ysel = 17: xsel = 62: selection = ""
RETURN

GetSel:
   sel = 0
   IF pnt = 2 THEN sel = source
   IF pnt = 3 THEN sel = store
   IF pnt = 6 THEN sel = difftimestep
   IF pnt = 8 THEN sel = monitor
   IF pnt = 9 THEN sel = gtrack
   IF pnt = 9 AND sel = -1 THEN sel = 3
   IF pnt = 10 THEN sel = doppler
   IF pnt = 11 AND pixelmode = 0 THEN sel = pixelstyle
   IF pnt = 11 AND pixelmode = 1 THEN sel = 2
RETURN

CheckSels:
   errcode = 0
   'overlap on particle menu:  errcode = 16
   IF source = 0 THEN
      IF rings(1) > 0 AND rmin(1) = 0 THEN errcode = 8
      IF rings(2) > 0 AND rmin(2) = 0 THEN errcode = 9
      IF rmax(1) < rmin(1) OR rings(1) > 1 AND rmax(1) = rmin(1) THEN errcode = 10
      IF rmax(2) < rmin(2) OR rings(2) > 1 AND rmax(2) = rmin(2) THEN errcode = 11
      IF rings(1) > 0 AND particles(1) = 0 THEN errcode = 12
      IF rings(2) > 0 AND particles(2) = 0 THEN errcode = 13
      IF gm(1) = 0 AND rings(1) > 0 THEN errcode = 6
      IF gm(2) = 0 AND rings(2) > 0 THEN errcode = 7
      IF rings(1) = 0 AND (tilt(1) <> 0 OR rmin(1) > 0 OR rmax(1) > 0 OR particles(1) > 0) THEN errcode = 14
      IF rings(2) = 0 AND (tilt(2) <> 0 OR rmin(2) > 0 OR rmax(2) > 0 OR particles(2) > 0) THEN errcode = 15
      IF gm(1) > 0 AND gm(2) > 0 AND x(1) = x(2) AND y(1) = y(2) AND z(1) = z(2) THEN errcode = 5
      IF numparticles = 0 THEN errcode = 17
   END IF
   IF source = 1 AND simname <> "" THEN
      filename = simname + ".SIM"
      fileindex = FREEFILE
      ON ERROR GOTO FileNotFound
      OPEN "" + filename FOR INPUT AS #fileindex
      IF errcode <> 18 THEN CLOSE #fileindex
   END IF
   IF (gtrack = 1 AND gm(1) = 0 OR gtrack = 2 AND gm(2) = 0 OR gtrack = -1 AND (gm(1) = 0 OR gm(2) = 0)) AND source = 0 THEN errcode = 4
   IF store > 0 AND recfreq = 0 THEN errcode = 3
   IF source = 1 AND store > 0 THEN errcode = 2
   IF (source = 1 OR store = 2) AND simname = "" THEN errcode = 1

   IF errcode > 0 THEN check = 0 ELSE check = 1
   IF errcode >= 1 AND errcode <= 4 OR errcode = 18 THEN skipinput = 1
   IF errcode >= 5 AND errcode <= 14 THEN nextif = 1

   IF source = 1 AND errcode = 0 THEN
       GOSUB GetLabel1
       COLOR 7, 0
       LOCATE ylab, xlab: PRINT label
       filename = simname + "I" + ".SET"
       fileindex = FREEFILE
       ON ERROR GOTO FileNotFound
       OPEN "" + filename FOR INPUT AS #fileindex
       INPUT #fileindex, simname, n, n, recfreq, timestep, difftimestep, n, n, n, gx(1), gy(1), gz(1), gvx(1), gvy(1), gvz(1), gm(1), rings(1), rmin(1), rmax(1), particles(1), tilt(1), gc(1), gx(2), gy(2), gz(2), gvx(2), gvy(2), gvz(2), gm(2), rings(2), rmin(2), rmax(2), particles(2), tilt(2), gc(2), IGdistance, IGvelocity, IGinclination, IGargument, IGt0
       CLOSE #fileindex

       store = 0
       source = 1
       FOR pnt = 1 TO 12
          GOSUB GetSel
          GOSUB GetSelection1
          LOCATE ysel, xsel: PRINT selection
          errcode = 0
       NEXT
       pnt = 1
   END IF
RETURN

IFUserInput:
   IF skipinput = 0 THEN
      DO
         inkeyv = INKEY$
      LOOP WHILE inkeyv = ""
      motion = 0
   END IF
  
   skipinput = 0

   IF LEN(inkeyv) > 1 THEN scancode = ASC(MID$(inkeyv, 2, 1)) ELSE scancode = 0
  
   IF scancode = 72 THEN pnt = pnt + -1: motion = 1
   IF scancode = 80 THEN pnt = pnt + 1: motion = 1
   IF scancode = 75 THEN sel = sel - 1: motion = 2
   IF scancode = 77 THEN sel = sel + 1: motion = 2
   IF inkeyv = CHR$(32) THEN sel = sel + 1: motion = 2
   IF scancode = 59 AND check = 0 THEN check = -1
   IF scancode = 59 AND check = 1 THEN runsim = 1
   IF scancode = 60 AND source = 0 THEN nextif = 1 ELSE nextif = 0
   IF scancode = 61 THEN exitflag = 1 ELSE exitflag = 0
   IF scancode = 63 THEN viewcat = 1 ELSE viewcat = 0
   inkeyv = ""
RETURN
END SUB

SUB ImageRotate
   'Rotates view of image about center of screen
 
   IF monitor = 1 THEN
      IF tracers = 1 THEN tracers = 0
      IF tracers = 2 THEN CLS
      pixelcolor = 0
      FOR particle = 1 TO numparticles
         CALL FindPixelStyle
         CALL DrawParticle
      NEXT particle
   END IF

   SELECT CASE optioncode
   CASE 1
      offa = offa + rotateinc
      IF CINT(offa * 10 / .174532) = 185 THEN offa = -175 * .174532 / 10
   CASE 2
      offa = offa - rotateinc
      IF CINT(offa * 10 / .174532) = -180 THEN offa = 180 * .174532 / 10
   CASE 3
      offb = offb + rotateinc
      IF CINT(offb * 10 / .174532) = 185 THEN offb = -175 * .174532 / 10
   CASE 4
      offb = offb - rotateinc
      IF CINT(offb * 10 / .174532) = -180 THEN offb = 180 * .174532 / 10
   CASE 5
      offc = offc + rotateinc
      IF CINT(offc * 10 / .174532) = 185 THEN offc = -175 * .174532 / 10
   CASE 6
      offc = offc - rotateinc
      IF CINT(offc * 10 / .174532) = -180 THEN offc = 180 * .174532 / 10
   END SELECT
   
   IF monitor = 1 THEN
      FOR particle = 1 TO numparticles
         CALL FindPixelStyle
         IF doppler = 1 THEN CALL DopplerColor
         pixelcolor = c(particle)
         CALL GraphicCalc
         CALL DrawParticle
      NEXT particle
   END IF
END SUB

SUB ImageSlide
   'Moves view of image along a plane
  
   IF tracers = 1 THEN tracers = 0
   IF tracers = 2 THEN CLS
   pixelcolor = 0
   FOR particle = 1 TO numparticles
      CALL FindPixelStyle
      CALL DrawParticle
   NEXT particle
  
   IF inkeyv = CHR$(50) THEN offy = offy - 2
   IF inkeyv = CHR$(52) THEN offx = offx + 2
   IF inkeyv = CHR$(54) THEN offx = offx - 2
   IF inkeyv = CHR$(56) THEN offy = offy + 2

   FOR particle = 1 TO numparticles
      CALL FindPixelStyle
      IF doppler = 1 THEN CALL DopplerColor
      pixelcolor = c(particle)
      CALL GraphicCalc
      CALL DrawParticle
   NEXT particle
END SUB

SUB ImageZoom
   'Increases/decreases apparent size of image

   IF inkeyv = CHR$(73) OR inkeyv = CHR$(105) THEN
      IF winx1 < 288 THEN
         winx1 = winx1 + 32
         winy1 = winy1 + 24
         winx2 = winx2 - 32
         winy2 = winy2 - 24
      END IF
   END IF
  
   IF inkeyv = CHR$(79) OR inkeyv = CHR$(111) THEN
      winx1 = winx1 - 32
      winy1 = winy1 - 24
      winx2 = winx2 + 32
      winy2 = winy2 + 24
   END IF

   WINDOW SCREEN (winx1, winy1)-(winx2, winy2)
   CLS

   FOR particle = 1 TO numparticles
      CALL FindPixelStyle
      IF doppler = 1 THEN CALL DopplerColor
      pixelcolor = c(particle)
      CALL GraphicCalc
      CALL DrawParticle
   NEXT particle
END SUB

SUB Initialize
   'Specifies defaults

   ' SIMULATION UNITS
   ' smu = mass that exerts 0.1 sdu/stu ^ 2 of force over 1 sdu       
   '                             = 10 ^ 39 kg
   ' sdu = 1 pixel at 100% zoom  = 10 ^ 19 m
   ' stu = 1 timestep            = 3.871 * 10 ^ 13 s
   '                             = 1.227 * 10 ^ 6 yr

   speed = 1
   tracers = 1
   store = 0
   source = 0
   difftimestep = 0
   autosnap = 0
   winx1 = 0
   winy1 = 0
   winx2 = 640
   winy2 = 480
   simname = "UNTITL"
   timestep = 1
   subtimestep = 1
   recfreq = 6
   monitor = 1
   gtrack = 0
   speed = 1
   contrast = 10

   doppler = 1
   pixelmode = 1

   gc(1) = 10
   gc(2) = 13
END SUB

SUB Initialize2
   'Sets up simulation and specifies unalterable initial conditions
  
   FOR particle = 1 TO maxparticles
      IF difftimestep = 0 THEN subtimestep(particle) = timestep * subtimestep
      IF difftimestep = 1 THEN subtimestep(particle) = .0625
   NEXT
   IF gtrack = 1 THEN track = 1
   IF gtrack = 2 THEN track = 2
   IF gtrack = -1 THEN track = -1

   CALL StoreInitial

   time = 0
   displaytime = 0
   dataentry = 0
   autoterminate = 0
   particle = 1
   difftimestep = 0
   pauseflag = 2
   display = 1
   IF track = -1 THEN track = 0

   IF c(1) = 15 THEN c2(1) = 8
   IF c(2) = 15 THEN c2(2) = 8

   SCREEN 12
   WINDOW SCREEN (winx1, winy1)-(winx2, winy2)

END SUB

SUB LimInput
   'Special user-defined input function disallowing input after specified char
  
   inputstring = LTRIM$(inputstring)
   LOCATE ypos, xpos + 1: PRINT inputstring
   n = LEN(inputstring)
   IF n > lim THEN n = lim

   DO
      IF n < lim THEN LOCATE ypos, xpos + n + 1: PRINT "_"

      DO
         inkeyv = INKEY$
         IF inkeyv = "" THEN inkeyv = CHR$(1)
         IF LEN(inkeyv) > 1 THEN scancode = ASC(MID$(inkeyv, 2, 1)) ELSE scancode = 0: inkeyv = UCASE$(inkeyv)
      LOOP UNTIL scancode = 72 OR scancode = 80 OR scancode >= 60 AND scancode <= 63 OR ASC(UCASE$(inkeyv)) > 64 AND ASC(UCASE$(inkeyv)) < 91 AND msy = 0 OR ASC(inkeyv) = 27 OR ASC(inkeyv) = 8 OR ASC(inkeyv) = 13 OR ASC(inkeyv) > 47 AND ASC(inkeyv) < 58 OR ASC(inkeyv) = 45 AND n = 0 AND msy = 1 OR ASC(inkeyv) = 46 AND msy = 1
    
      IF ASC(inkeyv) > 64 AND ASC(inkeyv) < 91 OR ASC(inkeyv) > 47 AND ASC(inkeyv) < 58 OR ASC(inkeyv) = 45 OR ASC(inkeyv) = 46 THEN
         IF n = lim THEN inputstring = LEFT$(inputstring, n - 1)
         inputstring = inputstring + inkeyv
      END IF
      IF ASC(inkeyv) = 8 AND n > 0 THEN
         inputstring = LEFT$(inputstring, n - 1)
         IF n < lim THEN LOCATE ypos, xpos + n + 1: PRINT " "
         n = n - 1
      END IF
      IF (ASC(inkeyv) = 13 OR scancode = 80 OR scancode = 72) AND n < lim THEN LOCATE ypos, xpos + n + 1: PRINT " "
      IF n < lim AND ASC(inkeyv) <> 8 THEN n = n + 1
      IF ASC(UCASE$(inkeyv)) > 64 AND ASC(UCASE$(inkeyv)) < 91 OR ASC(inkeyv) > 47 AND ASC(inkeyv) < 58 OR ASC(inkeyv) = 45 OR ASC(inkeyv) = 46 THEN LOCATE ypos, xpos + n: PRINT inkeyv
   LOOP UNTIL inkeyv = CHR$(27) OR scancode = 80 OR scancode = 72 OR inkeyv = CHR$(13) OR scancode >= 60 AND scancode <= 63
END SUB

SUB Pause
   'Waits for input here while paused with <space>
      
   pauseflag = 1
   DO
      IF display = 1 THEN CALL DisplayInfo
      CALL UserInput
   LOOP UNTIL inkeyv = CHR$(32) OR exitflag = 1
   pauseflag = 0
   IF display = 1 AND exitflag = 0 THEN CLS
   display = 0
END SUB

SUB RetroProject
   'Finds initial coordinates of nuclei by projecting orbit in reverse time
  
   IGinclination = IGinclination * .174532 / 10
   IGargument = IGargument * .174532 / 10

   ON ERROR GOTO 0

   xtrans0 = IGdistance
   ytrans0 = 0
   ztrans0 = 0
   GOSUB Transform
   x(1) = -xtrans2 / 2
   y(1) = -ytrans2 / 2
   z(1) = -ztrans2 / 2
   x(2) = xtrans2 / 2
   y(2) = ytrans2 / 2
   z(2) = ztrans2 / 2
   perir = (xtrans2 ^ 2 + ytrans2 ^ 2 + ztrans2 ^ 2) ^ (1 / 2)

   xtrans0 = 0
   ytrans0 = IGvelocity
   ztrans0 = 0
   GOSUB Transform
    
   vx(1) = -xtrans2 / 2
   vy(1) = -ytrans2 / 2
   vz(1) = -ztrans2 / 2
   vx(2) = xtrans2 / 2
   vy(2) = ytrans2 / 2
   vz(2) = ztrans2 / 2
   m(1) = gm(1)
   m(2) = gm(2)
  
   IGinclination = IGinclination * 10 / .174532
   IGargument = IGargument * 10 / .174532

   timestep = -timestep
   subtimestep(1) = timestep
   subtimestep(2) = timestep

   rx(1) = 0: ry(1) = 0: rz(1) = 0
   rx(2) = 0: ry(2) = 0: rz(2) = 0

   retro = 1
   numparticles = 2
   time = 0
   particle = 1
   R = 0
 
   DO UNTIL time <= IGt0
      CALL Core
      IF R < perir * .97 THEN EXIT DO
   LOOP
  
   retro = 0
  
   IF R < perir * .97 AND IGt0 < 0 THEN
      errcode = 19
      gx(1) = 0: gy(1) = 0: gz(1) = 0
      gvx(1) = 0: gvy(1) = 0: gvz(1) = 0
      gx(2) = 0: gy(2) = 0: gz(2) = 0
      gvx(2) = 0: gvy(2) = 0: gvz(2) = 0
      timestep = -timestep
      EXIT SUB
   END IF

   gx(1) = x(1): gy(1) = y(1): gz(1) = z(1)
   gx(2) = x(2): gy(2) = y(2): gz(2) = z(2)
   gvx(1) = vx(1): gvy(1) = vy(1): gvz(1) = vz(1)
   gvx(2) = vx(2): gvy(2) = vy(2): gvz(2) = vz(2)
 
   tilt(1) = 0
   timestep = -timestep
   subtimestep(2) = timestep
   subtimestep(2) = timestep
EXIT SUB

'*******************
Transform: 
   xtrans1 = xtrans0 * COS(IGargument) + ytrans0 * SIN(IGargument)
   ytrans1 = ytrans0 * COS(IGargument) - xtrans0 * SIN(IGargument)
   ztrans1 = ztrans0

   xtrans2 = xtrans1 * COS(IGinclination) + ztrans1 * SIN(IGinclination)
   ytrans2 = ytrans1
   ztrans2 = ztrans1 * COS(IGinclination) - xtrans1 * SIN(IGinclination)
RETURN
END SUB

SUB Snapshot
   'Records graphical data at moment
   CLS

   IF shotnum < 10 THEN filenum = "0" + RIGHT$(STR$(shotnum), 1)
   IF shotnum > 9 AND shotnum < 100 THEN filenum = RIGHT$(STR$(shotnum), 2)
   filename = simname + filenum + ".SNP"

   IF shotnum < 100 THEN shotnum = shotnum + 1 ELSE shotnum = shotnum - 100

   fileindex = FREEFILE
   OPEN "" + filename FOR OUTPUT AS #fileindex
   WRITE #fileindex, numparticles, winx1, winy1, winx2, winy2, time
   FOR particle = 1 TO numparticles
      WRITE #fileindex, xg(particle), yg(particle), c(particle)
   NEXT particle
   CLOSE #fileindex

   FOR n = 1 TO 260000
   NEXT n

   FOR particle = 1 TO numparticles
      CALL FindPixelStyle
      IF doppler = 1 THEN CALL DopplerColor
      pixelcolor = c(particle)
      CALL DrawParticle
   NEXT particle

   LOCATE 1, 40: PRINT "Snapshot image created as "; filename
END SUB

SUB StoreInitial
   'Either records initial statistics or loads previously recorded initstats

   IF store = 2 THEN
      filename = simname + ".INT"
      fileindex = FREEFILE
      OPEN "" + filename FOR OUTPUT AS #fileindex
      WRITE #fileindex, numparticles, recfreq, timestep, difftimestep
      FOR particle = 1 TO maxparticles
         WRITE #fileindex, c(particle)
      NEXT particle
      CLOSE #fileindex
   
      filename = simname + "I" + ".SET"
      fileindex = FREEFILE
      OPEN "" + filename FOR OUTPUT AS #fileindex
      WRITE #fileindex, simname, source, store, recfreq, timestep, difftimestep, speed, monitor, gtrack, gx(1), gy(1), gz(1), gvx(1), gvy(1), gvz(1), gm(1), rings(1), rmin(1), rmax(1), particles(1), tilt(1), gc(1), gx(2), gy(2), gz(2), gvx(2), gvy(2), gvz(2), gm(2), rings(2), rmin(2), rmax(2), particles(2), tilt(2), gc(2), IGdistance, IGvelocity, IGinclination, IGargument, IGt0
      CLOSE #fileindex
   END IF
 
   IF source = 1 THEN
      filename = simname + ".INT"
      fileindex = FREEFILE
      OPEN "" + filename FOR INPUT AS #fileindex
      INPUT #fileindex, numparticles, recfreq, timestep, difftimestep
      FOR particle = 1 TO maxparticles
         INPUT #fileindex, c(particle)
      NEXT particle
      CLOSE #fileindex
   
      filename = simname + ".SIM"
      fileindex = FREEFILE
      OPEN "" + filename FOR INPUT AS #fileindex
      FOR dataentry = 0 TO maxentry
         FOR particle = 1 TO numparticles
            INPUT #fileindex, sx(particle, dataentry), sy(particle, dataentry), sz(particle, dataentry)
         NEXT particle
      NEXT dataentry
      CLOSE #fileindex
   END IF
END SUB

SUB StoreSimulation
   'Records all particle positions on disk
  
   IF dataentry > maxentry THEN
      filename = simname + ".SIM"
      fileindex = FREEFILE
      OPEN "" + filename FOR OUTPUT AS #fileindex
      FOR dataentry = 0 TO maxentry
         FOR particle = 1 TO numparticles
            WRITE #fileindex, sx(particle, dataentry), sy(particle, dataentry), sz(particle, dataentry)
         NEXT particle
      NEXT dataentry
      CLOSE #fileindex
   END IF
END SUB

SUB TerminateSim
   'Prints final status, empties some variables

   IF store = 2 AND dataentry <= maxentry THEN
      ON ERROR GOTO FileNotFound
      KILL "" + simname + "I" + ".SET"
      KILL "" + simname + ".INT"
      ON ERROR GOTO 0
   END IF

   LOCATE 1, 30
   IF store = 2 AND autoterminate = 1 THEN PRINT "Simulation saved as "; filename; "at"; displaytime
   IF store = 2 AND autoterminate = 0 THEN PRINT "Simulation "; simname; "aborted by user at"; displaytime
   IF store <> 2 AND autoterminate = 1 THEN PRINT "Simulation terminated at end of data at"; displaytime
   IF store <> 2 AND autoterminate = 0 THEN PRINT "Simulation terminated by user at"; displaytime
  
   FOR particle = 1 TO numparticles
      FOR dataentry = 0 TO maxentry
         sx(particle, dataentry) = 0
         sy(particle, dataentry) = 0
         sz(particle, dataentry) = 0
      NEXT dataentry
      rx(particle) = 0
      ry(particle) = 0
      rz(particle) = 0
      x(particle) = 0
      y(particle) = 0
      z(particle) = 0
      vx(particle) = 0
      vy(particle) = 0
      vz(particle) = 0
      xg(particle) = 0
      yg(particle) = 0
   NEXT particle

   DO
   LOOP WHILE INKEY$ = ""
END SUB

SUB TrackParticle
   'Adjusts screen coordinates to follow particle or point

   IF track > 0 AND track = particle THEN
      centx = x(track)
      centy = y(track)
      centz = z(track)
   END IF
   IF track = -1 AND particle = 1 THEN
      centx = (x(1) * m(1) + x(2) * m(2)) / (m(1) + m(2))
      centy = (y(1) * m(1) + y(2) * m(2)) / (m(1) + m(2))
      centz = (z(1) * m(1) + z(2) * m(2)) / (m(1) + m(2))
   END IF
   IF track <> 0 THEN offx = -centx
   IF track <> 0 THEN offy = -centy
END SUB

SUB UserInput
   'Detects user input and routes program to requested subroutines

   inkeyv = INKEY$
   IF inkeyv = CHR$(27) THEN exitflag = 1
   IF inkeyv = CHR$(32) AND pauseflag = 0 OR pauseflag = 2 THEN CALL Pause
   IF inkeyv = CHR$(47) OR inkeyv = CHR$(63) THEN display = 1: CALL Pause
   IF inkeyv = CHR$(67) OR inkeyv = CHR$(99) THEN CLS
   IF inkeyv = CHR$(22) OR inkeyv = CHR$(45) THEN spdchg = -1: CALL ChgSpd
   IF inkeyv = CHR$(43) OR inkeyv = CHR$(61) THEN spdchg = 1: CALL ChgSpd
   IF inkeyv = CHR$(50) OR inkeyv = CHR$(52) OR inkeyv = CHR$(54) OR inkeyv = CHR$(56) THEN CALL ImageSlide
   IF inkeyv = CHR$(84) OR inkeyv = CHR$(116) THEN IF tracers = 0 THEN tracers = 1 ELSE tracers = 0
   IF inkeyv = CHR$(44) OR inkeyv = CHR$(60) THEN optioncode = 1: CALL ImageRotate
   IF inkeyv = CHR$(46) OR inkeyv = CHR$(62) THEN optioncode = 2: CALL ImageRotate
   IF inkeyv = CHR$(76) OR inkeyv = CHR$(108) THEN optioncode = 3: CALL ImageRotate
   IF inkeyv = CHR$(59) OR inkeyv = CHR$(58) THEN optioncode = 4: CALL ImageRotate
   IF inkeyv = CHR$(80) OR inkeyv = CHR$(112) THEN optioncode = 5: CALL ImageRotate
   IF inkeyv = CHR$(91) OR inkeyv = CHR$(123) THEN optioncode = 6: CALL ImageRotate
   IF inkeyv = CHR$(73) OR inkeyv = CHR$(105) OR inkeyv = CHR$(79) OR inkeyv = CHR$(111) THEN CALL ImageZoom
   IF inkeyv = CHR$(83) OR inkeyv = CHR$(115) THEN CALL Snapshot
   IF inkeyv = CHR$(78) OR inkeyv = CHR$(110) OR inkeyv = CHR$(77) OR inkeyv = CHR$(109) THEN CALL AdjustColors
   IF inkeyv = CHR$(74) OR inkeyv = CHR$(106) OR inkeyv = CHR$(75) OR inkeyv = CHR$(107) THEN CALL AdjustColors
END SUB