home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / cpm / basic / astrmenu.lbr / KEPLER.BZS / KEPLER.BAS
BASIC Source File  |  1987-04-26  |  2KB  |  91 lines
  1. 1 PRINT "From the August 1985 SKY & TELESCOPE, pp. 158-9."
  2. 2 PRINT
  3. 3 PRINT "This program solves Kepler's equation in orbit computation to compute"
  4. 4 PRINT "an ephemeris for a comet or asteroid from its orbital elements."
  5. 5 PRINT
  6. 6 PRINT "INPUT:  Perhelion distance (AU), eccentricity, days from perihelion."
  7. 7 PRINT
  8. 8 PRINT "OUTPUT:  True anomaly (degrees), distance from the sun (AU)." 
  9. 9 PRINT
  10. 10 REM     KEPLER'S EQUATION
  11. 11 REM
  12. 12 P1=3.14159265#: R1=180/P1
  13. 13 K=.01720209895#
  14. 14 INPUT "PERI DISTANCE Q ";Q
  15. 15 INPUT "ECCENTRICITY    ";E0
  16. 16 INPUT "DAYS FROM PERI  ";T
  17. 17 PRINT
  18. 18 IF E0>=.95 THEN 40
  19. 19 IF E0>=1 THEN 85
  20. 20 A1=Q/(1-E0): M=K*T*A1^(-1.5)
  21. 21 REM
  22. 22 REM     BINARY SEARCH
  23. 23 REM
  24. 24 F=SGN(M): M=ABS(M)/(2*P1)
  25. 25 M=(M-INT(M))*2*P1*F
  26. 26 IF M<0 THEN M=M+2*P1
  27. 27 F=1: IF M>P1 THEN F=-1
  28. 28 IF M>P1 THEN M=2*P1-M
  29. 29 E=P1/2: D=P1/4
  30. 30 FOR I1=1 TO 23
  31. 31 M1=E-E0*SIN(E)
  32. 32 E=E+SGN(M-M1)*D: D=D/2
  33. 33 NEXT I1
  34. 34 V=SQR((1+E0)/(1-E0)): E=E*F
  35. 35 V=2*ATN(V*SIN(E/2)/COS(E/2))
  36. 36 R=A1*(1-E0*COS(E))
  37. 37 PRINT "BINARY SEARCH -- "
  38. 38 GOTO 81
  39. 39 REM
  40. 40 REM      GAUSS METHOD
  41. 41 REM
  42. 42 GOSUB 69
  43. 43 A=SQR((1+9*E0)/10)
  44. 44 B=5*(1-E0)/(1+9*E0)
  45. 45 C=SQR(5*(1+E0)/(1+9*E0))
  46. 46 B1=3*A*K*T/SQR(2*Q*Q*Q)
  47. 47 B2=1: REM  INITIAL VALUE
  48. 48 W1=B2*B1: B3=ATN(2/W1)
  49. 49 T1=SIN(B3/2)/COS(B3/2)
  50. 50 S1=SGN(T1): T1=ABS(T1)
  51. 51 T2=T1^(1/3)*S1: G=ATN(T2)
  52. 52 S=2*COS(2*G)/SIN(2*G)
  53. 53 A2=B*S*S: B0=B2: B2=0
  54. 54 IF ABS(A2)>.3 THEN 19
  55. 55 FOR J=0 TO 7
  56. 56 B2=B2+B(J)*A2^J
  57. 57 NEXT J
  58. 58 IF ABS(B2-B0)>1E-08 THEN 48
  59. 59 C1=0
  60. 60 FOR J=0 TO 7
  61. 61 C1=C1+S(J)*A2^J
  62. 62 NEXT J
  63. 63 C1=SQR(1/C1)
  64. 64 V1=C*C1*S: D1=1/(1+A2*C1*C1)
  65. 65 V=2*ATN(V1): R=Q*D1*(1+V1*V1)
  66. 66 PRINT "GAUSS METHOD --"
  67. 67 GOTO 81
  68. 68 REM
  69. 69 REM   COEFFICIENTS
  70. 70 FOR J=0 TO 7: READ B(J): NEXT
  71. 71 FOR J=0 TO 7: READ S(J): NEXT
  72. 72 RETURN
  73. 73 DATA 1, 0, -0.017142857
  74. 74 DATA -0.003809524,-0.001104267
  75. 75 DATA -0.000367358,-0.000131675
  76. 76 DATA -0.000049577, 1, -0.8
  77. 77 DATA 0.04571429, 0.01523810
  78. 78 DATA 0.00562820, 0.00218783
  79. 79 DATA 0.00087905, 0.00036155
  80. 80 REM
  81. 81 IF V<0 THEN V=V+2*P1
  82. 82 PRINT "TRUE ANOMALY:  ";V*R1
  83. 83 PRINT "DISTANCE (AU): ";R
  84. 84 GOTO 86
  85. 85 PRINT "OUT OF RANGE"
  86. 86 RUN"ASTRMENU.BAS"
  87. *P1
  88. 82 PRINT "TRUE ANOMALY:  ";V*R1
  89. 83 PRINT "DISTANCE (AU): ";R
  90. 84 GOTO 86
  91. 85 PRINT "OUT OF RANG