/usr/share/gnudatalanguage/astrolib/helio.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 190 191 | PRO HELIO, JD, LIST, HRAD, HLONG, HLAT, RADIAN = radian
;+
; NAME:
; HELIO
; PURPOSE:
; Compute (low-precision) heliocentric coordinates for the planets.
; EXPLANATION:
; The mean orbital elements for epoch J2000 are used. These are derived
; from a 250 yr least squares fit of the DE 200 planetary ephemeris to a
; Keplerian orbit where each element is allowed to vary linearly with
; time. For dates between 1800 and 2050, this solution fits the
; terrestrial planet orbits to ~25" or better, but achieves only ~600"
; for Saturn.
;
; Use PLANET_COORDS (which calls HELIO) to get celestial (RA, Dec)
; coordinates of the planets
; CALLING SEQUENCE:
; HELIO, JD, LIST, HRAD, HLONG, HLAT, [/RADIAN]
; INPUTS:
; JD = Julian date, double precision scalar or vector
; LIST = List of planets array. May be a single number.
; 1 = merc, 2 = venus, ... 9 = pluto.
;
; OUTPUTS:
; HRAD = array of Heliocentric radii (A.U).
; HLONG = array of Heliocentric (ecliptic) longitudes (degrees).
; HLAT = array of Heliocentric latitudes (degrees).
; These output parameters will be dimensioned Nplanet by Ndate,
; where Nplanet is the number of elements of list, and Ndate is
; the number of elements of JD.
;
; OPTIONAL INPUT KEYWORD:
; /RADIAN - If set, then the output longitude and latitude are given in
; radians.
; EXAMPLE:
; (1) Find the current heliocentric positions of all the planets
;
; IDL> GET_JULDATE, jd ;Get current Julian date
; IDL> HELIO,jd,indgen(9)+1,hrad,hlong,hlat ;Get radius, long, and lat
;
; (2) Find heliocentric position of Mars on August 23, 2018
; IDL> JDCNV, 2018,08,23,0,jd
; IDL> HELIO,JD,4,HRAD,HLONG,HLAT
; ===> hrad = 1.6407 AU hlong = 124.3197 hlat = 1.7853
; For comparison, the JPL ephemeris gives
; hrad = 1.6407 AU hlong = 124.2985 hlat = 1.7845
; (3) Find the heliocentric positions of Mars and Venus for every day in
; November 2018
; IDL> JDCNV, 2018, 11, 1, 0, jd ;Julian date of November 1, 2000
; IDL> helio, jd+indgen(30), [4,2], hrad,hlong,hlat ;Mars=4, Venus=2
; hrad, hlong, and hlat will be dimensioned [2,30]
; first column contains Mars data, second column Venus
; COMMON BLOCKS:
; None
; ROUTINES USED:
; CIRRANGE - force angle between 0 and 2*!PI
; NOTES:
; (1) The calling sequence for this procedure was changed in August 2000
; (2) This program is based on the two-body model and thus neglects
; interactions between the planets. This is why the worst results
; are for Saturn. Use the procedure JPLEPHINTERP for more accurate
; positions using the JPL ephemeris. Also see JPL Horizons
; ( https://ssd.jpl.nasa.gov/horizons.cgi ) for a more accurate ephemeris
; generator online.
; (3) The coordinates are given for equinox 2000 and *not* the equinox
; of the supplied date(s)
; MODIFICATION HISTORY:
; R. Sterner. 20 Aug, 1986.
; Code cleaned up a bit W. Landsman December 1992
; Major rewrite, use modern orbital elements, vectorize, more accurate
; solution to Kepler's equation W. Landsman August 2000
; Wasn't working for planet vectors W. Landsman August 2000
; Work for more than 32767 positions S. Leach Jan 2009
; Use data from https://ssd.jpl.nasa.gov/txt/p_elem_t1.txt W. Landsman Jun 2017
;-
On_error,2
compile_opt idl2
if N_params() LT 3 then begin
print,'Syntax - Helio, jd, list, hrad, hlong, hlat, [/RADIAN]'
print,' jd - Scalar or vector Julian date'
print,' list - scalar or vector of planet numbers [1-9]'
print, $
' hrad, hlong, hlat - output heliocentric distance, longitude latitude'
return
endif
; Mean orbital elements taken from https://ssd.jpl.nasa.gov/txt/p_elem_t1.txt
; (1) semi-major axis in AU, (2) eccentricity, (3) inclination (degrees),
; (4) mean longitude (degrees), (5) longitude of perihelion (degrees)
; and (6) longitude of the ascending node (degrees)
;Mercury
PD = [ [0.38709927d, 0.20563593d, 7.00497902d, 252.25032350d, 77.45779628d , 48.33076593d], $
;Venus
[ 0.72333566d, 0.00677672d, 3.39467605d, 181.97909950d, 131.60246718d, 76.67984255d], $
;EMBary
[ 1.00000261d, 0.01671123d, -0.00001531d, 100.46457166d, 102.93768193d, 0.0], $
;Mars
[ 1.52371034d, 0.09339410d, 1.84969142d, -4.55343205d, -23.94362959d, 49.55953891d], $
;Jupiter
[ 5.20288700d, 0.04838624d, 1.30439695d, 34.39644051d, 14.72847983d, 100.47390909d], $
;Saturn
[ 9.53667594d, 0.05386179d, 2.48599187d, 49.95424423d, 92.59887831d, 113.66242448d], $
;Uranus
[ 19.18916464d, 0.04725744d, 0.77263783d, 313.23810451d, 170.95427630d, 74.01692503d], $
;Neptune
[ 30.06992276d, 0.00859048d, 1.77004347d, -55.12002969d, 44.96476227d, 131.78422574d], $
;Pluto
[ 39.48211675d, 0.24882730d, 17.14001206d, 238.92903833d, 224.06891629d, 110.30393684d]]
DPD = [[ 0.00000037d, 0.00001906d, -0.00594749d, 149472.67411175d, 0.16047689d, -0.12534081 ], $
[0.00000390d, -0.00004107d, -0.00078890d, 58517.81538729d, 0.00268329d, -0.27769418 ], $
[0.00000562d, -0.00004392d, -0.01294668d, 35999.37244981d, 0.32327364d, 0.0d], $
[ 0.00001847d, 0.00007882d, -0.00813131d, 19140.30268499d, 0.44441088d, -0.29257343d], $
[ -0.00011607d, -0.00013253d, -0.00183714d, 3034.74612775d, 0.21252668d, 0.20469106d],$
[ -0.00125060d, -0.00050991d, 0.00193609d, 1222.49362201d, -0.41897216d, -0.28867794d], $
[ -0.00196176d, -0.00004397d, -0.00242939d, 428.48202785d, 0.40805281d, 0.04240589d], $
[ 0.00026291d, 0.00005105d, 0.00035372d, 218.45945325d, -0.32241464d, -0.00508664d], $
[ -0.00031596d, 0.00005170d, 0.00004818d, 145.20780515d, -0.04062942d, -0.01183482d]]
JD0 = 2451545.0d ;Julian Date for Epoch 2000.0
radeg = 180/!DPI
;----------------- Days since Epoch ---------------
T = (JD - JD0)/36525.0d ;Time in centuries since 2000.0
ip = list-1
ntime = N_elements(t)
nplanet = N_elements(list)
hrad = fltarr(nplanet,ntime) & hlong = hrad & hlat = hrad
;----------------- Loop over dates --------------
for i =0L,ntime-1L do begin ;SML made longword
pd1 = pd[*,ip] + dpd[*,ip]*T[i]
a = pd1[0,*] ;semi-major axis
eccen = pd1[1,*] ;eccentricity
inc = pd1[2,*]/RADEG ;inclination in radians
L = pd1[3,*]/RADEG ;mean longitude
pi = pd1[4,*]/RADEG ;longitude of the perihelion
omega = pd1[5,*]/RADEG ;longitude of the ascending node
m = L - pi
cirrange,m,/RADIAN
e1 = m + (m + eccen*sin(m) - m)/(1 - eccen*cos(m) )
e = e1 + (m + eccen*sin(e1) - e1)/(1 - eccen*cos(e1) )
maxdif = max(abs(e-e1))
niter = 0
while (maxdif GT 1e-7) && (niter lt 10) do begin
e1 = e
e = e1 + (m + eccen*sin(e1) - e1)/(1 - eccen*cos(e1) )
maxdif = max(abs(e-e1))
niter++
endwhile
nu = 2*atan( sqrt( (1+eccen)/(1-eccen) )* tan(E/2)) ;true anomaly
hrad[0,i] = reform( a*(1 - eccen*cos(e) ) )
hlong[0,i] = reform (nu + pi)
hlat[0,i] = reform( asin(sin(hlong[*,i] - omega)*sin(inc) ) )
endfor
cirrange,hlong,/RADIAN
if ~keyword_set(RADIAN) then begin
hlong = hlong*RADEG
hlat = hlat*RADEG
endif
if N_elements(hrad) GT 1 then begin
hrad = reform(hrad,/over)
hlong = reform(hlong,/over)
hlat = reform(hlat,/over)
endif else begin
if N_elements(size(jd)) EQ 3 then begin ;scalar?
hrad = hrad[0]
hlong = hlong[0]
hlat = hlat[0]
endif
endelse
return
end
|