home *** CD-ROM | disk | FTP | other *** search
/ The CIA World Factbook 1992 / k3bimage.iso / sel / 04 / 0086 / j2000.bas < prev    next >
Encoding:
BASIC Source File  |  1991-12-02  |  9.4 KB  |  354 lines
  1. ' Program  "J2000"
  2.  
  3. ' copyright (C) 1986 by David Eagle
  4.  
  5. ' public domain on November 14, 1986
  6.  
  7. ' IBM-PC  << QuickBASIC Compiler Version 3.0 >>
  8.  
  9. ' conversion of stellar positions, proper motions, parallax and
  10. ' radial velocity from the standard epoch B1950 (FK4) to J2000 (FK5)
  11. ' reference: Astronomical Almanac, 1985, page X
  12.  
  13. '***************************************************************
  14.  
  15. OPTION BASE 1
  16. DEFDBL A-Z
  17.  
  18. DIM A1(3),A2(3),UPV(3),UVV(3),R1(3),R2(6),R3(6),V2(3),M(6,6)
  19.  
  20. COMMON SHARED RAS.HR.1950,RAS.MIN.1950,RAS.SEC.1950
  21. COMMON SHARED DEC.DEG.1950,DEC.MIN.1950,DEC.SEC.1950
  22. COMMON SHARED RAS.PMOTION.1950,DEC.PMOTION.1950
  23. COMMON SHARED PARALLAX.1950,RADIALV.1950
  24.  
  25. COMMON SHARED RAS.HR.2000,RAS.MIN.2000,RAS.SEC.2000
  26. COMMON SHARED DEC.DEG.2000,DEC.MIN.2000,DEC.SEC.2000
  27. COMMON SHARED RAS.PMOTION.2000,DEC.PMOTION.2000
  28. COMMON SHARED PARALLAX.2000,RADIALV.2000
  29.  
  30. CONST PI=3.141592653589793D0
  31. CONST PIDIV2=.5D0*PI
  32. CONST DTR=PI/180D0
  33. CONST RTD=180D0/PI
  34. CONST XE=PI/12D0
  35. CONST XER=12D0/PI
  36.  
  37. A1(1)=-.00000162557D0
  38. A1(2)=-.00000031919D0
  39. A1(3)=-.00000013843D0
  40. A2(1)=.001244D0
  41. A2(2)=-.001579D0
  42. A2(3)=-.00066D0
  43.  
  44. ' define elements of transformation matrix
  45.  
  46. M(1,1)=.9999256782D0
  47. M(1,2)=-.0111820611D0
  48. M(1,3)=-.0048579477D0
  49. M(1,4)=.00000242395018D0
  50. M(1,5)=-.00000002710663D0
  51. M(1,6)=-.00000001177656D0
  52. M(2,1)=.011182061D0
  53. M(2,2)=.9999374784D0
  54. M(2,3)=-.0000271765D0
  55. M(2,4)=.00000002710663D0
  56. M(2,5)=.00000242397878D0
  57. M(2,6)=-.00000000006587D0
  58. M(3,1)=.0048579479D0
  59. M(3,2)=-.0000271474D0
  60. M(3,3)=.9999881997D0
  61. M(3,4)=.00000001177656D0
  62. M(3,5)=-.00000000006582D0
  63. M(3,6)=.00000242410173D0
  64. M(4,1)=-.000551D0
  65. M(4,2)=-.238565D0
  66. M(4,3)=.435739D0
  67. M(4,4)=.99994704D0
  68. M(4,5)=-.01118251D0
  69. M(4,6)=-.00485767D0
  70. M(5,1)=.238514D0
  71. M(5,2)=-.002667D0
  72. M(5,3)=-.008541D0
  73. M(5,4)=.01118251D0
  74. M(5,5)=.99995883D0
  75. M(5,6)=-.00002718D0
  76. M(6,1)=-.435623D0
  77. M(6,2)=.012254D0
  78. M(6,3)=.002117D0
  79. M(6,4)=.00485767D0
  80. M(6,5)=-.00002714D0
  81. M(6,6)=1.00000956D0
  82.  
  83. ' initialization
  84.  
  85. CLS
  86. PRINT
  87. PRINT "Program J2000"
  88. PRINT "(C) Copyright 1986 by David Eagle"
  89. PRINT
  90. PRINT "Microsoft QuickBASIC Compiler"
  91. PRINT "(C) Copyright Microsoft Corp. 1982-1987"
  92. PRINT
  93. CALL KEYCHECK
  94.  
  95. CLS
  96. PRINT
  97. PRINT
  98. PRINT "Program introduction ( y = yes, n = no )"
  99. INPUT INTRO$
  100. IF INSTR("yY",INTRO$) THEN CALL INTRODUCTION
  101.  
  102. ' main program loop
  103.  
  104. DO
  105. CLS
  106. PRINT
  107. PRINT
  108. PRINT "Right ascension with respect to 1950"
  109. PRINT "( hours [ 0 - 23 ], minutes [ 0 - 59 ], seconds [ 0 - 59 ] )"
  110. INPUT RAS.HR.1950,RAS.MIN.1950,RAS.SEC.1950
  111. PRINT
  112. PRINT "Declination with respect to 1950"
  113. PRINT "( degrees [ -90 to +90 ], minutes [ 0 - 59 ], seconds [ 0 - 59 ] )"
  114. INPUT DEC.DEG.1950,DEC.MIN.1950,DEC.SEC.1950
  115. PRINT
  116. PRINT "Right ascension proper motion with respect to 1950"
  117. PRINT "( arc seconds per tropical century )"
  118. INPUT RAS.PMOTION.1950
  119. PRINT
  120. PRINT "Declination proper motion with respect to 1950"
  121. PRINT "( arc seconds per tropical century )"
  122. INPUT DEC.PMOTION.1950
  123. PRINT
  124. PRINT "Parallax with respect to 1950 ( arc seconds )"
  125. INPUT PARALLAX.1950
  126. PRINT
  127. PRINT "Radial velocity with respect to 1950 ( kilometers per second )"
  128. INPUT RADIALV.1950
  129.  
  130. ' convert 1950 right ascension and declination to radians
  131.  
  132. RAS.1950=XE*(RAS.HR.1950+RAS.MIN.1950/60D0+RAS.SEC.1950/3600D0)
  133. DEC.1950=DTR*SGN(DEC.DEG.1950)*(ABS(DEC.DEG.1950)+(DEC.MIN.1950*60D0+DEC.SEC.1950)/3600D0)
  134.  
  135. ' calculate unit position vector
  136.  
  137. UPV(1)=COS(RAS.1950)*COS(DEC.1950)
  138. UPV(2)=SIN(RAS.1950)*COS(DEC.1950)
  139. UPV(3)=SIN(DEC.1950)
  140.  
  141. K1=21.095D0*RADIALV.1950*PARALLAX.1950
  142.  
  143. ' calculate unit velocity vector
  144.  
  145. UVV(1)=-RAS.PMOTION.1950*SIN(RAS.1950)*COS(DEC.1950)-DEC.PMOTION.1950*COS(RAS.1950)*SIN(DEC.1950)+K1*UPV(1)
  146. UVV(2)=RAS.PMOTION.1950*COS(RAS.1950)*COS(DEC.1950)-DEC.PMOTION.1950*SIN(RAS.1950)*SIN(DEC.1950)+K1*UPV(2)
  147. UVV(3)=DEC.PMOTION.1950*COS(DEC.1950)+K1*UPV(3)
  148.  
  149. ' E-term constants
  150.  
  151. K2=A1(1)*UPV(1)+A1(2)*UPV(2)+A1(3)*UPV(3)
  152. K3=A2(1)*UPV(1)+A2(2)*UPV(2)+A2(3)*UPV(3)
  153.  
  154. FOR I%=1 TO 3
  155.     R1(I%)=UPV(I%)-A1(I%)+K2*UPV(I%)
  156.     R2(I%)=R1(I%)
  157.     V2(I%)=UVV(I%)-A2(I%)+K3*UPV(I%)
  158.     R2(I%+3)=V2(I%)
  159. NEXT I%
  160.  
  161. ' matrix multiplication
  162.  
  163. FOR I%=1 TO 6
  164.     A=0D0
  165.     FOR J%=1 TO 6
  166.         A=A+M(I%,J%)*R2(J%)
  167.     NEXT J%
  168.     R3(I%)=A
  169. NEXT I%
  170.  
  171. RM.2000=SQR(R3(1)^2+R3(2)^2+R3(3)^2)
  172.  
  173. ' compute declination wrt J2000
  174.  
  175. CALL ARCSIN(R3(3)/RM.2000,B)
  176. CALL CONVERT(B*RTD,DEC.DEG.2000,DEC.MIN.2000,DEC.SEC.2000)
  177.  
  178. ' compute right ascension wrt J2000
  179.  
  180. CALL ATAN3(R3(2),R3(1),C)
  181. CALL CONVERT(C*XER,RAS.HR.2000,RAS.MIN.2000,RAS.SEC.2000)
  182.  
  183. K4=R3(1)^2+R3(2)^2
  184.  
  185. ' compute proper motions wrt J2000
  186.  
  187. RAS.PMOTION.2000=(R3(1)*R3(5)-R3(2)*R3(4))/K4
  188. DEC.PMOTION.2000=(R3(6)*K4-R3(3)*(R3(1)*R3(4)+R3(2)*R3(5)))/(RM.2000^2*SQR(K4))
  189.  
  190. ' compute parallax and radial velocity wrt J2000
  191.  
  192. IF PARALLAX.1950=0D0 THEN
  193.   RADIALV.2000=RADIALV.1950
  194. ELSE
  195.   RADIALV.2000=(R3(1)*R3(4)+R3(2)*R3(5)+R3(3)*R3(6))/(21.095D0*PARALLAX.1950*RM.2000)
  196. END IF
  197.  
  198. PARALLAX.2000=PARALLAX.1950/RM.2000
  199.  
  200. CALL DISPLAY.DATA
  201.  
  202. CLS
  203. PRINT
  204. PRINT
  205. PRINT "Another selection ( y = yes, n = no )"
  206. INPUT SELECTION$
  207. LOOP UNTIL INSTR("nN",SELECTION$)
  208.  
  209. END
  210.  
  211.  
  212. '***************************************************************
  213.  
  214. SUB CONVERT(DEG,HR,MIN,SEC) STATIC
  215.  
  216.     ' convert degrees to hms or dms subroutine
  217.  
  218.     HR=FIX(DEG)
  219.     C=ABS(HR-DEG)*60D0
  220.     MIN=INT(C)
  221.     D=(C-MIN)*60D0
  222.     SEC=INT(D+.5D0)
  223.  
  224. END SUB
  225.  
  226. '***************************************************************
  227.  
  228. SUB DISPLAY.DATA STATIC
  229.  
  230.     ' display data subroutine
  231.  
  232.     FORMAT$="########.###"
  233.  
  234.     CLS
  235.     PRINT
  236.     PRINT TAB(24);"* 1950 < FK4 > *"
  237.     PRINT
  238.     a$=STR$(RAS.HR.1950)+" hours"+STR$(RAS.MIN.1950)+" minutes"+STR$(RAS.SEC.1950)+" seconds"
  239.     PRINT TAB(5);"Right ascension";TAB(60-LEN(a$));a$
  240.     a$=STR$(DEC.DEG.1950)+" degrees"+STR$(DEC.MIN.1950)+" minutes"+STR$(DEC.SEC.1950)+" seconds"
  241.     PRINT TAB(5);"Declination";TAB(60-LEN(a$));a$
  242.     PRINT
  243.     PRINT TAB(5);"Right ascension proper motion";
  244.     PRINT USING FORMAT$;TAB(48);RAS.PMOTION.1950
  245.     PRINT TAB(5);"Declination proper motion";
  246.     PRINT USING FORMAT$;TAB(48);DEC.PMOTION.1950
  247.     PRINT
  248.     PRINT TAB(5);"Parallax        (arc seconds)";
  249.     PRINT USING FORMAT$;TAB(48);PARALLAX.1950
  250.     PRINT TAB(5);"Radial velocity (kilometers per second)";
  251.     PRINT USING FORMAT$;TAB(48);RADIALV.1950
  252.     PRINT
  253.     PRINT TAB(24);"* 2000 < FK5 > *"
  254.     PRINT
  255.     a$=STR$(RAS.HR.2000)+" hour"+STR$(RAS.MIN.2000)+" minutes"+STR$(RAS.SEC.2000)+" seconds"
  256.     PRINT TAB(5);"Right ascension";TAB(60-LEN(a$));a$
  257.     a$=STR$(DEC.DEG.2000)+" degrees"+STR$(DEC.MIN.2000)+" minutes"+STR$(DEC.SEC.2000)+" seconds"
  258.     PRINT TAB(5);"Declination";TAB(60-LEN(a$));a$
  259.     PRINT
  260.     PRINT TAB(5);"Right ascension proper motion";
  261.     PRINT USING FORMAT$;TAB(48);RAS.PMOTION.2000
  262.     PRINT TAB(5);"Declination proper motion";
  263.     PRINT USING FORMAT$;TAB(48);DEC.PMOTION.2000
  264.     PRINT
  265.     PRINT TAB(5);"Parallax        (arc seconds)";
  266.     PRINT USING FORMAT$;TAB(48);PARALLAX.2000
  267.     PRINT TAB(5);"Radial velocity (kilometers per second)";
  268.     PRINT USING FORMAT$;TAB(48);RADIALV.2000
  269.     CALL KEYCHECK
  270.  
  271. END SUB
  272.  
  273. '***************************************************************
  274.  
  275. SUB ATAN3(SINANGLE,COSANGLE,ANGLE) STATIC
  276.  
  277.     ' four quadrant inverse tangent subroutine
  278.  
  279.     IF ABS(SINANGLE)<1D-08 THEN
  280.        ANGLE=(1-SGN(COSANGLE))*PIDIV2
  281.        EXIT SUB
  282.     END IF
  283.  
  284.        ANGLE=(2-SGN(SINANGLE))*PIDIV2
  285.  
  286.     IF ABS(COSANGLE)<1D-08 THEN
  287.       EXIT SUB
  288.     ELSE
  289.       ANGLE=ANGLE+SGN(SINANGLE)*SGN(COSANGLE)*(ABS(ATN(SINANGLE/COSANGLE))-PIDIV2)
  290.     END IF
  291.  
  292. END SUB
  293.  
  294. '***************************************************************
  295.  
  296. SUB ARCSIN(SINANGLE,ANGLE) STATIC
  297.  
  298.     ' inverse sine subroutine
  299.  
  300.     IF ABS(SINANGLE)>=1 THEN
  301.        ANGLE=SGN(SINANGLE)*PIDIV2
  302.     ELSE
  303.        ANGLE=ATN(SINANGLE/SQR(1-SINANGLE^2))
  304.     END IF
  305.  
  306. END SUB
  307.  
  308. '***************************************************************
  309.  
  310. SUB INTRODUCTION STATIC
  311.  
  312.     ' program introduction subroutine
  313.  
  314.     CLS
  315.     PRINT
  316.     PRINT TAB(10);"J2000 is an interactive QuickBASIC program which can be"
  317.     PRINT TAB(10);"used to convert stellar positions, proper motions,"
  318.     PRINT TAB(10);"parallax and radial velocity from the standard epoch"
  319.     PRINT TAB(10);"B1950 (FK4) to J2000 (FK5). The method used in this"
  320.     PRINT TAB(10);"program is based on the algorithm published in the 1985"
  321.     PRINT TAB(10);"edition of the Astronomical Almanac, Addendum to Section"
  322.     PRINT TAB(10);"B, page X. This technique is an approximation because"
  323.     PRINT TAB(10);"it does not account for systematic and individual star"
  324.     PRINT TAB(10);"corrections between the FK4 and FK5 epochs."
  325.     PRINT
  326.     PRINT TAB(10);"J2000 will first request the right ascension, declination,"
  327.     PRINT TAB(10);"right ascension proper motion, declination proper motion,"
  328.     PRINT TAB(10);"parallax and radial velocity of a stellar object with"
  329.     PRINT TAB(10);"respect to the 1950 (FK4) epoch. Each program prompt will"
  330.     PRINT TAB(10);"advise you of the proper units to input. If you do not"
  331.     PRINT TAB(10);"have a number for proper motion or parallax or radial"
  332.     PRINT TAB(10);"velocity, then input '0' for any of these values."
  333.     CALL KEYCHECK
  334.  
  335. END SUB
  336.  
  337. '***************************************************************
  338.  
  339. SUB KEYCHECK STATIC
  340.  
  341.     ' check user response subroutine
  342.  
  343.     PRINT
  344.     PRINT
  345.     PRINT TAB(20);"< press any key to continue >";
  346.  
  347.     A$=""
  348.     WHILE A$=""
  349.       A$=INKEY$
  350.     WEND
  351.  
  352. END SUB
  353.  
  354.