DEFDBL A-Z
 pi = 3.14159: pi2 = pi * 2: conv = 180 / pi
 red = 12: yellow = 14: green = 10
 SCREEN 12: WIDTH , 60: col = 15
 emr = 81.3: em = 384235  'e/m mass ratio & distance
 ce = em / (1 + emr): cm = em - ce 'e/centre/m dists
 gce = 5.17234 * 10 ^ 12        'gravity const earth
 gcm = gce / emr     'grav. const moon,6.36204*10^10
 gme = gce / em ^ 2            'moon accel. to earth
 wm = SQR(gme / cm)           'moon angular velocity
 xm = cm: ym = 0: xe = -ce: ye = 0'moon/earth coords
 a = 1.057785: ca = COS(a): sa = SIN(a)'L4 lead angle
 r4 = 381922.1                            'L4 radius
 v4 = wm * r4    'orbital velocity at L4 equilibrium
 'note: distances L4 to earth and moon both equal em
 'total energy of body at L4 equilibrium
 enL4 = v4 * v4 / 2 - gce / em - gcm / em 'energy L4
 pi = 3.14159: pi2 = pi * 2: conv = 180 / pi
 CLS : WINDOW (310000, -15000)-(410000, 60000)
 mp = 360 / (conv * wm)'moon period, sidereal, hours
      'Body launched at centre of moon face to earth
 CIRCLE (cm, 0), 1738, col                     'moon
 x = cm - 1740: y = 100: u = 0: v = wm * x + 7700: rd = x
 ti = .0001: t = 0: angm = 0: stage = 1: miny = 99
 LOCATE 2: PRINT " STAGE   FUEL  RADIUS  ENERGY  TIME"
 LOCATE 4: PRINT " Start     0  "; FIX(rd)
 LOCATE 44, 60: PRINT "1": LOCATE 46, 50: PRINT "MOON"
 FOR yy = -15000 TO 15000 STEP 150  '...contour plot
  FOR xx = 323000 TO 390000 STEP 150'...contour plot
   xe = xx + ce: de = SQR(yy * yy + xe * xe)
   xm = cm - xx: dm = SQR(yy * yy + xm * xm)
   IF dm > 12000 THEN                 'not near moon
    pe = (gce / de + gcm / dm) / 10 ^ 5  'PE deficit
    PSET (xx, yy), pe MOD 16             'plot pixel
   END IF
 NEXT: NEXT
 LINE (323000, -15000)-(390000, 15000), col, B'inset
DO              'main calculation loop
 angm = t * wm: c = COS(angm): s = SIN(angm)
 xm = cm * c: ym = cm * s: xe = -ce * c: ye = -ce * s
 dmx = xm - x: dmy = ym - y   'body/moon separations
 dm2 = dmx * dmx + dmy * dmy: dm = SQR(dm2)'distance
 gm = gcm / dm2           'acceleration towards moon
 mx = gm * dmx / dm: my = gm * dmy / dm  'components
 dex = xe - x: dey = ye - y  'body/earth separations
 de2 = dex * dex + dey * dey: de = SQR(de2)'distance
 eg = gce / de2          'acceleration towards earth
 ex = eg * dex / de: ey = eg * dey / de  'components
 gx = ex + mx: gy = ey + my   'total body acc. comps
 ang = ATN(y / x)               'orbit angle of body
 IF x < 0 THEN ang = ang + pi 'put in right quadrant
 IF x > 0 AND y < 0 THEN ang = ang + pi2      'ditto
 at = ang - angm                'angle ahead of moon
 IF at < 0 THEN at = at + pi2  'moon at 300-360 deg.
 xt = rd * COS(at): yt = rd * SIN(at)'relative coords
 PSET (xt, yt), col     'plot body relative position
 SELECT CASE stage   'actions specific to each stage
 CASE 1          'partial orbit of moon after launch
  IF yt < miny THEN 'stage 2 burn on far side of moon
   LOCATE 6: PRINT " End 1        "; FIX(rd)
   stage = 2: fuel = 0: ti = .001: col = red'red plot
   COLOR col: LOCATE 51, 63: PRINT "2"
   accx = 1000: accy = 3000       'thrust components
   acct = SQR(accx * accx + accy * accy)     'thrust
   burn = 2800 * acct / (accx + accy): accti = acct * ti
  'PRINT burn: DO: LOOP WHILE INKEY$ = ""
  END IF
 CASE 2 'far side burn
  IF fuel < burn THEN                   'during burn
   gx = gx - accx: gy = gy - accy: fuel = fuel + accti
  ELSE       'burn finished - start stage 3 coasting
   LOCATE 8: PRINT " End 2  "; INT(fuel); FIX(rd)
   stage = 3: ytc = 15000: col = 15
   COLOR col: LOCATE 54, 55: PRINT "3"
  END IF
 CASE 3 'coasting (in inset)
  IF yt > ytc THEN             'reached top of inset
   en = FIX((u * u + v * v) / 2 - gce / de - gcm / dm)
LOCATE10: PRINT " End 3  "; INT(fuel); FIX(rd); en; INT(t)
   LOCATE 41, 74: PRINT "MOON"
   stage = 4: WINDOW (-400000, -200000)-(400000, 400000)
   LINE (323000, -15000)-(390000, 15000), col, B'inset
   s60 = SIN(60 / conv)          'for 60 deg. radial
   CIRCLE (0, 0), cm, 15'moon orbit, white, & radial
   LINE (150000, 300000 * s60)-(cm / 2, cm * s60), 15
   CIRCLE (0, 0), r4, green'orbit of L4, green,& radial
   LINE (300000 * ca, 300000 * sa)-(r4 * ca, r4 * sa), green
   ti = .001: LOCATE 36, 71: PRINT "4": LOCATE 22, 65: PRINT "4"
  END IF
 CASE 4 'coasting (full screen)
  PSET (x, y), col      'plot body absolute position
  IF INT(t * 10) = 1080 THEN rdo = rd     '@ 108 hrs
  IF t > 120 THEN        'go to stage 5 at 120 hours
   en = FIX((u * u + v * v) / 2 - gce / de - gcm / dm)'TE
   LOCATE 16: PRINT " At L4 need   "; FIX(r4); FIX(enL4)
LOCATE12: PRINT " End 4  "; INT(fuel); FIX(rd); en; FIX(t)
   stage = 5: it = 9: ti = .004: col = red: COLOR col
   LOCATE 22, 71: PRINT "5": LOCATE 6, 37: PRINT "5"
   vr = .0357: va = .0004
  END IF
 CASE 5    'control approach to L4
  IF angm < pi THEN PSET (x, y), col  'abs. position
  IF INT(t / 12) > it THEN            'each 12 hours
   it = INT(t / 12)
   en = FIX((u * u + v * v) / 2 - gce / de - gcm / dm)'TE
   dvr = -vr * (rd - rdo)    'radial velocity adjust
   dva = va * (enL4 - en)   'orbital velocity adjust
   dv = SQR(dva * dva + dvr * dvr)'total velocity adj
   fuel = fuel + dv 'fuel used. Catch-up speed on L4
   vo = FIX(v * COS(ang) - u * SIN(ang) - v4)
'-------------------------------------
 'Gary, 110 lines of 52 characters to here.
  ' Going to 44 lines of 70 characters
'-------------------------------------
   LOCATE 14: PRINT " 5      "; INT(fuel); FIX(rd); en; INT(t); vo; "Delta v"
   vx = dvr * COS(ang) - dva * SIN(ang)                   'x component
   vy = dvr * SIN(ang) + dva * COS(ang)                   'y component
   u = u + vx: v = v + vy: rdo = rd
   IF angm < pi THEN
    LINE (x, y)-STEP(vx * 100, vy * 100), yellow           'thrustline
   END IF
    vxt = dvr * COS(at) - dva * SIN(at)                  'xt component
    vyt = dvr * SIN(at) + dva * COS(at)                  'yt component
    LINE (xt, yt)-STEP(vxt * 100, vyt * 100), yellow       'thrustline
  END IF
  IF at > a THEN  'reaches L4 radial
   LOCATE 14: PRINT " End 5  "; INT(fuel); FIX(rd); en; INT(t); vo; "Delta v"
   stage = 6: fgr = .002: fga = .01
   ti = .05: col = yellow: COLOR col
   LOCATE 20, 52: PRINT "L4": LOCATE 7, 61: PRINT "L4"
   LINE (r4 * ca - 2800, r4 * sa - 2000)-(r4 * ca + 2000, r4 * sa + 2100), green, B
   WINDOW (r4 * ca - 9000, r4 * sa - 7000)-(r4 * ca + 7000, r4 * sa + 5000)
   LINE ((r4 - 2000) * ca, (r4 - 2000) * sa)-(r4 * ca, r4 * sa), green
   LINE (r4 * ca - 3000 * sa, r4 * sa + 3000 * ca)-(r4 * ca + 2000 * sa, r4 * sa - 2000 * ca), green
   LINE (r4 * ca - 2800, r4 * sa - 2000)-(r4 * ca + 2000, r4 * sa + 2100), green, B
  END IF
 CASE 6 'final settling on L4
  en = FIX((u * u + v * v) / 2 - gce / de - gcm / dm)    'total energy
  gr = -fgr * (rd - rdo) / ti  'radial control acc.
  ga = (enL4 - en) * fga     'orbital control acc.
  gt = SQR(ga * ga + gr * gr)              'total control acceleration
  fuel = fuel + gt * ti: rdo = rd                           'fuel used
  LOCATE 18: PRINT " 6      "; INT(fuel); FIX(rd); en; INT(t)
  gx = gx + gr * COS(ang) - ga * SIN(ang)                 'x component
  gy = gy + gr * SIN(ang) + ga * COS(ang)                 'y component
 CASE ELSE: EXIT DO
 END SELECT
 u = u + ti * gx: v = v + ti * gy'new vel. components
 x = x + ti * u: y = y + ti * v      'new coordinates
 rd2 = x * x + y * y: rd = SQR(rd2)       'new radius
 t = t + ti          'increment time by time interval
 IF INKEY$ <> "" THEN
  DO: k$ = INKEY$: LOOP WHILE k$ = ""
  IF k$ = CHR$(27) THEN EXIT DO
 END IF
LOOP