;+
; NAME:
; UTC2BJD
; PURPOSE:
; Converts a Julian Date in Coordinated Universal Time (JD_UTC) to
; Barycentric Julian Date in Barycentric Dynamical Time
; (BJD_TDB) from anywhere on Earth or anywhere for which HORIZONS can
; generate an ephemeris. Uses the plane parallel approximation
; unless TARGET is specified and an ephemeris can be generated by
; HORIZONS for it (in general not a good approximation for objects
; in the Solar System). In the geocentric, plane parallel case, it
; agrees with BARYCEN to 50 ns, which agrees with AXBARY to ~1
; us. See below for notes on precision and dependencies.
;
; CALLING SEQUENCE
; BJD_TDB = UTC2BJD(jd_utc, ra, dec
; [,lat=lat, lon=lon, elevation=elevation,
; earthobs=obsname, spaceobs=spaceobs, tbase=tbase,
; /1950, /TIME_DIFF, /TT_IN])
;
; INPUTS:
; JD_UTC - A scalar or array of JDs (in UTC). Must be double
; precision.
; NOTE 1: If TBASE is used, the input must be JD_UTC - TBASE
; NOTE 2: If the keyword TT_IN is set, must be an
; array of JDs in TT.
; JD_TT = JD_UTC + 32.184 + N
; where N is the number of leap seconds since 1961.
; RA/DEC - A scalar specifying the Right Ascension and
; Declinatoin of the object, in decimal degrees. J2000
; is assumed unless /B1950 is set
;
; OPTIONAL INPUTS:
; SPACEOBS - A string input to HORIZONS specifying the space based
; observatory. Can be anything for which HORIZONS can
; generate an ephemeris. If set,
; EARTHOBS/LAT/LON/ELEVATION are ignored.
; http://www-int.stsci.edu/~sontag/spicedocs/req/naif_ids.html
; The name must be unique, or the telnet session will fail.
; NOTE: a file called SPACEOBS.bary.eph will be generated.
; EARTHOBS - A string input to OBSERVATORY.PRO that specifies an
; observatory. If set, LAT/LON/ELEVATION are ignored.
; LAT/LON - The latitude and longitude of the observatory, in
; decimal degrees. Positive longitude is *East*. If
; LAT/LON and OBSNAME are not specified, the geocenter
; will be assumed and the resultant accuracy will be
; ~20 ms. If OBSNAME is set, these are ignored.
; ELEVATION - The elevation of the observatory, in meters. If
; LAT/LON are not set, or if OBSNAME is set, this will
; be ignored. Sea level is assumed if not given.
; TARGET - A string input to HORIZONS specifying the Solar
; System target of observation. Can be anything for
; which HORIZONS can generate an ephemeris. If set,
; RA/DEC are ignored.
; The name must be unique, or the telnet session will fail.
; NOTE: a file called TARGET.bary.eph will be generated.
; PATH - The path to the ephemeris files generated by
; HORIZONS. Current directory assumed if not specified.
; TBASE - A scalar or n_elements(JD) vector that is the
; baseline already subtracted from the input JDs. For
; highest precision, TBASE = floor(jd) is recommended.
; BIPMFILE - The name of a file containing the BIPM-TAI
; offsets from ftp://tai.bipm.org/TFG/TT%28BIPM%29/
; Note that the previous year's file should be
; downloaded, and the .ext file should be concatenated
; at the end. As of 2010, this would be the TTBIPM.09
; and TTBIPM09.ext files. This is only necessary for
; 30 us precision.
; DISTANCE - The distance from the observer to the target, in
; Parsecs. If given, this will apply the second order
; term to the roemer delay which is accurate to at
; least 10^-14 seconds for any object outside the Solar
; System. The plane parallel approximation is accurate
; to at least 1 ms for any object outside the SS.
;
; OPTIONAL KEYWORDS:
; B1950 - If set, the input coordinates are assumed to be in equinox
; B1950 coordinates. Default is the J2000 equinox.
; TIME_DIFF - To return the difference in seconds instead (either
; to have the offset or for higher precision), use this
; keyword.
; TDB = jd_utc + utc2bjd(jd_utc,ra,dec,/time_diff)/86400.d0
; TT_IN - Set this keyword to use JD_TT as an input
; instead. This will skip the check for leap second
; updates and assume you have done this already
; USE_EOP - Use EOP data for high precision nutation angles. This
; is required for times more precise than 1 us. If set,
; you must have Craig Markwardt's EOPDATA, retrieve
; this file:
; ftp://maia.usno.navy.mil/ser7/finals.all
; and save it as $ASTRO_DATA/iers_final_a.dat
; This file must be kept up-to-date.
;
; OUTPUTS:
; BJD_TDB - The light travel time corrected, Barycentric
; Dynamical Time (TDB) in BJD for each given UTC. If
; /TIME_DIFF is set, then returns the time difference
; in seconds.
; TDB = utc2bjd(jd_utc,ra,dec)
; BJD_TDB-JD_UTC = utc2bjd(jd_utc,ra,dec,/time_diff)
;
; DESCRIPTION:
; TDB time is the current time standard (as of 2010). This routine
; follows the procedure here:
; http://lheawww.gsfc.nasa.gov/users/craigm/bary/
; but ignores the dispersion correction. BJD_TDB is defined:
; BJD_TDB = TOBS + GEOMETRIC + CLOCK + EINSTEIN - SHAPIRO
;
; where
;
; TOBS - the observed JD in UTC of the event on Earth
; GEOMETRIC - the light travel time from your position on earth to
; the barycenter of the solar system (~300 s)
; CLOCK - the correction from input time to TDB time
; EINSTEIN - the relativistic correction due to using the earth
; as your inertial frame (~1 ms)
; SHAPIRO - the time delay due to photon bending in the
; potential of the solar system (~1 us)
;
; *** PLEASE READ THE FOLLOWING NOTES ABOUT ACCURACY CAREFULLY ***
;
; NOTES (in order of accuracy of result):
;
; This routine transparently handles the conversion from UTC to TDB.
; Using TT will produce a systematic offset equal to the number of
; leapseconds required plus 32.184 seconds, unless the TT_IN keyword
; is set (66.184 seconds in 2009).
;
; HELIO_JD.PRO, BARYCEN.PRO, and many routines that calculate the
; position of astronomical objects assume the input "JD" is "JD_TT",
; a critical assumption that is not obvious to users unfamiliar with
; the complexities of precision timing. This has likely led to the
; fact that the term "BJD" has been used in the literature to mean
; both "BJD_UTC" and "BJD_TT" (similarly for "HJD"), which differ by
; ~1 minute! BJD_UTC and HJD_UTC are not continous, uniformly
; increasing timescales and should never be used in astronomy.
;
; Leap seconds are added unpredictably every ~1 year. A current file
; is vital for ~1 second accuracy. This program will wget an updated
; file the first time it runs after every January 1st and July 1st
; (when leap seconds are added). Therefore wget and an internet
; connection are required to use JD_UTC as an input. If wget fails,
; you can manually retrieve it from
; ftp://maia.usno.navy.mil/ser7/tai-utc.dat and put it in
; $ASTRO_DATA. You'll also have to update
; $ASTRO_DATA/exofast_lastupdate to have the current JD_UTC so this
; program knows the leap seconds are up to date. To skip this step
; altogether, input JD_TT and set the TT_IN keyword.
;
; UTC and UT1 may differ by as much as 0.9s (the difference between
; "delta t" and leap seconds + 32.184), and both may be called
; UT. UTC is returned by NTP servers, and is usually the value
; recorded in image headers, but this should be verified if UT is
; specified.
;
; For it to take into account your position on the earth, LAT and
; LON, or OBSNAME must be specified. If none of these are set, the
; assumed position is the center of the earth, and the results will
; be systematically offset by 8-21.3 ms (observation bias). This
; correction agrees with naive estimates (spherical earth at noon on
; the winter solstice at -23.5 degrees latitude and 0 longitude) to
; 2 ms. Based on this and a handful of similar tests, this program
; is believed to accurately implement Craig Marquardt's routines,
; which are reported to be limited by the GPS coordinates of your
; observatory.
;
; For 1 ms accuracy, double precision cannot hold the full
; JD. TBASE=floor(jd) should be used.
;
; For accuracies to 30 us, TT(BIPM) must be used instead of TT(TAI)
; = UTC + 32.184 + N. This is calculated ~1 month after the fact and
; is available here: ftp://tai.bipm.org/TFG/TT(BIPM)
; /TT_IN should be specified.
;
; ELEVATION or OBSNAME should be set to account for the light travel
; time to your elevation. If not set, the elevation is assumed to be
; at sea level and results will be accurate to ~10 us.
;
; USE_EOP data must be set (and the file must be kept up to date)
; for accuracies beyond 1 us.
;
; ****************************************************************
; ****The accuracy beyond this has not been thoroughly tested.****
; ***any accuracy beyond 1 us should be independently verified.***
; ****************************************************************
;
; For ~100 ns accuracy, days cannot hold the precision. The
; correction should be returned in seconds using the /TIME_DIFF
; keyword and applied carefully.
;
; For better than ~30 ns, the corrections are limited by the
; imperfect knowledge of your location with respect to the center of
; the earth (~10 m). A slight modification of this code to use ITRF
; measurements of your position on earth (ie from VLBI) will be
; required. Further verification of this routine at that level is
; strongly recommended, and heroic measures to ensure double
; precision is adequate may be required.
;
; Also note the TT to TDB correction uses the 791-term Fairhead &
; Bretagnon analytical approximation which has a maximum error of 23
; ns and rms error of 14 ns in the time range 1980-2000, compared to
; a numerical integration.
;
; DEPENDENCIES:
; IDL astronomy library
;
; wget and an internet connection (to update the leap second file)
; alternatively, the $ASTRO_DATA/tai-utc.dat can be manually updated
; and $ASTRO_DATA/exofast_lastupdate can be manually edited to
; contain the JD_UTC of the manual update. Note this only happens
; the first time it is run after every Jan 1st and Jul 1st.
;
; an environment variable "ASTRO_DATA", which specifies that path to
; JPLEPH.405, which can be found here:
; http://www.physics.wisc.edu/~craigm/idl/down/JPLEPH.405
; NOTE: tai-utc.dat and exofast_lastupdate will be placed in
; $ASTRO_DATA too.
; For manual updates, tai-utc.dat can be found here:
; ftp://maia.usno.navy.mil/ser7/tai-utc.dat
;
; If making use of HORIZONS ephemerides (for spacecraft observatories
; or SS targets), you need Expect, telnet, and horizons.exp in your path.
;
; Craig Markwardt's routines:
; http://www.physics.wisc.edu/~craigm/idl/down/
; TAI_UTC
; (for earth position correction)
; HPRSTATN HPRNUTANG QTEULER QTCOMPOSE QTMULT QTVROT EOPDATA
;
; REVISION HISTORY:
; 2011/07/18: Fixed bug in Einstein correction (2 us level);
; Tim Lister, LCOGT, J. Eastman, OSU.
; 2010/07/19:
; Now returns a scalar instead of a 1 element array when an
; observatory or coordinates are given with a scalar date
; 2010/05/12:
; Fixed bug in geodetic calculations. Added USE_EOP keyword. Now
; accurate to < 1 us for anywhere on Earth. Added PATH keyword.
; 2010/04/27:
; Automatic generation of ephemerides for non-earth-based observers/targets
; NOTE: Requires "Expect" and "horizons.exp" in your path.
; Computes spherical wavefront correction for SS targets
; 2010/04/08: Added TT_IN keyword
; 2009/12/10: Written by Jason Eastman (OSU)
; Based on BARYCEN.PRO by goehler -- major differences are:
; uses JD_UTC as input (not MJD_TT). Automatically converts to
; TT(TAI), and can convert to TT(BIPM) if given the table of
; BIPM-TAI offsets (see notes).
; Indexes JPL ephemeris with JD_TDB, not JD_TT
; automatically wgets of leap second file required to convert to TT
; allows user settable TBASE for higher precision
; uses your position on earth for higher precision (see notes)
function utc2bjd, jd_utc, ra, dec, B1950=b1950, DISTANCE=distance,$
SPACEOBS=spaceobs, EARTHOBS=earthobs, TARGET=target, $
LAT=lat, LON=lon, ELEVATION=elevation,STEPSIZE=stepsize,$
TBASE=tbase, TIME_DIFF=time_diff,TT_IN=TT_IN,$
bipmfile=bipmfile,USE_EOP=USE_EOP,PATH=path
;; speed of light in AU/sec
c = 0.00200398880422056596d0
au = 149597870.691d0 ;km/AU
if n_elements(path) eq 0 then path = '.'
if n_elements(stepsize) eq 0 then stepsize = 10
ntimes = n_elements(jd_utc)
if n_elements(tbase) eq 0 then tbase = 0
;; convert to UTC to TDB (same as JPL ephemeris time)
jd_tdb = jdutc2jdtdb(jd_utc, TT_IN=TT_IN, tbase=tbase, $
jd_tt=jd_tt, clock_corr=clock_corr,bipmfile=bipmfile)
;; --------------------------------------------------------
;; read the ephemeris file with +/- 1 day margin
;; --------------------------------------------------------
ephemfile = find_with_def('JPLEPH.405','ASTRO_DATA')
IF NOT file_test(ephemfile) THEN $
message,"Error: JPL Ephemeris file not found"
mindate = min(jd_tdb,max=maxdate) - 1.d0
maxdate += 1.d0
JPLEPHREAD,ephemfile,pinfo,pdata,[mindate,maxdate]+tbase, $
status=status, errmsg=errmsg
IF status EQ 0 THEN message,"Ephemeris file reading failed: " + errmsg
;; get the velocity of earth (required for the Einstein delay)
JPLEPHINTERP,pinfo,pdata,jd_tdb,x_earth,y_earth,z_earth,$
vx_earth, vy_earth,vz_earth,/earth, $
posunits="AU", tbase=tbase,velunits='AU/DAY',/velocity
;; Get the coordinates of the observatory
if n_elements(spaceobs) ne 0 then begin
;; space observatory
spacendx = where(pinfo.objname eq spaceobs)
if spacendx[0] ne -1 then begin
JPLEPHINTERP,pinfo,pdata,jd_tdb,x_obs,y_obs,z_obs,$
OBJECTNAME=spaceobs, posunits="AU", tbase=tbase
endif else begin
baryfile = path + '/bary.' + spaceobs + '.eph'
if not file_test(baryfile) then $
get_eph, jd_tdb+tbase, spaceobs, outfile=baryfile, stepsize=stepsize
;; read in the CSV HORIZONS emphemeris file
readcol, baryfile, jdlist, junk, xlist, ylist, zlist, $
format='a,a,d,d,d',delimiter=',',/silent
;; don't lose precision
base = double(strmid(jdlist,0,7)) - tbase
decimal = double(strmid(jdlist,7))
jdlist = base + decimal
;; make sure the ephemeris covers the time range
if min(jdlist) gt min(jd_tdb) or max(jdlist) lt max(jd_tdb) then $
message,'ERROR: ephemeris does not cover time range -- check ' $
+ baryfile
;; interpolate the target positions at each JD
x_obs = interpol(xlist, jdlist, jd_tdb, /quadratic)
y_obs = interpol(ylist, jdlist, jd_tdb, /quadratic)
z_obs = interpol(zlist, jdlist, jd_tdb, /quadratic)
endelse
endif else begin
;; geocenter (accurate to 21.3 ms)
x_obs = x_earth
y_obs = y_earth
z_obs = z_earth
;; correct for latitude/longitude
if n_elements(earthobs) ne 0 then begin
;; observing station is on earth and observatory is specified
observatory, earthobs, obs
lat1 = obs.latitude*!dpi/180.d0
lon1 = (360.d0-obs.longitude)*!dpi/180.d0 ;; EAST!
elevation = obs.altitude
endif else begin
;; observing station is on earth and lat/lon/elevation
if n_elements(lat) eq 1 and n_elements(lon) eq 1 then begin
lat1 = lat*!dpi/180.d0
lon1 = lon*!dpi/180.d0
if n_elements(elevation) eq 0 then elevation = 0
endif
endelse
;; calculate the cartesian coordinates of the observatory
;; described by the standard (WGS84) lat/lon
if n_elements(lat1) ne 0 then begin
;; formula proposed in HPRSTATN isn't accurate (CC = 1/C)!
F = 1.d0/298.257223563d0
CC = 1.d0/SQRT(COS(LAT1)^2 + (1.d0 - F)^2*SIN(LAT1)^2)
S = (1.d0 - F)^2 * CC
H = elevation/1000.d0 ; km
A = 6378.137d0 ; km
R_ITRF = [(A*CC + H)*COS(LAT1)*COS(LON1),$
(A*CC + H)*COS(LAT1)*SIN(LON1),$
(A*S + H)*SIN(LAT1)]
HPRSTATN, jd_tt, R_ITRF, R_ECI, V_ECI, TBASE=TBASE,USE_EOP=USE_EOP
sz = size(R_ECI)
if sz[2] eq 1 then begin
x_obs += R_ECI[0,0]/au
y_obs += R_ECI[1,0]/au
z_obs += R_ECI[2,0]/au
endif else begin
x_obs += R_ECI[0,*]/au
y_obs += R_ECI[1,*]/au
z_obs += R_ECI[2,*]/au
endelse
endif
endelse
;; calculate the Einstein delay relative to the geocenter
;; (TDB accounts for Einstein delay to geocenter)
einstein_corr = ((x_obs-x_earth)*vx_earth + $
(y_obs-y_earth)*vy_earth + $
(z_obs-z_earth)*vz_earth)/c^2/86400.d0
;; Get the coordinates of what we're looking at
if n_elements(target) ne 0 then begin
;; If a target has been specified, see if it's in the ephemeris
targetndx = where(pinfo.objname eq target)
if targetndx[0] ne -1 then begin
JPLEPHINTERP,pinfo,pdata,jd_tdb,x_obj,y_obj,z_obj,$
OBJECTNAME=target, posunits="AU", tbase=tbase
endif else begin
;; otherwise, telnet to HORIZONS and generate it
tarfile = path + '/bary.' + target + '.eph'
if not file_test(tarfile) then $
get_eph, jd_tdb+tbase, target, outfile=tarfile, stepsize=stepsize
;; read in the CSV HORIZONS emphemeris file
readcol, tarfile, jdlist, junk, xlist, ylist, zlist, $
format='a,a,d,d,d',delimiter=',',/silent
;; don't lose precision
base = double(strmid(jdlist,0,7)) - tbase
decimal = double(strmid(jdlist,7))
jdlist = base + decimal
;; make sure the ephemeris covers the time range
if min(jdlist) gt min(jd_tdb) or max(jdlist) lt max(jd_tdb) then $
message,'ERROR: ephemeris does not cover time range -- check ' $
+ baryfile
;; interpolate the target positions at each JD
x_obj = interpol(xlist, jdlist, jd_tdb, /quadratic)
y_obj = interpol(ylist, jdlist, jd_tdb, /quadratic)
z_obj = interpol(zlist, jdlist, jd_tdb, /quadratic)
endelse
;; get the geometric correction
d_bary2obj = sqrt(x_obj^2 + y_obj^2 + z_obj^2)
d_obs2obj = sqrt((x_obj-x_obs)^2+(y_obj-y_obs)^2+(z_obj-z_obs)^2)
x_hat = x_obj/d_obs2obj
y_hat = y_obj/d_obs2obj
z_hat = z_obj/d_obs2obj
geo_corr = (d_bary2obj - d_obs2obj)/c
endif else begin
;; precess the object coordinates if desired
if keyword_set(B1950) then jprecess,ra,dec,ra1,dec1 else begin
ra1 = ra
dec1 = dec
endelse
;; convert coordinates to radians:
ra1 = ra1*!dpi/180.d0
dec1 = dec1*!dpi/180.d0
x_hat = cos(dec1)*cos(ra1)
y_hat = cos(dec1)*sin(ra1)
z_hat = sin(dec1)
geo_corr = (x_obs*x_hat + y_obs*y_hat + z_obs*z_hat)/c
;; second order correction given the distance to the object
;; good for 10^(-14) seconds for any object outside the SS.
if n_elements(distance) eq 1 then begin
pc = 180.d0*3600.d0/!dpi ;; AU/pc
geo_corr += ((x_obs^2 + y_obs^2 + z_obs^2) - $
(x_obs*x_hat + y_obs*y_hat + z_obs*z_hat)^2)/$
(2*c*distance*pc)
endif
endelse
;; ------------------------------------------------------------
;; COMPUTE SHAPIRO CORRECTION (~20 us)
;; ------------------------------------------------------------
JPLEPHINTERP,pinfo,pdata,jd_tdb,x_sun,y_sun,z_sun,/sun,$
posunits="AU", tbase=tbase
;; distance of sun to observatory:
sun_dist = sqrt((x_sun - x_obs)^2 + (y_sun - y_obs)^2 + (z_sun - z_obs)^2)
;; cosine of angle between sun and observatory
costh = ((x_sun-x_obs)*x_hat+(y_sun-y_obs)*y_hat+(z_sun-z_obs)*z_hat)/sun_dist
const = 9.850981896617677d-6 ; 2*G*M_sun/c^3 (seconds)
shapiro_corr = const*alog(1-costh)
;; ------------------------------------------------------------
;; SUMMARIZE CORRECTIONS
;; ------------------------------------------------------------
delta_t = geo_corr + clock_corr + einstein_corr + shapiro_corr
;; return the desired value
if keyword_set(time_diff) then return, delta_t
return, jd_utc + delta_t/86400.d0
end