  ' NPOE.BAS               May 23, 2017
  
  ' Numerical Prediction of Orbital Events
  
  ' includes perturbations due to
  
  '   non-spherical earth gravity
  '   aerodynamic drag (us 76 model)
  '   solar radiation pressure
  '   point mass sun and moon
  
  ' MicroMite BASIC eXtreme version
  
  '''''''''''''''''''''''''''''''''
  
  option default float
  
  option base 1
  
  dim month$(12) as string
  
  dim xnut(11, 13)
  
  ' sun ephemeris
  
  dim xsl(50), xsr(50), xsa(50), xsb(50)
  
  ' RKF78 integration constants
  
  dim ch(13), alph(13), beta(13, 12)
  
  ' gravity model coefficients
  
  dim scoef(10, 10), ccoef(10, 10)
  
  ' state vector and orbital elements
  
  dim ri(3), vi(3), rf(3), vf(3)
  
  dim yold(6), ynew(6), oevi(32), oevf(32)
  
  ' initial gst, julian day and spacecraft properties
  
  dim gst0, jdate0, csrp0, bcoeff
  
  ' perturbation "flags"
  
  dim isun%, imoon%, idrag%, isrp%
  
  ' algorithm tolerances and working arrays
  
  dim tetol, rtol, tiwrk, yi(6), yfinal(6)
  
  ' orbital "targets"
  
  dim ttanom, talt, tdec, tfpa, targlat, tlong
  
  dim tlat, tspeed, trasc, ittype%
  
  ' astrodynamic and utility constants
  
  const dtr = pi / 180.0
  
  const rtd = 180.0 / pi
  
  const pi2 = 2.0 * pi
  
  const pidiv2 = pi / 2.0
  
  const atr = pi / 648000.0
  
  const seccon = 206264.8062470964
  
  ' earth gravitational constant (km**3/sec**2)
  
  const emu = 398600.436233
  
  ' moon gravitational constant (km**3/sec**2)
  
  const mmu = 4902.800076
  
  ' sun gravitational constant (km**3/sec**2)
  
  const smu = 132712440040.944
  
  ' earth inertial rotation rate (radians/second)
  
  const omega = 7.292115486e-5
  
  ' earth equatorial radius (kilometers)
  
  const req = 6378.1363
  
  ' earth flattening factor (non-dimensional)
  
  const flat = 1.0 / 298.257
  
  ' earth oblateness gravity coefficients (non-dimensional)
  
  const xj2 = 1.082626335439e-3
  
  const xj4 = -1.619974305782e-6
  
  ' astronomical unit (kilometers)
  
  const aunit = 149597870.691
  
  ' diameter of the sun (kilometers)
  
  const dsun = 696000.0
  
  ' diameter of the moon (kilometers)
  
  const dmoon = 1738.0
  
  ' solar constant
  
  const pss = 0.00456
  
  '''''''''''''''''''''''''''
  ' read solar ephemeris data
  '''''''''''''''''''''''''''
  
  for i% = 1 to 50
    
    read xsl(i%), xsr(i%), xsa(i%), xsb(i%)
    
  next i%
  
  data 403406,      0, 4.721964,      1.621043
  data 195207, -97597, 5.937458,  62830.348067
  data 119433, -59715, 1.115589,  62830.821524
  data 112392, -56188, 5.781616,  62829.634302
  data   3891,  -1556, 5.5474  , 125660.5691
  data   2819,  -1126, 1.5120  , 125660.9845
  data   1721,   -861, 4.1897  ,  62832.4766
  data      0,    941, 1.163   ,      0.813
  data    660,   -264, 5.415   , 125659.310
  data    350,   -163, 4.315   ,  57533.850
  data    334,      0, 4.553   ,    -33.931
  data    314,    309, 5.198   , 777137.715
  data    268,   -158, 5.989   ,  78604.191
  data    242,      0, 2.911   ,      5.412
  data    234,    -54, 1.423   ,  39302.098
  data    158,      0, 0.061   ,    -34.861
  data    132,    -93, 2.317   , 115067.698
  data    129,    -20, 3.193   ,  15774.337
  data    114,      0, 2.828   ,   5296.670
  data     99,    -47, 0.52    ,  58849.27
  data     93,      0, 4.65    ,   5296.11
  data     86,      0, 4.35    ,  -3980.70
  data     78,    -33, 2.75    ,  52237.69
  data     72,    -32, 4.50    ,  55076.47
  data     68,      0, 3.23    ,    261.08
  data     64,    -10, 1.22    ,  15773.85
  data     46,    -16, 0.14    , 188491.03
  data     38,      0, 3.44    ,  -7756.55
  data     37,      0, 4.37    ,    264.89
  data     32,    -24, 1.14    , 117906.27
  data     29,    -13, 2.84    ,  55075.75
  data     28,      0, 5.96    ,  -7961.39
  data     27,     -9, 5.09    , 188489.81
  data     27,      0, 1.72    ,   2132.19
  data     25,    -17, 2.56    , 109771.03
  data     24,    -11, 1.92    ,  54868.56
  data     21,      0, 0.09    ,  25443.93
  data     21,     31, 5.98    , -55731.43
  data     20,    -10, 4.03    ,  60697.74
  data     18,      0, 4.27    ,   2132.79
  data     17,    -12, 0.79    , 109771.63
  data     14,      0, 4.24    ,  -7752.82
  data     13,     -5, 2.01    , 188491.91
  data     13,      0, 2.65    ,    207.81
  data     13,      0, 4.98    ,  29424.63
  data     12,      0, 0.93    ,     -7.99
  data     10,      0, 2.21    ,  46941.14
  data     10,      0, 3.59    ,    -68.29
  data     10,      0, 1.50    ,  21463.25
  data     10,     -9, 2.55    , 157208.40
  
  ''''''''''''''''''''''''''''''''''''''
  ' read subset of IAU2000 nutation data
  ''''''''''''''''''''''''''''''''''''''
  
  for j% = 1 to 13
    
    for i% = 1 to 11
      
      read xnut(i%, j%)
      
    next i%
    
  next j%
  
  data  0,  0, 0,  0, 1, -172064161, -174666, 92052331,  9086,  33386, 15377
  data  0,  0, 2, -2, 2,  -13170906,   -1675,  5730336, -3015, -13696, -4587
  data  0,  0, 2,  0, 2,   -2276413,    -234,   978459,  -485,   2796,  1374
  data  0,  0, 0,  0, 2,    2074554,     207,  -897492,   470,   -698,  -291
  data  0,  1, 0,  0, 0,    1475877,   -3633,    73871,  -184,  11817, -1924
  data  0,  1, 2, -2, 2,    -516821,    1226,   224386,  -677,   -524,  -174
  data  1,  0, 0,  0, 0,     711159,      73,    -6750,     0,   -872,   358
  data  0,  0, 2,  0, 1,    -387298,    -367,   200728,    18,    380,   318
  data  1,  0, 2,  0, 2,    -301461,     -36,   129025,   -63,    816,   367
  data  0, -1, 2, -2, 2,     215829,    -494,   -95929,   299,    111,   132
  data  0,  0, 2, -2, 1,     128227,     137,   -68982,    -9,    181,    39
  data -1,  0, 2,  0, 2,     123457,      11,   -53311,    32,     19,    -4
  data -1,  0, 0,  2, 0,     156994,      10,    -1235,     0,   -168,    82
  
  '''''''''''''''''
  ' calendar months
  '''''''''''''''''
  
  month$(1) = "January"
  month$(2) = "February"
  month$(3) = "March"
  month$(4) = "April"
  month$(5) = "May"
  month$(6) = "June"
  month$(7) = "July"
  month$(8) = "August"
  month$(9) = "September"
  month$(10) = "October"
  month$(11) = "November"
  month$(12) = "December"
  
  '''''''''''''''''''''''''''''''
  ' US76 standard atmosphere data
  '''''''''''''''''''''''''''''''
  
  ' dimension atmospheric data arrays
  
  dim zs(49) as float, ps(49) as float, tms(49) as float, wms(49) as float
  
  ' read atmospheric data
  
  for i% = 1 to 49
    
    read zs(i%)
    
  next i%
  
  data 0.0, 11.019, 20.063, 32.162, 47.35, 51.413
  data 71.802, 86.0, 91.0, 94.0, 97.0, 100.0, 103.0
  data 106.0, 108.0, 110.0, 112.0, 115.0, 120.0, 125.0
  data 130.0, 135.0, 140.0, 145.0, 150.0, 155.0, 160.0
  data 165.0, 170.0, 180.0, 190.0, 210.0, 230.0, 265.0
  data 300.0, 350.0, 400.0, 450.0, 500.0, 550.0, 600.0
  data 650.0, 700.0, 750.0, 800.0, 850.0, 900.0, 950.0, 1000.0
  
  for i% = 1 to 49
    
    read tms(i%)
    
  next i%
  
  data 288.15,  216.65,  216.65,  228.65, 270.65
  data 270.65,  214.65,  186.95,  186.87, 187.74
  data 190.40,  195.08,  202.23,  212.89, 223.29
  data 240.00,  264.00,  300.00,  360.00, 417.23
  data 469.27,  516.59,  559.63,  598.78, 634.39
  data 666.80,  696.29,  723.13,  747.57, 790.07
  data 825.31,  878.84,  915.78,  955.20, 976.01
  data 990.06,  995.83,  998.22,  999.24, 999.67
  data 999.85,  999.93,  999.97,  999.99, 999.99
  data 1000.00, 1000.00, 1000.00, 1000.00
  
  for i% = 1 to 49
    
    read wms(i%)
    
  next i%
  
  data 28.9644, 28.9644, 28.9644, 28.9644, 28.9644
  data 28.9644, 28.9644, 28.9522, 28.889,  28.783
  data 28.620,  28.395,  28.104,  27.765,  27.521
  data 27.268,  27.020,  26.680,  26.205,  25.803
  data 25.436,  25.087,  24.749,  24.422,  24.103
  data 23.792,  23.488,  23.192,  22.902,  22.342
  data 21.809,  20.825,  19.952,  18.688,  17.726
  data 16.735,  15.984,  15.247,  14.330,  13.092
  data 11.505,   9.718,   7.998,   6.579,   5.543
  data  4.849,   4.404,   4.122,   3.940
  
  for i% = 1 to 49
    
    read ps(i%)
    
  next i%
  
  data 1013.25,  226.32, 54.7487, 8.68014, 1.10905
  data 0.66938,  0.039564, 3.7338e-3, 1.5381e-3
  data 9.0560e-4, 5.3571e-4, 3.2011e-4, 1.9742e-4
  data 1.2454e-4, 9.3188e-5, 7.1042e-5, 5.5547e-5
  data 4.0096e-5, 2.5382e-5, 1.7354e-5, 1.25054e-5
  data 9.3568e-6, 7.2028e-6, 5.6691e-6, 4.5422e-6
  data 3.6930e-6, 3.0395e-6, 2.5278e-6, 2.1210e-6
  data 1.5271e-6, 1.1266e-6, 6.4756e-7, 3.9276e-7
  data 1.7874e-7, 8.7704e-8, 3.4498e-8, 1.4518e-8
  data 6.4468e-9, 3.0236e-9, 1.5137e-9, 8.2130e-10
  data 4.8865e-10, 3.1908e-10, 2.2599e-10, 1.7036e-10
  data 1.3415e-10, 1.0873e-10, 8.9816e-11, 7.5138e-11
  
  ''''''''''''''''''''''''''''''''''''''''''
  ' read gravity model coefficients (GGM02C)
  ''''''''''''''''''''''''''''''''''''''''''
  
  for ii% = 1 to 42
    
    read i%, j%, ctmp, stmp
    
    ccoef(i% + 1, j% + 1) = ctmp
    
    scoef(i% + 1, j% + 1) = stmp
    
  next ii%
  
  data   2,    0,   -1.082626335439e-3,    0.000000000000
  data   3,    0,   +2.532473691333e-6,    0.000000000000
  data   4,    0,   +1.619974305782e-6,    0.000000000000
  data   5,    0,   +2.279051260821e-7,    0.000000000000
  data   6,    0,   -5.406167899402e-7,    0.000000000000
  data   7,    0,   +3.505229256321e-7,    0.000000000000
  data   8,    0,   +2.040167670105e-7,    0.000000000000
  data   2,    1,   -2.641160102687e-10,  +1.803286269386e-9
  data   3,    1,   +2.193163948488e-6,   +2.680587395880e-7
  data   4,    1,   -5.086610622286e-7,   -4.492660737921e-7
  data   5,    1,   -5.386781134151e-8,   -8.081643446849e-8
  data   6,    1,   -5.972058174582e-8,   +2.086344028373e-8
  data   7,    1,   +2.055896831162e-7,   +6.962045068478e-8
  data   8,    1,   +1.591517378410e-8,   +4.047286188441e-8
  data   2,    2,   +1.574576419563e-6,   -9.038679457258e-7
  data   3,    2,   +3.090482638647e-7,   -2.114302585649e-7
  data   4,    2,   +7.837677629613e-8,   +1.481259637805e-7
  data   5,    2,   +1.055330562324e-7,   -5.232899486269e-8
  data   6,    2,   +6.054894249564e-9,   -4.650080282823e-8
  data   7,    2,   +3.291018225110e-8,   +9.261551569961e-9
  data   8,    2,   +6.571982917625e-9,   +5.361626853052e-9
  data   3,    3,   +1.005789794582e-7,   +1.972224716900e-7
  data   4,    3,   +5.921491430169e-8,   -1.201057108321e-8
  data   5,    3,   -1.492756134350e-8,   -7.102436843719e-9
  data   6,    3,   +1.186709381073e-9,   +1.852695542356e-10
  data   7,    3,   +3.527838137946e-9,   -3.058718640148e-9
  data   8,    3,   -1.959204907845e-10,  -8.692343892569e-10
  data   4,    4,   -3.982102837836e-9,   +6.525061656322e-9
  data   5,    4,   -2.299570414087e-9,   +3.878954367218e-10
  data   6,    4,   -3.256426896345e-10,  -1.784566345807e-9
  data   7,    4,   -5.839481895654e-10,  -2.634600857631e-10
  data   8,    4,   -3.189480413465e-10,  +9.111352478877e-11
  data   5,    5,   +4.303834339282e-10,  -1.648158423212e-9
  data   6,    5,   -2.156240650756e-10,  -4.329815578525e-10
  data   7,    5,   +5.875565799080e-13,  +6.347191438055e-12
  data   8,    5,   -4.650824883957e-12,  +1.614425047070e-11
  data   6,    6,   +2.205569478246e-12,  -5.531070854282e-11
  data   7,    6,   -2.490475620563e-11,  +1.053602131659e-11
  data   8,    6,   -1.842244217832e-12,  +8.628417788820e-12
  data   7,    7,   +2.794245686804e-14,  +4.473695242495e-13
  data   8,    7,   +3.429697468438e-13,  +3.817951257292e-13
  data   8,    8,   -1.581200671882e-13,  +1.536769701558e-13
  
  '''''''''''''''''''''''''''''''''''''''
  ' define rkf78 integration coefficients
  '''''''''''''''''''''''''''''''''''''''
  
  ch(6) = 34.0 / 105.0
  ch(7) = 9.0 / 35.0
  ch(8) = ch(7)
  ch(9) = 9.0 / 280.0
  ch(10) = ch(9)
  ch(12) = 41.0 / 840.0
  ch(13) = ch(12)
  
  alph(2) = 2.0 / 27.0
  alph(3) = 1.0 / 9.0
  alph(4) = 1.0 / 6.0
  alph(5) = 5.0 / 12.0
  alph(6) = 0.5
  alph(7) = 5.0 / 6.0
  alph(8) = 1.0 / 6.0
  alph(9) = 2.0 / 3.0
  alph(10) = 1.0 / 3.0
  alph(11) = 1.0
  alph(13) = 1.0
  
  beta(2, 1) = 2.0 / 27.0
  beta(3, 1) = 1.0 / 36.0
  beta(4, 1) = 1.0 / 24.0
  beta(5, 1) = 5.0 / 12.0
  beta(6, 1) = 0.05
  beta(7, 1) = -25.0 / 108.0
  beta(8, 1) = 31.0 / 300.0
  beta(9, 1) = 2.0
  beta(10, 1) = -91.0 / 108.0
  beta(11, 1) = 2383.0 / 4100.0
  beta(12, 1) = 3.0 / 205.0
  beta(13, 1) = -1777.0 / 4100.0
  beta(3, 2) = 1.0 / 12.0
  beta(4, 3) = 1.0 / 8.0
  beta(5, 3) = -25.0 / 16.0
  beta(5, 4) = -beta(5, 3)
  beta(6, 4) = 0.25
  beta(7, 4) = 125.0 / 108.0
  beta(9, 4) = -53.0 / 6.0
  beta(10, 4) = 23.0 / 108.0
  beta(11, 4) = -341.0 / 164.0
  beta(13, 4) = beta(11, 4)
  beta(6, 5) = 0.2
  beta(7, 5) = -65.0 / 27.0
  beta(8, 5) = 61.0 / 225.0
  beta(9, 5) = 704.0 / 45.0
  beta(10, 5) = -976.0 / 135.0
  beta(11, 5) = 4496.0 / 1025.0
  beta(13, 5) = beta(11, 5)
  beta(7, 6) = 125.0 / 54.0
  beta(8, 6) = -2.0 / 9.0
  beta(9, 6) = -107.0 / 9.0
  beta(10, 6) = 311.0 / 54.0
  beta(11, 6) = -301.0 / 82
  beta(12, 6) = -6.0 / 41.0
  beta(13, 6) = -289.0 / 82.0
  beta(8, 7) = 13.0 / 900.0
  beta(9, 7) = 67.0 / 90.0
  beta(10, 7) = -19.0 / 60.0
  beta(11, 7) = 2133.0 / 4100.0
  beta(12, 7) = -3.0 / 205.0
  beta(13, 7) = 2193.0 / 4100.0
  beta(9, 8) = 3.0
  beta(10, 8) = 17.0 / 6.0
  beta(11, 8) = 45.0 / 82.0
  beta(12, 8) = -3.0 / 41.0
  beta(13, 8) = 51.0 / 82.0
  beta(10, 9) = -1.0 / 12.0
  beta(11, 9) = 45.0 / 164.0
  beta(12, 9) = 3.0 / 41.0
  beta(13, 9) = 33.0 / 164.0
  beta(11, 10) = 18.0 / 41.0
  beta(12, 10) = 6.0 / 41.0
  beta(13, 10) = 12.0 / 41.0
  beta(13, 12) = 1.0
  
  ''''''''''''''''''
  ' begin simulation
  ''''''''''''''''''
  
  sim_definition(ri(), vi(), target, fxnew, tfinal)
  
  ' compute orbital elements at initial time
  
  eci2orb2(0.0, ri(), vi(), oevi())
  
  print " "
  print "program npoe - numerical prediction of orbital events"
  print "-----------------------------------------------------"
  print " "
  
  print "degree of gravity model  ", lgrav%
  
  print "order of gravity model   ", mgrav%
  print " "
  print "integration tolerance    ", tetol
  
  print "root-finding tolerance   ", rtol
  
  print " "
  
  if (isun% = 1) then print "simulation includes solar perturbations"
  
  if (imoon% = 1) then print "simulation includes lunar perturbations"
  
  if (idrag% = 1) then print "simulation includes drag perturbations"
  
  if (isrp% = 1) then print "simulation includes srp perturbations"
  
  keycheck
  
  ' print initial conditions
  
  print " "
  print "initial conditions"
  print "------------------"
  print " "
  
  jd2str(jdate0)
  
  print " "
  
  svprint(ri(), vi())
  
  keycheck
  
  oeprint(oevi())
  
  keycheck
  
  print " "
  
  print "working ...."
  
  ' root-finding step size (seconds)
  
  dtstep = 300.0
  
  tf = 0.0
  
  ' load initial integration vector
  
  for i% = 1 to 3
    
    yold(i%) = ri(i%)
    
    yold(i% + 3) = vi(i%)
    
  next i%
  
  do
    
    if (ittype% <> 0) then
      
      ' initialize root-bracketed flag
      
      irflg% = 0
      
      ' initial step size guess
      
      h = 60.0
      
      ' save current value of objective function
      
      fxold = fxnew
      
      ' set initial time to final time
      
      ti = tf
      
      ' increment final time
      
      tf = ti + dtstep
      
      ' save initial time as left side of bracket
      
      tisaved = ti
      
      ' integrate equations of motion
      
      rkf78(ti, tf, h, yold(), ynew())
      
      ' compute current value of objective function
      ' as the state and orbital elements at tf
      
      for i% = 1 to 3
        
        ri(i%) = ynew(i%)
        
        vi(i%) = ynew(i% + 3)
        
      next i%
      
      eci2orb2(tf, ri(), vi(), oevi())
      
    end if
    
    ' save current value of objective function
    
    select case ittype%
        
      case 0
        
        ''''''''''''''''''''''''
        ' orbit propagation only
        ''''''''''''''''''''''''
        
        ' step size guess (seconds)
        
        h = 180.0
        
        rkf78(0.0, tfinal, h, yold(), ynew())
        
        for i% = 1 to 3
          
          rf(i%) = ynew(i%)
          
          vf(i%) = ynew(i% + 3)
          
        next i%
        
        ' compute orbital elements at final time
        
        eci2orb2(tfinal, rf(), vf(), oevf())
        
        ' compute current julian date
        
        jdate = jdate0 + tfinal / 86400.0
        
        ' print final conditions
        
        print " "
        print "final conditions"
        print "----------------"
        print " "
        
        jd2str(jdate)
        
        print " "
        
        svprint(rf(), vf())
        
        keycheck
        
        oeprint(oevf())
        
        keycheck
        
        exit do
        
      case 1
        
        ' geodetic altitude
        
        fxnew = oevi(16) - talt
        
      case 2
        
        ' geodetic latitude
        
        fxnew = oevi(14) - tlat
        
      case 3
        
        ' east longitude
        
        fxnew = oevi(15) - tlong
        
      case 4
        
        ' geocentric declination
        
        fxnew = oevi(13) - tdec
        
      case 5
        
        ' true anomaly
        
        fxnew = oevi(6) - ttanom
        
      case 6
        
        ' argument of latitude
        
        fxnew = oevi(8) - targlat
        
      case 7
        
        ' flight path angle
        
        fxnew = oevi(11) - tfpa
        
      case 8
        
        ' orbital speed
        
        fxnew = oevi(28) - tspeed
        
      case 9
        
        ' right ascension
        
        fxnew = oevi(12) - trasc
        
    end select
    
    ' reset integration vector
    
    for i% = 1 to 6
      
      yold(i%) = ynew(i%)
      
    next i%
    
    ' check to see if the user-defined objective
    ' function has been bracketed during this step
    
    select case ittype%
        
      case 3, 5, 6, 9
        
        ' true anomaly, argument of latitude,
        ' east longitude and right ascension
        
        ' check for sign change
        
        if (fxnew * fxold < 0.0) then
          
          ' check for no boundary crossing
          
          if (abs(fxnew - fxold) < pi) then irflg% = 1
          
        end if
        
        ' check for boundary crossing
        
        if (abs(fxnew - fxold) > pi) then
          
          ' check for target argument = 0
          
          if (target = 0.0) then
            
            irflg% = 1
            
          elseif (fxnew < 0.0) then
            
            ' null
            
          elseif (fxold > 0.0) then
            
            ' null
            
          end if
          
        end if
        
      case 1, 2, 4, 7, 8
        
        ' geodetic altitude and latitude, geocentric declination,
        ' flight path angle, orbital speed
        
        if (fxnew * fxold < 0.0) then irflg% = 1
        
    end select
    
    ' check if time to bail out
    
    if (tf >= tfinal) then exit do
    
    ' if bracketed, find root of objective function
    
    if (irflg% = 1) then
      
      ' load 'working' time and state vector array
      ' as values on right side of bracket
      
      tiwrk = tf
      
      for i% = 1 to 6
        
        yi(i%) = ynew(i%)
        
      next i%
      
      ' find root of objective function
      
      realroot(tisaved, tiwrk, troot, froot)
      
      ' evaluate converged solution
      
      npoe_func(troot, xroot)
      
      for i% = 1 to 3
        
        rf(i%) = yfinal(i%)
        
        vf(i%) = yfinal(i% + 3)
        
      next i%
      
      ' compute orbital elements at solution
      
      eci2orb2(troot, rf(), vf(), oevf())
      
      ' compute current julian date
      
      jdate = jdate0 + troot / 86400.0
      
      ' print results at each solution
      
      print " "
      print " "
      
      jd2str(jdate)
      
      print " "
      
      svprint(rf(), vf())
      
      keycheck
      
      oeprint(oevf())
      
      keycheck
      
      print " "
      
      print "working ...."
      
    end if
    
  loop
  
  print " "
  
end
  
  '''''''''''''''''''''''''''''''''''''''''''''''''''
  '''''''''''''''''''''''''''''''''''''''''''''''''''
  
sub sim_definition(ri(), vi(), target, fxnew, tfinal)
  
  ' define simulation inputs subroutine
  
  '''''''''''''''''''''''''''''''''''''
  
  local sma, ecc, xinc, argper, raan, tanom
  
  local cmonth, cday, cyear, uthr, utmin, utsec
  
  local jdate, oev(32)
  
  ''''''''''''''''''''''''''''''''''''
  ' initial classical orbital elements
  ''''''''''''''''''''''''''''''''''''
  
  ' semimajor axis (kilometers; semimajor axis > 0)
  
  sma = 8000.0
  
  ' orbital eccentricity (non-dimensional); 0 <= eccentricity < 1)
  
  ecc = 0.0265427
  
  ' orbital inclination (degrees); 0 <= inclination <= 180)
  
  xinc = 28.5
  
  ' argument of perigee (degrees); 0 <= argument of perigee <= 360)
  
  argper = 0.0
  
  ' right ascension of the ascending node (degrees); 0 <= RAAN <= 360)
  
  raan = 45.0
  
  ' true anomaly (degrees); 0 <= true anomaly <= 360)
  
  tanom = 0.0
  
  oev(1) = sma
  oev(2) = ecc
  oev(3) = xinc * dtr
  oev(4) = argper * dtr
  oev(5) = raan * dtr
  oev(6) = tanom * dtr
  
  ''''''''''''''''''''''''''''''''''''''''''''
  ' request initial calendar date and UTC time
  ''''''''''''''''''''''''''''''''''''''''''''
  
  cmonth = 4
  
  cday = 18
  
  cyear = 2017
  
  uthr = 0.0
  
  utmin = 0.0
  
  utsec = 0.0
  
  ' simulation period (days)
  
  ndays = 0.25
  
  ' simulation period (seconds)
  
  tfinal = 86400.0 * ndays
  
  ''''''''''''''''''''''''''''''
  ' algorithm control parameters
  ''''''''''''''''''''''''''''''
  
  ' integration error tolerance; a value of 1.0e-8 is recommended)
  
  tetol = 1.0e-8
  
  ' root-finding error tolerance; a value between 1.0e-4 and 1.0e-6 is recommended)
  
  rtol = 1.0e-4
  
  '''''''''''''''''''''''''''''''''''''''''''
  ' request degree and order of gravity model
  '''''''''''''''''''''''''''''''''''''''''''
  
  ' degree of the gravity model (zonals); 0 <= zonals <= 8)
  
  lgrav% = 8
  
  ' order of the gravity model (tesserals); 0 <= tesserals <= 8)
  
  mgrav% = 8
  
  '''''''''''''''''''''''
  ' orbital perturbations
  '''''''''''''''''''''''
  
  ' include solar perturbation (1 = yes, 0 = no)
  
  isun% = 1
  
  ' include lunar perturbation (1 = yes, 0 = no)
  
  imoon% = 1
  
  ' include drag perturbation (1 = yes, 0 = no)
  
  idrag% = 1
  
  ' include srp perturbation (1 = yes, 0 = no)
  
  isrp% = 1
  
  ' drag coefficient (non-dimensional)
  
  cd = 2.2
  
  ' cross-sectional area (square meters)
  
  areadrag = 40.0
  
  ' spacecraft mass (kilograms)
  
  scmass = 6078.0
  
  ' compute ballistic coefficient
  
  bcoeff = 1.0e-6 * areadrag * cd /scmass
  
  ' reflectivity constant (non-dimensional)
  
  reflect = 1.5
  
  ' cross-sectional area (square meters)
  
  areasrp = 10.0
  
  ' spacecraft mass (kilograms)
  
  scmass = 6078.0
  
  ' compute srp coefficient
  
  csrp0 = 0.000001 * reflect * pss * aunit^2 * areasrp / scmass
  
  ' determine initial eci state vector
  
  orb2eci(oev(), ri(), vi())
  
  ' compute julian date at initial time
  
  julian(cmonth, cday, cyear, jdate)
  
  jdate0 = jdate + uthr / 24.0 + utmin / 1440.0 + utsec / 86400.0
  
  ' compute greenwich sidereal time at initial utc (radians)
  
  gast2(jdate0, gst0)
  
  ' determine complete set of initial orbital elements
  
  eci2orb2(0, ri(), vi(), oev())
  
  '''''''''''''''''''''''''''''''''''''''
  ' select the orbital element to predict
  '''''''''''''''''''''''''''''''''''''''
  
  ' 1 = geodetic altitude
  
  ' 2 = geodetic latitude
  
  ' 3 = east longitude
  
  ' 4 = geocentric declination
  
  ' 5 = true anomaly
  
  ' 6 = argument of latitude
  
  ' 7 = flight path angle
  
  ' 8 = orbital speed
  
  ' 9 = right ascension
  
  ' NOTE: to propagate without searching
  '       for events, set ittype% = 0
  
  ittype% = 2
  
  ''''''''''''''''''''''''''''''''''
  ' define target value and compute
  ' initial value of objective function
  '''''''''''''''''''''''''''''''''''''
  
  select case ittype%
      
    case 0
      
      ' propagate without searching for events
      
    case 1
      
      ' geodetic altitude (kilometers)
      
      talt = 200.0
      
      fxnew = oev(16) - talt
      
    case 2
      
      ' geodetic latitude (degrees); -90 <= latitude <= +90)
      
      xlat = 20.0
      
      tlat = xlat * dtr
      
      fxnew = oev(14) - tlat
      
    case 3
      
      ' east longitude (degrees); 0 <= longitude <= 360)
      
      elong = 120.0
      
      tlong = elong * dtr
      
      fxnew = oev(15) - tlong
      
      target = tlong
      
    case 4
      
      ' geocentric declination (degrees); -90 <= declination <= +90)
      
      dec = 20.0
      
      tdec = dec * dtr
      
      fxnew = oev(13) - tdec
      
    case 5
      
      ' true anomaly (degrees); 0 <= true anomaly <= 360)
      
      tanom = 260.0
      
      ttanom = tanom * dtr
      
      fxnew = oev(6) - ttanom
      
      target = ttanom
      
    case 6
      
      ' argument of latitude (degrees); 0 <= arglat <= 360)
      
      arglat = 200.0
      
      targlat = arglat * dtr
      
      fxnew = oev(8) - targlat
      
      target = targlat
      
    case 7
      
      ' flight path angle (degrees); -90 <= fpa <= +90)
      
      fpa = 1.0
      
      tfpa = fpa * dtr
      
      fxnew = oev(11) - tfpa
      
    case 8
      
      ' orbital speed (kilometers/second)
      
      tspeed = 0.0
      
      fxnew = oev(28) - tspeed
      
    case 9
      
      ' right ascension (degrees); 0 <= right ascension <= 360)
      
      rasc = 100.0
      
      trasc = rasc * dtr
      
      fxnew = oev(12) - trasc
      
      target = trasc
      
  end select
  
end sub
  
  ''''''''''''''''''
  ''''''''''''''''''
  
sub npoe_func(x, fx)
  
  ' npoe objective function subroutine
  
  ' input
  
  '  x = current simulation time
  
  ' output
  
  '  fx = objective function at time = x
  
  ''''''''''''''''''''''''''''''''''''''
  
  local r(3), v(3), oev(32)
  
  ' step size guess (seconds)
  
  h = 180.0
  
  ' integrate from tiwrk to requested time = x
  
  rkf78(tiwrk, x, h, yi(), yfinal())
  
  ' compute orbital elements at current time
  
  for i% = 1 to 3
    
    r(i%) = yfinal(i%)
    
    v(i%) = yfinal(i% + 3)
    
  next i%
  
  eci2orb2(x, r(), v(), oev())
  
  ' compute value of objective function
  
  select case ittype%
      
    case 1
      
      ' geodetic altitude
      
      fx = oev(16) - talt
      
    case 4
      
      ' geocentric declination
      
      fx = oev(13) - tdec
      
    case 5
      
      ' true anomaly
      
      if (oev(6) > pi and ttanom = 0.0) then
        
        fx = oev(6) - pi2
        
      else
        
        fx = oev(6) - ttanom
        
      end if
      
    case 7
      
      ' flight path angle
      
      fx = oev(11) - tfpa
      
    case 6
      
      ' argument of latitude
      
      if (oev(8) > pi and targlat = 0.0) then
        
        fx = oev(8) - pi2
        
      else
        
        fx = oev(8) - targlat
        
      end if
      
    case 3
      
      ' east longitude
      
      if (oev(15) > pi and tlong = 0.0) then
        
        fx = oev(15) - pi2
        
      else
        
        fx = oev(15) - tlong
        
      end if
      
    case 2
      
      ' geodetic latitude
      
      fx = oev(14) - tlat
      
    case 8
      
      ' orbital speed
      
      fx = oev(28) - tspeed
      
    case 9
      
      ' right ascension
      
      if (oev(12) > pi and trasc = 0.0) then
        
        fx = oev(12) - pi2
        
      else
        
        fx = oev(12) - trasc
        
      end if
      
  end select
  
end sub
  
  ''''''''''''''''''''''
  ''''''''''''''''''''''
  
sub ceqm(t, y(), ydot())
  
  ' first order form of Cowell's
  ' equations of orbital motion
  
  ' includes perturbations due to
  
  '   non-spherical earth gravity
  '   aerodynamic drag (us 76 model)
  '   solar radiation pressure
  '   sun and moon
  
  ' input
  
  '  t   = current simulation time
  '  y() = current eci state vector
  
  ' output
  
  '  ydot() = eci acceleration vector
  
  '''''''''''''''''''''''''''''''''''
  
  local agrav(3), asun(3), amoon(3)
  
  local adrag(3), asrp(3), vtmp(3)
  
  local rsun(3), rmoon(3)
  
  local rs2sc(3), rm2sc(3)
  
  ' geopotential perturbations
  
  gravity(t, y(), agrav())
  
  if (isun% = 1) then
    
    '''''''''''''''''''''''''''''''
    ' point mass solar perturbation
    '''''''''''''''''''''''''''''''
    
    jdate = jdate0 + t / 86400.0
    
    sun(jdate, rasc, decl, rsun())
    
    ' compute heliocentric position vector of the spacecraft
    
    for i% = 1 to 3
      
      rs2sc(i%) = y(i%) - rsun(i%)
      
    next i%
    
    ' f(q) formulation
    
    for i% = 1 to 3
      
      vtmp(i%) = y(i%) - 2.0 * rsun(i%)
      
    next i%
    
    dot1 = y(1) * vtmp(1) + y(2) * vtmp(2) + y(3) * vtmp(3)
    
    dot2 = vdot(rsun(), rsun())
    
    qsun = dot1 / dot2
    
    fsun = qsun * ((3.0 + 3.0 * qsun + qsun * qsun) / (1.0 + (1.0 + qsun)^1.5))
    
    d3sun = vecmag(rs2sc())^3
    
    ' point-mass gravity of the sun
    
    for i% = 1 to 3
      
      asun(i%) = -smu * (y(i%) + fsun * rsun(i%)) / d3sun
      
    next i%
    
  end if
  
  if (imoon% = 1) then
    
    '''''''''''''''''''''''''''''''
    ' point mass lunar perturbation
    '''''''''''''''''''''''''''''''
    
    jdate = jdate0 + t / 86400.0
    
    moon(jdate, rasc, decl, rmoon())
    
    ' compute selenocentric position vector of the spacecraft
    
    for i% = 1 to 3
      
      rm2sc(i%) = y(i%) - rmoon(i%)
      
    next i%
    
    ' f(q) formulation
    
    for i% = 1 to 3
      
      vtmp(i%) = y(i%) - 2.0 * rmoon(i%)
      
    next i%
    
    dot1 = y(1) * vtmp(1) + y(2) * vtmp(2) + y(3) * vtmp(3)
    
    dot2 = vdot(rmoon(), rmoon())
    
    qmoon = dot1 / dot2
    
    fmoon = qmoon * ((3.0 + 3.0 * qmoon + qmoon * qmoon) / (1.0 + (1.0 + qmoon)^1.5))
    
    d3moon = vecmag(rm2sc())^3
    
    ' point-mass gravity of the moon
    
    for i% = 1 to 3
      
      amoon(i%) = -mmu * (y(i%) + fmoon * rmoon(i%)) / d3moon
      
    next i%
    
  end if
  
  if (idrag% = 1) then
    
    '''''''''''''''''''''''''''''''
    ' atmospheric drag perturbation
    '''''''''''''''''''''''''''''''
    
    aero_drag(y(), adrag())
    
  end if
  
  if (isrp% = 1) then
    
    '''''''''''''''''''''''''''''''''''''''
    ' solar radiation pressure perturbation
    '''''''''''''''''''''''''''''''''''''''
    
    srp(t, y(), asrp())
    
  end if
  
  ' total acceleration vector
  
  ydot(1) = y(4)
  
  ydot(2) = y(5)
  
  ydot(3) = y(6)
  
  for i% = 1 to 3
    
    ydot(i% + 3) = agrav(i%) + asun(i%) + amoon(i%) + adrag(i%) + asrp(i%)
    
  next i%
  
end sub
  
  ''''''''''''''''''''''''''''''
  ''''''''''''''''''''''''''''''
  
sub rkf78(ti, tf, h, xi(), xf())
  
  ' solve first order system of differential equations
  
  ' Runge-Kutta-Fehlberg 7(8) method
  
  ' input
  
  '  ti    = initial simulation time (seconds)
  '  tf    = final simulation time (seconds)
  '  h     = initial guess for integration step size
  '  xi()  = integration vector at time = ti
  
  ' output
  
  '  xf() = integration vector at time = tf
  
  '''''''''''''''''''''''''''''''''''''''''
  
  local xdot(6), xwrk(6), f(6, 13)
  
  local sdt, dt, temp, twrk
  
  ' compute integration "direction"
  
  sdt = sgn(tf - ti)
  
  dt = abs(h) * sdt
  
  do
    
    ' load "working" time and integration vector
    
    twrk = ti
    
    for i% = 1 to 6
      
      xwrk(i%) = xi(i%)
      
    next i%
    
    ' check for last delta-t
    
    if (abs(dt) > abs(tf - ti)) then dt = tf - ti
    
    ' check for end of integration period
    
    if (abs(ti - tf) < 1.0e-12) then
      
      for i% = 1 to 6
        
        xf(i%) = xi(i%)
        
      next i%
      
      ' finished - exit do loop
      
      exit do
      
    end if
    
    ' evaluate equations of motion
    
    ceqm(ti, xi(), xdot())
    
    for i% = 1 to 6
      
      f(i%, 1) = xdot(i%)
      
    next i%
    
    ' compute solution
    
    for k% = 2 to 13
      
      kk% = k% - 1
      
      for i% = 1 to 6
        
        temp = 0.0
        
        for j% = 1 to kk%
          
          temp = temp + beta(k%, j%) * f(i%, j%)
          
        next j%
        
        xi(i%) = xwrk(i%) + dt * temp
        
      next i%
      
      ti = twrk + alph(k%) * dt
      
      ceqm(ti, xi(), xdot())
      
      for j% = 1 to 6
        
        f(j%, k%) = xdot(j%)
        
      next j%
      
    next k%
    
    for i% = 1 to 6
      
      temp = 0.0
      
      for l% = 1 to 13
        
        temp = temp + ch(l%) * f(i%, l%)
        
      next l%
      
      xi(i%) = xwrk(i%) + dt * temp
      
    next i%
    
    ' truncation error calculations
    
    xerr = tetol
    
    for i% = 1 to 6
      
      ter = abs((f(i%, 1) + f(i%, 11) - f(i%, 12) - f(i%, 13)) * ch(12) * dt)
      
      tol = abs(xi(i%)) * tetol + tetol
      
      tconst = ter / tol
      
      if (tconst > xerr) then xerr = tconst
      
    next i%
    
    ' compute new step size
    
    dt = 0.8 * dt * (1.0 / xerr) ^ (1.0 / 8.0)
    
    if (xerr > 1.0) then
      
      ' reject current step and try again
      
      ti = twrk
      
      for i% = 1 to 6
        
        xi(i%) = xwrk(i%)
        
      next i%
      
    else
      
      ''''''''''''''''''''''''''''''''''''
      ' accept current step              '
      ''''''''''''''''''''''''''''''''''''
      ' perform graphics, additional     '
      ' calculations, etc. at this point '
      ''''''''''''''''''''''''''''''''''''
      
    end if
    
  loop
  
end sub
  
  ''''''''''''''''''''''''''
  ''''''''''''''''''''''''''
  
sub gravity(t, y(), agrav())
  
  ' N degree and M order
  ' gravitational acceleration
  
  ' input
  
  '  t   = simulation time (seconds)
  '  y() = state vector
  
  ' output
  
  '  agrav() = ECI gravitational acceleration
  '            vector (km/sec/sec)
  
  ''''''''''''''''''''''''''''''''
  
  local aj2(3), aj4(3), lamda, pmt, sr1, sr2
  
  local p(51, 51), cn(51), sn(51), tn(51)
  
  local r1, r2, r3, r5, d1, d2, e1, e2, e3
  
  local im%, i1%, i11%, im1%, im2%
  
  r2 = y(1) * y(1) + y(2) * y(2) + y(3) * y(3)
  
  r1 = sqr(r2)
  
  r3 = r2 * r1
  
  if (lgrav% = 0 and mgrav% = 0) then
    
    ' Keplerian motion
    
    agrav(1) = 0.0
    agrav(2) = 0.0
    agrav(3) = 0.0
    
  elseif (lgrav% = 2 and mgrav% = 0) then
    
    ' j2 only
    
    r5 = r2 * r3
    
    d1 = -1.5 * xj2 * req * req * emu / r5
    
    d2 = 1.0 - 5.0 * y(3) * y(3) / r2
    
    agrav(1) = agrav(1) + y(1) * d1 * d2
    
    agrav(2) = agrav(2) + y(2) * d1 * d2
    
    agrav(3) = agrav(3) + y(3) * d1 * (d2 + 2.0)
    
  elseif (lgrav% = 4 and mgrav% = 0) then
    
    ' j2 and j4 only
    
    r5 = r2 * r3
    
    d1 = -1.5 * xj2 * req * req * emu / r5
    
    d2 = 1.0 - 5.0 * y(3) * y(3) / r2
    
    aj2(1) = y(1) * d1 * d2
    
    aj2(2) = y(2) * d1 * d2
    
    aj2(3) = y(3) * d1 * (d2 + 2.0)
    
    r4 = r2 * r2
    
    r7 = r2 * r5
    
    r9 = r2 * r7
    
    d3 = 15.0 * xj4 * req^4 * emu / (8.0 * r7)
    
    aj4(1) = y(1) * d3 * (1.0 - 14.0 * y(3) * y(3) / r2 + 21.0 * y(3)^4 / r4)
    
    aj4(2) = y(2) * d3 * (1.0 - 14.0 * y(3) * y(3) / r2 + 21.0 * y(3)^4 / r4)
    
    aj4(3) = y(3) * d3 * (5.0 - 70.0 * y(3) * y(3) / (3.0 * r2) + 21.0 * y(3)^4 / r4)
    
    for i% = 1 to 3
      
      agrav(i%) = aj2(i%) + aj4(i%)
      
    next i%
    
  else
    
    ' user-defined degree and order gravity model
    
    sr2 = y(1) * y(1) + y(2) * y(2)
    
    sr1 = sqr(sr2)
    
    sphi = y(3) / r1
    
    phi = asin(sphi)
    
    ' right ascension of greenwich
    
    pmt = modulo(gst0 + omega * t)
    
    ' east longitude of the spacecraft
    
    lamda = modulo(atan3(y(2), y(1)) - pmt)
    
    im% = mgrav%
    
    if (mgrav% < lgrav%) then im% = mgrav% + 1
    
    legend(lgrav%, im%, sphi, p())
    
    angles(mgrav%, lamda, phi, cn(), sn(), tn())
    
    e1 = 0.0
    e2 = 0.0
    e3 = 0.0
    
    d1 = req / r1
    
    d2 = 1.0
    
    for i1% = 1 to lgrav%
      
      f1 = 0.0
      f2 = 0.0
      f3 = 0.0
      
      i11% = i1% + 1
      
      for im% = 0 to i1%
        
        if (im% > mgrav%) then exit for
        
        im2% = im% + 2
        im1% = im% + 1
        
        d3 = ccoef(i11%, im1%) * cn(im1%) + scoef(i11%, im1%) * sn(im1%)
        
        f1 = f1 + d3 * p(i11%, im1%)
        
        if (im2% <= i11%) then f2 = f2 + d3 * (p(i11%, im2%) - tn(im1%) * p(i11%, im1%))
        
        if (im2% > i11%) then f2 = f2 - d3 * tn(im1%) * p(i11%, im1%)
        
        if (im% <> 0) then f3 = f3 + im% * (scoef(i11%, im1%) * cn(im1%) - ccoef(i11%, im1%) * sn(im1%)) * p(i11%, im1%)
        
      next im%
      
      d2 = d2 * d1
      
      e1 = e1 + d2 * i11% * f1
      
      e2 = e2 + d2 * f2
      
      e3 = e3 + d2 * f3
      
    next i1%
    
    d3 = emu / r1
    
    e1 = -e1 * d3 / r1
    
    e2 = e2 * d3
    
    e3 = e3 * d3
    
    d1 = e1 / r1
    
    d2 = y(3) / (r2 * sr1) * e2
    
    d3 = e3 / sr2
    
    agrav(1) = agrav(1) + (d1 - d2) * y(1) - d3 * y(2)
    
    agrav(2) = agrav(2) + (d1 - d2) * y(2) + d3 * y(1)
    
    agrav(3) = agrav(3) + d1 * y(3) + sr1 * e2 / r2
    
  end if
  
  ' complete gravity acceleration vector
  
  for i% = 1 to 3
    
    agrav(i%) = agrav(i%) - emu * y(i%) / r3
    
  next i%
  
end sub
  
  ''''''''''''''''''''''''''''''''''''
  ''''''''''''''''''''''''''''''''''''
  
sub angles(m%, a, b, cn(), sn(), tn())
  
  ' multiple angles subroutine
  
  ''''''''''''''''''''''''''''
  
  cn(1) = 1.0
  sn(1) = 0.0
  tn(1) = 0.0
  
  if (m% = 0) then return
  
  cn(2) = cos(a)
  sn(2) = sin(a)
  tn(2) = tan(b)
  
  if (m% = 1) then return
  
  for i% = 2 to m%
    
    cn(i% + 1) = 2.0 * cn(2) * cn(i%) - cn(i% - 1)
    
    sn(i% + 1) = 2.0 * cn(2) * sn(i%) - sn(i% - 1)
    
    tn(i% + 1) = tn(2) + tn(i%)
    
  next i%
  
end sub
  
  ''''''''''''''''''''''''
  ''''''''''''''''''''''''
  
sub legend(n%, m%, x, p())
  
  ' Legendre polynomials subroutine
  
  '''''''''''''''''''''''''''''''''
  
  local y
  
  p(1, 1) = 1.0
  
  p(2, 1) = x
  
  for i% = 2 to n%
    
    p(i% + 1, 1) = ((2 * i% - 1) * x * p(i%, 1) -(i% - 1) * p(i% - 1, 1)) / i%
    
  next i%
  
  if (m% = 0) then return
  
  y = sqr(1.0 - x * x)
  
  p(2, 2) = y
  
  if (m% = 1) then
    
    ' null
    
  else
    
    for i% = 2 to m%
      
      p(i% + 1, i% + 1) = (2 * i% - 1) * y * p(i%, i%)
      
    next i%
    
  end if
  
  for i% = 2 to n%
    
    i1% = i% - 1
    
    for j% = 1 to i1%
      
      if (j% > m%) then exit for
      
      p(i% + 1, j% + 1) = (2 * i% - 1) * y * p(i%, j%)
      
      if ((i% - 2) >= j%) then
        
        p(i% + 1, j% + 1) = p(i% + 1, j% + 1) + p(i% - 1, j% + 1)
        
      end if
      
    next j%
    
  next i%
  
end sub
  
  ''''''''''''''''''''''
  ''''''''''''''''''''''
  
sub srp(t, sv(), asrp())
  
  ' ECI acceleration vector due to
  ' solar radiation pressure subroutine
  
  ' input
  
  '  t    = simulation time
  '  sv() = current state vector
  
  ' output
  
  '  asrp() = eci acceleration vector
  
  '''''''''''''''''''''''''''''''''''
  
  local rsun(3), usun(3), usat(3), rs2s(3)
  
  local csrp, ratm
  
  local jdate, rmsun, rscm, rs2sm, rs2sm3
  
  local a, b, c, d, e, p, q, u, v, w, x, y
  
  ' current julian date
  
  jdate = jdate0 + t / 86400.0
  
  ' compute solar ephemeris
  
  sun(jdate, rasc, decl, rsun())
  
  ' geocentric distance of the sun (kilometers)
  
  rmsun = vecmag(rsun())
  
  ' geocentric unit vector of the sun
  
  uvector(rsun(), usun())
  
  ' geocentric distance of the spacecraft
  
  rscm = sqr(sv(1) * sv(1) + sv(2) * sv(2) + sv(3) * sv(3))
  
  ' compute unit position vector of spacecraft
  
  for i% = 1 to 3
    
    usat(i%) = sv(i%) / rscm
    
  next i%
  
  ' determine shadow conditions
  
  a = usat(2) * usun(3) - usat(3) * usun(2)
  b = usat(3) * usun(1) - usat(1) * usun(3)
  c = usat(1) * usun(2) - usat(2) * usun(1)
  
  d = sqr(a * a + b * b + c * c)
  
  e = vdot(usat(), usun())
  
  u = asin(0.00460983743 / (rmsun / aunit))
  
  p = asin(0.0046951089 / (rmsun / aunit))
  
  if (e > 0.0) then
    
    q = -d
    
  else
    
    q = d
    
  end if
  
  ' increase the Earth's radius by 90 km
  ' to account for the atmosphere
  
  ratm = req + 90.0
  
  b = asin(ratm / rscm)
  
  v = b - u
  w = b + p
  
  x = sin(v)
  y = sin(w)
  
  ' determine shadow conditions
  
  if (q <= y and q > x) then
    
    ' penumbra
    
    csrp = 0.0
    
  elseif (q <= x and q >= 0.0) then
    
    ' umbra
    
    csrp = 0.0
    
  else
    
    ' sunlight
    
    csrp = csrp0
    
  end if
  
  ' compute vector from the satellite to the sun
  
  for i% = 1 to 3
    
    rs2s(i%) = sv(i%) - rsun(i%)
    
  next i%
  
  rs2sm = vecmag(rs2s())
  
  rs2sm3 = rs2sm * rs2sm * rs2sm
  
  ' compute acceleration vector
  
  for i% = 1 to 3
    
    asrp(i%) = csrp * rs2s(i%) / rs2sm3
    
  next i%
  
end sub
  
  ''''''''''''''''''''''''''
  ''''''''''''''''''''''''''
  
sub aero_drag(sv(), adrag())
  
  ' acceleration due to atmospheric drag subroutine
  
  ' us standard 1976 atmosphere model
  
  ' input
  
  '  sv() = ECI state vector
  
  ' output
  
  '  adrag() = ECI drag acceleration vector
  '            (kilometers/second/second)
  
  '''''''''''''''''''''''''''''''''''''''
  
  local v3(3), atmtmp, atmdns, atmprs, atmsos
  
  local rmag, dec, alt, xlat, vrel
  
  ' compute geodetic coordinates
  
  rmag = sqr(sv(1) * sv(1) + sv(2) * sv(2) + sv(3) * sv(3))
  
  dec = asin(sv(3) / rmag)
  
  geodet1(rmag, dec, alt, xlat)
  
  ' calculate atmospheric density (kg/km^3)
  
  us76(alt, atmtmp, atmdns, atmprs, atmsos)
  
  ' velocity vector of atmosphere relative to spacecraft
  
  v3(1) = sv(4) + omega * sv(2)
  v3(2) = sv(5) - omega * sv(1)
  v3(3) = sv(6)
  
  vrel = vecmag(v3())
  
  ' acceleration vector
  
  for i% = 1 to 3
    
    adrag(i%) = -0.5 * atmdns * bcoeff * vrel * v3(i%)
    
  next i%
  
end sub
  
  ''''''''''''''''''''''''''
  ''''''''''''''''''''''''''
  
sub orb2eci(oev(), r(), v())
  
  ' convert orbital elements to state vector subroutine
  
  ' input
  
  '  oev(1) = semimajor axis (kilometers)
  '  oev(2) = orbital eccentricity (non-dimensional)
  '  oev(3) = orbital inclination (radians)
  '           (0 <= oev(3) <= pi)
  '  oev(4) = argument of perigee (radians)
  '           (0 <= oev(4) <= 2 pi)
  '  oev(5) = right ascension of ascending node or
  '           east longitude of ascending node (radians)
  '           (0 <= oev(5) <= 2 pi)
  '  oev(6) = true anomaly (radians)
  '           (0 <= oev(6) <= 2 pi)
  
  ' output
  
  '  r() = geocentric, equatorial position vector (kilometers)
  '  v() = geocentric, equatorial velocity vector (kilometers/second)
  
  ' special notes
  
  '  (1) eci coordinates if oev(5) = raan
  
  '  (2) ecf coordinates if oev(5) = east longitude
  
  '  (3) valid for elliptic and hyperbolic orbits
  
  '''''''''''''''''''''''''''''''''''''''''''''''
  
  local sma, ecc, inc, argper, nangle, tanom
  
  local slr, rm, arglat, sarglat, carglat
  
  local sinc, cinc, snangle, cnangle, c1, c2, c3
  
  ' "unload" orbital elements array
  
  sma = oev(1)
  ecc = oev(2)
  inc = oev(3)
  argper = oev(4)
  nangle = oev(5)
  tanom = oev(6)
  
  ' semiparameter
  
  slr = sma * (1.0 - ecc * ecc)
  
  ' radius magnitude
  
  rm = slr / (1.0 + ecc * cos(tanom))
  
  ' argument of latitude
  
  arglat = argper + tanom
  
  sarglat = sin(arglat)
  
  carglat = cos(arglat)
  
  sinc = sin(inc)
  
  cinc = cos(inc)
  
  snangle = sin(nangle)
  
  cnangle = cos(nangle)
  
  c1 = sqr(emu / slr)
  
  c2 = ecc * cos(argper) + carglat
  
  c3 = ecc * sin(argper) + sarglat
  
  ' x, y and z components of the position vector
  
  r(1) = rm * (cnangle * carglat - snangle * cinc * sarglat)
  
  r(2) = rm * (snangle * carglat + cinc * sarglat * cnangle)
  
  r(3) = rm * sinc * sarglat
  
  ' x, y and z components of the velocity vector
  
  v(1) = -c1 * (cnangle * c3 + snangle * cinc * c2)
  
  v(2) = -c1 * (snangle * c3 - cnangle * cinc * c2)
  
  v(3) = c1 * c2 * sinc
  
end sub
  
  '''''''''''''''''''''''''''''''
  '''''''''''''''''''''''''''''''
  
sub eci2orb2(ut, r(), v(), oev())
  
  ' convert eci state vector to complete
  ' set of classical orbital elements
  
  ' input
  
  '  ut   = universal time (seconds)
  '         (0 <= ut <= 86400)
  '  r()  = eci position vector (kilometers)
  '  v()  = eci velocity vector (kilometers/second)
  
  ' output
  
  '  oev(1)  = semimajor axis (kilometers)
  '  oev(2)  = orbital eccentricity (non-dimensional)
  '            (0 <= eccentricity < 1)
  '  oev(3)  = orbital inclination (radians)
  '            (0 <= inclination <= pi)
  '  oev(4)  = argument of perigee (radians)
  '            (0 <= argument of perigee <= 2 pi)
  '  oev(5)  = right ascension of ascending node (radians)
  '            (0 <= raan <= 2 pi)
  '  oev(6)  = true anomaly (radians)
  '            (0 <= true anomaly <= 2 pi)
  '  oev(7)  = orbital period (seconds)
  '  oev(8)  = argument of latitude (radians)
  '            (0 <= argument of latitude <= 2 pi)
  '  oev(9)  = east longitude of ascending node (radians)
  '            (0 <= east longitude <= 2 pi)
  '  oev(10) = specific orbital energy (kilometer^2/second^2)
  '  oev(11) = flight path angle (radians)
  '            (-0.5 pi <= fpa <= 0.5 pi)
  '  oev(12) = right ascension (radians)
  '            (-2 pi <= right ascension <= 2 pi)
  '  oev(13) = declination (radians)
  '            (-0.5 pi <= declination <= 0.5 pi)
  '  oev(14) = geodetic latitude of subpoint (radians)
  '            (-0.5 pi <= latitude <= 0.5 pi)
  '  oev(15) = east longitude of subpoint (radians)
  '            (-2 pi <= latitude <= 2 pi)
  '  oev(16) = geodetic altitude (kilometers)
  '  oev(17) = geocentric radius of perigee (kilometers)
  '  oev(18) = geocentric radius of apogee (kilometers)
  '  oev(19) = perigee velocity (kilometers/second)
  '  oev(20) = apogee velocity (kilometers/second)
  '  oev(21) = geodetic altitude of perigee (kilometers)
  '  oev(22) = geodetic altitude of apogee (kilometers)
  '  oev(23) = geodetic latitude of perigee (radians)
  '            (-0.5 pi <= latitude <= 0.5 pi)
  '  oev(24) = geodetic latitude of apogee (radians)
  '            (-0.5 pi <= latitude <= 0.5 pi)
  '  oev(25) = mean motion (radians/second)
  '  oev(26) = mean anomaly (radians)
  '            (-2 pi <= mean anomaly <= 2 pi)
  '  oev(27) = eccentric anomaly (radians)
  '            (-2 pi <= eccentric anomaly <= 2 pi)
  '  oev(28) = velocity magnitude (kilometers/second)
  '  oev(29) = time until periapsis passage
  '  oev(30) = time until ascending node crossing
  '  oev(31) = nodal period
  '  oev(32) = anomalistic period
  
  '''''''''''''''''''''''''''''''
  
  local rmag, vmag, sma, p, q, const1
  
  local h, xk, x1, y1, eccm, inc, xlambdat
  
  local raan, argper, tanom
  
  local rhat(3), vhat(3), hhat(3), vtmp(3), ecc(3)
  
  local hv(3), fhat(3), ghat(3)
  
  ' position and velocity magnitude
  
  rmag = vecmag(r())
  
  vmag = vecmag(v())
  
  ' position and velocity unit vectors
  
  uvector(r(), rhat())
  
  uvector(v(), vhat())
  
  ' angular momentum vectors
  
  vcross(r(), v(), hv())
  
  uvector(hv(), hhat())
  
  ' eccentricity vector
  
  for i% = 1 to 3
    
    vtmp(i%) = v(i%) / emu
    
  next i%
  
  vcross(vtmp(), hv(), ecc())
  
  for i% = 1 to 3
    
    ecc(i%) = ecc(i%) - rhat(i%)
    
  next i%
  
  ' semimajor axis
  
  sma = 1.0 / (2.0 / rmag - vmag^2 / emu)
  
  ' keplerian period
  
  period = pi2 / sqr(emu / sma^3)
  
  p = hhat(1) / (1.0 + hhat(3))
  
  q = -hhat(2) / (1.0 + hhat(3))
  
  const1 = 1.0 / (1.0 + p^2 + q^2)
  
  fhat(1) = const1 * (1.0 - p^2 + q^2)
  fhat(2) = const1 * 2.0 * p * q
  fhat(3) = -const1 * 2.0 * p
  
  ghat(1) = const1 * 2.0 * p * q
  ghat(2) = const1 * (1.0 + p^2 - q^2)
  ghat(3) = const1 * 2.0 * q
  
  h = vdot(ecc(), ghat())
  
  xk = vdot(ecc(), fhat())
  
  x1 = vdot(r(), fhat())
  
  y1 = vdot(r(), ghat())
  
  ' orbital eccentricity
  
  eccm = sqr(h * h + xk * xk)
  
  ' orbital inclination
  
  inc = 2.0 * atn(sqr(p * p + q * q))
  
  ' true longitude
  
  xlambdat = atan3(y1, x1)
  
  ' check for equatorial orbit
  
  if (inc > 1.0e-10) then
    
    raan = atan3(p, q)
    
    xlan = modulo(raan - (gst0 + omega * ut))
    
  else
    
    raan = 0.0
    
    xlan = modulo(gst0 + omega * ut)
    
  end if
  
  ' check for circular orbit
  
  if (eccm > 1.0e-10) then
    
    argper = modulo(atan3(h, xk) - raan)
    
  else
    
    argper = 0.0
    
  end if
  
  ' true anomaly
  
  tanom = modulo(xlambdat - raan - argper)
  
  a = sin(tanom) * sqr(1.0 - eccm * eccm)
  
  b = eccm + cos(tanom)
  
  ' eccentric anomaly
  
  eanom = atan3(a, b)
  
  ' mean anomaly
  
  xmanom = eanom - eccm * sin(eanom)
  
  ' east longitude of subpoint
  
  a = modulo(gst0 + omega * ut)
  
  bx = sin(a)
  cx = cos(a)
  
  c1 = cx * r(1) + bx * r(2)
  
  c2 = cx * r(2) - bx * r(1)
  
  elong = atan3(c2, c1)
  
  energy = 0.5 * vmag * vmag - emu / rmag
  
  period = pi2 * sma * sqr(sma / emu)
  
  rperigee = sma * (1.0 - eccm)
  
  rapogee = sma * (1.0 + eccm)
  
  vperigee = sqr(2.0 * emu * rapogee / (rperigee * (rapogee + rperigee)))
  
  vapogee = sqr(2.0 * emu * rperigee / (rapogee * (rapogee + rperigee)))
  
  decper = asin(sin(inc) * sin(argper))
  
  geodet1(rperigee, decper, altper, xlatper)
  
  decapo = asin(sin(inc) * sin(pi + argper))
  
  geodet1(rapogee, decapo, altapo, xlatapo)
  
  rdotv = vdot(r(), v())
  
  fpa = asin(rdotv / rmag / vmag)
  
  dec = asin(r(3) / rmag)
  
  geodet1(rmag, dec, alt, xlat)
  
  ' mean motion
  
  xmm = sqr(emu / sma * sma * sma)
  
  ras = atan3(r(2), r(1))
  
  if (sma > 0.0) then
    
    ' trajectory event times and orbital periods
    
    tevent(sma, eccm, argper, tanom, ttpp, ttanc)
    
    tperiod(sma, eccm, inc, argper, tnodal, tanomal)
    
  else
    
    ttpp = 0.0
    ttanc = 0.0
    tnodal = 0.0
    tanomal = 0.0
    
  end if
  
  oev(1) = sma
  oev(2) = eccm
  oev(3) = inc
  oev(4) = argper
  oev(5) = raan
  oev(6) = tanom
  oev(7) = period
  oev(8) = modulo(tanom + argper)
  oev(9) = xlan
  oev(10) = energy
  oev(11) = fpa
  oev(12) = ras
  oev(13) = dec
  oev(14) = xlat
  oev(15) = elong
  oev(16) = alt
  oev(17) = rperigee
  oev(18) = rapogee
  oev(19) = vperigee
  oev(20) = vapogee
  oev(21) = altper
  oev(22) = altapo
  oev(23) = xlatper
  oev(24) = xlatapo
  oev(25) = xmm
  oev(26) = xmanom
  oev(27) = eanom
  oev(28) = vmag
  oev(29) = ttpp
  oev(30) = ttanc
  oev(31) = tnodal
  oev(32) = tanomal
  
end sub
  
  ''''''''''''''''''''
  ''''''''''''''''''''
  
sub gast2(jdate, gast)
  
  ' greenwich apparent sidereal time subroutine
  
  ' input
  
  '  jdate = julian date
  
  ' output
  
  '  gast = greenwich apparent sidereal time (radians)
  
  ''''''''''''''''''''''''''''''''''''''''''''''''''''
  
  LOCAL t, t2, t3, eqeq, dpsi, deps
  
  LOCAL th, tl, obm, obt, st, x
  
  local tjdh as integer, tjdl
  
  tjdh = int(jdate)
  
  tjdl = jdate - tjdh
  
  th = (tjdh - 2451545.0) / 36525
  
  tl = tjdl / 36525.0
  
  t = th + tl
  
  t2 = t * t
  
  t3 = t2 * t
  
  ' obtain equation of the equinoxes
  
  eqeq = 0.0
  
  ' obtain nutation parameters in seconds of arc
  
  nut2000_lp(jdate, dpsi, deps)
  
  ' compute mean obliquity of the ecliptic in seconds of arc
  
  obm = 84381.4480 - 46.8150 * t - 0.00059 * t2 + 0.001813 * t3
  
  ' compute true obliquity of the ecliptic in seconds of arc
  
  obliq = obm + deps
  
  ' compute equation of the equinoxes in seconds of time
  
  eqeq = (dpsi / 15.0) * cos(obliq / seccon)
  
  st = eqeq - 6.2e-6 * t3 + 0.093104 * t2 + 67310.54841 + 8640184.812866 * tl + 3155760000.0 * tl + 8640184.812866 * th + 3155760000.0 * th
  
  ' modulo 24 hours
  
  x = st / 3600.0
  
  gast = x - 24.0 * fix(x / 24.0)
  
  if (gast < 0.0) then
    
    gast = gast + 24.0
    
  end if
  
  ' convert to radians
  
  gast = pi2 * (gast / 24.0)
  
END sub
  
  ''''''''''''''''''''''''''''''''
  ''''''''''''''''''''''''''''''''
  
sub sun(jdate, rasc, decl, rsun())
  
  ' precision ephemeris of the Sun subroutine
  
  ' input
  
  '  jdate = julian ephemeris day
  
  ' output
  
  '  rsun = eci position of the sun
  
  '''''''''''''''''''''''''''''''''
  
  local u, a1, a2, psi, deps, meps, eps, seps, ceps
  
  local dl, dr, w, srl, crl, srb, crb, sra, cra
  
  u = (jdate - 2451545.0) / 3652500.0
  
  ' compute nutation in longitude
  
  a1 = 2.18 + u * (-3375.7 + u * 0.36)
  
  a2 = 3.51 + u * (125666.39 + u * 0.1)
  
  psi = 0.0000001 * (-834.0 * sin(a1) - 64.0 * sin(a2))
  
  ' compute nutation in obliquity
  
  deps = 0.0000001 * u * (-226938 + u * (-75 + u * (96926 + u * (-2491 - u * 12104))))
  
  meps = 0.0000001 * (4090928.0 + 446.0 * cos(a1) + 28.0 * cos(a2))
  
  eps = meps + deps
  
  obliq = eps
  
  seps = sin(eps)
  
  ceps = cos(eps)
  
  dl = 0.0
  
  dr = 0.0
  
  for i% = 1 to 50
    
    w = xsa(i%) + xsb(i%) * u
    
    dl = dl + xsl(i%) * sin(w)
    
    if (xsr(i%) <> 0.0) then
      
      dr = dr + xsr(i%) * cos(w)
      
    end if
    
  next i%
  
  dl = modulo(dl * 0.0000001 + 4.9353929 + 62833.196168 * u)
  
  dr = aunit * (dr * 0.0000001 + 1.0001026)
  
  rlsun = modulo(dl + 0.0000001 * (-993.0 + 17.0 * cos(3.1 + 62830.14 * u)) + psi)
  
  rb = 0.0
  
  ' compute geocentric declination and right ascension
  
  crl = cos(rlsun)
  srl = sin(rlsun)
  crb = cos(rb)
  srb = sin(rb)
  
  decl = asin(ceps * srb + seps * crb * srl)
  
  sra = -seps * srb + ceps * crb * srl
  
  cra = crb * crl
  
  rasc = atan3(sra, cra)
  
  ' geocentric equatorial position vector of the Sun (au)
  
  rsun(1) = dr * cos(rasc) * cos(decl)
  
  rsun(2) = dr * sin(rasc) * cos(decl)
  
  rsun(3) = dr * sin(decl)
  
end sub
  
  ''''''''''''''''''''''''''''''''''
  ''''''''''''''''''''''''''''''''''
  
sub moon(jdate, rasc, decl, rmoon())
  
  ' geocentric position of the moon subroutine
  
  ' input
  
  '  jdate = tdb julian day
  
  ' output
  
  '  rmoon = eci position vector of the moon (kilometers)
  
  '''''''''''''''''''''''''''''''''''''''''''''''''''''''
  
  local t1, t2, t3
  
  local t4, dpsi, deps
  
  LOCAL ll, d, lp
  
  local l, f, t
  
  LOCAL ve, ma, ju
  
  local rm, dv, pl
  
  local plat, plon
  
  LOCAL a, b
  
  local obliq
  
  ' get nutations and obliquity
  
  obliq_lp(jdate, dpsi, deps, obliq)
  
  ' evaluate lunar ephemeris
  
  t1 = (jdate - 2451545.0) / 36525.0
  
  t2 = t1 * t1
  t3 = t1 * t1 * t1
  t4 = t1 * t1 * t1 * t1
  
  ' calculate fundamental arguments (radians)
  
  ll = dtr * (218 + (18 * 60 + 59.95571) / 3600)
  
  ll = modulo(ll + atr * (1732564372.83264 * t1 - 4.7763 * t2 + .006681 * t3 - 0.00005522 * t4))
  
  d = dtr * (297 + (51 * 60 + .73512) / 3600)
  
  d = modulo(d + atr * (1602961601.4603 * t1 - 5.8681 * t2 + .006595 * t3 - 0.00003184 * t4))
  
  lp = dtr * (357 + (31 * 60 + 44.79306) / 3600)
  
  lp = modulo(lp + atr * (129596581.0474 * t1 - .5529 * t2 + 0.000147 * t3))
  
  l = dtr * (134 + (57 * 60 + 48.28096) / 3600)
  
  l = modulo(l + atr * (1717915923.4728 * t1 + 32.3893 * t2 + .051651 * t3 - 0.0002447 * t4))
  
  f = dtr * (93 + (16 * 60 + 19.55755) / 3600)
  
  f = modulo(f + atr * (1739527263.0983 * t1 - 12.2505 * t2 - .001021 * t3 + 0.00000417 * t4))
  
  t = dtr * (100 + (27 * 60 + 59.22059) / 3600)
  
  t = modulo(t + atr * (129597742.2758 * t1 - .0202 * t2 + .000009 * t3 + 0.00000015 * t4))
  
  ve = dtr * (181 + (58 * 60 + 47.28305) / 3600)
  
  ve = modulo(ve + atr * 210664136.43355 * t1)
  
  ma = dtr * (355 + (25 * 60 + 59.78866) / 3600)
  
  ma = modulo(ma + atr * 68905077.59284 * t1)
  
  ju = dtr * (34 + (21 * 60 + 5.34212) / 3600)
  
  ju = modulo(ju + atr * 10925660.42861 * t1)
  
  ' compute geocentric distance (kilometers)
  
  ' a(c,0,r) series
  
  rm = 385000.52899
  rm = rm - 20905.35504 * sin(l + pidiv2)
  rm = rm - 3699.11092 * sin(2 * d - l + pidiv2)
  rm = rm - 2955.96756 * sin(2 * d + pidiv2)
  rm = rm - 569.92512 * sin(2 * l + pidiv2)
  rm = rm + 246.15848 * sin(2 * d - 2 * l + pidiv2)
  rm = rm - 204.58598 * sin(2 * d - lp + pidiv2)
  rm = rm - 170.73308 * sin(2 * d + l + pidiv2)
  rm = rm - 152.13771 * sin(2 * d - lp - l + pidiv2)
  rm = rm - 129.62014 * sin(lp - l + pidiv2)
  rm = rm + 108.7427 * sin(d + pidiv2)
  rm = rm + 104.75523 * sin(lp + l + pidiv2)
  rm = rm + 79.66056 * sin(l - 2 * f + pidiv2)
  rm = rm + 48.8883 * sin(lp + pidiv2)
  rm = rm - 34.78252 * sin(4 * d - l + pidiv2)
  rm = rm + 30.82384 * sin(2 * d + lp + pidiv2)
  rm = rm + 24.20848 * sin(2 * d + lp - l + pidiv2)
  rm = rm - 23.21043 * sin(3 * l + pidiv2)
  rm = rm - 21.63634 * sin(4 * d - 2 * l + pidiv2)
  rm = rm - 16.67471 * sin(d + lp + pidiv2)
  rm = rm + 14.40269 * sin(2 * d - 3 * l + pidiv2)
  rm = rm - 12.8314 * sin(2 * d - lp + l + pidiv2)
  rm = rm - 11.64995 * sin(4 * d + pidiv2)
  rm = rm - 10.44476 * sin(2 * d + 2 * l + pidiv2)
  rm = rm + 10.32111 * sin(2 * d - 2 * f + pidiv2)
  rm = rm + 10.0562 * sin(2 * d - lp - 2 * l + pidiv2)
  rm = rm - 9.88445 * sin(2 * d - 2 * lp + pidiv2)
  rm = rm + 8.75156 * sin(2 * d - l - 2 * f + pidiv2)
  rm = rm - 8.37911 * sin(d - l + pidiv2)
  rm = rm - 7.00269 * sin(lp - 2 * l + pidiv2)
  rm = rm + 6.322 * sin(d + l + pidiv2)
  rm = rm + 5.75085 * sin(lp + 2 * l + pidiv2)
  rm = rm - 4.95013 * sin(2 * d - 2 * lp - l + pidiv2)
  rm = rm - 4.42118 * sin(2 * l - 2 * f + pidiv2)
  rm = rm + 4.13111 * sin(2 * d + l - 2 * f + pidiv2)
  rm = rm - 3.95798 * sin(4 * d - lp - l + pidiv2)
  rm = rm + 3.25824 * sin(3 * d - l + pidiv2)
  rm = rm - 3.1483 * sin(2 * f + pidiv2)
  rm = rm + 2.61641 * sin(2 * d + lp + l + pidiv2)
  rm = rm + 2.35363 * sin(2 * d + 2 * lp - l + pidiv2)
  rm = rm - 2.11713 * sin(2 * lp - l + pidiv2)
  rm = rm - 1.89704 * sin(4 * d - lp - 2 * l + pidiv2)
  rm = rm - 1.73853 * sin(d - 2 * l + pidiv2)
  rm = rm - 1.57139 * sin(4 * d - lp + pidiv2)
  rm = rm - 1.42255 * sin(4 * d + l + pidiv2)
  rm = rm - 1.41893 * sin(3 * d + pidiv2)
  rm = rm + 1.16553 * sin(2 * lp + l + pidiv2)
  rm = rm - 1.11694 * sin(4 * l + pidiv2)
  rm = rm + 1.06567 * sin(2 * lp + pidiv2)
  rm = rm - .93332 * sin(d + lp + l + pidiv2)
  rm = rm + .86243 * sin(3 * d - 2 * l + pidiv2)
  rm = rm + .85124 * sin(d + lp - l + pidiv2)
  rm = rm - .8488 * sin(2 * d - lp + 2 * l + pidiv2)
  rm = rm - .79563 * sin(d - 2 * f + pidiv2)
  rm = rm + .77854 * sin(2 * d - 4 * l + pidiv2)
  rm = rm + .77404 * sin(2 * d - 2 * l + 2 * f + pidiv2)
  rm = rm - .66968 * sin(2 * d + 3 * l + pidiv2)
  rm = rm - .65753 * sin(2 * d - 2 * lp + l + pidiv2)
  rm = rm + .65706 * sin(2 * d - lp - 2 * f + pidiv2)
  rm = rm + .59632 * sin(2 * d - l + 2 * f + pidiv2)
  rm = rm + .57879 * sin(4 * d + lp - l + pidiv2)
  rm = rm - .51423 * sin(4 * d - 3 * l + pidiv2)
  rm = rm - .50792 * sin(4 * d - 2 * f + pidiv2)
  rm = rm + .49755 * sin(d - lp + pidiv2)
  rm = rm + .49504 * sin(2 * d - lp - 3 * l + pidiv2)
  rm = rm + .47262 * sin(2 * d - 2 * l - 2 * f + pidiv2)
  rm = rm - .4225 * sin(6 * d - 2 * l + pidiv2)
  rm = rm - .42241 * sin(lp - 3 * l + pidiv2)
  rm = rm - .41071 * sin(2 * d - 3 * lp + pidiv2)
  rm = rm + .37852 * sin(d + 2 * l + pidiv2)
  rm = rm + .35508 * sin(lp + 3 * l + pidiv2)
  rm = rm + .34302 * sin(2 * d - 2 * lp - 2 * l + pidiv2)
  rm = rm + .33463 * sin(lp - l + 2 * f + pidiv2)
  rm = rm + .33225 * sin(d + lp - 2 * l + pidiv2)
  rm = rm + .32334 * sin(2 * d - lp - l - 2 * f + pidiv2)
  rm = rm - .32176 * sin(4 * d - l - 2 * f + pidiv2)
  rm = rm - .28663 * sin(6 * d - l + pidiv2)
  rm = rm + .28399 * sin(2 * d + 2 * l - 2 * f + pidiv2)
  rm = rm - .27904 * sin(4 * d - 2 * lp - l + pidiv2)
  rm = rm + .2556 * sin(3 * d - lp - l + pidiv2)
  rm = rm - .2481 * sin(lp + l - 2 * f + pidiv2)
  rm = rm + .24452 * sin(4 * d + lp + pidiv2)
  rm = rm + .23695 * sin(4 * d + lp - 2 * l + pidiv2)
  rm = rm - .21258 * sin(3 * d + lp - l + pidiv2)
  rm = rm + .21251 * sin(2 * d + lp + 2 * l + pidiv2)
  rm = rm + .20941 * sin(2 * d - lp + l - 2 * f + pidiv2)
  rm = rm - .20285 * sin(4 * d - lp + l + pidiv2)
  rm = rm + .20099 * sin(3 * d - 2 * f + pidiv2)
  rm = rm - .18567 * sin(lp - 2 * f + pidiv2)
  rm = rm - .18316 * sin(6 * d - 3 * l + pidiv2)
  rm = rm + .16857 * sin(2 * d + lp - 3 * l + pidiv2)
  rm = rm - .15802 * sin(lp + 2 * f + pidiv2)
  rm = rm - .15707 * sin(3 * d - lp + pidiv2)
  rm = rm - .14806 * sin(2 * d - 3 * lp - l + pidiv2)
  rm = rm + .14763 * sin(2 * d + 2 * lp + pidiv2)
  rm = rm + .14368 * sin(2 * d + lp - 2 * l + pidiv2)
  rm = rm - .13922 * sin(4 * d + 2 * l + pidiv2)
  rm = rm - .13617 * sin(2 * lp - 2 * l + pidiv2)
  rm = rm - .13571 * sin(2 * d + lp - 2 * f + pidiv2)
  rm = rm - .12805 * sin(4 * d - 2 * lp + pidiv2)
  rm = rm + .11411 * sin(d - lp - l + pidiv2)
  rm = rm + .10998 * sin(d - lp + l + pidiv2)
  rm = rm - .10887 * sin(2 * d + 2 * lp - 2 * l + pidiv2)
  rm = rm - .10833 * sin(4 * d - 2 * lp - 2 * l + pidiv2)
  rm = rm - .10766 * sin(3 * d + lp + pidiv2)
  rm = rm - .10326 * sin(l + 2 * f + pidiv2)
  rm = rm - .09938 * sin(d - 3 * l + pidiv2)
  rm = rm - .08587 * sin(6 * d + pidiv2)
  rm = rm - .07982 * sin(4 * d - 2 * l - 2 * f + pidiv2)
  rm = rm - 6.678e-02 * sin(6 * d - lp - 2 * l + pidiv2)
  rm = rm - 6.545e-02 * sin(3 * d + l + pidiv2)
  rm = rm + .06055 * sin(d + l - 2 * f + pidiv2)
  rm = rm - .05904 * sin(d + lp + 2 * l + pidiv2)
  rm = rm - .05888 * sin(5 * l + pidiv2)
  rm = rm - .0585 * sin(2 * d - lp + 3 * l + pidiv2)
  rm = rm - .05789 * sin(4 * d - lp - 2 * f + pidiv2)
  rm = rm - .05527 * sin(2 * d + lp + l - 2 * f + pidiv2)
  rm = rm + .05293 * sin(3 * d - lp - 2 * l + pidiv2)
  rm = rm - .05191 * sin(6 * d - lp - l + pidiv2)
  rm = rm + .05072 * sin(2 * lp + 2 * l + pidiv2)
  rm = rm - .0502 * sin(lp - 2 * l + 2 * f + pidiv2)
  rm = rm - .04843 * sin(3 * d - 3 * l + pidiv2)
  rm = rm + .0474 * sin(2 * d - 5 * l + pidiv2)
  rm = rm - .04736 * sin(2 * d + lp - l - 2 * f + pidiv2)
  rm = rm - .04608 * sin(2 * d - 2 * lp + 2 * l + pidiv2)
  rm = rm + .04591 * sin(5 * d - 2 * l + pidiv2)
  rm = rm - .04422 * sin(2 * d + 4 * l + pidiv2)
  rm = rm - .04316 * sin(4 * d - lp - 3 * l + pidiv2)
  rm = rm - .04232 * sin(d - l - 2 * f + pidiv2)
  rm = rm - .03894 * sin(3 * lp - l + pidiv2)
  rm = rm + .0381 * sin(3 * d + lp - 2 * l + pidiv2)
  rm = rm + .03734 * sin(2 * d - lp - l + 2 * f + pidiv2)
  rm = rm + .03729 * sin(d + 2 * lp + pidiv2)
  rm = rm + .03682 * sin(4 * d + lp + l + pidiv2)
  rm = rm + .03379 * sin(d + lp - 2 * f + pidiv2)
  rm = rm + .03265 * sin(lp + 2 * l - 2 * f + pidiv2)
  rm = rm + .03143 * sin(2 * d + 2 * f + pidiv2)
  rm = rm + .03024 * sin(2 * d - lp - 2 * l + 2 * f + pidiv2)
  rm = rm - .02948 * sin(d - 2 * lp + pidiv2)
  rm = rm - .02939 * sin(4 * d - 4 * l + pidiv2)
  rm = rm + .0291 * sin(2 * d - 3 * l - 2 * f + pidiv2)
  rm = rm - .02855 * sin(2 * d - 3 * lp + l + pidiv2)
  rm = rm + .02839 * sin(2 * d - 2 * lp - 2 * f + pidiv2)
  rm = rm - .02698 * sin(4 * d - lp - l - 2 * f + pidiv2)
  rm = rm - .02674 * sin(lp - 4 * l + pidiv2)
  rm = rm + .02658 * sin(4 * d + 2 * lp - 2 * l + pidiv2)
  rm = rm - .02471 * sin(d - l + 2 * f + pidiv2)
  rm = rm - .02436 * sin(6 * d - lp - 3 * l + pidiv2)
  rm = rm - .02399 * sin(4 * d + lp - 3 * l + pidiv2)
  rm = rm + .02368 * sin(d + 3 * l + pidiv2)
  rm = rm + .02334 * sin(2 * d - lp - 4 * l + pidiv2)
  rm = rm + .02304 * sin(lp + 4 * l + pidiv2)
  rm = rm + .02127 * sin(3 * lp + pidiv2)
  rm = rm - .02079 * sin(4 * d - lp + 2 * l + pidiv2)
  rm = rm - .02008 * sin(2 * d - 3 * l + 2 * f + pidiv2)
  
  ' a(p,0,r) series
  
  rm = rm + 1.0587 * sin(2 * t - 2 * ju + 2 * d - l + 90.11969000000001 * dtr)
  rm = rm + .72783 * sin(18 * ve - 16 * t - 2 * l + 116.54311 * dtr)
  rm = rm + .68256 * sin(18 * ve - 16 * t + 296.54574 * dtr)
  rm = rm + .59827 * sin(3 * ve - 3 * t + 2 * d - l + 89.98187 * dtr)
  rm = rm + .45648 * sin(ll + l - f + 270.00126 * dtr)
  rm = rm + .45276 * sin(ll - l - f + 90.00128 * dtr)
  rm = rm + .41011 * sin(2 * t - 3 * ju + 2 * d - l + 280.06924 * dtr)
  rm = rm + .20497 * sin(t - ju - 2 * d + 91.79862 * dtr)
  rm = rm + .20473 * sin(18 * ve - 16 * t - 2 * d - l + 116.54222 * dtr)
  rm = rm + .20367 * sin(18 * ve - 16 * t + 2 * d - l + 296.54299 * dtr)
  rm = rm + .16644 * sin(2 * ve - 2 * t - 2 * d + 90.36386 * dtr)
  rm = rm + .1578 * sin(4 * t - 8 * ma + 3 * ju + l + 194.98833 * dtr)
  rm = rm + .1578 * sin(4 * t - 8 * ma + 3 * ju - l + 14.98841 * dtr)
  rm = rm + .15751 * sin(t - ju - 2 * d + l + 91.74578 * dtr)
  rm = rm + .1445 * sin(2 * t - 2 * ju - l + 89.97863 * dtr)
  rm = rm + .13811 * sin(ve - t - l + 270.00993 * dtr)
  rm = rm + .13477 * sin(18 * ve - 16 * t - 2 * d + 116.53978 * dtr)
  rm = rm + .12671 * sin(18 * ve - 16 * t + 2 * d - 2 * l + 296.54238 * dtr)
  rm = rm + .12666 * sin(t - ju - l + 91.22751 * dtr)
  rm = rm + .12362 * sin(ve - t - 2 * d + 269.98523 * dtr)
  rm = rm + .12047 * sin(2 * ve - 2 * t + 2 * d - l + 269.99692 * dtr)
  rm = rm + .11998 * sin(ve - t + l + 90.01606 * dtr)
  rm = rm + .11617 * sin(2 * ve - 2 * t - 2 * d + l + 90.31081 * dtr)
  rm = rm + .11256 * sin(4 * t - 8 * ma + 3 * ju + 2 * d - l + 197.11421 * dtr)
  rm = rm + .11251 * sin(4 * t - 8 * ma + 3 * ju - 2 * d + l + 17.11263 * dtr)
  rm = rm + .11226 * sin(4 * t - 8 * ma + 3 * ju + 2 * d + 196.69224 * dtr)
  rm = rm + .11216 * sin(4 * t - 8 * ma + 3 * ju - 2 * d + 16.68897 * dtr)
  rm = rm + .10689 * sin(ll + 2 * d - f + 270.00092 * dtr)
  rm = rm + .10504 * sin(t - ju + l + 271.06726 * dtr)
  rm = rm + .1006 * sin(ve - t - 2 * d + l + 269.98452 * dtr)
  rm = rm + 9.932e-02 * sin(3 * ve - 3 * t - l + 90.1054 * dtr)
  rm = rm + .09554 * sin(ll - 2 * d - f + 90.00096 * dtr)
  rm = rm + .08508 * sin(ll - f + 270.00061 * dtr)
  rm = rm + 7.9450e-02 * sin(4 * ve - 4 * t + 2 * d - l + 89.99224 * dtr)
  rm = rm + .07725 * sin(2 * t - 3 * ju - l + 280.16516 * dtr)
  rm = rm + 7.054e-02 * sin(6 * ve - 8 * t + 2 * d - l + 77.22087 * dtr)
  rm = rm + 6.313e-02 * sin(2 * ju - 5 * l + l + 256.56163 * dtr)
  
  ' a(p,1,r) series
  
  rm = rm + .51395 * t1 * sin(2 * d - lp + pidiv2)
  rm = rm + .38245 * t1 * sin(2 * d - lp - l + pidiv2)
  rm = rm + .32654 * t1 * sin(lp - l + pidiv2)
  rm = rm + .26396 * t1 * sin(lp + l + 270 * dtr)
  rm = rm + .12302 * t1 * sin(lp + 270 * dtr)
  rm = rm + .07754 * t1 * sin(2 * d + lp + 270 * dtr)
  rm = rm + .06068 * t1 * sin(2 * d + lp - l + 270 * dtr)
  rm = rm + .0497 * t1 * sin(2 * d - 2 * lp + pidiv2)
  rm = rm + .04194 * t1 * sin(d + lp + pidiv2)
  rm = rm + .03222 * t1 * sin(2 * d - lp + l + pidiv2)
  rm = rm + .02529 * t1 * sin(2 * d - lp - 2 * l + 270 * dtr)
  rm = rm + .0249 * t1 * sin(2 * d - 2 * lp - l + pidiv2)
  rm = rm + .00149 * t2 * sin(2 * d - lp + pidiv2)
  rm = rm + .00111 * t2 * sin(2 * d - lp - l + pidiv2)
  
  rmm = rm
  
  ' compute geocentric ecliptic longitude (arc seconds)
  
  ' a(c,0,v) series
  
  dv = 22639.58578 * sin(l)
  dv = dv + 4586.4383 * sin(2 * d - l)
  dv = dv + 2369.91394 * sin(2 * d)
  dv = dv + 769.02571 * sin(2 * l)
  dv = dv - 666.4171 * sin(lp)
  dv = dv - 411.59567 * sin(2 * f)
  dv = dv + 211.65555 * sin(2 * d - 2 * l)
  dv = dv + 205.43582 * sin(2 * d - lp - l)
  dv = dv + 191.9562 * sin(2 * d + l)
  dv = dv + 164.72851 * sin(2 * d - lp)
  dv = dv - 147.32129 * sin(lp - l)
  dv = dv - 124.98812 * sin(d)
  dv = dv - 109.38029 * sin(lp + l)
  dv = dv + 55.17705 * sin(2 * d - 2 * f)
  dv = dv - 45.0996 * sin(l + 2 * f)
  dv = dv + 39.53329 * sin(l - 2 * f)
  dv = dv + 38.42983 * sin(4 * d - l)
  dv = dv + 36.12381 * sin(3 * l)
  dv = dv + 30.77257 * sin(4 * d - 2 * l)
  dv = dv - 28.39708 * sin(2 * d + lp - l)
  dv = dv - 24.35821 * sin(2 * d + lp)
  dv = dv - 18.58471 * sin(d - l)
  dv = dv + 17.95446 * sin(d + lp)
  dv = dv + 14.53027 * sin(2 * d - lp + l)
  dv = dv + 14.3797 * sin(2 * d + 2 * l)
  dv = dv + 13.89906 * sin(4 * d)
  dv = dv + 13.19406 * sin(2 * d - 3 * l)
  dv = dv - 9.67905 * sin(lp - 2 * l)
  dv = dv - 9.36586 * sin(2 * d - l + 2 * f)
  dv = dv + 8.60553 * sin(2 * d - lp - 2 * l)
  dv = dv - 8.45310 * sin(d + l)
  dv = dv + 8.05016 * sin(2 * d - 2 * lp)
  dv = dv - 7.63015 * sin(lp + 2 * l)
  dv = dv - 7.44749 * sin(2 * lp)
  dv = dv + 7.37119 * sin(2 * d - 2 * lp - l)
  dv = dv - 6.38315 * sin(2 * d + l - 2 * f)
  dv = dv - 5.74161 * sin(2 * d + 2 * f)
  dv = dv + 4.37401 * sin(4 * d - lp - l)
  dv = dv - 3.99761 * sin(2 * l + 2 * f)
  dv = dv - 3.20969 * sin(3 * d - l)
  dv = dv - 2.91454 * sin(2 * d + lp + l)
  dv = dv + 2.73189 * sin(4 * d - lp - 2 * l)
  dv = dv - 2.56794 * sin(2 * lp - l)
  dv = dv - 2.5212 * sin(2 * d + 2 * lp - l)
  dv = dv + 2.48889 * sin(2 * d + lp - 2 * l)
  dv = dv + 2.14607 * sin(2 * d - lp - 2 * f)
  dv = dv + 1.97773 * sin(4 * d + l)
  dv = dv + 1.93368 * sin(4 * l)
  dv = dv + 1.87076 * sin(4 * d - lp)
  dv = dv - 1.75297 * sin(d - 2 * l)
  dv = dv - 1.43716 * sin(2 * d + lp - 2 * f)
  dv = dv - 1.37257 * sin(2 * l - 2 * f)
  dv = dv + 1.26182 * sin(d + lp + l)
  dv = dv - 1.22412 * sin(3 * d - 2 * l)
  dv = dv + 1.18683 * sin(4 * d - 3 * l)
  dv = dv + 1.177 * sin(2 * d - lp + 2 * l)
  dv = dv - 1.16169 * sin(2 * lp + l)
  dv = dv + 1.07769 * sin(d + lp - l)
  dv = dv + 1.0595 * sin(2 * d + 3 * l)
  dv = dv - .99022 * sin(2 * d + l + 2 * f)
  dv = dv + .94828 * sin(2 * d - 4 * l)
  dv = dv + .75168 * sin(2 * d - 2 * lp + l)
  dv = dv - .66938 * sin(lp - 3 * l)
  dv = dv - .63521 * sin(4 * d + lp - l)
  dv = dv - .58399 * sin(d + 2 * l)
  dv = dv - .58331 * sin(d - 2 * f)
  dv = dv + .57156 * sin(6 * d - 2 * l)
  dv = dv - .56064 * sin(2 * d - 2 * l - 2 * f)
  dv = dv - .55692 * sin(d - lp)
  dv = dv - .54592 * sin(lp + 3 * l)
  dv = dv - .53571 * sin(2 * d - 2 * l + 2 * f)
  dv = dv + .4784 * sin(2 * d - lp - 3 * f)
  dv = dv - .45379 * sin(2 * d + 2 * l - 2 * f)
  dv = dv - .42622 * sin(2 * d - lp - l + 2 * f)
  dv = dv + .42033 * sin(4 * f)
  dv = dv + .4134 * sin(lp + 2 * f)
  dv = dv + .40423 * sin(3 * d)
  dv = dv + .39451 * sin(6 * d - l)
  dv = dv - .38213 * sin(2 * d - lp + 2 * f)
  dv = dv - .37451 * sin(2 * d - lp + l - 2 * f)
  dv = dv - .35758 * sin(4 * d + lp - 2 * l)
  dv = dv + .34965 * sin(d + lp - 2 * l)
  dv = dv + .33979 * sin(2 * d - 3 * lp)
  dv = dv - .32866 * sin(3 * l + 2 * f)
  dv = dv + .30872 * sin(4 * d - 2 * lp - l)
  dv = dv + .30155 * sin(lp - l - 2 * f)
  dv = dv + .30086 * sin(4 * d - l - 2 * f)
  dv = dv + .2942 * sin(2 * d - 2 * lp - 2 * l)
  dv = dv + .29255 * sin(6 * d - 3 * l)
  dv = dv - .29022 * sin(2 * d + lp + 2 * l)
  
  ' a(p,2,v) series
  
  dv = dv + .00487 * t2 * sin(lp)
  dv = dv - .0015 * t2 * sin(2 * d - lp - l + pi)
  dv = dv - .0012 * t2 * sin(2 * d - lp + pi)
  dv = dv + .00108 * t2 * sin(lp - l)
  dv = dv + .0008 * t2 * sin(lp + l)
  
  ' a(p,0,v) series
  
  dv = dv + 14.24883 * sin(18 * ve - 16 * t - l + dtr * 26.54261)
  dv = dv + 7.06304 * sin(ll - f + dtr * .00094)
  dv = dv + 1.14307 * sin(2 * t - 2 * ju + 2 * d - l + dtr * 180.11977)
  dv = dv + .901140 * sin(4 * t - 8 * ma + 3 * ju + dtr * 285.98707)
  dv = dv + .82155 * sin(ve - t + dtr * 180.00988)
  dv = dv + .78811 * sin(18 * ve - 16 * t - 2 * l + dtr * 26.54324)
  dv = dv + .7393 * sin(18 * ve - 16 * t + dtr * 26.5456)
  dv = dv + .64371 * sin(3 * ve - 3 * t + 2 * d - l + dtr * 179.98144)
  dv = dv + .6388 * sin(t - ju + dtr * 1.2289)
  dv = dv + .56341 * sin(10 * ve - 3 * t - l + dtr * 333.30551)
  dv = dv + .49331 * sin(ll + l - f + .00127 * dtr)
  dv = dv + .49141 * sin(ll - l - f + .00127 * dtr)
  dv = dv + .44532 * sin(2 * t - 3 * ju + 2 * d - l + 10.07001 * dtr)
  dv = dv + .36061 * sin(ll + f + .00071 * dtr)
  dv = dv + .34355 * sin(2 * ve - 3 * t + 269.95393 * dtr)
  dv = dv + .32455 * sin(t - 2 * ma + 318.13776 * dtr)
  dv = dv + .30155 * sin(2 * ve - 2 * t + .20448 * dtr)
  dv = dv + .28938 * sin(t + d - f + 95.13523 * dtr)
  dv = dv + .28281 * sin(2 * t - 3 * ju + 2 * d - 2 * l + 10.03835 * dtr)
  dv = dv + .24515 * sin(2 * t - 2 * ju + 2 * d - 2 * l + .08642 * dtr)
  
  ' a(p,1,v) series
  
  dv = dv + 1.6768 * t1 * sin(lp)
  dv = dv + .51642 * t1 * sin(2 * d - lp - l + pi)
  dv = dv + .41383 * t1 * sin(2 * d - lp + pi)
  dv = dv + .37115 * t1 * sin(lp - l)
  dv = dv + .2756 * t1 * sin(lp + l)
  dv = dv + .25425 * t1 * sin(18 * ve - 16 * t - l + 114.5655 * dtr)
  dv = dv + 7.1178e-02 * t1 * sin(2 * d + lp - l)
  dv = dv + .06128 * t1 * sin(2 * d + lp)
  dv = dv + .04516 * t1 * sin(d + lp + pi)
  dv = dv + .04048 * t1 * sin(2 * d - 2 * lp + pi)
  dv = dv + .03747 * t1 * sin(2 * lp)
  dv = dv + .03707 * t1 * sin(2 * d - 2 * lp - l + pi)
  dv = dv + .03649 * t1 * sin(2 * d - lp + l + pi)
  dv = dv + .02438 * t1 * sin(lp - 2 * l)
  dv = dv + .02165 * t1 * sin(2 * d - lp - 2 * l + pi)
  dv = dv + .01923 * t1 * sin(lp + 2 * l)
  
  plon = modulo(ll + atr * dv + dpsi)
  
  ' compute geocentric ecliptic latitude (arc seconds)
  
  ' a(c,0,u) series
  
  pl = 18461.23868 * sin(f)
  pl = pl + 1010.16707 * sin(l + f)
  pl = pl + 999.69358 * sin(l - f)
  pl = pl + 623.65243 * sin(2 * d - f)
  pl = pl + 199.48374 * sin(2 * d - l + f)
  pl = pl + 166.5741 * sin(2 * d - l - f)
  pl = pl + 117.26069 * sin(2 * d + f)
  pl = pl + 61.91195 * sin(2 * l + f)
  pl = pl + 33.3572 * sin(2 * d + l - f)
  pl = pl + 31.75967 * sin(2 * l - f)
  pl = pl + 29.57658 * sin(2 * d - lp - f)
  pl = pl + 15.56626 * sin(2 * d - 2 * l - f)
  pl = pl + 15.12155 * sin(2 * d + l + f)
  pl = pl - 12.09414 * sin(2 * d + lp - f)
  pl = pl + 8.86814 * sin(2 * d - lp - l + f)
  pl = pl + 7.95855 * sin(2 * d - lp + f)
  pl = pl + 7.43455 * sin(2 * d - lp - l - f)
  pl = pl - 6.73143 * sin(lp - l - f)
  pl = pl + 6.57957 * sin(4 * d - l - f)
  pl = pl - 6.46007 * sin(lp + f)
  pl = pl - 6.29648 * sin(3 * f)
  pl = pl - 5.63235 * sin(lp - l + f)
  pl = pl - 5.3684 * sin(d + f)
  pl = pl - 5.31127 * sin(lp + l + f)
  pl = pl - 5.07591 * sin(lp + l - f)
  pl = pl - 4.83961 * sin(lp - f)
  pl = pl - 4.80574 * sin(d - f)
  pl = pl + 3.98405 * sin(3 * l + f)
  pl = pl + 3.67446 + sin(4 * d - f)
  pl = pl + 2.99848 * sin(4 * d - l + f)
  pl = pl + 2.79864 * sin(l - 3 * f)
  pl = pl + 2.41388 * sin(4 * d - 2 * l + f)
  pl = pl + 2.18631 * sin(2 * d - 3 * f)
  pl = pl + 2.14617 * sin(2 * d + 2 * l - f)
  pl = pl + 1.76598 * sin(2 * d - lp + l - f)
  pl = pl - 1.62442 * sin(2 * d - 2 * l + f)
  pl = pl + 1.5813 * sin(3 * l - f)
  pl = pl + 1.51975 * sin(2 * d + 2 * l + f)
  pl = pl - 1.51563 * sin(2 * d - 3 * l - f)
  pl = pl - 1.31782 * sin(2 * d + lp - l + f)
  pl = pl - 1.26427 * sin(2 * d + lp + f)
  pl = pl + 1.19187 * sin(4 * d + f)
  pl = pl + 1.13461 * sin(2 * d - lp + l + f)
  pl = pl + 1.08578 * sin(2 * d - 2 * lp - f)
  pl = pl - 1.01938 * sin(l + 3 * f)
  pl = pl - .822710 * sin(2 * d + lp + l - f)
  pl = pl + .80422 * sin(d + lp - f)
  pl = pl + .80259 * sin(d + lp + f)
  pl = pl - .79319 * sin(lp - 2 * l - f)
  pl = pl - .79101 * sin(2 * d + lp - l - f)
  pl = pl - .66741 * sin(d + l + f)
  pl = pl + .65022 * sin(2 * d - lp - 2 * l - f)
  pl = pl - .63881 * sin(lp + 2 * l + f)
  pl = pl + .63371 * sin(4 * d - 2 * l - f)
  pl = pl + .59577 * sin(4 * d - lp - l - f)
  pl = pl - .58893 * sin(d + l - f)
  pl = pl + .47338 * sin(4 * d + l - f)
  pl = pl - .42989 * sin(d - l - f)
  pl = pl + .41494 * sin(4 * d - lp - f)
  pl = pl + .3835 * sin(2 * d - 2 * lp + f)
  pl = pl - .35183 * sin(3 * d - f)
  pl = pl + .33881 * sin(4 * d - lp - l + f)
  pl = pl + .32906 * sin(2 * d - l - 3 * f)
  
  ' a(p,0,u) series
  
  pl = pl + 8.04508 * sin(ll + 180.00071 * dtr)
  pl = pl + 1.51021 * sin(t + d + 276.68007 * dtr)
  pl = pl + .63037 * sin(18 * ve - 16 * t - l + f + 26.54287 * dtr)
  pl = pl + .63014 * sin(18 * ve - 16 * t - l - f + 26.54272 * dtr)
  pl = pl + .45586 * sin(ll - l + .00075 * dtr)
  pl = pl + .41571 * sin(ll + l + 180.00069 * dtr)
  pl = pl + .32622 * sin(ll - 2 * f + .00086 * dtr)
  pl = pl + .29854 * sin(ll - 2 * d + .00072 * dtr)
  
  ' a(p,1,u) series
  
  pl = pl + .0743 * t1 * sin(2 * d - lp - f + pi)
  pl = pl + .03043 * t1 * sin(2 * d + lp - f)
  pl = pl + .02229 * t1 * sin(2 * d - lp - l + f + pi)
  pl = pl + .01999 * t1 * sin(2 * d - lp + f + pi)
  pl = pl + .01869 * t1 * sin(2 * d - lp - l - f + pi)
  pl = pl + .01696 * t1 * sin(lp - l - f)
  pl = pl + .01623 * t1 * sin(lp + f)
  
  plat = atr * pl
  
  ' compute geocentric right ascension and declination
  
  a = sin(plon) * cos(obliq) - tan(plat) * sin(obliq)
  
  b = cos(plon)
  
  rasc = atan3(a, b)
  
  decl = asin(sin(plat) * cos(obliq) + cos(plat) * sin(obliq) * sin(plon))
  
  ' compute the geocentric position vector of the moon (au)
  
  rmoon(1) = rmm * cos(rasc) * cos(decl)
  
  rmoon(2) = rmm * sin(rasc) * cos(decl)
  
  rmoon(3) = rmm * sin(decl)
  
end sub
  
  '''''''''''''''''''''''''''''''''''''''''''
  '''''''''''''''''''''''''''''''''''''''''''
  
sub us76(alt, atmtmp, atmdns, atmprs, atmsos)
  
  ' U.S. 1976 standard atmosphere subroutine
  
  ' input
  
  '  alt = altitude (kilometers)
  
  ' output
  
  '  atmtmp = temperature (degrees Kelvin)
  '  atmdns = density (kilograms/cubic kilometer)
  '  atmprs = pressure (newton/square meter)
  '  atmsos = speed of sound (meters/second)
  
  ''''''''''''''''''''''''''''''''''''''''''
  
  local i%, j%
  
  z = alt
  
  ' altitude must be less than 1000 kilometers
  
  if (z > 1000.0) then
    
    z = 999.0
    
  end if
  
  ' determine positive altitude values
  
  if (z >= 0.0) then
    
    ' set initial values
    
    ro = 6356.766
    
    go = 9.80665
    
    wmo = 28.9644
    
    rs = 8314.32
    
    zm = z * 1000.0
    
    rom = ro * 1000.0
    
    ' ------------------------------------
    ' altitudes greater than 86 kilometers
    ' ------------------------------------
    
    if (z >= 86.0) then
      
      ' find pointer location
      
      for i% = 8 to 48
        
        if (zs(i%) <= z and z < zs(i% + 1)) then
          
          exit for
          
        end if
        
      next i%
      
      ' calculate temperature
      
      if (i% = 8) then
        
        t = tms(9)
        
      elseif (not (i% < 16 or i% >= 19)) then
        
        t = 240.0 + 12.0 * (z - 110.0)
        
      elseif (i% >= 19) then
        
        xi = (z - 120.0) * (ro + 120.0) / (ro + z)
        
        t = 1000.0 - 640.0 * exp(-0.01875 * xi)
        
      else
        
        t = 263.1905 - 76.3232 * sqr(1.0 - ((z - 91.0) / 19.9429)^2)
        
      end if
      
      ' set pointers
      
      j% = i%
      
      if (i% = 48) then
        
        j% = i% - 1
        
      end if
      
      ' calculate altitude indices values for molecular weight calculation
      
      z0 = zs(j%)
      
      z1 = zs(j% + 1)
      
      z2 = zs(j% + 2)
      
      ' calculate molecular weight
      
      wma = wms(j%) * (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2))
      
      wma = wma + wms(j% + 1) * (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2))
      
      wma = wma + wms(j% + 2) * (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1))
      
      ' calculate atmospheric pressure indices
      
      alp0 = log(ps(j%))
      
      alp1 = log(ps(j% + 1))
      
      alp2 = log(ps(j% + 2))
      
      ' calculate pressure
      
      alpa = alp0 * (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2))
      
      alpa = alpa + alp1 * (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2))
      
      alpa = alpa + alp2 * (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1))
      
      ' set second side of molecular weight and pressure value
      
      alpb = alpa
      
      wmb = wma
      
      ' perform same calculations for far side values for smoothing
      
      if (not (i% = 8 or i% = 48)) then
        
        j% = j% - 1
        
        ' calculate altitude indices values for molecular weight calculation
        
        z0 = zs(j%)
        
        z1 = zs(j% + 1)
        
        z2 = zs(j% + 2)
        
        ' calculate pressure indices values for pressure calculation
        
        alp0 = log(ps(j%))
        
        alp1 = log(ps(j% + 1))
        
        alp2 = log(ps(j% + 2))
        
        ' calculate atmospheric pressure
        
        alpb = alp0 * (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2))
        
        alpb = alpb + alp1 * (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2))
        
        alpb = alpb + alp2 * (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1))
        
        ' calculate molecular weight
        
        wmb = wms(j%) * (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2))
        
        wmb = wmb + wms(j% + 1) * (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2))
        
        wmb = wmb + wms(j% + 2) * (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1))
        
      end if
      
      ' calculate pressure and molecular weight from derived constants
      
      p = 100.0 * exp((alpa + alpb) / 2.0)
      
      wm = (wma + wmb) / 2.0
      
    else
      
      ' find pointer location
      
      for i% = 1 to 7
        
        if (zs(i%) <= z and z < zs(i% + 1)) then
          
          exit for
          
        end if
        
      next i%
      
      ' calculate altitude breakpoints
      
      zl = ro * zs(i%) / (ro + zs(i%))
      
      zu = ro * zs(i% + 1) / (ro + zs(i% + 1))
      
      ' convert to meters
      
      zlm = zl * 1000.0
      
      zum = zu * 1000.0
      
      wm = wmo
      
      ht = (ro * z) / (ro + z)
      
      hm = ht * 1000.0
      
      g = (tms(i% + 1) - tms(i%)) / (zu - zl)
      
      gm = g * 0.001
      
      ' calculate pressure and temperature
      
      if (g < 0.0 or g > 0.0) then
        
        p = ps(i%) * ((tms(i%) / (tms(i%) + g * (ht - zl)))^((go * wmo) / (rs * gm))) * 100.0
        
      else
        
        p = ps(i%) * exp(-(go * wmo * (hm - zlm)) / (rs * tms(i%))) * 100.0
        
      end if
      
      t = tms(i%) + g * (ht - zl)
      
    end if
    
    ' calculate density
    
    d = (wm * p) / (rs * t)
    
  else
    
    ' set default values for negative altitudes
    
    t = 288.15
    
    p = 101325.0
    
    d = 1.224999155888
    
  end if
  
  ' calculate speed of sound
  
  if (z < 91.0) then
    
    ss = 20.046707 * sqr(t)
    
  else
    
    ss = 274.0392912
    
  end if
  
  atmdns = d * (1000.0)^3
  
  atmtmp = t
  
  atmsos = ss
  
  atmprs = p
  
end sub
  
  ''''''''''''''''''''''''''''''''
  ''''''''''''''''''''''''''''''''
  
sub realroot(x1, x2, xroot, froot)
  
  ' real root of a single non-linear function subroutine
  
  ' input
  
  '  x1  = lower bound of search interval
  '  x2  = upper bound of search interval
  
  ' output
  
  '  xroot = real root of f(x) = 0
  '  froot = function value
  
  ' note: requires sub npoe_func
  
  ''''''''''''''''''''''''''''''
  
  local eps, a, b, c, d, e, fa, fb, fcc, tol1
  
  local xm, p, q, r, s, xmin, tmp
  
  eps = 2.23e-16
  
  e = 0.0
  
  a = x1
  
  b = x2
  
  npoe_func(a, fa)
  
  npoe_func(b, fb)
  
  fcc = fb
  
  for iter% = 1 to 50
    
    if (fb * fcc > 0.0) then
      
      c = a
      
      fcc = fa
      
      d = b - a
      
      e = d
      
    end if
    
    if (abs(fcc) < abs(fb)) then
      
      a = b
      
      b = c
      
      c = a
      
      fa = fb
      
      fb = fcc
      
      fcc = fa
      
    end if
    
    tol1 = 2.0 * eps * abs(b) + 0.5 * rtol
    
    xm = 0.5 * (c - b)
    
    if (abs(xm) <= tol1 or fb = 0.0) then exit for
    
    if (abs(e) >= tol1 and abs(fa) > abs(fb)) then
      
      s = fb / fa
      
      if (a = c) then
        
        p = 2.0 * xm * s
        
        q = 1.0 - s
        
      else
        
        q = fa / fcc
        
        r = fb / fcc
        
        p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0))
        
        q = (q - 1.0) * (r - 1.0) * (s - 1.0)
        
      end if
      
      if (p > 0.0) then q = -q
      
      p = abs(p)
      
      min = abs(e * q)
      
      tmp = 3.0 * xm * q - abs(tol1 * q)
      
      if (min < tmp) then min = tmp
      
      if (2.0 * p < min) then
        
        e = d
        
        d = p / q
        
      else
        
        d = xm
        
        e = d
        
      end if
      
    else
      
      d = xm
      
      e = d
      
    end if
    
    a = b
    
    fa = fb
    
    if (abs(d) > tol1) then
      
      b = b + d
      
    else
      
      b = b + sgn(xm) * tol1
      
    end if
    
    npoe_func(b, fb)
    
  next iter%
  
  froot = fb
  
  xroot = b
  
end sub
  
  ''''''''''''''''''''''''''''''''
  ''''''''''''''''''''''''''''''''
  
sub julian(month, day, year, jday)
  
  ' Gregorian date to julian day subroutine
  
  ' input
  
  '  month = calendar month
  '  day   = calendar day
  '  year  = calendar year (all four digits)
  
  ' output
  
  '  jday = julian day
  
  ' special notes
  
  '  (1) calendar year must include all digits
  
  '  (2) will report October 5, 1582 to October 14, 1582
  '      as invalid calendar dates and exit
  
  '''''''''''''''''''''''''''''''''''''''''
  
  local a, b, c, m, y
  
  y = year
  
  m = month
  
  b = 0.0
  
  c = 0.0
  
  if (m <= 2.0) then
    
    y = y - 1.0
    
    m = m + 12.0
    
  end if
  
  if (y < 0.0) then c = -0.75
  
  if (year < 1582.0) then
    
    ' null
    
  elseif (year > 1582.0) then
    
    a = fix(y / 100.0)
    
    b = 2.0 - a + fix(a / 4.0)
    
  elseif (month < 10.0) then
    
    ' null
    
  elseif (month > 10.0) then
    
    a = fix(y / 100.0)
    
    b = 2.0 - a + fix(a / 4.0)
    
  elseif (day <= 4.0) then
    
    ' null
    
  elseif (day > 14.0) then
    
    a = fix(y / 100.0)
    
    b = 2.0 - a + fix(a / 4.0)
    
  else
    
    print "this date does not exist!!"
    
    exit
    
  end if
  
  jday = fix(365.25 * y + c) + fix(30.6001 * (m + 1.0)) + day + b + 1720994.5
  
end sub
  
  ''''''''''''''''''''''''''''''''
  ''''''''''''''''''''''''''''''''
  
sub gdate (jday, month, day, year)
  
  ' Julian day to Gregorian date subroutine
  
  ' input
  
  '  jday = julian day
  
  ' output
  
  '  month = calendar month
  '  day   = calendar day
  '  year  = calendar year
  
  ''''''''''''''''''''''''
  
  local a, b, c, d, e, f, z, alpha
  
  z = fix(jday + 0.5)
  
  f = jday + 0.5 - z
  
  if (z < 2299161.0) then
    
    a = z
    
  else
    
    alpha = fix((z - 1867216.25) / 36524.25)
    
    a = z + 1.0 + alpha - fix(alpha / 4.0)
    
  end if
  
  b = a + 1524.0
  
  c = fix((b - 122.1) / 365.25)
  
  d = fix(365.25 * c)
  
  e = fix((b - d) / 30.6001)
  
  day = b - d - fix(30.6001 * e) + f
  
  if (e < 13.5) then
    
    month = e - 1.0
    
  else
    
    month = e - 13.0
    
  end if
  
  if (month > 2.5) then
    
    year = c - 4716.0
    
  else
    
    year = c - 4715.0
    
  end if
  
end sub
  
  '''''''''''''''
  '''''''''''''''
  
sub jd2str(jdutc)
  
  ' convert julian day to calendar date and UTC time
  
  ''''''''''''''''''''''''''''''''''''''''''''''''''
  
  gdate (jdutc, cmonth, day, year)
  
  print "calendar date      ", month$(cmonth); " "; STR$(int(day)); " "; str$(year)
  
  print " "
  
  thr0 = 24.0 * (day - int(day))
  
  thr = int(thr0)
  
  tmin0 = 60.0 * (thr0 - thr)
  
  tmin = int(tmin0)
  
  tsec = 60.0 * (tmin0 - tmin)
  
  ' fix seconds and minutes for rollover
  
  if (tsec >= 59.9) then
    
    tsec = 0.0
    
    tmin = tmin + 1.0
    
  end if
  
  ' fix minutes for rollover
  
  if (tmin >= 60.0) then
    
    tmin = 0.0
    
    thr = thr + 1.0
    
  end if
  
  print "UTC                ", str$(thr) + " hours " + str$(tmin) + " minutes " + str$(tsec, 0, 2) + " seconds"
  
end sub
  
  '''''''''''''''''''''''''''''''
  '''''''''''''''''''''''''''''''
  
sub nut2000_lp(jdate, dpsi, deps)
  
  ' low precison nutation based on iau 2000a subroutine
  
  ' this function evaluates a short nutation series and returns approximate
  ' values for nutation in longitude and nutation in obliquity for a given
  ' tdb julian date. in this mode, only the largest 13 terms of the iau 2000a
  ' nutation series are evaluated.
  
  ' input
  
  '  jdate = tdb julian date
  
  ' output
  
  '  dpsi = nutation in longitude in arcseconds
  
  '  deps = nutation in obliquity in arcseconds
  
  '''''''''''''''''''''''''''''''''''''''''''''
  
  LOCAL rev = 1296000.0
  
  LOCAL el, elp, f, d, omega
  
  LOCAL i as integer, arg
  
  LOCAL t = (jdate - 2451545.0) / 36525.0
  
  ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  ' computation of fundamental (delaunay) arguments from simon et al. (1994)
  ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
  
  ' mean anomaly of the moon
  
  el = (485868.249036 + t * (1717915923.2178 + t * (31.8792 + t * (0.051635 + t * (-0.00024470)))) mod rev) / seccon
  
  ' mean anomaly of the sun
  
  elp = (1287104.79305 + t * (129596581.0481 + t * (-0.5532 + t * (0.000136 + t * (- 0.00001149)))) mod rev) / seccon
  
  ' mean argument of the latitude of the moon
  
  f = (335779.526232 + t * (1739527262.8478 + t * (-12.7512 + t * (-0.001037 + t * (0.00000417)))) mod rev) / seccon
  
  ' mean elongation of the moon from the sun
  
  d = (1072260.70369 + t * (1602961601.2090 + t * (- 6.3706 + t * (0.006593 + t * (- 0.00003169)))) mod rev) / seccon
  
  ' mean longitude of the ascending node of the moon (from simon section 3.4(b.3), precession = 5028.8200 arcsec/cy)
  
  omega = (450160.398036 + t * (- 6962890.5431 + t * (7.4722 + t * (0.007702 + t * (- 0.00005939)))) mod rev) / seccon
  
  dpsi = 0.0
  
  deps = 0.0
  
  ' sum nutation series terms
  
  for i = 13 to 1 step -1
    
    arg = xnut(1, i) * el + xnut(2, i) * elp + xnut(3, i) * f + xnut(4, i) * d + xnut(5, i) * omega
    
    dpsi = (xnut(6, i) + xnut(7, i) * t) * sin(arg) + xnut(10, i) * cos(arg) + dpsi
    
    deps = (xnut(8, i) + xnut(9, i) * t) * cos(arg) + xnut(11, i) * sin(arg) + deps
    
  next i
  
  dpsi = 1.0e-7 * dpsi
  
  deps = 1.0e-7 * deps
  
  ' add in out-of-phase component of principal (18.6-year) term
  ' (to avoid small but long-term bias in results)
  
  dpsi = dpsi + 0.0033 * cos(omega)
  
  deps = deps + 0.0015 * sin(omega)
  
END sub
  
  ''''''''''''''''''''''''''''''''''''
  ''''''''''''''''''''''''''''''''''''
  
sub obliq_lp(jdate, dpsi, deps, obliq)
  
  ' nutations and true obliquity subroutine
  
  ' input
  
  '  jdate = julian date
  
  ' output
  
  '  dpsi = nutation in longitude in radians
  
  '  deps = nutation in obliquity in radians
  
  '  obliq = true obliquity of the ecliptic in radians
  
  ''''''''''''''''''''''''''''''''''''''''''''''''''''
  
  LOCAL t, t2, t3, eqeq
  
  LOCAL th, tl, obm, obt, st
  
  LOCAL tjdh, tjdl
  
  ' split the julian date
  
  tjdh = int(jdate)
  
  tjdl = jdate - tjdh
  
  ' fundamental time units
  
  th = (tjdh - 2451545.0) / 36525.0
  
  tl = tjdl / 36525.0
  
  t = th + tl
  
  t2 = t * t
  
  t3 = t2 * t
  
  ' obtain equation of the equinoxes
  
  eqeq = 0.0
  
  ' obtain nutations
  
  nut2000_lp(jdate, dpsi, deps)
  
  ' compute mean obliquity of the ecliptic in seconds of arc
  
  obm = 84381.4480 - 46.8150 * t - 0.00059 * t2 + 0.001813 * t3
  
  ' compute true obliquity of the ecliptic in seconds of arc
  
  obt = obm + deps
  
  ' return elements in radians
  
  deps = atr * deps
  
  dpsi = atr * dpsi
  
  obliq = atr * obt
  
END sub
  
  ''''''''''''''''''''''''''''''
  ''''''''''''''''''''''''''''''
  
sub geodet1(rmag, dec, alt, lat)
  
  ' geodetic latitude and altitude subroutine
  
  ' series solution
  
  ' input
  
  '  rmag = geocentric radius (kilometers)
  '  dec  = geocentric declination (radians)
  '         (+north, -south -pi/2 <= dec <= +pi/2)
  
  ' output
  
  '  alt = geodetic altitude (kilometers)
  '  lat = geodetic latitude (radians)
  '        (+north, -south -pi/2 <= lat <= +pi/2)
  
  ' global constants
  
  '  req  = equatorial radius (kilometers)
  '  flat = flattening factor (non-dimensional)
  
  '  reference: NASA TN D-7522
  
  ''''''''''''''''''''''''''''
  
  local a, n, o, p, q, r, s
  
  n = req / rmag
  o = flat * flat
  
  a = 2.0 * dec
  p = sin(a)
  q = cos(a)
  
  a = 4.0 * dec
  r = sin(a)
  s = cos(a)
  
  ' geodetic latitude (radians)
  
  lat = dec + flat * n * p + o * n * r * (n - 0.25)
  
  ' geodetic altitude (kilometers)
  
  alt = rmag + req * (flat * 0.5 * (1.0 - q) + o * (0.25 * n - 0.0625) * (1.0 - s) - 1.0)
  
end sub
  
  ''''''''''''''''''''''''''''''''''''''''''''''
  ''''''''''''''''''''''''''''''''''''''''''''''
  
sub tevent(sma, ecc, argper, tanom, ttpp, ttanc)
  
  ' trajectory event times subroutine
  
  ' input
  
  '  sma    = semimajor axis (kilometers)
  '  ecc    = orbital eccentricity
  '  argper = argument of perigee (radians)
  '  tanom  = true anomaly (radians)
  
  ' output
  
  '  ttpp  = time until perigee passage (seconds)
  '  ttanc = time until ascending node crossing (seconds)
  
  '''''''''''''''''''''''''''''''''''''''''''''''''''''''
  
  ' Keplerian period
  
  period = pi2 * sma * sqr(sma / emu)
  
  a = sin(tanom) * sqr(1.0 - ecc * ecc)
  b = ecc + cos(tanom)
  
  ' eccentric anomaly
  
  eanom2 = atan3(a, b)
  
  c = -eanom2 + ecc * sin(eanom2)
  
  ' time until perigee passage
  
  ttpp = c * period / pi2
  
  if (ttpp < 0.0) then
    
    ttpp = period + ttpp
    
  end if
  
  ' time until ascending node crossing
  
  d = pi2 - argper
  
  a = sin(d) * sqr(1.0 - ecc * ecc)
  b = ecc + cos(d)
  
  eanom1 = atan3(a, b)
  
  c = eanom1 - eanom2 - ecc * (sin(eanom1) - sin(eanom2))
  
  ttanc = c * period / pi2
  
  if (ttanc < 0.0) then ttanc = period + ttanc
  
end sub
  
  '''''''''''''''''''''''''''''''''''''''''''''''''
  '''''''''''''''''''''''''''''''''''''''''''''''''
  
sub tperiod(sma, ecc, inc, argper, tnodal, tanomal)
  
  ' orbital periods subroutine
  
  ' input
  
  '  sma    = semimajor axis
  '  ecc    = orbital eccentricity
  '  inc    = orbital inclination (radians)
  '  argper = argument of perigee (radians)
  
  ' output
  
  '  tnodal  = nodal period (seconds)
  '  tanomal = anomalistic period (seconds)
  
  '''''''''''''''''''''''''''''''''''''''''
  
  local a, b, e, ar, ep, sw, pp
  
  tkepler = pi2 * sma * sqr(sma / emu)
  
  e = sin(inc) * sin(inc)
  
  ar = (sma / req) * (sma / req)
  
  ep = 1.0 - ecc * ecc
  
  sw = sin(argper)
  
  pp = 1.0 + ecc * cos(argper)
  
  a = -0.75 * xj2 * (4.0 - 5.0 * e) / (ar * sqr(ep) * pp * pp)
  
  b = -1.5 * xj2 * pp * pp * pp / (ar * ep * ep * ep)
  
  ' nodal period
  
  tnodal = tkepler * (1.0 + a + b)
  
  ' anomalistic period
  
  tanomal = tkepler * (1.0 - 1.5 * xj2 * (1.0 - 3.0 * e * sw * sw) / (ar * (1.0 - ecc) ^ 3))
  
end sub
  
  ''''''''''''''''
  ''''''''''''''''
  
sub oeprint(oev())
  
  ' print classical orbital elements subroutine
  
  if (oev(1) > 0.0) then
    
    tau = pi2 * sqr(oev(1) * oev(1) * oev(1) / emu)
    
  else
    
    tau = 0.0
    
  end if
  
  print "semimajor axis          ", str$(oev(1), 0, 4) + " kilometers"
  print "eccentricity            ", str$(oev(2), 0, 10)
  print "orbital inclination     ", str$(rtd * oev(3), 0, 4) + " degrees"
  print "argument of perigee     ", str$(rtd * oev(4), 0, 4) + " degrees"
  print "argument of latitude    ", str$(rtd * oev(8), 0, 4) + " degrees"
  print "RAAN                    ", str$(rtd * oev(5), 0, 4) + " degrees"
  print "true anomaly            ", str$(rtd * oev(6), 0, 4) + " degrees"
  print " "
  print "Keplerian period        ", str$(tau / 60.0, 0, 4) + " minutes"
  print "nodal period            ", str$(oev(31) / 60.0, 0, 4) + " minutes"
  print "anomalistic period      ", str$(oev(32) / 60.0, 0, 4) + " minutes"
  print " "
  print "geodetic latitude       ", str$(rtd * oev(14), 0, 4) + " degrees"
  print "geodetic altitude       ", str$(oev(16), 0, 4) + " kilometers"
  print " "
  print "flight path angle       ", str$(rtd * oev(11), 0, 4) + " degrees"
  print " "
  print "right ascension         ", str$(rtd * oev(12), 0, 4) + " degrees"
  print "declination             ", str$(rtd * oev(13), 0, 4) + " degrees"
  print "east longitude          ", str$(rtd * oev(15), 0, 4) + " degrees"
  
end sub
  
  '''''''''''''''''''
  '''''''''''''''''''
  
sub svprint(r(), v())
  
  ' print position and velocity vectors subroutine
  
  local rmag, vmag
  
  print "position vector x-component  ", str$(r(1), 0, -12) + " kilometers"
  print "position vector y-component  ", str$(r(2), 0, -12) + " kilometers"
  print "position vector z-component  ", str$(r(3), 0, -12) + " kilometers"
  
  rmag = vecmag(r())
  
  print " "
  
  print "position vector magnitude    ", str$(rmag, 0, -12) + " kilometers"
  
  print " "
  
  print "velocity vector x-component  ", str$(v(1), 0, -12) + " kilometers/second"
  print "velocity vector y-component  ", str$(v(2), 0, -12) + " kilometers/second"
  print "velocity vector z-component  ", str$(v(3), 0, -12) + " kilometers/second"
  
  vmag = vecmag(v())
  
  print " "
  
  print "velocity vector magnitude    ", str$(vmag, 0, -12) + " kilometers/second"
  
  
end sub
  
  '''''''''''''''''''''''''
  '''''''''''''''''''''''''
  
function modulo(x) as float
  
  ' modulo 2 pi function
  
  ''''''''''''''''''''''
  
  local a
  
  a = x - pi2 * fix(x / pi2)
  
  if (a < 0.0) then
    
    a = a + pi2
    
  end if
  
  modulo = a
  
end function
  
  '''''''''''''''''''''''''''
  '''''''''''''''''''''''''''
  
function atan3(a, b) as float
  
  ' four quadrant inverse tangent function
  
  ' input
  
  '  a = sine of angle
  '  b = cosine of angle
  
  ' output
  
  '  atan3 = angle (0 =< atan3 <= 2 * pi; radians)
  
  ''''''''''''''''''''''''''''''''''''''''''''''''
  
  local c
  
  if (abs(a) < 1.0e-10) then
    
    atan3 = (1.0 - sgn(b)) * pidiv2
    
    exit function
    
  else
    
    c = (2.0 - sgn(a)) * pidiv2
    
  endif
  
  if (abs(b) < 1.0e-10) then
    
    atan3 = c
    
    exit function
    
  else
    
    atan3 = c + sgn(a) * sgn(b) * (abs(atn(a / b)) - pidiv2)
    
  endif
  
end function
  
  '''''''''''''''''''''''''''
  '''''''''''''''''''''''''''
  
function vecmag(a()) as float
  
  ' vector magnitude function
  
  ' input
  
  '  { a } = column vector ( 3 rows by 1 column )
  
  ' output
  
  '  vecmag = scalar magnitude of vector { a }
  
  vecmag = sqr(a(1) * a(1) + a(2) * a(2) + a(3) * a(3))
  
end function
  
  ''''''''''''''''''''
  ''''''''''''''''''''
  
sub uvector (a(), b())
  
  ' unit vector subroutine
  
  ' input
  
  '  a = column vector (3 rows by 1 column)
  
  ' output
  
  '  b = unit vector (3 rows by 1 column)
  
  '''''''''''''''''''''''''''''''''''''''
  
  local amag, i%
  
  amag = vecmag(a())
  
  for i% = 1 to 3
    
    if (amag <> 0.0) then
      
      b(i%) = a(i%) / amag
      
    else
      
      b(i%) = 0.0
      
    end if
    
  next i%
  
end sub
  
  ''''''''''''''''''''''''''''''
  ''''''''''''''''''''''''''''''
  
function vdot(a(), b()) as float
  
  ' vector dot product function
  
  ' c = { a } dot { b }
  
  ' input
  
  '  n%    = number of rows
  '  { a } = column vector with n rows
  '  { b } = column vector with n rows
  
  ' output
  
  '  vdot = dot product of { a } and { b }
  
  ''''''''''''''''''''''''''''''''''''''''
  
  local c = 0.0, i%
  
  for i% = 1 to 3
    
    c = c + a(i%) * b(i%)
    
  next i%
  
  vdot = c
  
end function
  
  '''''''''''''''''''''''
  '''''''''''''''''''''''
  
sub vcross(a(), b(), c())
  
  ' vector cross product subroutine
  
  ' { c } = { a } x { b }
  
  ' input
  
  '  { a } = vector a ( 3 rows by 1 column )
  '  { b } = vector b ( 3 rows by 1 column )
  
  ' output
  
  '  { c } = { a } x { b } ( 3 rows by 1 column )
  
  c(1) = a(2) * b(3) - a(3) * b(2)
  c(2) = a(3) * b(1) - a(1) * b(3)
  c(3) = a(1) * b(2) - a(2) * b(1)
  
end sub
  
  ''''''''''''''''''''''''
  ''''''''''''''''''''''''
  
sub matxvec(a(), b(), c())
  
  ' matrix/vector multiplication subroutine
  
  ' { c } = [ a ] * { b }
  
  ' input
  
  '  a  = matrix a ( 3 rows by 3 columns )
  '  b  = vector b ( 3 rows )
  
  ' output
  
  '  c = vector c ( 3 rows )
  
  ''''''''''''''''''''''''''
  
  local s, i%, j%
  
  for i% = 1 to 3
    
    s = 0.0
    
    for j% = 1 to 3
      
      s = s + a(i%, j%) * b(j%)
      
    next j%
    
    c(i%) = s
    
  next i%
  
end sub

  ''''''''''
  ''''''''''
  
sub keycheck

    ' check user response subroutine
    
    ''''''''''''''''''''''''''''''''
    
    local a$

    print " "
    
    print "< press any key to continue >"
    
    a$ = ""
    
    while a$ = ""
    
       a$ = inkey$
       
    wend

end sub
