'lag45.bas
 SCREEN 12: VIEW: VIEW PRINT: WIDTH , 60
 WINDOW (80000, 285000)-(260000, 420000)
 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: c = COS(a): s = SIN(a) 'L4 lead angle
 r4 = 381922.1                            'L4 radius
 v4 = wm * r4    'velocity of body at L4 equilibrium
 'note: distances L4 to earth and moon both equal em
 enL4 = v4 * v4 / 2 - gce / em - gcm / em 'energy L4
 CIRCLE (0, 0), cm, 15, .86, 1.25        'moon orbit
 CIRCLE (0, 0), r4, 11, .96, 1.15          'L4 orbit
 LINE (0, 0)-(r4 * c, r4 * s), 11         'L4 radial
 LINE (0, 0)-(189865, 328855), 15    '60 deg. radial
 LINE (184097, 332119)-(178571, 322156), 15 '61 deg.
 LINE (178272, 335282)-(172924, 325223), 15 '62 deg.
 LINE (129000, 322000)-STEP(10000, 10000), 15, B
 LOCATE 3, 15: PRINT "LEAD ANGLE": LOCATE 10, 4
 PRINT "0": LOCATE 12,34:PRINT "ORBIT": LOCATE 19,13
 PRINT "RADIUS": LOCATE 27, 3: PRINT "ORBIT  FUEL "
 LOCATE 46, 19: PRINT "10,000 km sq"
 LOCATE 26, 36: PRINT "Start at L4, +2000 km radial"
 'Body 2000 km outside L4 with moon angular velocity
 rd = r4 + 2000: x = rd * c: y = rd * s 'body rd/x/y
 vel = wm * rd: u = -vel * s: v = vel * c'vel. comps
 col = 10: colf = 12
'col = 13: colf = 11               'plotting colours
 xp = 85000: fam = 50000: faa = 355000
 fa = 5 * 10 ^ 4: fa2 = 407890
 pi = 3.14159: pi2 = 6.28318
 fgr = .0018: fga = 6 * 10 ^ -7: sx10 = 140
 ti = .02: t = 0: t10 = 0: orbit = 0
 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 accel comps
  IF cont THEN                     'if under control
  'radial control acc. - for radial velocity damping
   drd = rd - rdo: gr = -fgr * drd / ti
   LINE (xp, r4)-STEP(0, gr * fa), colf 'plot thrust
  'control along orbit - correcting total energy
   ga = fga *(enL4 - (u*u + v*v)/2 +gce/de + gcm/dm)
   LINE (xp, fa2)-STEP(0, -ga*fa), colf 'plot thrust
   gt = SQR(ga * ga + gr * gr) 'total control accel.
   fuel = fuel + ABS(gt * ti)             'fuel used
   gx = gx + gr * COS(ang) - ga * SIN(ang)   'rel. x
   gy = gy + gr * SIN(ang) + ga * COS(ang)   'rel. y
  END IF
  rdo = rd                        'record old radius
  ang = ATN(y / x)              'orbit angle of body
  IF x < 0 THEN ang = ang + pi'put in right quadrant
  at = ang - angm      'body/moon angular separation
  xt = rd *COS(at): yt = rd *SIN(at) 'posn rel. moon
  IF at < 0 THEN at = at + pi2  'moon at 300-360 deg
  PSET (xt, yt), col    'plot body relative position
  PSET (xp, rd), col            'plot radius vs time
  PSET (xp, at * fam + faa), col 'plot angle vs time
  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: t10 = t10 + ti        'increment times
  IF t10 >10 THEN xp =xp +sx10 'incr. screen x coord
  IF t10 > 10 THEN t10 = ti               'reset t10
  IF t > 654 THEN                        'next orbit
   t = ti: orbit = orbit + 1  'reset t for new orbit
   LOCATE 10, (orbit * 4 + 3) MOD 80: PRINT orbit
   IF cont THEN              'print fuel consumption
   LOCATE orbit*2+8, 3:PRINT orbit-1;"   ";INT(fuel)
   END IF
   IF orbit = 10 THEN cont = 1    'switch on control
   IF orbit = 10 THEN col = 15   'change plot colour
   LINE (xp, 367000)-(xp, 417000), 15 'vertical line
  END IF
 LOOP UNTIL orbit = 18
 DO:  LOOP WHILE INKEY$ = ""