;+ ; 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