The CIA World Factbook 1992/k3bimage.iso/sel/04/0086/kepler.bas

BASIC Source File | 1991-12-02 | 5.2 KB | 211 lines

' PROGRAM "KEPLER"

' NUMERICAL SOLUTION OF KEPLER'S EQUATION FOR
' ELLIPTIC, PARABOLIC AND HYPERBOLIC ORBITS

' PUBLIC DOMAIN

' IBM-PC  << QUICKBASIC COMPILER VERSION 3.0 >>

DEFDBL A-Z

CONST PI=3.14159265D0
CONST RTD=180/PI
CONST DTR=PI/180

COMMON SHARED ECC,MANOM.RAD,EANOM.RAD

'**********************************************************

CLS
PRINT
PRINT
PRINT "Program KEPLER"
PRINT
PRINT "Microsoft QuickBASIC Compiler"
PRINT "(C) Copyright Microsoft Corp. 1982-1987"
PRINT
CALL KEYCHECK

CLS
PRINT
PRINT
PRINT "Program introduction ( y = yes, n = no )"
INPUT A$
IF INSTR("yY",A$) THEN CALL INTRODUCTION

DO
   CLS
   PRINT
   PRINT
   PRINT TAB(15); "Kepler's Equation"
   PRINT
   PRINT
   PRINT "Mean anomaly ( degrees [ 0 - 360 ] )"
   INPUT MANOM.DEG
   MANOM.RAD=DTR*MANOM.DEG
   PRINT
   PRINT "Orbital eccentricity ( non-dimensional )"
   INPUT ECC

   IF ECC<>0D0 THEN
      CALL KEPLER
   ELSE
      EANOM.RAD=MANOM.RAD
   END IF

   ' PRINT RESULTS

   FORMAT$="####.######"

   CLS
   PRINT
   PRINT
   PRINT TAB(22);"Kepler's Equation"
   PRINT
   PRINT
   PRINT TAB(5);"Eccentricity       ( non-dimensional )";TAB(45);
   PRINT USING FORMAT$;ECC
   PRINT
   PRINT TAB(5);"Mean anomaly       ( degrees )";TAB(45);
   PRINT USING FORMAT$;MANOM.DEG
   PRINT
   IF ECC<1 THEN PRINT TAB(5);"Eccentric anomaly  ( degrees )";TAB(45);
   IF ECC=1 THEN PRINT TAB(5);"Parabolic anomaly  ( degrees )";TAB(45);
   IF ECC>1 THEN PRINT TAB(5);"Hyperbolic anomaly ( degrees )";TAB(45);
   PRINT USING FORMAT$;EANOM.RAD*RTD
   CALL KEYCHECK

   CLS
   PRINT
   PRINT
   PRINT "Another selection ( y = yes, n = no )"
   INPUT SELECTION$
LOOP UNTIL INSTR("nN",SELECTION$)

END

'**********************************************************

SUB KEPLER STATIC

    ' KEPLER'S EQUATION SUBROUTINE

    EANOM.RAD=LOG(2*MANOM.RAD/ECC+1.8D0)
    IF ECC<1D0 THEN CALL ESTIMATE

    DO
       IF ECC<1D0 THEN L=ECC*SIN(EANOM.RAD)
       IF ECC>=1D0 THEN L=.5D0*ECC*(EXP(EANOM.RAD)-EXP(-EANOM.RAD))

       IF ECC<1D0 THEN M=EANOM.RAD-L-MANOM.RAD
       IF ECC>=1D0 THEN M=L-EANOM.RAD-MANOM.RAD

       IF ABS(M)<=1D-8 THEN EXIT DO

       IF ECC<1D0 THEN P=ECC*COS(EANOM.RAD)
       IF ECC>=1D0 THEN P=.5D0*ECC*(EXP(EANOM.RAD)+EXP(-EANOM.RAD))
       IF ECC<1D0 THEN Q=1D0-P
       IF ECC>=1D0 THEN Q=P-1D0

       R=-M/Q
       R=-M/(Q+.5D0*R*L)
       R=-M/(Q+.5D0*R*L+R*R*P/6D0)
       R=-M/(Q+.5D0*R*L+R*R*P/6D0-R*R*R*L/24D0)
       EANOM.RAD=EANOM.RAD+R
    LOOP

END SUB

'**********************************************************

SUB ESTIMATE STATIC

    ' EMPIRICAL ESTIMATE SUBROUTINE

    A=ECC^2
    B=MANOM.RAD^2

    EANOM.RAD=A*(-.0579437D0*B+.02665324D0*MANOM.RAD+.05868658D0)
    EANOM.RAD=EANOM.RAD+ECC*(-.19174923D0*B+.46905672D0*MANOM.RAD+.600725329D0)
    EANOM.RAD=EANOM.RAD-.04987627D0*B+1.17552263D0*MANOM.RAD-.12324871D0

END SUB

'**********************************************************

SUB INTRODUCTION STATIC

    ' INTRODUCTION SUBROUTINE

    CLS
    PRINT
    PRINT TAB(10);"PROGRAM 'KEPLER' IS AN INTERACTIVE BASIC"
    PRINT TAB(10);"PROGRAM WHICH CAN BE USED TO SOLVE"
    PRINT TAB(10);"KEPLER'S EQUATION. THE ORBITAL MOTION"
    PRINT TAB(10);"OF SATELLITES, COMETS AND ALL CELESTIAL"
    PRINT TAB(10);"BODIES CAN BE DETERMINED FROM KEPLER'S"
    PRINT TAB(10);"EQUATION. THIS EQUATION IS WRITTEN AS"
    PRINT
    PRINT TAB(10);"       M = E - e SIN E "
    PRINT
    PRINT TAB(10);"FOR ELLIPTIC ORBITS OR"
    PRINT
    PRINT TAB(10);"       M = e SINH F - F"
    PRINT
    PRINT TAB(10);"FOR PARABOLIC OR HYPERBOLIC ORBITS."
    PRINT
    PRINT TAB(10);"KEPLER'S EQUATION IS A TRANSCENDENTAL"
    PRINT TAB(10);"EQUATION WHICH IS USUALLY SOLVED WITH"
    PRINT TAB(10);"AN ITERATIVE NUMERICAL METHOD."
    CALL KEYCHECK

    CLS
    PRINT
    PRINT TAB(10);"   IN BOTH OF THESE EQUATIONS, 'M'"
    PRINT TAB(10);"IS CALLED THE 'MEAN ANOMALY' AND 'E'"
    PRINT TAB(10);"IS CALLED THE 'ORBITAL ECCENTRICITY'."
    PRINT TAB(10);"IN THE FIRST EQUATION, 'E' IS CALLED"
    PRINT TAB(10);"THE 'ECCENTRIC ANOMALY' WHILE IN THE"
    PRINT TAB(10);"SECOND EQUATION 'F' IS CALLED THE"
    PRINT TAB(10);"'PARABOLIC' OR 'HYPERBOLIC' ANOMALY."
    PRINT
    PRINT TAB(10);"THE TRIG FUNCTION 'SINH' IN THE SECOND"
    PRINT TAB(10);"EQUATION IS THE HYPERBOLIC SINE."
    PRINT
    PRINT TAB(10);"   THE TYPE OF ORBIT IS DETERMINED"
    PRINT TAB(10);"BY THE ECCENTRICITY."
    PRINT
    PRINT TAB(10);"       0 < e < 1   ELLIPTIC"
    PRINT TAB(10);"           e = 1   PARABOLIC"
    PRINT TAB(10);"           e > 1   HYPERBOLIC"
    CALL KEYCHECK

    CLS
    PRINT
    PRINT TAB(10);"THE ACTUAL POSITION OF A BODY IN ITS"
    PRINT TAB(10);"ORBIT IS CALLED 'TRUE ANOMALY' AND"
    PRINT TAB(10);"CAN BE DETERMINED FROM ECCENTRICITY"
    PRINT TAB(10);"AND THE ECCENTRIC OR PARABOLIC OR"
    PRINT TAB(10);"HYPERBOLIC ANOMALY."
    CALL KEYCHECK

END SUB

'**********************************************************

SUB KEYCHECK STATIC

    ' CHECK USER RESPONSE SUBROUTINE

    PRINT
    PRINT
    PRINT TAB(15);"< press any key to continue >"

    A$=""
    WHILE A$=""
      A$=INKEY$
    WEND

END SUB