/usr/share/gnudatalanguage/astrolib/xyz.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 | pro xyz,date,x,y,z,xvel,yvel,zvel,equinox=equinox
;+
; NAME:
; XYZ
; PURPOSE:
; Calculate geocentric X,Y, and Z and velocity coordinates of the Sun
; EXPLANATION:
; Calculates geocentric X,Y, and Z vectors and velocity coordinates
; (dx, dy and dz) of the Sun. (The positive X axis is directed towards
; the equinox, the y-axis, towards the point on the equator at right
; ascension 6h, and the z axis toward the north pole of the equator).
; Typical position accuracy is <1e-4 AU (15000 km).
;
; CALLING SEQUENCE:
; XYZ, date, x, y, z, [ xvel, yvel, zvel, EQUINOX = ]
;
; INPUT:
; date: reduced julian date (=JD - 2400000), scalar or vector
;
; OUTPUT:
; x,y,z: scalars or vectors giving heliocentric rectangular coordinates
; (in A.U) for each date supplied. Note that sqrt(x^2 + y^2
; + z^2) gives the Earth-Sun distance for the given date.
; xvel, yvel, zvel: velocity vectors corresponding to X, Y and Z.
;
; OPTIONAL KEYWORD INPUT:
; EQUINOX: equinox of output. Default is 1950.
;
; EXAMPLE:
; What were the rectangular coordinates and velocities of the Sun on
; Jan 22, 1999 0h UT (= JD 2451200.5) in J2000 coords? NOTE:
; Astronomical Almanac (AA) is in TDT, so add 64 seconds to
; UT to convert.
;
; IDL> xyz,51200.5+64.d/86400.d,x,y,z,xv,yv,zv,equinox = 2000
;
; Compare to Astronomical Almanac (1999 page C20)
; X (AU) Y (AU) Z (AU)
; XYZ: 0.51456871 -0.76963263 -0.33376880
; AA: 0.51453130 -0.7697110 -0.3337152
; abs(err): 0.00003739 0.00007839 0.00005360
; abs(err)
; (km): 5609 11759 8040
;
; NOTE: Velocities in AA are for Earth/Moon barycenter
; (a very minor offset) see AA 1999 page E3
; X VEL (AU/DAY) YVEL (AU/DAY) Z VEL (AU/DAY)
; XYZ: -0.014947268 -0.0083148382 -0.0036068577
; AA: -0.01494574 -0.00831185 -0.00360365
; abs(err): 0.000001583 0.0000029886 0.0000032077
; abs(err)
; (km/sec): 0.00265 0.00519 0.00557
;
; PROCEDURE CALLS:
; PRECESS_XYZ
; REVISION HISTORY
; Original algorithm from Almanac for Computers, Doggett et al. USNO 1978
; Adapted from the book Astronomical Photometry by A. Henden
; Written W. Landsman STX June 1989
; Correct error in X coefficient W. Landsman HSTX January 1995
; Added velocities, more terms to positions and EQUINOX keyword,
; some minor adjustments to calculations
; P. Plait/ACC March 24, 1999
;-
On_error,2
if (n_params() eq 0) then begin
print,'Syntax - XYZ, date, x, y, z, [ xvel, yvel, zvel, EQUINOX= ]'
print,' (date is REDUCED Julian date (JD - 2400000.0) )'
return
endif
picon = !DPI/180.0d
t = (date - 15020.0d0)/36525.0d0 ;Relative Julian century from 1900
; NOTE: longitude arguments below are given in *equinox* of date.
; Precess these to equinox 1950 to give everything an even footing.
; Compute argument of precession from equinox of date back to 1950
pp = (1.396041d + 0.000308d*(t + 0.5d))*(t-0.499998d)
; Compute mean solar longitude, precessed back to 1950
el = 279.696678D + 36000.76892D*t + 0.000303d*t*t - pp
; Compute Mean longitude of the Moon
c = 270.434164d + 480960.d*t + 307.883142d*t - 0.001133d*t*t - pp
; Compute longitude of Moon's ascending node
n = 259.183275d - 1800.d*t - 134.142008d*t + 0.002078d*t*t - pp
; Compute mean solar anomaly
g = 358.475833d + 35999.04975d*t - 0.00015d*t*t
; Compute the mean jupiter anomaly
j = 225.444651d + 2880.0d*t + 154.906654d*t*t
; Compute mean anomaly of Venus
v = 212.603219d + 58320.d*t + 197.803875d*t + 0.001286d*t*t
; Compute mean anomaly of Mars
m = 319.529425d + 19080.d*t + 59.8585d*t + 0.000181d*t*t
; Convert degrees to radians for trig functions
el = el*picon
g = g*picon
j = j*picon
c = c*picon
v = v*picon
n = n*picon
m = m*picon
; Calculate X,Y,Z using trigonometric series
X = 0.999860d*cos(el) $
- 0.025127d*cos(g - el) $
+ 0.008374d*cos(g + el) $
+ 0.000105d*cos(g + g + el) $
+ 0.000063d*t*cos(g - el) $
+ 0.000035d*cos(g + g - el) $
- 0.000026d*sin(g - el - j) $
- 0.000021d*t*cos(g + el) $
+ 0.000018d*sin(2.d*g + el - 2.d*v) $
+ 0.000017d*cos(c) $
- 0.000014d*cos(c - 2.d*el) $
+ 0.000012d*cos(4.d*g + el - 8.d*m + 3.d*j) $
- 0.000012d*cos(4.d*g - el - 8.d*m + 3.d*j) $
- 0.000012d*cos(g + el - v) $
+ 0.000011d*cos(2.d*g + el - 2.d*v) $
+ 0.000011d*cos(2.d*g - el - 2.d*j)
Y = 0.917308d*sin(el) $
+ 0.023053d*sin(g - el) $
+ 0.007683d*sin(g + el) $
+ 0.000097d*sin(g + g + el) $
- 0.000057d*t*sin(g - el) $
- 0.000032d*sin(g + g - el) $
- 0.000024d*cos(g - el - j) $
- 0.000019d*t*sin(g + el) $
- 0.000017d*cos(2.d0*g + el - 2.d0*v) $
+ 0.000016d*sin(c) $
+ 0.000013d*sin(c - 2.d0*el ) $
+ 0.000011d*sin(4.d0*g + el - 8.d0*m + 3.d0*j) $
+ 0.000011d*sin(4.d0*g - el - 8.d0*m + 3.d0*j) $
- 0.000011d*sin(g + el - v) $
+ 0.000010d*sin(2.d0*g + el - 2.d0*v ) $
- 0.000010d*sin(2.d0*g - el - 2.d0*j )
Z = 0.397825d*sin(el) $
+ 0.009998d*sin(g-el) $
+ 0.003332d*sin(g+el) $
+ 0.000042d*sin(g+g+el) $
- 0.000025d*t*sin(g-el) $
- 0.000014d*sin(g+g-el) $
- 0.000010d*cos(g-el-j)
;Precess_to new equator?
if keyword_set(equinox) then precess_xyz, x, y, z, 1950, equinox
if N_params() LE 3 then return
XVEL = -0.017200d * sin(el) $
-0.000288d * sin(g + el) $
-0.000005d * sin(2.d0*g + el) $
-0.000004d * sin(c) $
+0.000003d * sin(c - 2.d0*el) $
+0.000001d *t * sin(g+el) $
-0.000001d * sin(2.d0*g-el)
YVEL = 0.015780 * cos(el) $
+0.000264 * cos(g + el) $
+0.000005 * cos(2.d0*g + el) $
+0.000004 * cos(c) $
+0.000003 * cos(c - 2.d0*el) $
-0.000001 * t * cos(g + el)
ZVEL = 0.006843 * cos(el) $
+0.000115 * cos(g + el) $
+0.000002 * cos(2.d0*g + el) $
+0.000002 * cos(c) $
+0.000001 * cos(c - 2.d0*el)
;Precess to new equator?
if keyword_set(equinox) then precess_xyz, xvel, yvel, zvel, 1950, equinox
return
end
|