;+ ; NAME: ; JPLEPHREAD ; ; AUTHOR: ; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 ; craigm@lheamail.gsfc.nasa.gov ; UPDATED VERSIONs can be found on my WEB PAGE: ; http://cow.physics.wisc.edu/~craigm/idl/idl.html ; ; PURPOSE: ; Open and read JPL DE200 or DE405 Ephemeride FITS File ; ; MAJOR TOPICS: ; Planetary Orbits, Interpolation ; ; CALLING SEQUENCE: ; JPLEPHREAD, FILENAME, INFO, RAWDATA, JDLIMITS, STATUS=, ERRMSG= ; ; DESCRIPTION: ; ; JPLEPHREAD opens and reads the JPL DE200 or DE405 planetary ; ephemerides, as available in FITS format. The user must have the ; IDL Astronomy Library installed to use this routine. ; ; This routine is the initialization stage of a two-stage process to ; interpolate the JPL ephemeris. In this first stage, the file is ; opened, and the relevant portions of the table are read and stored ; into the two variables INFO and RAWDATA. In the second stage, the ; user actually interpolates the ephemeris for the desired bodies ; and to the desired ephemeris time using JPLEPHINTERP. ; ; Users must decide ahead of time the approximate dates of interest, ; and pass this range in the JDLIMITS parameter. Any date covered ; by the ephemeris is valid. ; ; JPLEPHREAD is able to read files of the following format: ; DE200 - Chebyshev - FITS format - Note 1 ; DE405 - Chebyshev - FITS format - Note 1 ; DE200 - Taylor - FITS format - Note 2 ; ; Note 1 - Chebyshev formatted FITS files are available in the ; AXBARY package by Arnold Rots, found here: ; ftp://heasarc.gsfc.nasa.gov/xte/calib_data/clock/bary/ ; or at the Markwardt FTP site: ; ftp://cow.physics.wisc.edu/pub/craigm/bary/ ; ; Note 2 - Taylor-series based ephemerides have been available for ; years in the FTOOLS / LHEASOFT package produced by NASA's ; Goddard Space Flight Center. The original file is ; de200_new.fits, which covers the years 1959-2000, ; inclusive. A newer file is named ; de200_1950-2050_v2.fits, and covers the years 1959-2050. ; See Markwardt FTP site for these files. ; ; PARAMETERS: ; ; FILENAME - name of ephemeris file (scalar string). ; ; INFO - upon completion, information about the ephemeris data is ; returned in this parameter in the form of a structure. ; Users must not modify INFO, although several fields are ; useful and may be accessed read-only: ; TSTART/TSTOP (start and stop time of data in Julian ; days); ; C (speed of light in km/s); ; DENUM (development ephemeris number [200 or 405]) ; AU (1 astronomical unit, in units of light-seconds) ; ; RAWDATA - upon completion, raw ephemeris data is returned in this ; parameter. Users are not meant to access this data ; directly, but rather to pass it to JPLEPHINTERP. ; ; JDLIMITS - a two-element vector (optional), describing the desired ; time range of interest. The vector should have the ; form [TSTART, TSTOP], where TSTART and TSTOP are the ; beginning and ending times of the range, expressed in ; Julian days. ; Default: entire table is read (note, this can be ; several megabytes) ; ; ; KEYWORD PARAMETERS: ; ; STATUS - upon completion, a value of 1 indicates success, and 0 ; indicates failure. ; ; ERRMSG - upon completion, an error message is returned in this ; keyword. If there were no errors, then the returned ; value is the empty string, ''. ; ; ; EXAMPLE: ; ; Find position of earth at ephemeris time 2451544.5 JD. Units are ; in Astronomical Units. ; ; JPLEPHREAD, 'JPLEPH.200', pinfo, pdata, [2451544D, 2451545D] ; ; JPLEPHINTERP, pinfo, pdata, 2451544.5D, xearth, yearth, zearth, $ ; /EARTH, posunits='AU' ; ; ; REFERENCES: ; ; AXBARY, Arnold Rots. ; ftp://heasarc.gsfc.nasa.gov/xte/calib_data/clock/bary/ ; ; HORIZONS, JPL Web-based ephermis calculator (Ephemeris DE406) ; http://ssd.jpl.nasa.gov/horizons.html ; ; JPL Export Ephmeris FTP Site ; ftp://navigator.jpl.nasa.gov/pub/ephem/export/ ; (ephemeris files are available here, however, they must be ; converted to FITS format using the "bin2eph" utility found in ; AXBARY) ; ; JPL Export Ephemeris CD-ROM - Ordering Information ; http://www.willbell.com/software/jpl.htm ; ; Standish, E.M. 1982, "Orientation of the JPL Ephemerides, ; DE200/LE200, to the Dynamical Equinox of J2000", Astronomy & ; Astrophysics, vol. 114, pp. 297-302. ; ; Standish, E.M.: 1990, "The Observational Basis for JPL's DE200, ; the planetary ephemeris of the Astronomical Almanac", Astronomy ; & Astrophysics, vol. 233, pp. 252-271. ; ; SEE ALSO ; JPLEPHREAD, JPLEPHINTERP, JPLEPHTEST ; ; MODIFICATION HISTORY: ; Written and Documented, CM, Jun 2001 ; ; $Id: jplephread.pro,v 1.6 2001/07/01 03:32:02 craigm Exp $ ; ;- ; Copyright (C) 2001, Craig Markwardt ; This software is provided as is without any warranty whatsoever. ; Permission to use, copy and distribute unmodified copies for ; non-commercial purposes, and to modify and use for personal or ; internal use, is granted. All other rights are reserved. ;- function jplephpar, header, parname, default=default, fatal=fatal forward_function fxpar ; '$Id: jplephread.pro,v 1.6 2001/07/01 03:32:02 craigm Exp $' value = fxpar(header, parname, Count = N_value) if N_value EQ 0 then begin if keyword_set(fatal) then $ message, 'ERROR: keyword '+strupcase(parname)+' was not found' return, default endif return, value end function jplephval, names, values, name, default=default, fatal=fatal wh = where(names EQ strupcase(name), ct) if ct EQ 0 then begin if keyword_set(fatal) then $ message, 'ERROR: value '+strupcase(name)+' was not found in file' return, default endif return, values[wh[0]] end pro jplephread, filename, info, raw, jdlimits, $ status=status, errmsg=errmsg status = 0 printerror = 1 - arg_present(errmsg) errmsg = '' if n_params() EQ 0 then begin message, 'USAGE: JPLEPHREAD, filename, info, rawdata, jdlimits', /info return endif ; if n_elements(jdlimits) LT 2 then begin ; errmsg = 'ERROR: You must specify JDLIMITS' ; return ; endif fxbopen, unit, filename, 1, ephhead, errmsg=errmsg if errmsg NE '' then $ if printerror then message,errmsg else return extname = strtrim(fxpar(ephhead, 'EXTNAME'),2) ttype1 = strtrim(fxpar(ephhead, 'TTYPE1'),2) if (extname EQ 'EPHEM' AND ttype1 EQ 'EARTH') then begin ;; This is the DE200_NEW format (standard FTOOLS) nrows = fxpar(ephhead, 'NAXIS2') tstart = fxpar(ephhead, 'TSTART') tstop = fxpar(ephhead, 'TSTOP') timedel = jplephpar(ephhead, 'TIMEDEL', default=1D) ;; 1-day default ;; Constants from XTEBARYCEN.F C=2.99792458D+8 TWOPI=6.28318530717958648D0 DAYSEC=1.D0/86400.D0 AULTSC=499.004782D0 GAUSS=0.01720209895D0 RSCHW=(GAUSS^2)*(AULTSC^3)*(DAYSEC^2) SUNRAD=2.315D0 if n_elements(jdlimits) GE 2 then begin if (min(jdlimits) LT tstart OR $ max(jdlimits) GT tstop) then begin errmsg = 'ERROR: '+filename+$ ' does not cover the time of interest' fxbclose, unit return endif ;; Expand by one row either side rowlimits = floor((jdlimits-tstart)/timedel) + [-2,2] rowlimits = rowlimits > 1 < nrows endif else begin jdlimits = [tstart, tstop] rowlimits = [1L, nrows] endelse ;; Read raw data fxbread, unit, cearth, 'EARTH', rowlimits, errmsg=errmsg if errmsg EQ '' then $ fxbread, unit, csun, 'SUN', rowlimits, errmsg=errmsg if errmsg EQ '' then $ fxbread, unit, ctdb2tdt, 'TIMEDIFF', rowlimits, errmsg=errmsg fxbclose, unit if errmsg NE '' then $ if printerror then message,errmsg else return nr = rowlimits[1]-rowlimits[0]+1 t0 = dindgen(nr)*timedel - (jdlimits[1]-jdlimits[0])/2D dtt = spl_init(t0, ctdb2tdt) raw = reform(dblarr(18, nr), 18, nr, /overwrite) raw[0 :11,*] = cearth * c/1000D ;; units of lt-s raw[12:14,*] = csun * c/1000D ;; units of lt-s/day raw[15, *] = t0 raw[16 ,*] = ctdb2tdt raw[17 ,*] = dtt jdlimits1 = (rowlimits+[-1,0])*timedel + tstart info = {nrows: nrows, tstart: tstart, tstop: tstop, $ timedel: timedel, format: 'DENEW', $ denum: 200L, c: c, emrat: 0.012150586D, $ au: aultsc, msol: rschw, sunrad: sunrad, $ jdlimits: jdlimits1, jdrows: nr } endif else if (extname EQ 'DE1' AND ttype1 EQ 'Cname') then begin ;; This is the BINEPH2FITS format (either DE200 or DE405) ;; --------------------------------------------- ;; First extension contains parameter data fxbread, unit, cname, 'Cname' fxbread, unit, cvalue, 'Cvalue' cname = strtrim(cname,2) denum = 0L & clight = 0D & emrat = 0D & au = 0D msol = 0D & radsol = 0D denum = round(jplephval(cname, cvalue, 'DENUM', /fatal)) clight = jplephval(cname, cvalue, 'CLIGHT', /fatal) emrat = jplephval(cname, cvalue, 'EMRAT', /fatal) au = jplephval(cname, cvalue, 'AU', /fatal) msol = jplephval(cname, cvalue, 'GMS', /fatal) radsol = jplephval(cname, cvalue, 'RADS', default=-1D) if radsol EQ -1D then $ radsol = jplephval(cname, cvalue, 'ASUN', default=-1D) emrat = 1D / (1D + emrat) if clight EQ 0 then begin errmsg = 'ERROR: Could not load physical constants from '+filename fxbclose, unit return endif x = au / clight msol = msol * x * x * x / 86400D^2 ;; GM_sun (in lt sec) radsol = radsol / clight ;; Solar radius (lt sec) clight = clight * 1000 ;; Speed of light (m/s) fxbclose, unit ;; --------------------------------------------- ;; Second extension contains accounting data fxbopen, unit, filename, 2, ephhead, errmsg=errmsg if errmsg NE '' then $ if printerror then message,errmsg else return extname = strtrim(fxpar(ephhead, 'EXTNAME'),2) if extname NE 'DE2' then begin errmsg = 'ERROR: '+filename+' is not a JPL ephemeris file' fxbclose, unit return endif fxbread, unit, ephobj, 'Object', errmsg=errmsg if errmsg EQ '' then $ fxbread, unit, ephptr, 'Pointer', errmsg=errmsg if errmsg EQ '' then $ fxbread, unit, ephncoeff, 'NumCoeff', errmsg=errmsg if errmsg EQ '' then $ fxbread, unit, ephnsub, 'NumSubIntv', errmsg=errmsg fxbclose, unit if errmsg NE '' then begin errmsg = 'ERROR: could not read '+filename+' extension 2' if printerror then message,errmsg else return endif ;; Trim each object name to first word only for i = 0, n_elements(ephobj)-1 do begin ephobj[i] = strupcase((str_sep(ephobj[i], ' '))[0]) endfor ;; --------------------------------------------- ;; Third extension contains Chebyshev coefficients fxbopen, unit, filename, 3, ephhead, errmsg=errmsg if errmsg NE '' then return extname = strtrim(fxpar(ephhead, 'EXTNAME'),2) if extname NE 'DE3' then begin errmsg = 'ERROR: '+filename+' is not a JPL ephemeris file' fxbclose, unit if printerror then message,errmsg else return endif nrows = fxpar(ephhead, 'NAXIS2') tstart = fxpar(ephhead, 'TSTART') tstop = fxpar(ephhead, 'TSTOP') timedel = jplephpar(ephhead, 'TIMEDEL', default=32D) ;; 32-day default if floor((tstop-tstart + 0.5)/timedel) NE nrows then begin errmsg = 'ERROR: Incorrect number of rows in '+filename fxbclose, unit if printerror then message,errmsg else return endif if n_elements(jdlimits) GE 2 then begin if (min(jdlimits) LT tstart OR $ max(jdlimits) GT tstop) then begin errmsg = 'ERROR: '+filename+$ ' does not cover the time of interest' fxbclose, unit if printerror then message,errmsg else return endif ;; Expand by two rows either side rowlimits = floor((jdlimits-tstart)/timedel) + [-2,2] rowlimits = rowlimits > 1 < nrows endif else begin jdlimits = [tstart, tstop] rowlimits = [1L, nrows] endelse ;; Read raw data dims = fxbdimen(unit, 'ChebCoeffs') fxbread, unit, coeffs, 'ChebCoeffs', rowlimits, errmsg=errmsg fxbclose, unit if errmsg NE '' then $ if printerror then message,errmsg else return raw = reform(coeffs, [dims, rowlimits[1]-rowlimits[0]+1], /overwrite) jdlimits1 = (rowlimits+[-1,0])*timedel + tstart if (abs(min(raw[0,*]) - jdlimits1[0]) GT 1d-6 OR $ abs(max(raw[1,*]) - jdlimits1[1]) GT 1d-6) then begin errmsg = 'ERROR: JDLIMITS and time column do not match' if printerror then message,errmsg else return endif nr = rowlimits[1]-rowlimits[0]+1 info = {nrows: nrows, tstart: tstart, tstop: tstop, $ timedel: timedel, format: 'BINEPH2FITS', $ denum: denum, c: clight, emrat: emrat, $ au: au*1000/clight, msol: msol, sunrad: radsol, $ jdlimits: jdlimits1, jdrows: nr, $ objname: ephobj, ptr: ephptr, ncoeff: ephncoeff, $ nsub: ephnsub} ; aufac: 1D/clight, velfac: 2D/(timedel*86400D), $ endif else begin errmsg = 'ERROR: '+filename+' was not in a recognized format' fxbclose, unit if printerror then message,errmsg else return endelse errmsg = '' status = 1 return end