;+ ; NAME: ; EOPDATA ; ; 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: ; Read and interpolate tabulated earth orientation parameters ; ; MAJOR TOPICS: ; Geometry ; ; CALLING SEQUENCE: ; EOPDATA, JDUTC, PMX, PMY, UT1_UTC, DPSI, DEPS, $ ; /RESET, FILENAME=, ANGUNITS=, TBASE= ; ; DESCRIPTION: ; ; The procedure EOPDATA reads, interpolates and returns Earth ; orientation parameters used for precision earth-base astronomy ; applications. ; ; ** NOTE: The user is responsible for downloading and maintaining an ; up-to-date file of earth orientation parameters from the ; International Earth Rotation Service. See below. ** ; ; This interface is somewhat provisional. See OPEN QUESTIONS below. ; ; The values returned are described below. These descriptions are ; taken from the Explanatory Supplement to IERS Bulletins A and B. ; ; * PMX and PMY, the coordinates of the Celestial Ephemeris Pole ; (CEP) relative to the earth-fixed International Reference Pole ; (IRP). The x-axis is in the direction of the IERS Reference ; Meridian (IRM), the y-axis is in the direction 90 degrees West ; longitude. The time series of PMX and PMY is referred to as ; "polar motion." ; ; These are the coordinates of the earth rotation pole, as seen ; in an *earth-fixed* coordinate system. A station whose ; coordinates are given in earth-fixed coordinates referred to ; the ITRS can be transformed to the earth-fixed coordinates ; referred to the true rotation pole of date using the following ; matrix transformation: ; ; R_TRUE = RX(PMY) ## RY(PMX) ## R_ITRS ; ; where the matrices RX and RY are defined below. ; ; * UT1, the the rotation angle about the pole. It is related to ; the Greenwich mean sidereal time (GMST) by a conventional ; relationship (Aoki et al., 1982). It gives access to the ; direction of the International Reference Meridian IRM in the ; ICRS, reckoned around the CEP axis. It is expressed as the ; difference UT1-UTC. Thus, the value of UT1 is computed as: ; ; UT1 = UT1_UTC + UTC ; ; where UTC is the UTC time, expressed in seconds. ; ; * DPSI and DEPS are the offsets in longitude and obliquity of the ; celestial pole with respect to its direction defined using the ; conventional IAU precession/nutation theory. An a priori ; correction model is available in the IERS Conventions (1996), ; (McCarthy, 1996). The expressions to compute the nutation ; angles are: ; ; DEPS_TRUE = DEPS_1980 + DEPS ;; Nutation in obliquity ; DPSI_TRUE = DPSI_1980 + DPSI ;; Nutation in longitude ; ; where DPSI_1980 and DEPS_1980 are the nutation values ; determined from the IAU 1980 Nutation Theory; and DPSI_TRUE and ; DEPS_TRUE are the nutations to be used as arguments to further ; precession and nutation computations. ; ; For requested times which are between tabular values, a linear ; interpolation is performed. This is not exactly the correct ; procedure, and can result in errors of +/- 0.1 mas in the earth ; polar motion and 1 usec in UT1 (see McCarthy & Gambis 1997). ; ; ; DATA FILES and MAINTENANCE ; ; The user is responsible for downloading and maintaining the earth ; orientation parameters file as supplied by the IERS. The format ; of the files is the "Final" EOP data ASCII format. They can be ; downloaded here: ; ; ftp://maia.usno.navy.mil/ser7/finals.all ;; from May 1976-present ; ftp://maia.usno.navy.mil/ser7/finals.data ;; from Jan 1992-present ; ; The user must place this file in a known location, and in *at ; least the first call*, this filename must be passed using the ; FILENAME keyword. ; ; EOPDATA will load the data once on the first call, and keep a ; cached copy for subsequent calls. On a daily basis the file will ; be reloaded in case the quantities have been updated from the ; server. A reload of data can be forced using the RESET keyword. ; ; ROTATION MATRICES ; ; The rotation matrices RX(T) and RY(T) mentioned above in relation ; to polar motion are: ; ; RX(T) =EQ= [[1,0,0], [0,cos(T),sin(T)], [0,-sin(T),cos(T)]] ; RY(T) =EQ= [[cos(T),0,-sin(T)], [0,1,0], [sin(T),0,cos(T)]] ; RZ(T) =EQ= [[cos(T),sin(T),0], [-sin(T),cos(T),0], [0,0,1]] ; ; and are meant to be applied to a vector R as, RX(T) ## R. ; ; ; OPEN QUESTIONS ; ; How will the transition to a new IERS EOP series be accomplished? ; Using a keyword? ; ; Should there be a quality flag? The EOP file contains a ; "predicted" flag, and also there are rows which contain no value ; at all. These should probably be flagged somehow. ; ; ; INPUTS: ; ; JDUTC - a vector or scalar, the UTC time for which earth ; orientation parameters are requested, expressed in Julian ; Days. The value of the keyword TBASE is added to this ; quantity to arrive at the actual Julian date. ; ; OUTPUTS: ; ; PMX, PMY - the earth-fixed angular coordinates of the celestial ; ephemeris pole, measured in ANGUNITS units. ; ; UT1_UTC - the value of UT1 - UTC, expressed in seconds. ; ; DPSI, DEPS - the corrections to the IAU 1980 theory of Nutation, ; for nutation in longitude and obliquity, expressed in ; ANGUNITS units. ; ; KEYWORD PARAMETERS: ; ; FILENAME - scalar string, on the first call, the name of the file ; from which earth orientation parameters will be read. ; Default value: (none) ; ; TBASE - a fixed epoch time (Julian days) to be added to each value ; of JDUTC. Since subtraction of large numbers occurs with ; TBASE first, the greatest precision is achieved when TBASE ; is expressed as a nearby julian epoch, JDUTC is expressed ; as a small offset from the fixed epoch. ; Default: 0 ; ; ANGUNITS - scalar string, output units of angular parameters. ; Possible values are 'ARCSEC' or 'RADIAN'. ; Default value: 'RADIAN' ; ; RESET - if set, forces EOP file to be re-read. ; ; ; EXAMPLE: ; ; ; SEE ALSO: ; ; HPRNUTANG, TAI_UTC (Markwardt Library) ; PRECESS, PREMAT, JPRECESS, BPRECESS (IDL Astronomy Library) ; ; ; REFERENCES: ; ; Aoki, S., Guinot, B., Kaplan, G.H., Kinoshita, H., McCarthy, D.D., ; Seidelmann, P.K., 1982: Astron. Astrophys., 105, 359-361. ; ; McCarthy, D. D. (ed.) 1996: IERS Conventions, IERS T.N. 21. ; http://maia.usno.navy.mil/conventions.html ; ; McCarthy, D. \& Gambis, D. 1997, "Interpolating the IERS Earth ; Orientation Data," IERS Gazette No. 13, ; http://maia.usno.navy.mil/iers-gaz13 ; Instructions for high precision EOP data interpolation, not done ; in this procedure. ; ; Ray, J. & Gambis, D. 2001, "Explanatory Supplement to IERS ; Bulletins A and B," ; http://hpiers.obspm.fr/iers/bul/bulb/explanatory.html ; ; Explains meanings of earth orientation parameters used and ; returned by this procedure. ; ; Definition of Final EOP data format ; ftp://maia.usno.navy.mil/ser7/readme.finals ; ; MODIFICATION HISTORY: ; Written, 30 Jan 2002, CM ; Documented, 14 Feb 2002, CM ; Add default message, 01 Mar 2002, CM ; More robust handling of input file, 10 Mar 2002, CM ; Fix bug in interpolation of nutation and polar motion adjustments, ; thanks to Tim Lister, 2014-10-09, CM ; ; $Id: eopdata.pro,v 1.7 2014/10/20 21:36:16 cmarkwar Exp $ ; ;- ; Copyright (C) 2002, 2014, Craig Markwardt ; This software is provided as is without any warranty whatsoever. ; Permission to use, copy, modify, and distribute modified or ; unmodified copies is granted, provided this copyright and disclaimer ; are included unchanged. ;- pro eopdata_read, file, jd, pmx, pmy, ut1, dpsi, deps, status=status status = 0 get_lun, unit openr, unit, file, error=err if err NE 0 then begin free_lun, unit message, 'ERROR: could not open '+file endif buffersize = 2048 ngood = 0L jd = 0d & pmx = 0D & pmy = 0D & ut1 = 0D & dpsi = 0d & deps = 0d ss = strarr(buffersize) while NOT eof(unit) do begin ;; Read data from disk on_ioerror, TRIM readf, unit, ss ;; On the last pass we may get less than a full buffer's worth TRIM: cc = (fstat(unit)).transfer_count if cc EQ 0 then goto, DONE if cc LT buffersize then ss = ss(0:cc-1) ;; Parse the parameters mjd1 = double(strcompress(strmid(ss,7,9),/remove_all)) NEXT0: on_ioerror, NEXT1 pmx1 = double(strcompress(strmid(ss,18,10),/remove_all)) NEXT1: on_ioerror, NEXT2 pmy1 = double(strcompress(strmid(ss,37,10),/remove_all)) NEXT2: on_ioerror, NEXT3 ut11 = double(strcompress(strmid(ss,58,11),/remove_all)) NEXT3: on_ioerror, NEXT4 dps1 = double(strcompress(strmid(ss,97,11),/remove_all)) NEXT4: on_ioerror, NEXT5 dep1 = double(strcompress(strmid(ss,116,11),/remove_all)) NEXT5: jd = [jd, mjd1] pmx = [pmx, pmx1] & pmy = [pmy, pmy1] ut1 = [ut1, ut11] dpsi = [dpsi, dps1] deps = [deps, dep1] ss(*) = '' endwhile DONE: free_lun, unit if n_elements(jd) EQ 1 then begin message, 'ERROR: could not read data from '+file endif AS2R = !dpi/180d/3600d ;; Arcsec to radians jd = jd(1:*) pmx = pmx(1:*) * AS2R ;; Convert to radians pmy = pmy(1:*) * AS2R ;; Convert to radians ut1 = ut1(1:*) ;; Already in seconds dpsi = dpsi(1:*) * (0.001D * AS2R) ;; Convert from mas to radians deps = deps(1:*) * (0.001D * AS2R) ;; Convert from mas to radians status = 1 return end pro eopdata, jdutc, pmx, pmy, ut1_utc, dpsi, deps, reset=reset, $ filename=filename0, angunits=angunits0, tbase=tbase0 common eopdata_table, mjd0, pmx0, pmy0, ut10, dpsi0, deps0, $ ntable, mjdmin, mjdmax, leap0, $ timestamp, oldfile if n_params() EQ 0 AND n_elements(filename0) EQ 0 then begin USAGE: message, 'USAGE:', /info message, 'EOPDATA, JDUTC, PMX, PMY, UT1_UTC, DPSI, DEPS, '+$ '[FILENAME=filename, ANGUNITS=angunits, TBASE=tbase', /info message, " ANGUNITS is one of 'ARCSEC' or 'RADIAN'", /info return endif if n_elements(mjd0) EQ 0 OR keyword_set(reset) OR $ (n_elements(jdutc) EQ 0 AND n_elements(filename0) GT 0) then begin RELOAD_COMMON: forward_function get_xtecal ;; Find filename if n_elements(filename0) EQ 0 then begin filename = '' ;; First try: the old file if n_elements(oldfile) NE 0 then $ filename = oldfile ;; Second try: use standard IDL Astronomy Library data directory if filename EQ '' then $ filename = find_with_def('iers_final_a.dat','ASTRO_DATA') ;; Third try: Markwardt-specific if filename EQ '' then begin catch, catcherr if catcherr EQ 0 then $ filename = get_xtecal() + 'clock/finals.data' catch, /cancel endif ;; Could not find it, so trigger a fatal error if filename EQ '' then $ message, 'ERROR: Could not find EOP data file' endif else begin filename = strtrim(filename0(0),2) endelse eopdata_read, filename, mjd0, pmx0, pmy0, ut10, dpsi0, deps0, $ status=status ;; Make a series of pseudo leap seconds so that we can ;; interpolate linearly below. wh = where(abs(ut10(1:*) - ut10) GT 0.8 AND ut10 NE 0, ct) leap0 = lonarr(n_elements(ut10)) if ct GT 0 then begin wh = [wh, n_elements(ut10)-1] & ct = ct + 1 for i = 0, ct-2 do $ leap0(wh(i)+1:wh(i+1)) = (leap0(wh(i)) $ - round(ut10(wh(i)+1)-ut10(wh(i)))) endif if status NE 1 then begin message, 'ERROR: could not read EOP data from '+filename endif timestamp = systime(1) oldfile = filename if n_elements(jdutc) EQ 0 then return endif if n_params() EQ 0 then goto, USAGE if systime(1) - timestamp GT 86400d then goto, RELOAD_COMMON if n_elements(tbase0) EQ 0 then $ tbase = 0d $ else $ tbase = double(tbase0) if n_elements(angunits0) EQ 0 then $ angunits = 'RADIAN' $ else $ angunits = strtrim(strupcase(strcompress(angunits0(0))),2) ;; Convert from Julian days to MJD t = (jdutc(*) + (tbase - 2400000.5d)) tmin = min(t, max=tmax) wh = where(mjd0 GE tmin-5 AND mjd0 LE tmax+5 AND ut10 NE 0, ct) if ct EQ 0 then begin OUT_OF_BOUNDS: message, 'ERROR: input time was out of bounds' endif mjd1 = mjd0(wh) ii = value_locate(mjd1, t) jj = wh(ii) ;; This is UT1 - UTC(TSTART), so it should be continuous. Thus, we ;; can interpolate it. The normal UT1-UTC series has discontinuities. ut1 = (ut10+leap0)(jj) ;; Linear interpolation dt = (t-mjd1(ii))/(mjd1(ii+1)-mjd1(ii)) ut1_utc = ut10(jj) + dt * (ut1(ii+1)-ut1(ii)) ;; Interpolate DPSI and DEPS, the adjustments to the nutation angles dpsi = dpsi0(jj) + dt * (dpsi0(jj+1)-dpsi0(jj)) deps = deps0(jj) + dt * (deps0(jj+1)-deps0(jj)) ;; Polar motion parameters pmx = pmx0(jj) + dt * (pmx0(jj+1)-pmx0(jj)) pmy = pmy0(jj) + dt * (pmy0(jj+1)-pmy0(jj)) ;; Units conversions case angunits of 'ARCSEC': begin R2AS = 3600d*180d/!dpi ;; Radian to arcsec dpsi = dpsi * R2AS deps = deps * R2AS pmx = pmx * R2AS pmy = pmy * R2AS end 'RADIAN': begin dummy = 1 end else: begin message, 'ERROR: angular unit '+angunits+$ ' was not recognized' end end return end