'  **************************     SOLEX 5.0     ****************************
'
'  Numerical Integration for Real Time Production of Fundamental Ephemerides
'        by a Fast Extrapolation Method Devised for the Solar System
'            A. Vitagliano, Celestial Mechanics 1997, 66, 293-308
'
'  ..............   Copyright (C) 1996-1999  Aldo Vitagliano    ............
'  ........ now incorporating an extended Minor Planets' data library ......
'
'  *************************************************************************
'    If you recompile and link this code with a Microsoft Basic compiler
'  (Quick Basic or PDS) you will loose a lot of speed. You can restore the
'  full speed to the modified EXE file by running the utility SPEEDUP.EXE
'  included in the package. This code is portable to a PowerBasic 3.x (DOS)
'  compiler by removing the DECLARE statements and activating the initial
'  line OPTION BINARY BASE 1.  Mail the author for further information.
'                        E-mail  alvitagl@unina.it
'  *************************************************************************
'
'  ----------------------  Last Revised 1999/06/13  ------------------------
'
'  *************************************************************************
'
DECLARE SUB AddBody (NEX%, N%, DT#, X#(), Y#(), Z#(), VX#(), VY#(), VZ#(), PLN() AS STRING * 8, Rj0#(), Rj00#(), M0#(), ABSM#())
DECLARE SUB CAcc (N%, X#(), Y#(), Z#(), AX#(), AY#(), AZ#(), M#())
DECLARE FUNCTION DeltaT# (TJ#)
DECLARE SUB DispToggles (N%, TimUnName() AS STRING * 4, PLN() AS STRING * 8, FLGEO%, FLECL%, FLPR%, FLNUTHIGH%, FLABP%, FLABS%, FLWRITE%, FL1ST%, FLBACK%, FLINV%, TIMUN%, SCAL%)
DECLARE FUNCTION Dj2000# (YR#, MN#, DY#)
DECLARE SUB DrawFrame (CL%, BKCL%)
DECLARE SUB DrawTitle (CL%, RO%, CO%)
DECLARE SUB FindAster (FLFOUND%)
DECLARE SUB FixBaric (N%, M0#(), X#(), Y#(), Z#(), VX#(), VY#(), VZ#(), XC0#, YC0#, ZC0#)
DECLARE SUB FixV (N%, DT#, DT0#, VX#(), VY#(), VZ#())
DECLARE SUB GetStart (TJ#, Nday#, NUM%, NST%, DjLib() AS LONG, M0#(), Rj0#(), ABSM#(), PLN() AS STRING * 8, X#(), Y#(), Z#(), VX#(), VY#(), VZ#())
DECLARE SUB Header ()
DECLARE SUB InputFile (NF$)
DECLARE SUB Integrate (Njump&, N%, T#, X#(), Y#(), Z#(), VX#(), VY#(), VZ#())
DECLARE SUB Nutat (TJ#, DPSI#, DEPS#)
DECLARE SUB OutHeadings (FLDE%, FLORT%, FLGEO%, FLTOP%, FLHOR%, FLECL%, FLPR%, FLNUTHIGH%, FLABP%, FLABS%, GLAT#, GLON#)
DECLARE SUB PlanMag (N%, XO#(), YO#(), ZO#(), ABSM#(), PH#(), MAGP#())
DECLARE SUB PrecAng (DJP#, FLNUTLOW%, FLNUTHIGH%, FLDE%)
DECLARE SUB PrintAngOut (FLGEO%, FL1ST%, N%, T#, ABSM#(), PLN() AS STRING * 8, XO#(), YO#(), ZO#(), R#(), YY%, MM%, DD%, HH%, MN%, SS%)
DECLARE SUB PrintOut (Djump#, TimUnName() AS STRING * 4, ABSM#(), PLN() AS STRING * 8, FL1ST%, FLGEO%, N%, T#, TC#, XO#(), YO#(), ZO#(), R#(), SCAL%)
DECLARE SUB RadToDMS (X#, SG$, DG%, HMIN%, SEC#, NPREC%, FL360%)
DECLARE SUB RadToHMS (X#, HR%, HMIN%, SEC#, NPREC%)
DECLARE SUB ReadFile (NF$, TJ#, NEX%, X#(), Y#(), Z#(), VX#(), VY#(), VZ#(), Rj0#(), Rj00#(), PLN() AS STRING * 8, M0#(), ABSM#())
DECLARE FUNCTION RelPar# (X#, Y#, Z#, VX#, VY#, VZ#)
DECLARE SUB RotAx (N%, X#(), Y#(), Z#(), FLECL%)
DECLARE SUB RotInv (N%, X#(), Y#(), Z#())
DECLARE SUB ScanLib (FLDE%, NST%, YMIN#, YMAX#, DjLib() AS LONG)
DECLARE SUB SelOption (Outp%)
DECLARE SUB SetDeltaT (DELT0#)
DECLARE SUB SetGeoCoord (GLAT#, GLON#, ZH#)
DECLARE SUB SetOrder (DT#, NO%)
DECLARE SUB SitePos (FLTOP%, FLGEO%, T#, DT#, GLAT#, GLON#, ZH#, XSIT#, YSIT#, ZSIT#, VXSIT#, VYSIT#, VZSIT#)
DECLARE SUB STORM8 (N%, X#(), Y#(), Z#(), VX#(), VY#(), VZ#())
DECLARE FUNCTION Tsid# (TU#, GLON#)
DECLARE SUB TToDate (TP#, Y%, MO%, D%, HR%, MINUT%, SEC%)

'---------------------------------------------------------------------------
' OPTION BINARY BASE 1 '... Must be activated if compiled with PowerBasic
'---------------------------------------------------------------------------

DEFINT F, H-L, N
DEFDBL A-E, G, M, O-Z

DIM M0(20), M(20), ABSM(20), Rj0(20), X(20), Y(20), Z(20), VX(20), VY(20), VZ(20)
DIM R(20), XO(20), YO(20), ZO(20), TIT$(5), AK(7, 8), DRj0(20), Rj00(20)
DIM DR1(20), DR2(20), DR3(20), DR4(20)
DIM DjLib(4095) AS LONG, PLN(20) AS STRING * 8
DIM TimUnFactor(4), TimUnName(4) AS STRING * 4

DEF FNATN2 (YW, XW)
IF XW = 0 THEN ZZZ = PIMEZ * SGN(YW) ELSE ZZZ = ATN(YW / XW)
   IF XW < 0 THEN ZZZ = ZZZ + PI * SGN(YW)
   FNATN2 = ZZZ
END DEF

' ............... Coefficients of the Extrapolation Series .............
DATA 1, 1, 0, 0, 0, 0, 0, 0, 0
DATA 3, -1, 4, 0, 0, 0, 0, 0, 0
DATA 240, 10, -256, 486, 0, 0, 0, 0, 0
DATA 2520, -7, 896, -6561, 8192, 0, 0, 0, 0
DATA 362880, 42, -24576, 531441, -2097152, 1953125, 0, 0, 0
DATA 19958400, -66, 168960, -9743085, 92274688, -244140625, 181398528, 0, 0

' the following combined coefficients are specific for the Moon orbit
' with n=7 and H=3d and are recommended for long integrations.
 DATA 922442320800, 33462, -480397632, 71157646989, -1394145099776, 7143798828125, -10686833516736, 5788944826368, 0

' the following are specific for n=7 and H=3.75d. They ensure the best
' compromise between speed and accuracy in short integrations (up to 50 yrs)
'DATA 758788430400,33462,-411126144,59537158083,-1157698420736,5911865234375,-8827555291776,4773050843136,0
'
DATA 49037788800, -26,1153152,-387420489,14394851328,-128173828125,322333846848,-549755813888,390625000000
FOR I = 0 TO 7: FOR J = 0 TO 8: READ AK(I, J): NEXT: NEXT
' ............ adopted sequence of stepsizes
DATA 1, 2, 3, 4, 5, 6, 8, 10
FOR K = 0 TO 7: READ QK(K): NEXT
' .... Time conversion factors to day, hr, min, sec
DATA 1, "day ", 24, "hr  ", 1440, "min ", 86400, "sec "
FOR K = 0 TO 3: READ TimUnFactor(K), TimUnName(K): NEXT
PI = 4 * ATN(1): RD = 180 / PI
GRAV0 = .0004980621312# '... Gravitational constant (time unit is 1 d)
VL = 25902.06837#  ' ....... Speed of light (length unit is 1 Gm)
PHAS = .040700012# '....................... adopted tidal phase angle
CV = .0864#        '............. converts velocity from Km/s to Gm/d
FLDE = 0
OPEN "TITLE.S50" FOR INPUT AS #1
FOR I = 1 TO 4: LINE INPUT #1, TIT$(I): NEXT
CLOSE #1
CALL Header

DO
5    '..................... Start
  CALL SelOption(Outp%): FLFILE = 0
6 SELECT CASE Outp%
     CASE 6: CLS : EXIT DO
     CASE 5: COLOR 7, 0, 0: CLS : SHELL: COLOR 14, 1, 1: GOTO 5
     CASE 3
        FLFILE = -1
        Outp% = 2: GOTO 6
     CASE 4
        FLASTER = -1
        Outp% = 3: GOTO 6
     CASE 1 TO 2
        IF Outp% <> FLDE THEN
           FLDE = Outp%
           CLS : LOCATE 10, 30: PRINT "Scanning Library ... ";
           CALL ScanLib(FLDE, NST, YMIN, YMAX, DjLib())


7          DATA -5.3412578D-008,6.48446607D-008,-3.30323978D-007,-6.3974422D-009,4.755770D-015,2.64932795D-005
           READ QJ0, QJ1, QJ2, QJ3, QJ4, QJ5  '... Quadrupole and related parameters
           DATA -5.39096D-008,6.48514D-008,-3.303186D-007,-6.42605D-009,5.1352373D-015,2.6448D-005

           '... empirical coefficients for secular longitude adjustements
           '... of Planets Mercury ... Jupiter in DE406

           DATA +5.48746e-013, -1.65320e-014, +5.02699e-016, -5.04797e-018
           DATA +7.18056e-013, +2.56972e-014, +1.88422e-015, -4.97720e-017
           DATA +1.95101e-012, -8.71138e-014, +1.76489e-015, -6.59207e-018
           DATA -1.22642e-010, +1.20640e-011, -3.94359e-013, +4.03462e-015
           DATA -9.14339e-012, -5.51662e-013, +1.58014e-014, -9.56279e-017

           '... parameters are different if option is DE406

           ERASE DR1, DR2, DR3, DR4
           IF FLDE = 2 THEN
              READ QJ0, QJ1, QJ2, QJ3, QJ4, QJ5
              FOR I = 1 TO 5: READ DR1(I), DR2(I), DR3(I), DR4(I): NEXT
           END IF
           RESTORE 7
        END IF
  END SELECT

  QLOV = QJ4 / PHAS: Qtid = QJ4
  IF FLDE = 2 THEN '... Oblate Sun and Williams Obliquity if DE406
     AU = 149.597870691# '... Astronomic Unit (Gm)
     Qs2 = .0000001453#: OBL0 = 84381.40617#
  ELSE
     AU = 149.59787066# '... Astronomic Unit (Gm)
     Qs2 = 0: OBL0 = 84381.448#
  END IF
  OBL0 = OBL0 * PI / 648000#
  SEQ0 = SIN(OBL0): CEQ0 = COS(OBL0)

  NTER = 3: NLUN = 13  '... Earth is body n. 3, Moon is body n. 13
  FLGEO = -1: FLPR = -1: FLECL = 0: FLABS = -1: FLABP = -1: FLORT = 0: FLDEC = 0
  FLNUTHIGH = -1: FLWRITE = 0: FLHOR = 0

  OutRange% = 1
  IF FLFILE THEN OutRange% = 0
  WHILE OutRange%
     CLS
     CALL DrawFrame(14, 1)
     CALL DrawTitle(11, 2, 14)
     COLOR 14: LOCATE 10, 23
     INPUT ; "  Enter Date (YYYY,MM,DD)   ", YY, MM, DD
     IF YY >= YMIN AND YY <= YMAX THEN OutRange% = 0
     IF OutRange% THEN
        SOUND 900, 4: CLS
        PRINT "  Allowed Range is  "; YMIN; " to "; YMAX
        PRINT "  Press a key to continue  ";
        M$ = INPUT$(1)
     END IF
  WEND
  IF FLFILE THEN
     IF FLASTER THEN
        CALL FindAster(FLFOUND)
        IF FLFOUND THEN FLFOUND = 0: NF$ = "DAT.TMP" ELSE GOTO 5
        FLASTER = 0
     ELSE
        CALL InputFile(NF$)
        IF LEN(NF$) < 3 THEN 5
     END IF
     CALL ReadFile(NF$, TJ, NEX, X(), Y(), Z(), VX(), VY(), VZ(), Rj0(), Rj00(), PLN(), M0(), ABSM())
     CALL TToDate(TJ + .5, A%, M%, D%, H%, MINUT%, SEC%)
     YY = A%: MM = M%: DD = D% + H% / 24 + MINUT% / 1440 + SEC% / 86400
  ELSE
     TJ = Dj2000(YY, MM, DD) - .5
  END IF
  LOCATE 11, 23: INPUT ; "  Enter Tab. Interval (d)   ", TabStep0
  LOCATE 12, 23: INPUT ; "  Enter # of Output Records ", Nout
  IF Nout <= 0 THEN Nout = 1
  IF TabStep0 = 0 THEN TabStep0 = 1
  D% = SGN(TabStep0): TabStep0 = ABS(TabStep0): FLBACK = (D% < 0)
  CALL GetStart(TJ, Nday#, N, NST, DjLib(), M0(), Rj0(), ABSM(), PLN(), X(), Y(), Z(), VX(), VY(), VZ())
  FOR I = 0 TO N: Rj00(I) = Rj0(I): NEXT
  RJ1 = 6 * M0(0) * GRAV0 / (VL * VL)
  MS% = QK(0): DT0 = 1
  T = TJ + .5
  IF FLDE = 2 THEN
     IF TJ > 24515098 THEN DRj0(5) = .000000000008# ELSE DRj0(5) = -.000000000012#
  END IF

  ' ..........................................................................
8 S% = 1
  IF Nday# < 0 THEN S% = -1  ' ......  S% is negative for backwards integration
 ' DFSTEP = 3.75              ' ......  Default Stepsize
  DFSTEP = 3        ' Optimal stepsize to be used in long integrations after
                     ' activating the corresponding DATA statement.
  DT = DFSTEP: NO = 6         ' ......  Initial Stepsize and Number of stages
  Nday# = ABS(Nday#): Djump = Nday#
  Njump& = INT(Nday# / DFSTEP) ' ......  Number of steps to reach target epoch
  SStep = Nday# - DFSTEP * Njump& '...  Size of last step to reach target epoch
  LOCATE 15, 25: PRINT "Jumping to Epoch ";
  PRINT USING "#####/##/##"; YY; MM; DD;
  COLOR 23: PRINT " ...";
  COLOR 14: TC = TIMER

  FOR I = 0 TO N
     VX(I) = S% * VX(I) * CV ' ......  Velocity is converted to Gm/d
     VY(I) = S% * VY(I) * CV
     VZ(I) = S% * VZ(I) * CV
  NEXT
  '.......................................... Baricenter is set to 0,0,0
  CALL FixBaric(N, M0(), X(), Y(), Z(), VX(), VY(), VZ(), XC0, YC0, ZC0)
  CALL FixV(N, DT, DT0, VX(), VY(), VZ()) '.. Time unit is set equal to stepsize
  DJP = T + 30 * S% - S% / 4096  '.... Time at which reference frame is adjusted
  CALL PrecAng(DJP, 1, 0, FLDE)           '.... Precession angles are computed
  '                               .... and frame is rotated accordingly
  CALL RotAx(N, X(), Y(), Z(), 0): CALL RotAx(N, VX(), VY(), VZ(), 0)
  '   .................... Integration up to one step before target epoch
  CALL Integrate(Njump&, N, T, X(), Y(), Z(), VX(), VY(), VZ())

  IF SStep THEN  ' ............................ Last step to target epoch
     DT0 = DT: DT = SStep: NO = 7
     CALL FixV(N, DT, DT0, VX(), VY(), VZ())
     CALL STORM8(N, X(), Y(), Z(), VX(), VY(), VZ()) ' ... Integration Routine
     T = T + MS% * S% * DT
  END IF
  IF FLFILE THEN
     CALL RotInv(N, X(), Y(), Z()): CALL RotInv(N, VX(), VY(), VZ())
     CALL AddBody(NEX, N, DT, X(), Y(), Z(), VX(), VY(), VZ(), PLN(), Rj0(), Rj00(), M0(), ABSM())
     CALL RotAx(N, X(), Y(), Z(), 0): CALL RotAx(N, VX(), VY(), VZ(), 0)
  END IF
  FLFILE = 0
  TIMUN% = 0: TabStep = TabStep0: SCAL% = 0
  TabSpan = Nout * TabStep0: TT0 = T

  P% = S% * D%
  IF P% < 0 THEN ' ...... Runs backwards if Tabulation interval is negative
     S% = D%
     FOR I = 0 TO N: VX(I) = -VX(I): VY(I) = -VY(I): VZ(I) = -VZ(I): NEXT
  END IF
  CLS

  DO  ' ............. Ephemerides are computed and tabulated
     FLSTEP = -1: NDUMMY = N
     DO

        IF FLGEO THEN NC0 = 3 ELSE NC0 = 0
        '..... Further computations require absolute coordinates for the Moon
        XL = X(NLUN): YL = Y(NLUN): ZL = Z(NLUN)
        X(NLUN) = XL + X(3): Y(NLUN) = YL + Y(3): Z(NLUN) = ZL + Z(3)
        NC = NC0
        DELTA = DeltaT(T)
        CALL SitePos(FLTOP, FLGEO, T - DELTA, DT, GLAT, GLON, ZH, XSIT, YSIT, ZSIT, VXSIT, VYSIT, VZSIT)
        FOR I = 0 TO N
           IF I = N THEN NC = 3
           DX = X(I) - X(NC) - XSIT: DY = Y(I) - Y(NC) - YSIT: DZ = Z(I) - Z(NC) - ZSIT
           R(I) = SQR(DX * DX + DY * DY + DZ * DZ)
        NEXT I

        IF (FLABS OR FLABP) THEN   '....... Computes Aberration

           DLIGHT = VL * DT   '....... Light Distance in 1 step
           NC = NC0
           FOR I = 0 TO N
              IF I = N THEN FLMOO = 1: NC = 3 ELSE FLMOO = 0'.. Vmoon is still relative
              DL = R(I) / DLIGHT '.... Light-Time from Planet to Central body (in steps)

              XO(I) = X(I) + (VX(I) + (VX(NC) + VXSIT) * (FLABP + FLMOO)) * DL * S% * FLABS
              YO(I) = Y(I) + (VY(I) + (VY(NC) + VYSIT) * (FLABP + FLMOO)) * DL * S% * FLABS
              ZO(I) = Z(I) + (VZ(I) + (VZ(NC) + VZSIT) * (FLABP + FLMOO)) * DL * S% * FLABS
           NEXT I
        ELSE
           FOR I = 0 TO N: XO(I) = X(I): YO(I) = Y(I): ZO(I) = Z(I): NEXT
        END IF
        X(NLUN) = XL: Y(NLUN) = YL: Z(NLUN) = ZL '.. Restore Geocentric Coordinates
        IF FLGEO THEN XO(3) = XO(3) + XSIT: YO(3) = YO(3) + YSIT: ZO(3) = ZO(3) + ZSIT
        '... Frame rotations performed according to the option selected

        CALL RotInv(N, XO(), YO(), ZO())   '.... Rotation back to J2000

        IF FLPR THEN TP = T ELSE TP = 0 '.... T=2000.0 if Precession = Off
                                      '.... Rotation to mean equinox of date
        CALL PrecAng(TP, 0, FLNUTHIGH, FLDE)  '.... or True equinox if FLNUTHIGH=1
        CALL RotAx(N, XO(), YO(), ZO(), FLECL)
        CALL PrecAng(DJP, 1, 0, FLDE)

        CALL PrintOut(Djump, TimUnName(), ABSM(), PLN(), FL1ST%, FLGEO, N, T, TC, XO(), YO(), ZO(), R(), SCAL%)
        IF FLWRITE = 0 THEN
           CALL OutHeadings(FLDE, FLORT, FLGEO, FLTOP, FLHOR, FLECL, FLPR, FLNUTHIGH, FLABP, FLABS, GLAT, GLON)
           CALL DispToggles(N, TimUnName(), PLN(), FLGEO, FLECL, FLPR, FLNUTHIGH, FLABP, FLABS, FLWRITE, FL1ST, FLBACK, FLINV, TIMUN%, SCAL%)
        END IF
        IF FLEXIT THEN 5 ' ............ Back to start
        TC = TIMER: Djump = 0
        IF FLSTEP THEN EXIT DO
     LOOP
     IF FLGO THEN
        FLGO = 0: CALL DrawFrame(14, 1)
        LOCATE 10, 22
        INPUT ; "  Enter Date (YYYY,MM,DD)   ", YY, MM, DD
        TJF = Dj2000(YY, MM, DD)
        Nday# = TJF - T
        DT0 = DT: DT = DFSTEP: NO = 6: SCV = S% * CV * DT
        CALL FixV(N, DT, DT0, VX(), VY(), VZ()): DT0 = 1
        FOR I = 0 TO N
           VX(I) = VX(I) / SCV: VY(I) = VY(I) / SCV: VZ(I) = VZ(I) / SCV
        NEXT I
        CALL RotInv(N, X(), Y(), Z()): CALL RotInv(N, VX(), VY(), VZ())
        GOTO 8
     END IF
     IF FLINV THEN  ' ... Time direction has been inverted
        S% = -S%
        FOR I = 0 TO N: VX(I) = -VX(I): VY(I) = -VY(I): VZ(I) = -VZ(I): NEXT
     END IF
     IF FLBACK THEN D% = -1 ELSE D% = 1 ' ... Check time direction

     '... Integration param. are set according to Tabulation Step
     TabStep = TabStep0 * 10 ^ SCAL% / TimUnFactor(TIMUN%)
     D1 = TabStep / 3: ND& = INT(D1)
     IF ND& < D1 THEN ND& = ND& + 1
     DT0 = DT
     DT = TabStep / ND&  ' ....... Optimal stepsize for given tabulation interval
     CALL SetOrder(DT, NO) '..... Optimal order (no. of stages) for given stepsize
     CALL FixV(N, DT, DT0, VX(), VY(), VZ()) '... Time unit scaled to stepsize

     CALL Integrate(ND&, N, T, X(), Y(), Z(), VX(), VY(), VZ())
     Djump = TabStep

  LOOP WHILE ABS(T - TT0) < TabSpan
  CLOSE

LOOP
COLOR 7, 0, 0: CLS

END

SUB AddBody (NEX, N, DT, X(), Y(), Z(), VX(), VY(), VZ(), PLN() AS STRING * 8, Rj0(), Rj00(), M0(), ABSM())
   SHARED NLUN, CV, S%
   FOR I = N + 1 TO N + NEX
      X(I) = X(0) + X(I): Y(I) = Y(0) + Y(I): Z(I) = Z(0) + Z(I)
      VX(I) = VX(0) + S% * VX(I) * CV * DT
      VY(I) = VY(0) + S% * VY(I) * CV * DT
      VZ(I) = VZ(0) + S% * VZ(I) * CV * DT
   NEXT
   N = N + NEX
   SWAP X(N), X(NLUN): SWAP Y(N), Y(NLUN): SWAP Z(N), Z(NLUN)
   SWAP VX(N), VX(NLUN): SWAP VY(N), VY(NLUN): SWAP VZ(N), VZ(NLUN)
   SWAP PLN(N), PLN(NLUN): SWAP M0(N), M0(NLUN): SWAP Rj0(N), Rj0(NLUN)
   SWAP ABSM(N), ABSM(NLUN)
   NLUN = N
END SUB

'
'..... Computes all the Accelerations in the system of bodies 0 to N
'
SUB CAcc (N, X(), Y(), Z(), AX(), AY(), AZ(), M()) STATIC
     SHARED Qs2, QJ0, QJ1, QJ2, QJ3, Qtid, QJ5, QLOV, Rj0(), RJ1, NTER, NLUN, SEQ, CEQ
     FOR I = 0 TO N: AX(I) = 0: AY(I) = 0: AZ(I) = 0: NEXT
' ..... Lunar Geocentric coordinates are stored and converted to baricentric
     XL = X(NLUN): YL = Y(NLUN): ZL = Z(NLUN): X(NLUN) = XL + X(NTER): Y(NLUN) = YL + Y(NTER): Z(NLUN) = ZL + Z(NTER)
     AX(0) = 0: AY(0) = 0: AZ(0) = 0
     FOR J = 1 TO N  ' ...... Computes accelerations by the Sun
        DX = X(J) - X(0): DY = Y(J) - Y(0): DZ = Z(J) - Z(0)
        R2INV = 1 / (DX * DX + DY * DY + DZ * DZ): R3INV = R2INV * SQR(R2INV)
        AFX = (1 + Rj0(J) + RJ1 * R3INV / R2INV + Qs2 * R2INV) * R3INV

        '... RJ0 and RJ1 are Relativistic perturbations
        '... Qs2 is the radial part of the Sun Oblateness J2 Perturbation

        IF J = NTER THEN '... If Planet is Earth then Oblateness is considered
           RZR = DZ * DZ * R2INV
           AFZ = AFX + QJ2 * (RZR - .6#) * R3INV * R2INV: AFX = AFX + QJ2 * (RZR - .2#) * R3INV * R2INV
           ARJ = AFX * M(J): ARJZ = AFZ * M(J): ARI = -AFX * M(0): ARIZ = -AFZ * M(0)
        ELSE
           IF J = NLUN THEN RL2 = R2INV '...... Sun-Moon r^2 is stored
           ARJ = AFX * M(J): ARJZ = ARJ: ARI = -AFX * M(0): ARIZ = ARI
        END IF
        AX(J) = DX * ARI: AY(J) = DY * ARI: AZ(J) = DZ * ARIZ
        AX(0) = AX(0) + DX * ARJ: AY(0) = AY(0) + DY * ARJ: AZ(0) = AZ(0) + DZ * ARJZ
     NEXT J

     FOR I = 1 TO N - 1
        FOR J = I + 1 TO N
           IF I = NTER AND J = NLUN THEN '........ Earth-Moon acceleration
              R2INV = 1 / (XL * XL + YL * YL + ZL * ZL): R3INV = R2INV * SQR(R2INV): ARJ = M(J) * R3INV: ARI = -M(I) * R3INV
              RZR = ZL * ZL * R2INV
              DZ1 = ZL * CEQ - YL * SEQ: RZR1 = DZ1 * DZ1 * R2INV
              AFX = 1 + QJ0 + (QJ1 + QJ2 * RZR + QJ3 * RZR1 + QJ5 * RL2) * R2INV
              AFY = AFX * YL + (.4# * QJ3 * DZ1 * SEQ + (QLOV * YL - Qtid * XL) * R3INV) * R2INV
              AFZ = AFX * ZL + (-.4# * (QJ2 * ZL + QJ3 * DZ1 * CEQ) + QLOV * ZL * R3INV) * R2INV
              AFX = AFX * XL + (QLOV * XL + Qtid * YL) * R3INV * R2INV
              AX(J) = AX(J) + ARI * AFX: AX(I) = AX(I) + ARJ * AFX
              AY(J) = AY(J) + ARI * AFY: AY(I) = AY(I) + ARJ * AFY
              AZ(J) = AZ(J) + ARI * AFZ: AZ(I) = AZ(I) + ARJ * AFZ
           ELSE     '......................... Planet-Planet accelerations
              DX = X(J) - X(I): DY = Y(J) - Y(I): DZ = Z(J) - Z(I)
              R2INV = 1 / (DX * DX + DY * DY + DZ * DZ): R3INV = R2INV * SQR(R2INV): ARJ = M(J) * R3INV: ARI = -M(I) * R3INV
              AX(J) = AX(J) + DX * ARI: AY(J) = AY(J) + DY * ARI: AZ(J) = AZ(J) + DZ * ARI
              AX(I) = AX(I) + DX * ARJ: AY(I) = AY(I) + DY * ARJ: AZ(I) = AZ(I) + DZ * ARJ
           END IF
        NEXT J
     NEXT I
     '.............................  Restores Lunar Geocentric Coordinates
     '......................... and computes Lunar Geocentric Acceleration
     X(NLUN) = XL: Y(NLUN) = YL: Z(NLUN) = ZL
     AX(NLUN) = AX(NLUN) - AX(NTER): AY(NLUN) = AY(NLUN) - AY(NTER): AZ(NLUN) = AZ(NLUN) - AZ(NTER)
END SUB

FUNCTION DeltaT (TJ)
   SHARED FLDT, DELT0
   '.... TJ is the number of Julian Days from J2000

   IF FLDT = 0 THEN DeltaT = 0: EXIT FUNCTION
   IF FLDT = 2 THEN DeltaT = DELT0: EXIT FUNCTION

   T = TJ / 365.25: AJ& = CLNG(2000 + T)
   IF AJ& < 948 THEN
      DELTA = (2177! + T * (4.97 + T * .00441))
   ELSEIF AJ& < 1600 THEN
      DELTA = 50.6 + T * (.675 + T * .00225)
   ELSEIF AJ& < 1700 THEN
      DELTA = 1791 + T * (11.672 + T * .0191234)
   ELSEIF AJ& < 1900 THEN
      DELTA = 6
   ELSE
      DELTA = 67 + .67# * T
   END IF
   DeltaT = DELTA / 86400
END FUNCTION

SUB DispToggles (N, TimUnName() AS STRING * 4, PLN() AS STRING * 8, FLGEO, FLECL, FLPR, FLNUTHIGH, FLABP, FLABS, FLWRITE, FL1ST, FLBACK, FLINV, TIMUN%, SCAL%)
   SHARED FLDEC, FLJD, FLDT, FLTOP, DELT0, GLAT, GLON, ZH, FLGO, FLHOR, FLEXIT, FLSTEP, FLORT
   DIM STATUS(-1 TO 2) AS STRING * 4
   DIM CL%(-1 TO 2)
   DIM StepScal(-2 TO 2) AS STRING * 4
   STATUS(0) = "Off ": STATUS(-1) = "On  ": STATUS(1) = "Aut ": STATUS(2) = "Man "
   StepScal(-2) = "/100": StepScal(-1) = "/10 ": StepScal(0) = "*1  ": StepScal(1) = "*10 ": StepScal(2) = "*100"
   CL%(0) = 12: CL%(-1) = 10: CL%(1) = 10: CL%(2) = 15
   V$ = "": FLEXIT = 0: FLINV = 0: L = 5: C% = 58
   COLOR 15, 0
   LOCATE L, C%
   PRINT "Key    Action    Status"; : L = L + 1
   COLOR CL%(FLORT): LOCATE L, C%: PRINT " O   Orthogonal    "; STATUS(FLORT);
   COLOR CL%(FLGEO): LOCATE L + 1, C%: PRINT " G   Geocentric    "; STATUS(FLGEO);
   COLOR CL%(FLTOP): LOCATE L + 2, C%: PRINT " T   Topocentric   "; STATUS(FLTOP);
   COLOR CL%(FLECL): LOCATE L + 3, C%: PRINT " E   Ecliptic      "; STATUS(FLECL);
   COLOR CL%(FLHOR): LOCATE L + 4, C%: PRINT " H   Horizontal    "; STATUS(FLHOR);
   COLOR CL%(FLPR): LOCATE L + 5, C%: PRINT " P   Precession    "; STATUS(FLPR);
   COLOR CL%(FLNUTHIGH): LOCATE L + 6, C%: PRINT " N   Nutation      "; STATUS(FLNUTHIGH);
   COLOR CL%(FLABP): LOCATE L + 7, C%: PRINT " A   Aberration    "; STATUS(FLABP);
   COLOR CL%(FLABS): LOCATE L + 8, C%: PRINT " L   Light-time    "; STATUS(FLABS);
   COLOR CL%(FLDT): LOCATE L + 9, C%: PRINT " D   Delta-T       "; STATUS(FLDT);
   COLOR CL%(FLBACK): LOCATE L + 10, C%: PRINT " B   Back-step     "; STATUS(FLBACK);
   COLOR CL%(FLDEC): LOCATE L + 11, C%: PRINT " M   Metric Units  "; STATUS(FLDEC);
   COLOR 3: LOCATE L + 12, C%: PRINT " F   File Output   "; TAB(80); " ";
   COLOR 3: LOCATE L + 13, C%: PRINT " U   Time Unit     "; TimUnName(TIMUN%);
   COLOR 3: LOCATE L + 14, C%: PRINT " S   Step Scaling  "; StepScal(SCAL%);
   COLOR 3: LOCATE L + 15, C%: PRINT " C   Geogr. Coord."; TAB(80); " ";
   COLOR 3: LOCATE L + 16, C%: PRINT " J   Jump to Date"; TAB(80); " ";
   COLOR 7: LOCATE L + 17, C%: PRINT "Esc  Main Menu"; TAB(80); " ";
   LOCATE 25, 11: COLOR 12, 0
   PRINT " Press a Key to select Option or Spacebar for Next Step ... ";
   COLOR 14, 1
   WHILE V$ = ""
      V$ = UCASE$(INKEY$)
      FLSTEP = 0
      SELECT CASE V$
         CASE "O": FLORT = NOT FLORT: IF FLORT THEN FLHOR = 0
         CASE "G": FLGEO = NOT FLGEO: IF FLGEO = 0 THEN FLTOP = 0: FLHOR = 0
         CASE "T": FLTOP = NOT FLTOP
            IF FLTOP THEN
               FLGEO = -1: FLPR = -1
               IF FLDT = 0 THEN FLDT = 1
            ELSE
               FLHOR = 0
            END IF
         CASE "E": FLECL = NOT FLECL: IF FLECL THEN FLHOR = 0
         CASE "H": FLHOR = NOT FLHOR
            IF FLHOR THEN
               FLTOP = -1: FLGEO = -1: FLECL = 0: FLPR = -1
               IF FLDT = 0 THEN FLDT = 1
            END IF
         CASE "P": FLPR = NOT FLPR: IF FLPR = 0 THEN FLNUTHIGH = 0
         CASE "N": FLNUTHIGH = NOT FLNUTHIGH: IF FLNUTHIGH THEN FLPR = -1
         CASE "A": FLABP = NOT FLABP: IF FLABP THEN FLABS = -1
         CASE "L": FLABS = NOT FLABS: IF FLABS = 0 THEN FLABP = 0
         CASE "F": FLWRITE = NOT FLWRITE
         CASE "B": FLBACK = NOT FLBACK: FLSTEP = -1: FLINV = -1
         CASE "M": FLDEC = NOT FLDEC
         CASE "U": TIMUN% = TIMUN% + 1: IF TIMUN% > 3 THEN TIMUN% = 0
         CASE "S": SCAL% = SCAL% + 1: IF SCAL% > 2 THEN SCAL% = -2
         CASE "D": FLDT = FLDT + 1: IF FLDT > 2 THEN FLDT = 0
                   IF FLDT = 2 THEN CALL SetDeltaT(DELT0)
         CASE "C": CALL SetGeoCoord(GLAT, GLON, ZH)
         CASE "J": FLGO = -1: FLSTEP = -1
         CASE CHR$(27): FLEXIT = -1
         CASE " ": FLSTEP = -1
         CASE CHR$(13): FLSTEP = -1
      END SELECT
   WEND
   COLOR 14
   IF FLGEO = 0 THEN FLABP = 0: FLABS = 0
   IF FLWRITE THEN
      CALL DrawFrame(14, 1)
      LOCATE 3, 22
      PRINT "Screen Output will be suspended and";
      LOCATE 4, 18
      PRINT "Records will be written to the following files: ";
      IF N > 13 THEN I0 = N - 13 ELSE I0 = 0
      LOCATE 6, 34: PRINT "OUTPUT  .LOG";
      FOR I = I0 TO N: LOCATE 7 + I, 34: PRINT PLN(I); ".OUT"; : NEXT
      LOCATE 23, 32: INPUT ; "O. K.  (Y or N) "; OK$
      CLS
      IF UCASE$(OK$) <> "Y" THEN FLWRITE = 0 ELSE FL1ST = -1
      IF FLWRITE THEN
         CLS : LOCATE 10, 22: PRINT "Do you want Time to be recorded as:";
         LOCATE 12, 25: PRINT "Date   (D) ";
         LOCATE 13, 25: PRINT "or JD  (J) ?"; : VD$ = ""
         WHILE VD$ <> "D" AND VD$ <> "J": VD$ = UCASE$(INKEY$): WEND
         IF VD$ = "J" THEN FLJD = -1 ELSE FLJD = 0
      END IF
   END IF

END SUB

'.... Computes Julian Day from Date (YR, MN, DY)
'.... Dj2000 is no. of Julian Days from 2000, January 1, 12h TDT (2000.0)
'
FUNCTION Dj2000 (YR, MN, DY)

    IF MN = 0 THEN MN = 1
    IF DY = 0 THEN DY = 1

    MN1 = MN: YR1 = YR: B = 0
    IF MN < 3 THEN MN1 = MN + 12: YR1 = YR1 - 1
    FLJUL = (YR < 1582) OR (YR = 1582 AND MN < 10) OR (YR = 1582 AND MN = 10 AND DY < 15)
    IF FLJUL = 0 THEN A = INT(YR1 / 100): B = 2 - A + INT(A / 4)
    Dj1 = (INT(365.25 * (YR1 + 4716)) + INT(30.6001# * (MN1 + 1)) + B - 1524.5 - 2451545)
    Dj2000 = Dj1 + DY

END FUNCTION

SUB DrawFrame (CL%, BKCL%)
    CLS : COLOR CL%, BKCL%
    PRINT CHR$(201); STRING$(78, 205); CHR$(187);
    FOR I = 2 TO 24
      LOCATE I, 1: PRINT CHR$(186);
      LOCATE I, 80: PRINT CHR$(186);
    NEXT
    LOCATE 25, 1
    PRINT CHR$(200); STRING$(78, 205); CHR$(188);

END SUB

SUB DrawTitle (CL%, RO%, CO%)
    SHARED TIT$()
    LOCATE RO%, CO%: COLOR CL%
    FOR I = 1 TO 4
       PRINT TIT$(I); : LOCATE RO% + I, CO%
    NEXT
END SUB

SUB FindAster (FLFOUND)
    DIM NOME AS STRING * 12, UCAS AS STRING * 12
    N = 13: DJ = 2451100.5#
    OPEN "ASTERB.BIN" FOR BINARY AS #3: OPEN "DAT.TMP" FOR OUTPUT AS #4

    PRINT #4, DJ
    A$ = "00": NREC& = LOF(3) / 64 - 1: FLFOUND = -1: NADD = 0
    CLS
    DO
       LOCATE 9, 1: PRINT SPACE$(79): PRINT SPACE$(79): LOCATE 9, 13
       IF FLFOUND = 0 THEN PRINT "Asteroid "; A$; " not in Library"
       INPUT ; "Enter Asteroid Number or Name (0 when finished): ", A$
       IF A$ = "" OR A$ = "0" THEN EXIT DO
       NA = VAL(A$): UC$ = LTRIM$(RTRIM$(UCASE$(A$)))
       FLFOUND = 0: IF NA THEN NC = 1 ELSE NC = 2
       FOR I = 1 TO NC
          FOR K& = 0 TO NREC&
             SEEK 3, K& * 64 + 1
             GET #3, , NA0: GET #3, , NOME
             IF NA = 0 THEN
                UC0$ = RTRIM$(UCASE$(NOME))
                IF I = 1 THEN
                   IF UC$ = UC0$ THEN FLFOUND = -1: EXIT FOR
                ELSEIF INSTR(UC0$, UC$) THEN
                   FLFOUND = -1: EXIT FOR
                END IF
             ELSEIF NA = NA0 THEN
                FLFOUND = -1: EXIT FOR
             END IF
          NEXT K&
          IF FLFOUND THEN EXIT FOR
       NEXT I
       IF FLFOUND THEN
          IF NA0 = 3 OR NA0 > 4 OR NA0 = 0 THEN
             NADD = NADD + 1
             GET #3, , MG%: MAGP = MG% / 100#
             IF NA0 THEN STB% = 1 ELSE STB% = 3
             PRINT #4, "0 ";
             PRINT #4, USING "###.## "; MAGP;
             PRINT #4, "0  "; MID$(NOME, STB%, 8)
             FOR J = 1 TO 6
               GET #3, , X: PRINT #4, X; : IF J MOD 3 = 0 THEN PRINT #4, ""
             NEXT
             IF NADD > 4 THEN EXIT DO
          END IF
       END IF
    LOOP
    IF NADD > 0 THEN FLFOUND = -1 ELSE FLFOUND = 0
    CLOSE #3: CLOSE #4

END SUB

SUB FixBaric (N, M0(), X(), Y(), Z(), VX(), VY(), VZ(), XC0, YC0, ZC0)
     SHARED NTER, NLUN
     SM = 0: SMX = 0: SMY = 0: SMZ = 0: SVX = 0: SVY = 0: SVZ = 0
     FOR I = 0 TO N
        SM = SM + M0(I)
        IF I = NLUN THEN
           XL = X(I): YL = Y(I): ZL = Z(I): VXL = VX(I): VYL = VY(I): VZL = VZ(I)
           X(I) = XL + X(NTER): Y(I) = YL + Y(NTER): Z(I) = ZL + Z(NTER): VX(I) = VXL + VX(NTER): VY(I) = VYL + VY(NTER): VZ(I) = VZL + VZ(NTER)
        END IF
        SMX = SMX + X(I) * M0(I): SMY = SMY + Y(I) * M0(I): SMZ = SMZ + Z(I) * M0(I): SVX = SVX + VX(I) * M0(I): SVY = SVY + VY(I) * M0(I): SVZ = SVZ + VZ(I) * M0(I)
     NEXT I
     IF SM = 0 THEN SM = 1
     XC0 = SMX / SM: YC0 = SMY / SM: ZC0 = SMZ / SM: VCX = SVX / SM: VCY = SVY / SM: VCZ = SVZ / SM
     FOR I = 0 TO N
        IF I = NLUN THEN
           X(I) = XL: Y(I) = YL: Z(I) = ZL: VX(I) = VXL: VY(I) = VYL: VZ(I) = VZL
        ELSE
           VX(I) = (VX(I) - VCX): VY(I) = (VY(I) - VCY): VZ(I) = (VZ(I) - VCZ)
        END IF
     NEXT I

END SUB

SUB FixV (N, DT, DT0, VX(), VY(), VZ()) STATIC
    SHARED CV, GRAV0, GRAV, M(), M0()
    GRAV = GRAV0 * DT * DT
    FOR I = 0 TO N
       M(I) = M0(I) * GRAV
       VX(I) = VX(I) * DT / DT0
       VY(I) = VY(I) * DT / DT0
       VZ(I) = VZ(I) * DT / DT0
    NEXT

END SUB

'
'... This subroutine reads from the file PL???.BIN the starting conditions
'... closest in time to the target epoch and returns the number of days
'... needed to reach the target epoch.
'... Within PL???.BIN, starting positions and velocities are stored in
'... compressed form. Each coordinate is stored in a total of 6 bytes
'... (one INT + one LONG)
'... Scale factors are used to convert them to their final values.
'
SUB GetStart (TJ, Nday#, NUM, NST, DjLib() AS LONG, M0(), Rj0(), ABSM(), PLN() AS STRING * 8, X(), Y(), Z(), VX(), VY(), VZ())
    SHARED FLDE, QJ6, QJ7, DR1(), DR2(), DR3(), DR4()
    DIM XX(2, 3)
    maxlong = 4294967295#
    IF FLDE = 1 THEN SSL$ = "PLANS.BIN": PLD$ = "PLANSB.DAT"
    IF FLDE = 2 THEN SSL$ = "PL406B.BIN": PLD$ = "PL406B.DAT"
    OPEN SSL$ FOR BINARY AS #4
    GET #4, , NUM: LBLOCK = NUM * 36 + 4
    T = 999999999: P& = 3
    CLOSE #4
    FOR J = 1 TO NST
       DJ& = DjLib(J)
       IF TJ = DJ& THEN PF& = P&: T = TJ: EXIT FOR
       IF ABS(TJ - DJ&) < ABS(TJ - T) THEN T = CDBL(DJ&): PF& = P& '...T is nearest day
       P& = P& + LBLOCK
    NEXT J
    Nday# = TJ - T
    N$ = SSL$
    OPEN PLD$ FOR INPUT AS #1
    OPEN N$ FOR BINARY AS #2: SEEK #2, PF&
    GET #2, , DJ&: TJ = DJ&  ' .... DJ& is the number of Jdays from J2000
    IF FLDE = 2 THEN
       IF DJ& > -450 THEN
          QJ6 = -1.17889D-19
          QJ7 = -2.02562D-20
          DR1(4) = .000000000982415#: DR2(4) = -.000000000498614#
          DR3(4) = 8.50072D-11: DR4(4) = -4.54639D-12
          DR1(5) = 1.23706D-11
       ELSE
          QJ6 = -2.79292D-19
          QJ7 = 1.49958D-21
          DR1(4) = -.000000000122642#: DR2(4) = .000000000012064#
          DR3(4) = -3.94359D-13: DR4(4) = 4.03462D-15
          DR1(5) = -9.14339D-12
       END IF
    END IF
    INPUT #1, M0(0), Rj0(0), PLN(0)
    X(0) = 0: Y(0) = 0: Z(0) = 0: VX(0) = 0: VY(0) = 0: VZ(0) = 0
    FOR I = 1 TO NUM
       INPUT #1, M0(I), Rj0(I), ABSM(I), SF, SFV, PLN(I)
       ABSM(I) = ABSM(I) / 100
       FOR K = 1 TO 2
          FOR J = 1 TO 3
             GET #2, , Xint%: GET #2, , Xlong&
             IF Xlong& < 0 THEN xfr = Xlong& + maxlong ELSE xfr = Xlong&
             XX(K, J) = (Xint% + xfr / maxlong) / SF
          NEXT J
          SF = SFV
       NEXT K
       X(I) = XX(1, 1): Y(I) = XX(1, 2): Z(I) = XX(1, 3)
       VX(I) = XX(2, 1): VY(I) = XX(2, 2): VZ(I) = XX(2, 3)
    NEXT I
    CLOSE

END SUB

SUB Header
   SCREEN 0: COLOR 14, 1, 1
   CALL DrawFrame(14, 1)
   CALL DrawTitle(11, 2, 14)
   V$ = ""
   COLOR 14, 1: LOCATE 7, 12
   PRINT "Computation of Fundamental and Minor Planets' Ephemerides";
   LOCATE 8, 18: PRINT "by Numerical Integration of the Solar System";
   LOCATE 9, 24: PRINT "using a Fast Extrapolation Method";
   COLOR 11
   LOCATE 11, 18: PRINT "(A. Vitagliano, Cel. Mech. 1997, 66, 293-308)";
   COLOR 14
   LOCATE 13, 17: PRINT "adapted to both DE200 and DE406 JPL Ephemerides";
   COLOR 11, 4
   LOCATE 15, 14: PRINT TAB(20); "Copyright (C) 1996-99 by Aldo Vitagliano"; TAB(68);
   LOCATE 16, 14: PRINT " Dip. di Chimica, Universita' di Napoli 'Federico II'"; TAB(68);
   LOCATE 17, 14: PRINT TAB(20); "via Mezzocannone 4, 80134 Napoli, Italy"; TAB(68);
   LOCATE 18, 14: PRINT TAB(27); " E-mail: alvitagl@unina.it"; TAB(68);
   COLOR 10, 1
   LOCATE 20, 25: PRINT "This program is of PUBLIC DOMAIN ";
   LOCATE 21, 8: PRINT "It can be freely copied and distributed, subject to the condition";
   LOCATE 22, 7: PRINT "that the program is not modified and this header is left unchanged.";
   COLOR 12
   LOCATE 24, 29: PRINT "Press a key to continue ... ";
   WHILE V$ = "": V$ = INKEY$: WEND
   COLOR 14
END SUB

SUB InputFile (NF$)

29 CLS
   LOCATE 10, 24: INPUT ; " Enter Data FileName: ", NF$
   IF NF$ = "" THEN EXIT SUB
   OPEN NF$ FOR APPEND AS #1: L = LOF(1): CLOSE #1
   IF L < 3 THEN KILL NF$: GOTO 29

END SUB

'........... Adjusts reference frame if time is right for
'........... and calls the integration routine STORM8
SUB Integrate (Njump&, N, T, X(), Y(), Z(), VX(), VY(), VZ())
    SHARED MS%, S%, DT, DJP, FLDE
    IF Njump& > 100 THEN COLOR 12, 0: LOCATE 14, 14: PRINT SPACE$(54);
    TINC = MS% * S% * DT: T0 = TIMER: T1 = T0
    FOR JJ& = 1 TO Njump&
       IF ABS(T - DJP) > 30 THEN
          DJP = DJP + 60 * S%
          CALL RotInv(N, X(), Y(), Z()): CALL RotInv(N, VX(), VY(), VZ())
          CALL PrecAng(DJP, 1, 0, FLDE)
          CALL RotAx(N, X(), Y(), Z(), 0): CALL RotAx(N, VX(), VY(), VZ(), 0)
          IF Njump& > 100 THEN
             BAR% = 54 * (CSNG(JJ&) / CSNG(Njump&))
             COLOR 12: LOCATE 14, 14: PRINT STRING$(BAR%, 178);
          END IF
       END IF
       CALL STORM8(N, X(), Y(), Z(), VX(), VY(), VZ())
       T = T + TINC
    NEXT JJ&
    COLOR 14, 1
END SUB

SUB Nutat (TJ, DPSI, DEPS)

      ' .... Computes Nutation in Longitude and Obliquity according to the
      ' .... IAU series for Nutation, truncated to the terms larger than
      ' .... 0.25 milliarcsecs. Final values of DPSI and DEPS are in radians.

      SHARED PI
      DR = PI / 180: SJ = TJ / 36525#
      SJ2 = SJ * SJ: SJ3 = SJ2 * SJ

      ' .... Fundamental Arguments (Meeus, Astronomical Algorithms, p. 132)

      D = (297.85036# + 445267.11148# * SJ - .0019142# * SJ2 + SJ3 / 189474#) * DR
      M = (357.52772# + 35999.05034# * SJ - .0001603# * SJ2 - SJ3 / 300000#) * DR
      M1 = (134.96298# + 477198.867398# * SJ + .0086972# * SJ2 + SJ3 / 56250#) * DR
      DF = (93.27191# + 483202.017538# * SJ - .0036825# * SJ2 + SJ3 / 327270#) * DR
      OM = (125.04452# - 1934.136261# * SJ + .0020708# * SJ2 + SJ3 / 450000#) * DR

      '............. Series for Longitude

      DPSI = -(171996 + 174.2# * SJ) * SIN(OM) - (13187 + 1.6 * SJ) * SIN(-2 * D + 2 * DF + 2 * OM) - (2274 + .2 * SJ) * SIN(2 * DF + 2 * OM) + (2062 + .2 * SJ) * SIN(2 * OM) + (1426 - 3.4 * SJ) * SIN(M) + (712 + .1 * SJ) * SIN(M1)
      DPSI = DPSI + (-517 + 1.2 * SJ) * SIN(-2 * D + M + 2 * DF + 2 * OM) - (386 + .4 * SJ) * SIN(2 * DF + OM) - 301 * SIN(M1 + 2 * DF + 2 * OM) + (217 - .5 * SJ) * SIN(-2 * D - M + 2 * DF + 2 * OM) - 158 * SIN(-2 * D + M1)
      DPSI = DPSI + (129 + .1 * SJ) * SIN(-2 * D + 2 * DF + OM) + 123 * SIN(-M1 + 2 * DF + 2 * OM) + 63 * SIN(2 * D) + (63 + .1 * SJ) * SIN(M1 + OM) - 59 * SIN(2 * D - M1 + 2 * DF + 2 * OM) - (58 + .1 * SJ) * SIN(-M1 + OM)
      DPSI = DPSI - 51 * SIN(M1 + 2 * DF + OM)
      DPSI = DPSI + 48 * SIN(-2 * D + 2 * M1) + 46 * SIN(-2 * M1 + 2 * DF + OM) - 38 * SIN(2 * D + 2 * DF + 2 * OM) - 31 * SIN(2 * M1 + 2 * DF + 2 * OM) + 29 * SIN(2 * M1) + 29 * SIN(-2 * D + M1 + 2 * DF + 2 * OM) + 26 * SIN(2 * DF)
      DPSI = DPSI - 22 * SIN(2 * DF - 2 * D) + 21 * SIN(2 * DF - M1) + (17 - .1# * SJ) * SIN(2 * M) + 16 * SIN(2 * D - M1 + OM) - (16 - .1# * SJ) * SIN(2 * (OM + DF + M - D)) - 15 * SIN(M + OM) - 13 * SIN(OM + M1 - 2 * D) - 12 * SIN(OM - M)
      DPSI = DPSI + 11 * SIN(2 * (M1 - DF)) - 10 * SIN(2 * D - M1 + 2 * DF) - 8 * SIN(2 * D + M1 + 2 * DF + 2 * OM) + 7 * SIN(M + 2 * DF + 2 * OM) - 7 * SIN(M + M1 - 2 * D) - 7 * SIN(2 * DF + 2 * OM - M) - 7 * SIN(2 * D + 2 * DF + OM)
      DPSI = DPSI + 6 * SIN(2 * D + M1)
      DPSI = DPSI + 6 * SIN(2 * (OM + DF + M1 - D)) + 6 * SIN(OM + 2 * DF + M1 - 2 * D) - 6 * SIN(2 * D - 2 * M1 + OM) - 6 * SIN(2 * D + OM) + 5 * SIN(M1 - M) - 5 * SIN(OM + 2 * DF - M - 2 * D) - 5 * SIN(OM - 2 * D) - 5 * SIN(OM + 2 * DF + 2 * M1)
      DPSI = DPSI + 4 * SIN(OM - 2 * M1 - 2 * D) + 4 * SIN(OM + 2 * DF + M - 2 * D) + 4 * SIN(M1 - 2 * DF) - 4 * SIN(M1 - D) - 4 * SIN(M - 2 * D) - 4 * SIN(D) + 3 * SIN(2 * DF + M1) - 3 * SIN(2 * (OM + DF - M1)) - 3 * SIN(M1 - M - D)
      DPSI = DPSI - 3 * SIN(M1 + M)
      DPSI = DPSI - 3 * SIN(2 * OM + 2 * DF + M1 - M) - 3 * SIN(2 * OM + 2 * DF - M1 - M + 2 * D) - 3 * SIN(2 * OM + 2 * DF + 3 * M1) - 3 * SIN(2 * OM + 2 * DF - M + 2 * D)
      DPSI = DPSI * DR / 36000000#

      '............. Series for Obliquity

      DEPS = (92025 + 8.9 * SJ) * COS(OM) + (5736 - 3.1 * SJ) * COS(-2 * D + 2 * DF + 2 * OM) + (977 - .5 * SJ) * COS(2 * DF + 2 * OM) + (-895 + .5 * SJ) * COS(2 * OM) + (54 - .1 * SJ) * COS(M) - 7 * COS(M1)
      DEPS = DEPS + (224 - .6 * SJ) * COS(-2 * D + M + 2 * DF + 2 * OM) + 200 * COS(2 * DF + OM) + (129 - .1 * SJ) * COS(M1 + 2 * DF + 2 * OM) + (-95 + .3 * SJ) * COS(-2 * D - M + 2 * DF + 2 * OM) - 70 * COS(-2 * D + 2 * DF + OM)
      DEPS = DEPS - 53 * COS(-M1 + 2 * DF + 2 * OM) - 33 * COS(M1 + OM) + 26 * COS(2 * D - M1 + 2 * DF + 2 * OM) + 32 * COS(-M1 + OM) + 27 * COS(M1 + 2 * DF + OM) - 24 * COS(-2 * M1 + 2 * DF + OM)
      DEPS = DEPS + 16 * COS(2 * (D + DF + OM)) + 13 * COS(2 * (M1 + DF + OM)) - 12 * COS(2 * OM + 2 * DF + M1 - 2 * D) - 10 * COS(OM + 2 * DF - M1) - 8 * COS(2 * D - M1 + OM) + 7 * COS(2 * (OM + DF + M - D)) + 9 * COS(M + OM)
      DEPS = DEPS + 7 * COS(OM + M1 - 2 * D) + 6 * COS(OM - M) + 5 * COS(OM + 2 * DF - M1 + 2 * D) + 3 * COS(2 * OM + 2 * DF + M1 + 2 * D) - 3 * COS(2 * OM + 2 * DF + M) + 3 * COS(2 * OM + 2 * DF - M) + 3 * COS(OM + 2 * DF + 2 * D)
      DEPS = DEPS - 3 * COS(2 * (OM + DF + M1 - D)) - 3 * COS(OM + 2 * DF + M1 - 2 * D) + 3 * COS(OM - 2 * M1 + 2 * D) + 3 * COS(OM + 2 * D) + 3 * COS(OM + 2 * DF - M - 2 * D) + 3 * COS(OM - 2 * D) + 3 * COS(OM + 2 * DF + 2 * M1)
      DEPS = DEPS * DR / 36000000#

END SUB

SUB OutHeadings (FLDE, FLORT, FLGEO, FLTOP, FLHOR, FLECL, FLPR, FLNUTHIGH, FLABP, FLABS, GLAT, GLON)
SHARED RD, COORD$, Equ$, CoordType$, Origin$, Frame$, Ref$, Nut$, Aberr$

   COORD$ = "Coordinates "
   Equ$ = "Equinox "
   IF FLDE = 1 THEN DE$ = "DE200": CLDE% = 10
   IF FLDE = 2 THEN DE$ = "DE406": CLDE% = 13
   IF FLORT THEN CoordType$ = "Orthogonal " ELSE CoordType$ = "  Angular  "
   IF FLGEO THEN Origin$ = "Geocentric  " ELSE Origin$ = "Heliocentric "
   IF FLTOP THEN Origin$ = "Topocentric "
   IF FLECL THEN Frame$ = " Ecliptic  " ELSE IF FLHOR = 0 THEN Frame$ = "Equatorial " ELSE Frame$ = "Horizontal "
   IF FLPR THEN Ref$ = "of Date " ELSE Ref$ = "of J2000"
   IF FLNUTHIGH THEN Nut$ = "True " ELSE Nut$ = "Mean "
   IF FLABS THEN
      Aberr$ = " (Corrected for Light-time) "
      IF FLABP THEN Aberr$ = " (Corrected for Aberr. & L-time)"
   ELSE
      Aberr$ = " (Geometric Positions)  "
   END IF

   COLOR 15, 4
   LOCATE 1, 1: PRINT TAB(5); Origin$ + CoordType$ + Frame$ + COORD$; TAB(54);
   LOCATE 2, 1: PRINT Nut$ + Equ$ + Ref$ + Aberr$; TAB(54);
   LOCATE 3, 1: COLOR CLDE%, 1: PRINT DE$; TAB(49); DE$;
   COLOR 15, 1

   FLAPP = FLGEO AND FLNUTHIGH AND FLABP AND (FLECL = 0)
   FLAST = FLGEO AND (FLPR = 0) AND FLABS AND (FLABP = 0) AND (FLECL = 0) AND (FLTOP = 0)
   DefCoord$ = SPACE$(29)
   IF FLAPP AND FLTOP = 0 THEN DefCoord$ = " - Apparent Coordinates -  "
   IF FLAST THEN DefCoord$ = "- Astrometric Coordinates -"
   IF FLGEO = 0 THEN DefCoord$ = "  - Moon is Geocentric -   "
   LOCATE 3, 11
   IF FLTOP THEN
      IF SGN(GLON) THEN EW$ = "E" ELSE EW$ = "W"
      IF SGN(GLAT) THEN SN$ = "N" ELSE SN$ = "S"
      PRINT USING "Lat =###.##" + CHR$(248) + " " + SN$; ABS(GLAT * RD);
      PRINT USING "   Lon =####.##" + CHR$(248) + " " + EW$ + SPACE$(2); ABS(GLON * RD);
   ELSE
      PRINT "   "; DefCoord$; SPACE$(5);
   END IF
   COLOR 14

END SUB

SUB PlanMag (N, XO(), YO(), ZO(), ABSM(), PH(), MAGP())
   SHARED AU

   '.... XO, YO and ZO are in units of Gm

   XE = XO(3) - XO(0): YE = YO(3) - YO(0): ZE = ZO(3) - ZO(0)
   RE2 = XE * XE + YE * YE + ZE * ZE
   FOR I = 1 TO N
      IF I <> 3 THEN
         DXS = XO(I) - XO(0): DYS = YO(I) - YO(0): DZS = ZO(I) - ZO(0)
         DXE = XO(I) - XO(3): DYE = YO(I) - YO(3): DZE = ZO(I) - ZO(3)
         R2 = DXS * DXS + DYS * DYS + DZS * DZS
         DELTA2 = DXE * DXE + DYE * DYE + DZE * DZE
         R = SQR(R2): DELTA = SQR(DELTA2): RDEL = R * DELTA
         COSI = (R2 + DELTA2 - RE2) / (2 * RDEL)
         PH(I) = (1 + COSI) / 2  '..... Illuminated Fraction
         XMAG = LOG(.0000446# * RDEL / SQR(PH(I)))
         MAGP(I) = 2.1715 * XMAG + ABSM(I)
      END IF
   NEXT I
END SUB

SUB PrecAng (DJP, FLNUTLOW, FLNUTHIGH, FLDE)

'... Precession angles according to Laskar, Astr. Astrophys. 1986, 157, 59-70.
'... in the case of DE406 according to Williams, Astron. J. 1994, 108, 711.

'... During the integration a nutation correction is applied limited to the
'... first term in obliquity.
'... Full nutation is called if the Nutation toggle is active.

  SHARED PI, SEQ, CEQ, SEI, CEI, SPA, CPA, SEN, CEN, QJ4, QJ6, QJ7, QJ8, Qtid
  SHARED Rj00(), Rj0(), DR1(), DR2(), DR3(), DR4(), NLUN, NTER

  DEPS = 0: DPSI = 0
  MJ = DJP / 365250#: SJ = MJ * 10
  ome = (125.0445222# + SJ * (-1934.1362608# + SJ * (.00207833# + SJ * 2.22E-06))) * PI / 180
  come = COS(ome): some = SIN(ome)
  SJR = SJ: IF SJR < -51 THEN SJR = -51
  IF SJR > 11 THEN SJR = 11
  SJR2 = SJR * SJR
  FOR I = 0 TO 20: Rj0(I) = Rj00(I): NEXT

  '... Empirical series of small corrections to adjust on DE406
  IF FLDE = 2 THEN FOR I = 0 TO 20: Rj0(I) = Rj0(I) + DR1(I) * ABS(SJR) + DR2(I) * SJR2 + DR3(I) * SJR2 * ABS(SJR) + DR4(I) * SJR2 * SJR2: NEXT

  Rj0(NLUN) = Rj0(NTER)
  IF FLNUTLOW THEN
     DEPS = 9.2 * come * PI / 648000 ' Nutation change in Obliquity
     DPSI = -17.2 * some * PI / 648000 'Nutation change in Longitude
  END IF
  IF FLNUTHIGH THEN CALL Nutat(DJP, DPSI, DEPS) '... Full Nutation Correction

  PA = (((((-.000000000866# * MJ - .00000004759#) * MJ + .0000002424#) * MJ + .000013095#) * MJ + .00017451#) * MJ - .0018055#) * MJ - .235316#
  EN = (((6.6402D-16 * MJ - 2.69151D-15) * MJ - 1.547021D-12) * MJ + 7.521313D-12) * MJ
  EI = (((((1.2147D-16 * MJ + 7.3759D-17) * MJ - 8.26287D-14) * MJ + 2.50341D-13) * MJ + 2.4650839D-11) * MJ - 5.4000441D-11) * MJ + 1.32115526D-09
  EQ0 = (((((.000000000245# * MJ + .00000000579#) * MJ + .0000002787#) * MJ + .000000712#) * MJ - .00003905#) * MJ - .0024967#) * MJ - .005138#

  IF FLDE = 2 THEN
     PA = (((PA * MJ + .076#) * MJ + 110.5414#) * MJ + 50287.91959#) * MJ * PI / 648000#
     EN = ((((((EN + .00000000019#) * MJ - .00000000354#) * MJ - .00000018103#) * MJ + .000000126#) * MJ + .00007436169#) * MJ - .04207794833#) * MJ + 3.052115282424#
     EI = (((EI * MJ - .0000006012#) * MJ - .0000162442#) * MJ + .00227850649#) * MJ
     EQANG = ((((EQ0 * MJ + 1.9989) * MJ - .0175) * MJ - 468.3396#) * MJ + 84381.406173#) * PI / 648000#
  ELSE
     PA = (((PA * MJ + .07732#) * MJ + 111.1971) * MJ + 50290.966#) * MJ * PI / 648000#
     EN = ((((((EN + 6.3190131D-10) * MJ - 3.48388152D-09) * MJ - 1.813065896D-07) * MJ + 2.75036225D-08) * MJ + .000074394531426#) * MJ - .042078604317#) * MJ + 3.052112654975#
     EI = (((EI * MJ - 5.998737027D-07) * MJ - .000016242797091#) * MJ + .002278495537#) * MJ
     EQANG = ((((EQ0 * MJ + 1.99925) * MJ - .0155) * MJ - 468.093) * MJ + 84381.448#) * PI / 648000#
  END IF

  EQANG = EQANG + DEPS
  PAS = PA + EN + DPSI

  SEQ = SIN(EQANG): CEQ = COS(EQANG): SEI = SIN(EI): CEI = COS(EI)
  SPA = SIN(PAS): CPA = COS(PAS): SEN = SIN(EN): CEN = COS(EN)

  Qtid = QJ4 '.... If DE406 tidal parameters changes with time
  IF FLDE = 2 THEN Qtid = QJ4 + SJR * (QJ6 + SJR * QJ7)

END SUB

SUB PrintAngOut (FLGEO, FL1ST, N, T, ABSM(), PLN() AS STRING * 8, XO(), YO(), ZO(), R(), YY%, MM%, DD%, HH%, MN%, SS%)

SHARED AU, RD, FLECL, FLDEC, FLHOR, COORD$, Equ$, CoordType$, Origin$, Frame$, Ref$, Nut$, Aberr$
SHARED PI, GLAT, GLON, TSD, FLJD, FLWRITE

DIM MAGP(20), PH(20)

DJW = T + 2451545#
FORMOON$ = "###.########"
IF FLDEC THEN FORM1$ = "#####.####" ELSE FORM1$ = "###.######"
FORM2$ = "### ## ##"
FORM3$ = "## ## ##."
IF FLDEC = 0 THEN
   DMS$ = "(" + CHR$(248) + " ' " + CHR$(34) + ")"
   HHH$ = "(h m s)": DUN$ = "(AU)"
ELSE
   DMS$ = " ( " + CHR$(248) + " ) "
   HHH$ = " ( " + CHR$(248) + " ) "
   DUN$ = "(Gm)"
END IF

FORM10$ = "####.######    "
FORM11$ = " ####.###      "

HEAD$ = "   RA " + HHH$ + "    Decl " + DMS$ + "    True D " + DUN$
HEADECL$ = "  long " + DMS$ + "    lat " + DMS$ + "    True D " + DUN$
HEADHOR$ = "   Azim (" + CHR$(248) + ")        alt (" + CHR$(248) + ")      Phase   Magn"
LOGF = N + 2
IF N > 13 THEN I0 = N - 13 ELSE I0 = 0
IF FLWRITE = 0 THEN
   LOCATE 5, 1: COLOR 15, 4: PRINT "Body"; TAB(10);
   IF FLECL THEN PRINT HEADECL$ ELSE IF FLHOR = 0 THEN PRINT HEAD$ ELSE PRINT HEADHOR$
END IF
IF FLGEO THEN NC0 = 3 ELSE NC0 = 0
IF FLHOR THEN
   SGLON = SIN(GLON): SGLAT = SIN(GLAT)
   CGLON = COS(GLON): CGLAT = COS(GLAT)
   CALL PlanMag(N, XO(), YO(), ZO(), ABSM(), PH(), MAGP())
END IF
L = 5: NC = NC0
COLOR 10, 1
FOR I = 0 TO N
    NF = I + 1
    L = L + 1: '... Set printing line
    IF I = NC0 THEN L = L - 1: GOTO 10 '... if Origin body then go to next planet

    '..................... Computes angular coordinates

    X = XO(I) - XO(NC): Y = YO(I) - YO(NC): Z = ZO(I) - ZO(NC)
          '... Moon has always geocentric coordinates
    IF I = N THEN X = XO(I) - XO(3): Y = YO(I) - YO(3): Z = ZO(I) - ZO(3)

    RXY = SQR(X * X + Y * Y)
    RA = FNATN2(Y, X)
    DEC = ATN(Z / RXY)
    IF (I < N AND FLDEC = 0) THEN RHO = R(I) / AU ELSE RHO = R(I)'.. if Moon unit is Gm

    IF FLHOR THEN  '... Compute Horizontal Coordinates
       SDEC = SIN(DEC): CDEC = COS(DEC)
       ANGOR = TSD - RA: SINH = SIN(ANGOR): COSH = COS(ANGOR)
       SINALT = SGLAT * SDEC + CGLAT * CDEC * COSH: IF ABS(SINALT) > 1 THEN SINALT = SGN(SINALT)
       COSALT = SQR(1 - SINALT * SINALT)
       SINAZ = -CDEC * SINH / COSALT
       COSAZ = -(SDEC * CGLAT - CDEC * SGLAT * COSH) / COSALT
       RA = -FNATN2(SINAZ, COSAZ)
       IF COSALT <> 0 THEN DEC = ATN(SINALT / COSALT) ELSE DEC = PI * SGN(SINALT)
       
    END IF

    '... Output Precision is set to 0.1 arcsec for Asteroids,
    '... to 0.01 arcsecs for Planets and Moon

    IF I > 9 AND I < N THEN DF% = 1: HF% = 2 ELSE DF% = 2: HF% = 3
    IF FLHOR THEN DF% = 0: HF% = 1
    NPRECD = 10 ^ DF%: NPRECH = 10 ^ HF%

    IF FLECL = 0 THEN
       CALL RadToHMS(RA, HR, HMIN, SEC, NPRECH)
    ELSE
       CALL RadToDMS(RA, SG$, EDEG%, EMIN%, ESEC, NPRECD, 1)
    END IF

    CALL RadToDMS(DEC, SG$, DEG%, MINUT%, ARCSEC, NPRECD, 0)

    FORMH$ = FORM3$ + STRING$(HF%, "#") + " "
    IF DF% > 0 THEN
       FORMD$ = FORM2$ + "." + STRING$(DF%, "#")
    ELSE
       FORMD$ = FORM2$
    END IF
    IF I > 9 AND I < N THEN
       FORMR$ = FORM1$ + "  ": CLB% = 2: CLP% = 11
    ELSEIF I = N THEN
       FORMR$ = FORMOON$: CLB% = 1: CLP% = 10
    ELSE
       FORMR$ = FORM1$ + "##": CLB% = 1: CLP% = 10
    END IF
    IF FLWRITE = 0 THEN
       LOCATE L, 1
       COLOR CLP%, CLB%: PRINT PLN(I); "   "; : COLOR 14, CLB%
       IF FLDEC = 0 AND FLHOR = 0 THEN
          IF FLECL = 0 THEN PRINT USING FORMH$; HR; HMIN; SEC;
          IF FLECL THEN PRINT USING FORMD$; EDEG%; EMIN%; ESEC; '... Ecl Long
          PRINT TAB(27); : PRINT SG$;
          FORMDEC$ = MID$(FORMD$, 2)
          PRINT USING FORMDEC$; DEG%; MINUT%; ARCSEC;
       ELSE
          IF FLHOR = 0 THEN FMDG$ = FORM10$ ELSE FMDG$ = FORM11$
          PRINT USING FMDG$; RA * RD; DEC * RD;
       END IF
       PRINT TAB(42);
       IF FLHOR = 0 THEN
          PRINT USING FORMR$; RHO;
       ELSE
          IF I > 0 THEN PRINT USING "##.##  ###.#"; PH(I); MAGP(I);  ELSE PRINT SPACE$(12)
       END IF
       IF I = N THEN
          IF FLHOR = 0 THEN COLOR 15, 1: PRINT " Gm "; : COLOR 14, 1 ELSE PRINT SPACE$(4);
       END IF

    ELSE
       PRINT #LOGF, USING "## "; I; : PRINT #LOGF, PLN(I); "  ";
       IF FLDEC = 0 AND FLHOR = 0 THEN
          IF FLECL = 0 THEN PRINT #LOGF, USING FORMH$; HR; HMIN; SEC;
          IF FLECL THEN PRINT #LOGF, USING FORMD$; EDEG%; EMIN%; ESEC; '... Ecl Long
          PRINT #LOGF, TAB(32); : PRINT #LOGF, SG$;
          PRINT #LOGF, USING FORMD$; DEG%; MINUT%; ARCSEC;
       ELSE
          IF FLHOR = 0 THEN FMDG$ = FORM10$ ELSE FMDG$ = FORM11$
          PRINT #LOGF, USING FMDG$; RA * RD; DEC * RD;
       END IF
          PRINT #LOGF, TAB(51);
          IF FLHOR = 0 THEN
             PRINT #LOGF, USING FORMR$; RHO;
             IF I = N THEN PRINT #LOGF, " (Gm)"
          ELSE
             IF I > 0 THEN PRINT #LOGF, USING "##.##  ###.#"; PH(I); MAGP(I);
          END IF
       PRINT #LOGF, ""
       IF I >= I0 THEN
          IF FL1ST THEN
             PRINT #NF, Origin$ + CoordType$ + Frame$ + COORD$ + SPACE$(10)
             PRINT #NF, Nut$ + Equ$ + Ref$ + Aberr$
          END IF
          IF FLJD = 0 THEN
             PRINT #NF, USING "#####/##/##"; YY%; MM%; DD%;
             PRINT #NF, USING "  ##:##:##  "; HH%; MN%; SS%;
          ELSE
             PRINT #NF, USING "########.######  "; DJW;
          END IF
          IF FLDEC = 0 AND FLHOR = 0 THEN
             IF FLECL = 0 THEN PRINT #NF, USING FORMH$; HR; HMIN; SEC;
             IF FLECL THEN PRINT #NF, USING FORMD$; EDEG%; EMIN%; ESEC; '... Ecl Long
             PRINT #NF, "    "; : PRINT #NF, SG$;
             PRINT #NF, USING FORMD$; DEG%; MINUT%; ARCSEC;
          ELSE
             IF FLHOR = 0 THEN FMDG$ = FORM10$ ELSE FMDG$ = FORM11$
             PRINT #NF, USING FMDG$; RA * RD; DEC * RD;
          END IF
          IF FLHOR = 0 THEN PRINT #NF, "   "; : PRINT #NF, USING FORMR$; RHO
          IF FLHOR THEN PRINT #NF, USING "##.##  ###.#"; PH(I); MAGP(I)
       END IF

    END IF
10 '
NEXT I

END SUB

SUB PrintOut (Djump, TimUnName() AS STRING * 4, ABSM(), PLN() AS STRING * 8, FL1ST, FLGEO, N, T, TC, XO(), YO(), ZO(), R(), SCAL%)
   SHARED AU, FLDEC, FLWRITE, FLSTEP, FLDT, FLJD, FLHOR, FLORT, FLECL, COORD$, Equ$, CoordType$, Origin$, Frame$, Ref$, Nut$, Aberr$, TIMUN%, TimUnFactor(), TabStep0

   LAU$ = "   ###.########"
   IF FLDEC = 0 THEN
      L$ = LAU$: S$ = "   ###.######  "
      HEADPL$ = "    X (AU)         Y (AU)         Z (AU)"
   ELSE
      L$ = "  ######.######": S$ = "  ######.####  "
      HEADPL$ = "    X (Gm)         Y (Gm)         Z (Gm)"
   END IF
   IF FLDT = 0 THEN TUT$ = "TDT" ELSE TUT$ = "UT "
   IF FLHOR THEN TUT$ = "LT "
   LOGF = N + 2
   DELTA = DeltaT(T)
   DJW = T + 2451545#
   CALL TToDate(T + .5 - DELTA, A%, M%, D%, H%, MINUT%, SEC%)
   IF FLGEO THEN NC0 = 3 ELSE NC0 = 0
   X0 = XO(NC0): Y0 = YO(NC0): Z0 = ZO(NC0)

   IF FLWRITE THEN
      IF N > 13 THEN I0 = N - 13 ELSE I0 = 0
      IF FL1ST THEN
         FLSTEP = -1: CLOSE
         FOR I = I0 TO N
            NF = I + 1: FILNAM$ = PLN(I) + ".OUT"
            OPEN FILNAM$ FOR OUTPUT AS NF
         NEXT I
         LOGF = N + 2
         OPEN "OUTPUT.LOG" FOR OUTPUT AS LOGF

         PRINT #LOGF, Origin$ + CoordType$ + Frame$ + COORD$ + SPACE$(10)
         PRINT #LOGF, Nut$ + Equ$ + Ref$ + Aberr$
      END IF
      PRINT #LOGF, ""
      PRINT #LOGF, USING "Date:#####/##/##"; A%; M%; D%;
      PRINT #LOGF, USING "   " + TUT$ + " :   ##:##:##"; H%; MINUT%; SEC%
      PRINT #LOGF, ""
   END IF

   LOCATE 1, 64: PRINT USING "Date:#####/##/##"; A%; M%; D%;
   LOCATE 2, 64: PRINT USING TUT$ + " :   ##:##:##"; H%; MINUT%; SEC%;
   LOCATE 3, 64: PRINT USING "Step :#####.## "; TabStep0 * 10 ^ SCAL%;
   PRINT LEFT$(TimUnName(TIMUN%), 1);
   LOCATE 4, 64: IF FLDT THEN PRINT USING "DeltaT =###### s"; DELTA * 86400;  ELSE PRINT SPACE$(16);

   LOCATE N + 6, 1: COLOR 15, 4: PRINT TAB(15);
   TJMP = INT((TIMER - TC) * 10 + .5) / 10
   IF Djump < 3653 THEN
      PRINT USING "#### "; Djump * TimUnFactor(TIMUN%);
      PRINT LEFT$(TimUnName(TIMUN%), 1); " jumped in";
      PRINT TJMP; "s";
   ELSE
      Yjump = Djump / 365.25: MJMP = TJMP / 60
      PRINT USING "####.# y jumped in###.## min"; Yjump; MJMP;
   END IF
   PRINT TAB(54);
   COLOR 14, 1
   IF FLORT THEN
      IF FLWRITE = 0 THEN LOCATE 5, 1: COLOR 15, 4: PRINT "Body"; TAB(12); HEADPL$; "  ";
      COLOR 14, 1
      K = 0
      FOR I = 0 TO N - 1
         NF = I + 1: K = K + 1
         IF I < 10 THEN
           FM$ = L$: CLB% = 1: CLP% = 10
           IF I = NC0 THEN I = I + 1: NF = I + 1
         ELSE
           FM$ = S$: K = I: CLB% = 2: CLP% = 11
         END IF
         X = XO(I) - X0: Y = YO(I) - Y0: Z = ZO(I) - Z0
         IF FLDEC = 0 THEN XP = X / AU: YP = Y / AU: ZP = Z / AU ELSE XP = X: YP = Y: ZP = Z
         IF FLWRITE = 0 THEN
            LOCATE K + 5, 1
            COLOR CLP%, CLB%: PRINT LEFT$(PLN(I), 6); TAB(9); : COLOR 14
            PRINT USING FM$; XP; YP; ZP;
         ELSE
            PRINT #LOGF, USING "##  "; I;
            PRINT #LOGF, PLN(I); "  ";
            PRINT #LOGF, USING FM$; XP; YP; ZP
            IF I >= I0 THEN
               IF FL1ST THEN
                  Org$ = Origin$: IF I = N THEN Org$ = "Geocentric "
                  PRINT #NF, Org$ + CoordType$ + Frame$ + COORD$ + SPACE$(10)
                  PRINT #NF, Nut$ + Equ$ + Ref$ + Aberr$
               END IF
               IF FLJD = 0 THEN
                  PRINT #NF, USING "#####/##/##"; A%; M%; D%;
                  PRINT #NF, USING "  ##:##:##  "; H%; MINUT%; SEC%;
               ELSE
                  PRINT #NF, USING "########.######  "; DJW;
               END IF
               PRINT #NF, USING FM$; XP; YP; ZP
            END IF
         END IF
      NEXT I
      NF = N + 1
      IF FLWRITE = 0 THEN
         LOCATE I + 5, 1: CLP% = 10: CLB% = 1
         COLOR CLP%, CLB%: PRINT LEFT$(PLN(I), 6); TAB(9); : COLOR 14
         PRINT USING LAU$; XO(I) - XO(3); YO(I) - YO(3); ZO(I) - ZO(3);
         COLOR 15, 1: PRINT " Gm";
         COLOR 14, 1
      ELSE
         PRINT #LOGF, ""
         PRINT #LOGF, "Geocentric Coordinates"; TAB(24); "(Gm)"; TAB(40); "(Gm)"; TAB(56); "(Gm)"
         PRINT #LOGF, USING "##  "; I; : PRINT #LOGF, PLN(I); "  ";
         PRINT #LOGF, USING LAU$; XO(I) - XO(3); YO(I) - YO(3); ZO(I) - ZO(3)
         PRINT #LOGF, ""
         IF I >= I0 THEN
            IF FL1ST THEN
               PRINT #NF, Origin$ + CoordType$ + Frame$ + COORD$ + SPACE$(10)
               PRINT #NF, Nut$ + Equ$ + Ref$ + Aberr$
            END IF
            IF FLJD = 0 THEN
               PRINT #NF, USING "#####/##/##"; A%; M%; D%;
               PRINT #NF, USING "  ##:##:##  "; H%; MINUT%; SEC%;
            ELSE
               PRINT #NF, USING "########.######  "; DJW;
            END IF
            PRINT #NF, USING LAU$; XO(I) - XO(3); YO(I) - YO(3); ZO(I) - ZO(3)
         END IF
      END IF

'      LOCATE I + 10, 24
   ELSE
      CALL PrintAngOut(FLGEO, FL1ST, N, T, ABSM(), PLN(), XO(), YO(), ZO(), R(), A%, M%, D%, H%, MINUT%, SEC%)
   END IF
   FL1ST = 0
END SUB

SUB RadToDMS (X, SG$, DG%, HMIN, SEC, NPREC, FL360)
      '.......... NPREC MUST be an integer power of 10
      SHARED RD
      XS = X * RD
      XS = X * RD
      IF FL360 THEN SG$ = "": G1 = 0: G2 = 360: GA = 360 ELSE G1 = -90: G2 = 90: GA = 180
      WHILE XS < G1: XS = XS + GA: WEND
      WHILE XS > G2: XS = XS - GA: WEND
      IF FL360 = 0 THEN SS% = SGN(XS): IF SS% >= 0 THEN SG$ = "+" ELSE SG$ = "-"
      DG% = ABS(FIX(XS))
      XMN = (ABS(XS) - DG%) * 60: HMIN = INT(XMN): XSEC = (XMN - HMIN) * 60
      HSEC = INT(XSEC): HCENT = INT((XSEC - HSEC) * NPREC + .5)
      IF HCENT = NPREC THEN HCENT = 0: HSEC = HSEC + 1
      IF HSEC = 60 THEN HSEC = 0: HMIN = HMIN + 1
      IF HMIN = 60 THEN HMIN = 0: DG% = DG% + 1
      SEC = HSEC + CDBL(HCENT / NPREC)
END SUB

SUB RadToHMS (X, HR, HMIN, SEC, NPREC)
      '...... NPREC MUST be an integer power of 10
      SHARED PI
      XS = X * 12 / PI
      WHILE XS < 0: XS = XS + 24: WEND
      WHILE XS > 24: XS = XS - 24: WEND
      HR = INT(XS): XMN = (XS - HR) * 60: HMIN = INT(XMN)
      XSEC = (XMN - HMIN) * 60: HSEC = INT(XSEC)
      HMILL = INT((XSEC - HSEC) * NPREC + .5)
      IF HMILL = NPREC THEN HMILL = 0: HSEC = HSEC + 1
      IF HSEC = 60 THEN HSEC = 0: HMIN = HMIN + 1
      IF HMIN = 60 THEN HMIN = 0: HR = HR + 1
      SEC = HSEC + CDBL(HMILL) / NPREC
END SUB

SUB ReadFile (NF$, TJ, NEX, X(), Y(), Z(), VX(), VY(), VZ(), Rj0(), Rj00(), PLN() AS STRING * 8, M0(), ABSM())
   SHARED CV

   OPEN NF$ FOR INPUT AS #1
   INPUT #1, DJ: TJ = DJ - 2451545 - .5
   I = 13: NEX = 0
   DO WHILE NOT EOF(1)
     I = I + 1
     INPUT #1, M0(I), ABSM(I), Rj, PLN(I): IF EOF(1) THEN EXIT DO
     INPUT #1, X(I), Y(I), Z(I)
     INPUT #1, VX(I), VY(I), VZ(I)
     Rj0(I) = Rj: Rj00(I) = Rj
     IF Rj = 0 THEN Rj0(I) = RelPar(X(I), Y(I), Z(I), VZ(I) * CV, VY(I) * CV, VZ(I) * CV)
     NEX = NEX + 1: IF NEX > 4 THEN EXIT DO
   LOOP
   CLOSE #1
END SUB

FUNCTION RelPar (X, Y, Z, VX, VY, VZ)
   SHARED VL, GRAV0
   M0 = 1989095.324#
   MU = M0 * GRAV0
   R = SQR(X * X + Y * Y + Z * Z)
   V2 = VX * VX + VY * VY + VZ * VZ
   A = 1 / (2 / R - V2 / MU)
   E0 = V2 / 2 - MU / R
'   RDV = X * VX + Y * VY + Z * VZ
'   VR = RDV / R  '... Radial speed
'   D = R * VR / SQR(MU)
'   E2 = (1 - R / A) ^ 2 + D * D / A
    RPAR = 18 * E0 / (VL * VL)
    IF E0 > 0 THEN RPAR = 0
    RelPar = RPAR

END FUNCTION

' ...... Rotates coordinates to account for precession from J2000 to date
' ...... Five successive rotations are performed, according to
' ...... Bretagnon & Francou, (1988) Astron. Astrophys, 202, 309-315.
'
SUB RotAx (N, X(), Y(), Z(), FLECL)
    SHARED SEQ0, CEQ0, SEQ, CEQ, SEI, CEI, SEN, CEN, SPA, CPA
       FOR I = 0 TO N
          Y = Y(I)
          Y(I) = Y * CEQ0 + Z(I) * SEQ0: Z(I) = Z(I) * CEQ0 - Y * SEQ0
          X = X(I)
          X(I) = X * CEN + Y(I) * SEN: Y(I) = Y(I) * CEN - X * SEN
          Y = Y(I)
          Y(I) = Y * CEI + Z(I) * SEI: Z(I) = Z(I) * CEI - Y * SEI
          X = X(I)
          X(I) = X * CPA - Y(I) * SPA: Y(I) = Y(I) * CPA + X * SPA
          IF FLECL = 0 THEN
             Y = Y(I)
             Y(I) = Y * CEQ - Z(I) * SEQ: Z(I) = Z(I) * CEQ + Y * SEQ
          END IF
       NEXT I
END SUB

'
' ............ Rotates back the coordinates from date to J2000
'
SUB RotInv (N, X(), Y(), Z())
    SHARED SEQ0, CEQ0, SEQ, CEQ, SEI, CEI, SEN, CEN, SPA, CPA
       FOR I = 0 TO N
          Y = Y(I)
          Y(I) = Y * CEQ + Z(I) * SEQ: Z(I) = Z(I) * CEQ - Y * SEQ
          X = X(I)
          X(I) = X * CPA + Y(I) * SPA: Y(I) = Y(I) * CPA - X * SPA
          Y = Y(I)
          Y(I) = Y * CEI - Z(I) * SEI: Z(I) = Z(I) * CEI + Y * SEI
          X = X(I)
          X(I) = X * CEN - Y(I) * SEN: Y(I) = Y(I) * CEN + X * SEN
          Y = Y(I)
          Y(I) = Y * CEQ0 - Z(I) * SEQ0: Z(I) = Z(I) * CEQ0 + Y * SEQ0
       NEXT I

END SUB

'
'... This subroutine scans the Library file PLANS.BIN, storing the Julian day
'... of each set of starting conditions (measured from J2000) in the array
'... DjLib() and returning the upper and lower date boundary of the Library
'
SUB ScanLib (FLDE, NST, YMIN, YMAX, DjLib() AS LONG)
    IF FLDE = 1 THEN SSL$ = "PLANS.BIN"
    IF FLDE = 2 THEN SSL$ = "PL406B.BIN"
    OPEN SSL$ FOR INPUT AS #3: CLOSE #3
    OPEN SSL$ FOR BINARY AS #4
    GET #4, , NUM: LBLOCK = NUM * 36 + 4
    T = 999999999: P& = 3: J = 0: DMIN = 10000000: DMAX = -10000000
    DO WHILE NOT EOF(4)
       J = J + 1
       GET #4, P&, DJ&: IF EOF(4) THEN EXIT DO
       DjLib(J) = DJ&  ' .... no of Julian days from J2000
       IF DJ& < DMIN THEN DMIN = DJ&
       IF DJ& > DMAX THEN DMAX = DJ&
       P& = P& + LBLOCK
    LOOP
    NST = J - 1
    YMIN = INT(DMIN / 365.25 + 1986): YMAX = INT(DMAX / 365.25 + 2014)
    CLOSE #4

END SUB

SUB SelOption (Outp%)

   CLS
   CALL DrawFrame(14, 1)
   CALL DrawTitle(11, 2, 14)
   COLOR 14: L = 8: C% = 32
   LOCATE L, C%: PRINT "Select Option";
   LOCATE L + 2, C%: COLOR 10: PRINT "1.     DE200";
   LOCATE L + 3, C%: COLOR 13: PRINT "2.     DE406";
   LOCATE L + 4, C%: COLOR 15: PRINT "3.     File ";
   LOCATE L + 5, C%: COLOR 15: PRINT "4.     AstLib";
   COLOR 14
   LOCATE L + 7, C%: PRINT "5.     Shell to DOS";
   LOCATE L + 8, C%: PRINT "6.     Quit";

   M$ = ""
   WHILE M$ < "1" OR M$ > "6" OR LEN(M$) > 1: M$ = INKEY$: WEND
   Outp% = VAL(M$)

END SUB

SUB SetDeltaT (DELT0)
   COLOR 0, 3
   FOR I = 1 TO 5
      LOCATE 19 + I, 63
      PRINT SPACE$(17);
   NEXT
   LOCATE 21, 65
   PRINT "Enter Delta-T ";
   LOCATE 23, 65
   INPUT ; DELT0
   DELT0 = DELT0 / 86400
   COLOR 14, 1
   CLS
END SUB

SUB SetGeoCoord (GLAT, GLON, ZH)
   SHARED RD
   COLOR 0, 3
   FOR I = 1 TO 7
      LOCATE 17 + I, 58
      PRINT SPACE$(22);
   NEXT
   SS% = SGN(GLON): IF SS% > 0 THEN EW$ = " (E)" ELSE EW$ = " (W)"
   LOCATE 19, 59
   PRINT USING "Lat = +###.###" + CHR$(248); GLAT * RD;
   LOCATE 20, 59
   PRINT USING "Lon = +###.###" + CHR$(248); GLON * RD;
   IF SS% THEN PRINT EW$;  ELSE PRINT "    ";
   LOCATE 21, 59
   PRINT USING " h = ###### m"; ZH * 1000000000#;
   LOCATE 23, 59: PRINT " Accept ? (Y or N)"
   V$ = ""
   WHILE V$ <> "Y" AND V$ <> "N": V$ = UCASE$(INKEY$): WEND
   IF V$ = "N" THEN
      FOR I = 1 TO 7
         LOCATE 17 + I, 58
         PRINT SPACE$(22);
      NEXT
      LOCATE 19, 59: PRINT "East Long. positive";
      LOCATE 21, 59: INPUT ; "Lat (deg) "; GLAT: GLAT = GLAT / RD
      LOCATE 22, 59: INPUT ; "Lon (deg) "; GLON: GLON = GLON / RD
      LOCATE 23, 59: INPUT ; " h (m)  "; ZH: ZH = ZH / 1000000000#
   END IF
   COLOR 14, 1
   CLS

END SUB

'..... Select optimal order (no. of stages in STORM8) according to stepsize
SUB SetOrder (DT, NO)

    IF DT = 3 THEN
       NO = 6
    ELSEIF DT > 1.5 THEN
       NO = 7
    ELSEIF DT > 1 THEN
       NO = 6
    ELSEIF DT > .5 THEN
       NO = 5
    ELSEIF DT > .25 THEN
       NO = 4
    ELSEIF DT > .125 THEN
       NO = 3
    ELSEIF DT > .03125 THEN
       NO = 2
    ELSE
       NO = 1
    END IF

END SUB

SUB SitePos (FLTOP, FLGEO, T, DT, GLAT, GLON, ZH, XSIT, YSIT, ZSIT, VXSIT, VYSIT, VZSIT)
   SHARED RD, S%, CV, TSD
   IF FLTOP = 0 OR FLGEO = 0 THEN
      XSIT = 0: YSIT = 0: ZSIT = 0
      VXSIT = 0: VYSIT = 0: VZSIT = 0
      EXIT SUB
   END IF
   E2 = .006694378#: R0 = .006378137#
   SGLAT = SIN(GLAT): CGLAT = COS(GLAT)

   TSD = Tsid(T, GLON)
   STSD = SIN(TSD): CTSD = COS(TSD)
   SQE = SQR(1 - E2 * SGLAT * SGLAT)
   XSIT = CGLAT * CTSD * (ZH + R0 / SQE)
   YSIT = CGLAT * STSD * (ZH + R0 / SQE)
   ZSIT = SGLAT * (ZH + R0 * (1 - E2) / SQE)
   V0 = .4638312# * COS(GLAT) * S% * CV * DT  '... Earth rotation speed
   '... Sign of components of rotational speed
   IF XSIT >= 0 THEN
      IF YSIT > 0 THEN PX% = -1: PY% = 1 ELSE PX% = 1: PY% = 1
   ELSE
      IF YSIT > 0 THEN PX% = -1: PY% = -1 ELSE PX% = 1: PY% = -1
   END IF
   VX = SQR(V0 * V0 / (1 + XSIT * XSIT / (YSIT * YSIT)))
   IF VX > V0 THEN VX = V0
   VY = SQR(V0 * V0 - VX * VX)
   ' .... Components of rotational speed
   VXSIT = PX% * VX: VYSIT = PY% * VY: VZSIT = 0

END SUB

'
'. Integration partly based on Extrapol. method for second order conservative
'. conservative equations (Numerical Recipes 2 Ed Cambridge Univ. Press, 1992)
'
SUB STORM8 (N, X(), Y(), Z(), VX(), VY(), VZ()) STATIC
     SHARED AK(), M(), M0(), GRAV, MS%, NO, QK(), CV, DT
     DIM XN(20), YN(20), ZN(20), DELTX(20), DELTY(20), DELTZ(20)
     DIM AX0(20), AY0(20), AZ0(20), AX(20), AY(20), AZ(20)
     DIM X1(7, 20), Y1(7, 20), Z1(7, 20), VX1(7, 20), VY1(7, 20), VZ1(7, 20)
     GR0 = GRAV
     CALL CAcc(N, X(), Y(), Z(), AX0(), AY0(), AZ0(), M()) '... Accelerations
     '......................... Computes NO+1 trial values of the variables
     FOR L = 0 TO NO '.......... NO+1 is the number of stages of the method
         ev = QK(L) / MS%: NSTEP = QK(L) - 1: ev2 = 2 * ev
         GRAV = GR0 / (ev * ev) '... Time unit is scaled so that always h=1
         FOR I = 0 TO N
            M(I) = M0(I) * GRAV
            DELTX(I) = (VX(I) + AX0(I) / ev2) / ev: DELTY(I) = (VY(I) + AY0(I) / ev2) / ev: DELTZ(I) = (VZ(I) + AZ0(I) / ev2) / ev
            XN(I) = X(I) + DELTX(I): YN(I) = Y(I) + DELTY(I): ZN(I) = Z(I) + DELTZ(I)
         NEXT I
         FOR K = 1 TO NSTEP
            CALL CAcc(N, XN(), YN(), ZN(), AX(), AY(), AZ(), M())
            FOR I = 0 TO N
               DELTX(I) = DELTX(I) + AX(I): DELTY(I) = DELTY(I) + AY(I): DELTZ(I) = DELTZ(I) + AZ(I)
               XN(I) = XN(I) + DELTX(I): YN(I) = YN(I) + DELTY(I): ZN(I) = ZN(I) + DELTZ(I)
            NEXT I
         NEXT K
         CALL CAcc(N, XN(), YN(), ZN(), AX(), AY(), AZ(), M())
         FOR I = 0 TO N '....... Stores current trial values of the variables
            X1(L, I) = XN(I): Y1(L, I) = YN(I): Z1(L, I) = ZN(I)
            VX1(L, I) = (DELTX(I) + AX(I) / 2) * ev: VY1(L, I) = (DELTY(I) + AY(I) / 2) * ev: VZ1(L, I) = (DELTZ(I) + AZ(I) / 2) * ev
         NEXT I
     NEXT L
     FOR I = 0 TO N '....... Performs the Extrapolation to h=0 ...............
        X(I) = (AK(NO, 8) * X1(7, I) + AK(NO, 7) * X1(6, I) + AK(NO, 6) * X1(5, I) + AK(NO, 5) * X1(4, I) + AK(NO, 4) * X1(3, I) + AK(NO, 3) * X1(2, I) + AK(NO, 2) * X1(1, I) + AK(NO, 1) * X1(0, I)) / AK(NO, 0)
        Y(I) = (AK(NO, 8) * Y1(7, I) + AK(NO, 7) * Y1(6, I) + AK(NO, 6) * Y1(5, I) + AK(NO, 5) * Y1(4, I) + AK(NO, 4) * Y1(3, I) + AK(NO, 3) * Y1(2, I) + AK(NO, 2) * Y1(1, I) + AK(NO, 1) * Y1(0, I)) / AK(NO, 0)
        Z(I) = (AK(NO, 8) * Z1(7, I) + AK(NO, 7) * Z1(6, I) + AK(NO, 6) * Z1(5, I) + AK(NO, 5) * Z1(4, I) + AK(NO, 4) * Z1(3, I) + AK(NO, 3) * Z1(2, I) + AK(NO, 2) * Z1(1, I) + AK(NO, 1) * Z1(0, I)) / AK(NO, 0)
        VX(I) = (AK(NO, 8) * VX1(7, I) + AK(NO, 7) * VX1(6, I) + AK(NO, 6) * VX1(5, I) + AK(NO, 5) * VX1(4, I) + AK(NO, 4) * VX1(3, I) + AK(NO, 3) * VX1(2, I) + AK(NO, 2) * VX1(1, I) + AK(NO, 1) * VX1(0, I)) / AK(NO, 0)
        VY(I) = (AK(NO, 8) * VY1(7, I) + AK(NO, 7) * VY1(6, I) + AK(NO, 6) * VY1(5, I) + AK(NO, 5) * VY1(4, I) + AK(NO, 4) * VY1(3, I) + AK(NO, 3) * VY1(2, I) + AK(NO, 2) * VY1(1, I) + AK(NO, 1) * VY1(0, I)) / AK(NO, 0)
        VZ(I) = (AK(NO, 8) * VZ1(7, I) + AK(NO, 7) * VZ1(6, I) + AK(NO, 6) * VZ1(5, I) + AK(NO, 5) * VZ1(4, I) + AK(NO, 4) * VZ1(3, I) + AK(NO, 3) * VZ1(2, I) + AK(NO, 2) * VZ1(1, I) + AK(NO, 1) * VZ1(0, I)) / AK(NO, 0)
     NEXT
     GRAV = GR0 '.......... Restores MG values for time unit = stepsize H
     FOR I = 0 TO N: M(I) = M0(I) * GRAV: NEXT
END SUB

FUNCTION Tsid (TU, GLON)
   SHARED RD
   '... Local Sidereal Time
   Tsid = (280.46061837# + TU * (360.98564736629# + TU * (2.907879D-13 - TU * 5.3016D-22))) / RD + GLON
END FUNCTION

SUB TToDate (TP#, Y%, MO%, D%, HR%, MINUT%, SEC%)
     SHARED FLHOR, GLON, PI
     T# = TP# + 2451545#
     IF FLHOR THEN TZ = (INT(GLON * 12 / PI + .5)) / 24 ELSE TZ = 0
     GGI = INT(T# + TZ): GGF = T# - GGI + TZ
     IF GGI > 2299160 THEN AGG = INT((GGI - 1867216.25#) / 36524.25#): BGG = GGI + 1 + AGG - INT(AGG / 4) ELSE BGG = GGI
     CGG = BGG + 1524: DGG = INT((CGG - 122.1) / 365.25): EGG = INT(365.25 * DGG): GGH = INT((CGG - EGG) / 30.6001#): GG = CGG - EGG + GGF - INT(30.6001 * GGH)
     IF GGH < 14 THEN MO% = GGH - 1 ELSE MO% = GGH - 13
     IF MO% > 2 THEN Y% = DGG - 4716 ELSE Y% = DGG - 4715
     D% = INT(GG): TH = (GG - D%) * 24: HR% = INT(TH): TMIN = (TH - HR%) * 60: MINUT% = INT(TMIN): SEC% = (TMIN - MINUT%) * 60
     IF SEC% = 60 THEN SEC% = 0: MINUT% = MINUT% + 1: IF MINUT% = 60 THEN MINUT% = 0: HR% = HR% + 1: IF HR% = 24 THEN HR% = 0: D% = D% + 1

END SUB

