'APOLLO7.BAS          'from moon orbit to earth atmosphere
SCREEN 12: WINDOW (-10000, -25000)-(390000, 275000)   'VGA
pi2 = 8 * ATN(1)                     '2pi = 8 * 45 degrees
emmr = 81.3:emd = 384235 'earth/moon mass ratio & distance
ce = emd / (1 + emmr): cm = emd - ce 'e/centre/m distances
Erad = 6378: gce = .0098 * Erad * Erad 'earth & moon radii
Mrad = 1738: gcm = gce / emmr       '...and gravity consts
mpx = 350000: mpy = 210000: mmag = 8       'magnified moon
CIRCLE (mpx, mpy), Mrad * mmag, 10
dx = Mrad * mmag * COS(.1): dy = Mrad * mmag * SIN(.1)
LINE (mpx - dx, mpy - dy)-(mpx + dx, mpy + dy), 10
PAINT (mpx, mpy - 1000), 10, 10         'sunlit hemisphere
epx = 45000: epy = 90000: emag = 8        'magnified earth
CIRCLE (epx, epy), Erad * emag, 9
dx = Erad * emag * COS(.12): dy = Erad * emag * SIN(.12)
LINE (epx - dx, epy - dy)-(epx + dx, epy + dy), 9
PAINT (epx, epy - 1000), 9, 9           'sunlit hemisphere
LINE (epx - dx, epy - dy)-STEP(30000, 0),15   'splash-down
angm0 = 0: angm = 0: c = 1: s = 0  'angle of moon at start
dangm = pi2 /(654.2 *3600) 'moon angular vel., radians/sec
mv = dangm * cm  'moon velocity in orbit = 1.012639 km/sec
dm1 = Mrad + 112: mg = gcm / (dm1 * dm1)  'orbit gravity &
cv = SQR(mg * dm1): period = pi2 * dm1 / (cv * 60)    'vel
'PRINT cv; period            '1.628045 km/s  118.9964 secs
sector = 1: t = 0     'start at 0 sec (on x axis, 119 hrs)
ts = 57120: ti = 1'start TLI burn @ ts secs, time interval
xe = -ce: ye = 0: xm = cm: ym = 0  'earth & moon positions
x = xm + dm1: y = 0: u = 0: v = mv -cv 'initial pos. & vel
mc = 17112: tl = 9318: thr = tl * .0098     'mass & thrust
SIl = 310: ff = tl / SIl   'specific impulse and fuel flow
LOCATE 28, 50: PRINT "(a)Moon & CM at 119 hrs  (a)"
DO: t = t + ti              'THIS IS THE MAIN PROGRAM LOOP
               'FIRST WE HAVE THE CODE USED IN ALL SECTORS
 angm = angm0 + dangm * t: c = COS(angm): s = SIN(angm)
 xm = cm * c: ym = cm * s: xe = -ce * c: ye = -ce * s
 dmx = xm - x: dmy = ym - y 'Apollo/moon x & y separations
 dm2 = dmx * dmx + dmy * dmy: dm = SQR(dm2) '  "  distance
 gm = gcm / dm2                 'acceleration towards moon
 mx = gm * dmx / dm: my = gm * dmy / dm 'accel. components
 dex = xe - x: dey = ye - y'Apollo/earth x & y separations
 de2 = dex * dex + dey * dey: de = SQR(de2)  ' "  distance
 eg = gce / de2                'acceleration towards earth
 ex = eg * dex / de: ey = eg * dey / de 'accel. components
 gx = ex + mx + atx: gy = ey + my + aty  'total acc. comps
 u = u + ti * gx: v = v + ti * gy      'increment velocity
 x = x + ti * u: y = y + ti * v        'increment position
 w2 = u * u + v * v: w = SQR(w2)    'magnitude of velocity
 IF t MOD 60 = 0 THEN       'each minute adjust min and hr
  tm = t / 60: hr = INT(tm / 60): min = tm MOD 60
 END IF
 PSET (x, y), 15  'plot Apollo position (bottom of screen)
 SELECT CASE sector'FOLLOWING CODE SPECIFIC TO EACH SECTOR
  CASE 1:    'COMMAND MODULE/SERVICE MODULE ORBITTING MOON
   IF t = ts THEN      'start TEI burn to escape from moon
    vrelmr=2.63:sector=2:COLOR 12:LOCATE 7,76: PRINT "(b)"
    LOCATE 20, 47: PRINT "(b)Trans Earth Injection burn"
   END IF         'Sets velocity relative to moon required
  CASE 2:       'TRANS EARTH INJECTION BURN
   um = -mv * s: vm = mv * c  'components of moon velocity
   vrelm = SQR((u - um) * (u - um) + (v - vm) * (v - vm))
   IF vrelm < vrelmr THEN 'until reaches required velocity
    acc = thr / (mc - ff / 2): mc = mc - ff 'accel. & mass
    ang = ATN(dmy / dmx)  'thrust parallel to moon surface
    atx = acc * SIN(ang): aty = -acc * COS(ang)'components
    LINE (mpx,mpy)-STEP((x -xm) *mmag, (y -ym) * mmag), 12
   ELSE                   'when exceeded required velocity
    sector = 3: atx = 0: aty = 0: f = CINT((t - ts) * ff)
    LOCATE 21,50: PRINT t-ts;"secs,"; f;"kg fuel"; "  (b)"
    angm0 = angm: t = 0              'reset time and angm0
   END IF
  CASE 3:           'COASTING OUT FROM MOON AFTER TEI BURN
   IF t MOD 60 = 0 THEN  'each minute, plot posn near moon
    PSET (mpx + (x - xm) * mmag, mpy + (y - ym) * mmag),15
    IF t = 7200 THEN         'end Sector 3 at 2 hrs and...
     LOCATE 22, 25: COLOR 15
     PRINT"50            40         30        20       10"
     sector = 4: ti = 5  '...change to 5 sec time interval
    END IF
   END IF
  CASE 4:     'TRIP TO EARTH (until final correction burn)
   IF tm MOD 60 = 0 AND min = 0 THEN           'each hour,
    LINE (x, y)-STEP(0, 3000), 15   '1 & 10 hr time scales
    IF hr MOD 10 = 0 THEN LINE (x, y)-STEP(0, 8000), 15
   END IF
   IF hr > 54 THEN   'position plot near (magnified) earth
    PSET (epx + (x - xe) * emag, epy + (y - ye) * emag),15
    IF de < 24000 THEN        'correction burn at 24000 km
     LOCATE 3, 40: COLOR 12: PRINT "(c)Correction burn"
     LOCATE 26, 12: PRINT "(c)Correction burn": COLOR 15
     LINE -STEP(30000, 0),15:LINE (x, y)-STEP(30000, 0),15
     sector = 5: v = v + .143 'adjust y velocity component
     td = t + 1800         'set up time to mark SM discard
    END IF      'end IF de < 24000
   END IF       'end IF hr > 55
  CASE 5:                     'FROM CORRECTION TO RE-ENTRY
    PSET (epx + (x - xe) * emag, epy + (y - ye) * emag),15
    IF t = td THEN
     LINE -STEP(30000, 0), 15
     LOCATE 8, 25: PRINT "(d)Service Module discarded"
    END IF
    IF de < 6490 THEN     'at 112 km altitude, end program
     LINE -STEP(30000, 0),15: LINE (x,y)-STEP(30000, 0),15
     CIRCLE (xe, ye), Erad, 9:finish = 1 'draw small earth
     ange = ATN((y - ye) / (xe - x)): angv = ATN(u / v)
     reang = (ange - angv) * 57.3: LOCATE 18, 10
     PRINT "(e)Re-entry angle";CINT(reang); "deg @ 112 km"
     LOCATE 28, 9: PRINT "(e)Re-entry"
     LOCATE 20, 10: PRINT "(f)Splash-down at dawn"
    END IF
 END SELECT
LOOP UNTIL finish = 1
DO: LOOP WHILE INKEY$ = ""            'Well! Look at that.