'Code Listing 1	 
 DEFDBL A-Z   'Use double precision, all variables
 pi = 3.14159: pi2 = pi * 2: conv = 180 / pi
 limit = 50000: col = 15: col2 = 13
 emr = 81.3: em =384235 'e/m mass ratio & distance
 ce = em / 82.3: cm = em - ce  'earth/c/moon dists
 gce = 5.17234 * 10 ^ 12        'grav. const earth
 gcm = gce / emr 'grav. const moon, 6.362041*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  'e/m positions
 r1 = 321571.4: r2 = 444053.6: r3 = 386180.3
 SCREEN 12:WINDOW(-450000,-375000)-(550000,375000)
 LOCATE 16, 62:PRINT "L1": LOCATE 16,72:PRINT "L2"
 LOCATE 15, 36: PRINT "E": LOCATE 16,37: PRINT "C"
 LOCATE 16, 5: PRINT "L3": LOCATE 15,67: PRINT "M"
 LOCATE 3, 51: PRINT "L4": LOCATE 28,51:PRINT "L5"
 LOCATE 5, 59: PRINT "Moon orbit No 0"
 LOCATE 6, 59: PRINT "Moon angle"
 PSET (-ce, 0), 15: PSET (cm,0), 15: PSET (0,0),15
 rd = r1 - 10 : x = rd            '10 km inside L1
 'rd = r2 - 10 : x = rd           '10 km inside L2
 'rd = r3 + 100 : x = -rd       '100 km outside L3
 y = 0: u = 0: v = x * wm     'position & velocity
 ti = .05: ti2 = .005: orbit = 0: dm = 60000
DO
 IF dm < limit THEN dt = ti2 ELSE dt = ti
 angm = angm + dt * wm:c = COS(angm):s = SIN(angm)
 cang = CINT(angm * conv): LOCATE 6, 72:PRINT cang
 IF orbit = 4 AND cang = 130 THEN
   col = 11: col2 = 14        'change plot colours
 END IF           'meant for use with L2 body only
 xm = cm *c: ym = cm *s: xe = -ce *c: ye = -ce * s
 dmx = xm - x: dmy = ym - y      'coordinate diffs
 dm2 = dmx * dmx + dmy * dmy: dm = SQR(dm2)  'sep.
 gm = gcm /dm2: mx = gm *dmx /dm: my =gm *dmy / dm
 dex = xe - x: dey = ye - y      'coordinate diffs
 de2 = dex * dex + dey * dey: de = SQR(de2)
 eg = gce /de2: ex = eg *dex /de: ey =eg *dey / de
 gx = ex + mx: gy = ey + my
 ang = ATN(y / x)   'arctan gives ambiguous result
 IF x < 0 THEN ang = ang + pi    'remove ambiguity
 at = ang -angm: xt =rd *COS(at): yt =rd * SIN(at)
 PSET (xt, yt), col       'plot relative and/or...
 PSET (x, y), col2        '...absolute positions
 u = u + dt * gx: v = v + dt * gy    'new velocity
 x = x + dt * u: y = y + dt * v      'new position
 rd2 = x * x + y * y: rd = SQR(rd2)  'new radius
 IF angm > pi2 THEN        'end of each moon orbit
  angm = 0: orbit = orbit + 1
  LOCATE 5, 72: PRINT orbit
 END IF
 IF INKEY$ <> "" THEN EXIT DO     'any key escapes
LOOP
DO: LOOP WHILE INKEY$ = ""             'look at it