Index: trunk/psLib/src/astronomy/Makefile.am =================================================================== --- trunk/psLib/src/astronomy/Makefile.am (revision 3114) +++ trunk/psLib/src/astronomy/Makefile.am (revision 3114) @@ -0,0 +1,38 @@ +#Makefile for astronomy functions of psLib +# +INCLUDES = \ + -I$(top_srcdir)/src/collections \ + -I$(top_srcdir)/src/dataManip \ + -I$(top_srcdir)/src/fileUtils \ + -I$(top_srcdir)/src/image \ + -I$(top_srcdir)/src/sysUtils \ + $(all_includes) + +noinst_LTLIBRARIES = libpslibastronomy.la + +libpslibastronomy_la_SOURCES = \ + psTime.c \ + psMetadata.c \ + psMetadataIO.c \ + psCoord.c \ + psAstrometry.c \ + aoppa.f aopqk.f oapqk.f airmas.f eqeqx.f geoc.f refco.f aoppat.f \ + dranrm.f dcs2c.f refz.f refro.f dcc2s.f gmst.f atms.f atmt.f nutc.f drange.f + +BUILT_SOURCES = psAstronomyErrors.h + +EXTRA_DIST = psAstronomyErrors.dat psAstronomyErrors.h + +psAstronomyErrors.h: psAstronomyErrors.dat + perl $(top_srcdir)/src/parseErrorCodes.pl --data=$? $@ + +pslibincludedir = $(includedir)/pslib +pslibinclude_HEADERS = \ + psTime.h \ + psMetadata.h \ + psMetadataIO.h \ + psCoord.h \ + psAstrometry.h \ + psPhotometry.h \ + slalib.h + Index: trunk/psLib/src/astronomy/airmas.f =================================================================== --- trunk/psLib/src/astronomy/airmas.f (revision 3114) +++ trunk/psLib/src/astronomy/airmas.f (revision 3114) @@ -0,0 +1,75 @@ + DOUBLE PRECISION FUNCTION sla_AIRMAS (ZD) +*+ +* - - - - - - - +* A I R M A S +* - - - - - - - +* +* Air mass at given zenith distance (double precision) +* +* Given: +* ZD d Observed zenith distance (radians) +* +* The result is an estimate of the air mass, in units of that +* at the zenith. +* +* Notes: +* +* 1) The "observed" zenith distance referred to above means "as +* affected by refraction". +* +* 2) Uses Hardie's (1962) polynomial fit to Bemporad's data for +* the relative air mass, X, in units of thickness at the zenith +* as tabulated by Schoenberg (1929). This is adequate for all +* normal needs as it is accurate to better than 0.1% up to X = +* 6.8 and better than 1% up to X = 10. Bemporad's tabulated +* values are unlikely to be trustworthy to such accuracy +* because of variations in density, pressure and other +* conditions in the atmosphere from those assumed in his work. +* +* 3) The sign of the ZD is ignored. +* +* 4) At zenith distances greater than about ZD = 87 degrees the +* air mass is held constant to avoid arithmetic overflows. +* +* References: +* Hardie, R.H., 1962, in "Astronomical Techniques" +* ed. W.A. Hiltner, University of Chicago Press, p180. +* Schoenberg, E., 1929, Hdb. d. Ap., +* Berlin, Julius Springer, 2, 268. +* +* Original code by P.W.Hill, St Andrews +* +* P.T.Wallace Starlink 18 March 1999 +* +* Copyright (C) 1999 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION ZD + + DOUBLE PRECISION SECZM1 + + + SECZM1 = 1D0/(COS(MIN(1.52D0,ABS(ZD))))-1D0 + sla_AIRMAS = 1D0 + SECZM1*(0.9981833D0 + : - SECZM1*(0.002875D0 + 0.0008083D0*SECZM1)) + + END Index: trunk/psLib/src/astronomy/aoppa.f =================================================================== --- trunk/psLib/src/astronomy/aoppa.f (revision 3114) +++ trunk/psLib/src/astronomy/aoppa.f (revision 3114) @@ -0,0 +1,193 @@ + SUBROUTINE sla_AOPPA (DATE, DUT, ELONGM, PHIM, HM, + : XP, YP, TDK, PMB, RH, WL, TLR, AOPRMS) +*+ +* - - - - - - +* A O P P A +* - - - - - - +* +* Precompute apparent to observed place parameters required by +* sla_AOPQK and sla_OAPQK. +* +* Given: +* DATE d UTC date/time (modified Julian Date, JD-2400000.5) +* DUT d delta UT: UT1-UTC (UTC seconds) +* ELONGM d mean longitude of the observer (radians, east +ve) +* PHIM d mean geodetic latitude of the observer (radians) +* HM d observer's height above sea level (metres) +* XP d polar motion x-coordinate (radians) +* YP d polar motion y-coordinate (radians) +* TDK d local ambient temperature (DegK; std=273.15D0) +* PMB d local atmospheric pressure (mB; std=1013.25D0) +* RH d local relative humidity (in the range 0D0-1D0) +* WL d effective wavelength (micron, e.g. 0.55D0) +* TLR d tropospheric lapse rate (DegK/metre, e.g. 0.0065D0) +* +* Returned: +* AOPRMS d(14) star-independent apparent-to-observed parameters: +* +* (1) geodetic latitude (radians) +* (2,3) sine and cosine of geodetic latitude +* (4) magnitude of diurnal aberration vector +* (5) height (HM) +* (6) ambient temperature (TDK) +* (7) pressure (PMB) +* (8) relative humidity (RH) +* (9) wavelength (WL) +* (10) lapse rate (TLR) +* (11,12) refraction constants A and B (radians) +* (13) longitude + eqn of equinoxes + sidereal DUT (radians) +* (14) local apparent sidereal time (radians) +* +* Notes: +* +* 1) It is advisable to take great care with units, as even +* unlikely values of the input parameters are accepted and +* processed in accordance with the models used. +* +* 2) The DATE argument is UTC expressed as an MJD. This is, +* strictly speaking, improper, because of leap seconds. However, +* as long as the delta UT and the UTC are consistent there +* are no difficulties, except during a leap second. In this +* case, the start of the 61st second of the final minute should +* begin a new MJD day and the old pre-leap delta UT should +* continue to be used. As the 61st second completes, the MJD +* should revert to the start of the day as, simultaneously, +* the delta UTC changes by one second to its post-leap new value. +* +* 3) The delta UT (UT1-UTC) is tabulated in IERS circulars and +* elsewhere. It increases by exactly one second at the end of +* each UTC leap second, introduced in order to keep delta UT +* within +/- 0.9 seconds. +* +* 4) IMPORTANT -- TAKE CARE WITH THE LONGITUDE SIGN CONVENTION. +* The longitude required by the present routine is east-positive, +* in accordance with geographical convention (and right-handed). +* In particular, note that the longitudes returned by the +* sla_OBS routine are west-positive, following astronomical +* usage, and must be reversed in sign before use in the present +* routine. +* +* 5) The polar coordinates XP,YP can be obtained from IERS +* circulars and equivalent publications. The maximum amplitude +* is about 0.3 arcseconds. If XP,YP values are unavailable, +* use XP=YP=0D0. See page B60 of the 1988 Astronomical Almanac +* for a definition of the two angles. +* +* 6) The height above sea level of the observing station, HM, +* can be obtained from the Astronomical Almanac (Section J +* in the 1988 edition), or via the routine sla_OBS. If P, +* the pressure in millibars, is available, an adequate +* estimate of HM can be obtained from the expression +* +* HM ~ -29.3D0*TSL*LOG(P/1013.25D0). +* +* where TSL is the approximate sea-level air temperature in +* deg K (see Astrophysical Quantities, C.W.Allen, 3rd edition, +* section 52). Similarly, if the pressure P is not known, +* it can be estimated from the height of the observing +* station, HM as follows: +* +* P ~ 1013.25D0*EXP(-HM/(29.3D0*TSL)). +* +* Note, however, that the refraction is proportional to the +* pressure and that an accurate P value is important for +* precise work. +* +* 7) Repeated, computationally-expensive, calls to sla_AOPPA for +* times that are very close together can be avoided by calling +* sla_AOPPA just once and then using sla_AOPPAT for the subsequent +* times. Fresh calls to sla_AOPPA will be needed only when +* changes in the precession have grown to unacceptable levels or +* when anything affecting the refraction has changed. +* +* Called: sla_GEOC, sla_REFCO, sla_EQEQX, sla_AOPPAT +* +* P.T.Wallace Starlink 24 October 2003 +* +* Copyright (C) 2003 P.T.Wallace and CCLRC +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE,DUT,ELONGM,PHIM,HM,XP,YP,TDK,PMB, + : RH,WL,TLR,AOPRMS(14) + + DOUBLE PRECISION sla_EQEQX + +* 2Pi + DOUBLE PRECISION D2PI + PARAMETER (D2PI=6.283185307179586476925287D0) + +* Seconds of time to radians + DOUBLE PRECISION S2R + PARAMETER (S2R=7.272205216643039903848712D-5) + +* Speed of light (AU per day) + DOUBLE PRECISION C + PARAMETER (C=173.14463331D0) + +* Ratio between solar and sidereal time + DOUBLE PRECISION SOLSID + PARAMETER (SOLSID=1.00273790935D0) + + DOUBLE PRECISION CPHIM,XT,YT,ZT,XC,YC,ZC,ELONG,PHI,UAU,VAU + + + +* Observer's location corrected for polar motion + CPHIM = COS(PHIM) + XT = COS(ELONGM)*CPHIM + YT = SIN(ELONGM)*CPHIM + ZT = SIN(PHIM) + XC = XT-XP*ZT + YC = YT+YP*ZT + ZC = XP*XT-YP*YT+ZT + IF (XC.EQ.0D0.AND.YC.EQ.0D0) THEN + ELONG = 0D0 + ELSE + ELONG = ATAN2(YC,XC) + END IF + PHI = ATAN2(ZC,SQRT(XC*XC+YC*YC)) + AOPRMS(1) = PHI + AOPRMS(2) = SIN(PHI) + AOPRMS(3) = COS(PHI) + +* Magnitude of the diurnal aberration vector + CALL sla_GEOC(PHI,HM,UAU,VAU) + AOPRMS(4) = D2PI*UAU*SOLSID/C + +* Copy the refraction parameters and compute the A & B constants + AOPRMS(5) = HM + AOPRMS(6) = TDK + AOPRMS(7) = PMB + AOPRMS(8) = RH + AOPRMS(9) = WL + AOPRMS(10) = TLR + CALL sla_REFCO(HM,TDK,PMB,RH,WL,PHI,TLR,1D-10, + : AOPRMS(11),AOPRMS(12)) + +* Longitude + equation of the equinoxes + sidereal equivalent of DUT +* (ignoring change in equation of the equinoxes between UTC and TDB) + AOPRMS(13) = ELONG+sla_EQEQX(DATE)+DUT*SOLSID*S2R + +* Sidereal time + CALL sla_AOPPAT(DATE,AOPRMS) + + END Index: trunk/psLib/src/astronomy/aoppat.f =================================================================== --- trunk/psLib/src/astronomy/aoppat.f (revision 3114) +++ trunk/psLib/src/astronomy/aoppat.f (revision 3114) @@ -0,0 +1,62 @@ + SUBROUTINE sla_AOPPAT (DATE, AOPRMS) +*+ +* - - - - - - - +* A O P P A T +* - - - - - - - +* +* Recompute the sidereal time in the apparent to observed place +* star-independent parameter block. +* +* Given: +* DATE d UTC date/time (modified Julian Date, JD-2400000.5) +* (see AOPPA source for comments on leap seconds) +* +* AOPRMS d(14) star-independent apparent-to-observed parameters +* +* (1-12) not required +* (13) longitude + eqn of equinoxes + sidereal DUT +* (14) not required +* +* Returned: +* AOPRMS d(14) star-independent apparent-to-observed parameters: +* +* (1-13) not changed +* (14) local apparent sidereal time (radians) +* +* For more information, see sla_AOPPA. +* +* Called: sla_GMST +* +* P.T.Wallace Starlink 1 July 1993 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE,AOPRMS(14) + + DOUBLE PRECISION sla_GMST + + + + AOPRMS(14) = sla_GMST(DATE)+AOPRMS(13) + + END Index: trunk/psLib/src/astronomy/aopqk.f =================================================================== --- trunk/psLib/src/astronomy/aopqk.f (revision 3114) +++ trunk/psLib/src/astronomy/aopqk.f (revision 3114) @@ -0,0 +1,259 @@ + SUBROUTINE sla_AOPQK (RAP, DAP, AOPRMS, AOB, ZOB, HOB, DOB, ROB) +*+ +* - - - - - - +* A O P Q K +* - - - - - - +* +* Quick apparent to observed place (but see note 8, below, for +* remarks about speed). +* +* Given: +* RAP d geocentric apparent right ascension +* DAP d geocentric apparent declination +* AOPRMS d(14) star-independent apparent-to-observed parameters: +* +* (1) geodetic latitude (radians) +* (2,3) sine and cosine of geodetic latitude +* (4) magnitude of diurnal aberration vector +* (5) height (HM) +* (6) ambient temperature (T) +* (7) pressure (P) +* (8) relative humidity (RH) +* (9) wavelength (WL) +* (10) lapse rate (TLR) +* (11,12) refraction constants A and B (radians) +* (13) longitude + eqn of equinoxes + sidereal DUT (radians) +* (14) local apparent sidereal time (radians) +* +* Returned: +* AOB d observed azimuth (radians: N=0,E=90) +* ZOB d observed zenith distance (radians) +* HOB d observed Hour Angle (radians) +* DOB d observed Declination (radians) +* ROB d observed Right Ascension (radians) +* +* Notes: +* +* 1) This routine returns zenith distance rather than elevation +* in order to reflect the fact that no allowance is made for +* depression of the horizon. +* +* 2) The accuracy of the result is limited by the corrections for +* refraction. Providing the meteorological parameters are +* known accurately and there are no gross local effects, the +* observed RA,Dec predicted by this routine should be within +* about 0.1 arcsec for a zenith distance of less than 70 degrees. +* Even at a topocentric zenith distance of 90 degrees, the +* accuracy in elevation should be better than 1 arcmin; useful +* results are available for a further 3 degrees, beyond which +* the sla_REFRO routine returns a fixed value of the refraction. +* The complementary routines sla_AOP (or sla_AOPQK) and sla_OaAP +* (or sla_OAPQK) are self-consistent to better than 1 micro- +* arcsecond all over the celestial sphere. +* +* 3) It is advisable to take great care with units, as even +* unlikely values of the input parameters are accepted and +* processed in accordance with the models used. +* +* 4) "Apparent" place means the geocentric apparent right ascension +* and declination, which is obtained from a catalogue mean place +* by allowing for space motion, parallax, precession, nutation, +* annual aberration, and the Sun's gravitational lens effect. For +* star positions in the FK5 system (i.e. J2000), these effects can +* be applied by means of the sla_MAP etc routines. Starting from +* other mean place systems, additional transformations will be +* needed; for example, FK4 (i.e. B1950) mean places would first +* have to be converted to FK5, which can be done with the +* sla_FK425 etc routines. +* +* 5) "Observed" Az,El means the position that would be seen by a +* perfect theodolite located at the observer. This is obtained +* from the geocentric apparent RA,Dec by allowing for Earth +* orientation and diurnal aberration, rotating from equator +* to horizon coordinates, and then adjusting for refraction. +* The HA,Dec is obtained by rotating back into equatorial +* coordinates, using the geodetic latitude corrected for polar +* motion, and is the position that would be seen by a perfect +* equatorial located at the observer and with its polar axis +* aligned to the Earth's axis of rotation (n.b. not to the +* refracted pole). Finally, the RA is obtained by subtracting +* the HA from the local apparent ST. +* +* 6) To predict the required setting of a real telescope, the +* observed place produced by this routine would have to be +* adjusted for the tilt of the azimuth or polar axis of the +* mounting (with appropriate corrections for mount flexures), +* for non-perpendicularity between the mounting axes, for the +* position of the rotator axis and the pointing axis relative +* to it, for tube flexure, for gear and encoder errors, and +* finally for encoder zero points. Some telescopes would, of +* course, exhibit other properties which would need to be +* accounted for at the appropriate point in the sequence. +* +* 7) The star-independent apparent-to-observed-place parameters +* in AOPRMS may be computed by means of the sla_AOPPA routine. +* If nothing has changed significantly except the time, the +* sla_AOPPAT routine may be used to perform the requisite +* partial recomputation of AOPRMS. +* +* 8) At zenith distances beyond about 76 degrees, the need for +* special care with the corrections for refraction causes a +* marked increase in execution time. Moreover, the effect +* gets worse with increasing zenith distance. Adroit +* programming in the calling application may allow the +* problem to be reduced. Prepare an alternative AOPRMS array, +* computed for zero air-pressure; this will disable the +* refraction corrections and cause rapid execution. Using +* this AOPRMS array, a preliminary call to the present routine +* will, depending on the application, produce a rough position +* which may be enough to establish whether the full, slow +* calculation (using the real AOPRMS array) is worthwhile. +* For example, there would be no need for the full calculation +* if the preliminary call had already established that the +* source was well below the elevation limits for a particular +* telescope. +* +* 9) The azimuths etc produced by the present routine are with +* respect to the celestial pole. Corrections to the terrestrial +* pole can be computed using sla_POLMO. +* +* Called: sla_DCS2C, sla_REFZ, sla_REFRO, sla_DCC2S, sla_DRANRM +* +* P.T.Wallace Starlink 24 October 2003 +* +* Copyright (C) 2003 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION RAP,DAP,AOPRMS(14),AOB,ZOB,HOB,DOB,ROB + +* Breakpoint for fast/slow refraction algorithm: +* ZD greater than arctan(4), (see sla_REFCO routine) +* or vector Z less than cosine(arctan(Z)) = 1/sqrt(17) + DOUBLE PRECISION ZBREAK + PARAMETER (ZBREAK=0.242535625D0) + + INTEGER I + + DOUBLE PRECISION SPHI,CPHI,ST,V(3),XHD,YHD,ZHD,DIURAB,F, + : XHDT,YHDT,ZHDT,XAET,YAET,ZAET,AZOBS, + : ZDT,REFA,REFB,ZDOBS,DZD,DREF,CE, + : XAEO,YAEO,ZAEO,HMOBS,DCOBS,RAOBS + + DOUBLE PRECISION sla_DRANRM + + + +* Sin, cos of latitude + SPHI = AOPRMS(2) + CPHI = AOPRMS(3) + +* Local apparent sidereal time + ST = AOPRMS(14) + +* Apparent RA,Dec to Cartesian -HA,Dec + CALL sla_DCS2C(RAP-ST,DAP,V) + XHD = V(1) + YHD = V(2) + ZHD = V(3) + +* Diurnal aberration + DIURAB = AOPRMS(4) + F = (1D0-DIURAB*YHD) + XHDT = F*XHD + YHDT = F*(YHD+DIURAB) + ZHDT = F*ZHD + +* Cartesian -HA,Dec to Cartesian Az,El (S=0,E=90) + XAET = SPHI*XHDT-CPHI*ZHDT + YAET = YHDT + ZAET = CPHI*XHDT+SPHI*ZHDT + +* Azimuth (N=0,E=90) + IF (XAET.EQ.0D0.AND.YAET.EQ.0D0) THEN + AZOBS = 0D0 + ELSE + AZOBS = ATAN2(YAET,-XAET) + END IF + +* Topocentric zenith distance + ZDT = ATAN2(SQRT(XAET*XAET+YAET*YAET),ZAET) + +* +* Refraction +* ---------- + +* Fast algorithm using two constant model + REFA = AOPRMS(11) + REFB = AOPRMS(12) + CALL sla_REFZ(ZDT,REFA,REFB,ZDOBS) + +* Large zenith distance? + IF (COS(ZDOBS).LT.ZBREAK) THEN + +* Yes: use rigorous algorithm + +* Initialize loop (maximum of 10 iterations) + I = 1 + DZD = 1D1 + DO WHILE (ABS(DZD).GT.1D-10.AND.I.LE.10) + +* Compute refraction using current estimate of observed ZD + CALL sla_REFRO(ZDOBS,AOPRMS(5),AOPRMS(6),AOPRMS(7), + : AOPRMS(8),AOPRMS(9),AOPRMS(1), + : AOPRMS(10),1D-8,DREF) + +* Remaining discrepancy + DZD = ZDOBS+DREF-ZDT + +* Update the estimate + ZDOBS = ZDOBS-DZD + +* Increment the iteration counter + I = I+1 + END DO + END IF + +* To Cartesian Az/ZD + CE = SIN(ZDOBS) + XAEO = -COS(AZOBS)*CE + YAEO = SIN(AZOBS)*CE + ZAEO = COS(ZDOBS) + +* Cartesian Az/ZD to Cartesian -HA,Dec + V(1) = SPHI*XAEO+CPHI*ZAEO + V(2) = YAEO + V(3) = -CPHI*XAEO+SPHI*ZAEO + +* To spherical -HA,Dec + CALL sla_DCC2S(V,HMOBS,DCOBS) + +* Right Ascension + RAOBS = sla_DRANRM(ST+HMOBS) + +* Return the results + AOB = AZOBS + ZOB = ZDOBS + HOB = -HMOBS + DOB = DCOBS + ROB = RAOBS + + END Index: trunk/psLib/src/astronomy/astronomy.i =================================================================== --- trunk/psLib/src/astronomy/astronomy.i (revision 3114) +++ trunk/psLib/src/astronomy/astronomy.i (revision 3114) @@ -0,0 +1,9 @@ +/* astronomy headers */ +%include "psAstrometry.h" +%include "psAstronomyErrors.h" +%include "psCoord.h" +%include "psMetadata.h" +%include "psMetadataIO.h" +%include "psPhotometry.h" +%include "psTime.h" + Index: trunk/psLib/src/astronomy/atms.f =================================================================== --- trunk/psLib/src/astronomy/atms.f (revision 3114) +++ trunk/psLib/src/astronomy/atms.f (revision 3114) @@ -0,0 +1,57 @@ + SUBROUTINE sla__ATMS (RT, TT, DNT, GAMAL, R, DN, RDNDR) +*+ +* - - - - - +* A T M S +* - - - - - +* +* Internal routine used by REFRO +* +* Refractive index and derivative with respect to height for the +* stratosphere. +* +* Given: +* RT d height of tropopause from centre of the Earth (metre) +* TT d temperature at the tropopause (deg K) +* DNT d refractive index at the tropopause +* GAMAL d constant of the atmospheric model = G*MD/R +* R d current distance from the centre of the Earth (metre) +* +* Returned: +* DN d refractive index at R +* RDNDR d R * rate the refractive index is changing at R +* +* P.T.Wallace Starlink 14 July 1995 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION RT,TT,DNT,GAMAL,R,DN,RDNDR + + DOUBLE PRECISION B,W + + + B = GAMAL/TT + W = (DNT-1D0)*EXP(-B*(R-RT)) + DN = 1D0+W + RDNDR = -R*B*W + + END Index: trunk/psLib/src/astronomy/atmt.f =================================================================== --- trunk/psLib/src/astronomy/atmt.f (revision 3114) +++ trunk/psLib/src/astronomy/atmt.f (revision 3114) @@ -0,0 +1,71 @@ + SUBROUTINE sla__ATMT (R0, T0, ALPHA, GAMM2, DELM2, + : C1, C2, C3, C4, C5, C6, R, T, DN, RDNDR) +*+ +* - - - - - +* A T M T +* - - - - - +* +* Internal routine used by REFRO +* +* Refractive index and derivative with respect to height for the +* troposphere. +* +* Given: +* R0 d height of observer from centre of the Earth (metre) +* T0 d temperature at the observer (deg K) +* ALPHA d alpha ) +* GAMM2 d gamma minus 2 ) see HMNAO paper +* DELM2 d delta minus 2 ) +* C1 d useful term ) +* C2 d useful term ) +* C3 d useful term ) see source +* C4 d useful term ) of sla_REFRO +* C5 d useful term ) +* C6 d useful term ) +* R d current distance from the centre of the Earth (metre) +* +* Returned: +* T d temperature at R (deg K) +* DN d refractive index at R +* RDNDR d R * rate the refractive index is changing at R +* +* Note that in the optical case C5 and C6 are zero. +* +* P.T.Wallace Starlink 30 May 1997 +* +* Copyright (C) 1997 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION R0,T0,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6, + : R,T,DN,RDNDR + + DOUBLE PRECISION TT0,TT0GM2,TT0DM2 + + + T = MAX(MIN(T0-ALPHA*(R-R0),320D0),100D0) + TT0 = T/T0 + TT0GM2 = TT0**GAMM2 + TT0DM2 = TT0**DELM2 + DN = 1D0+(C1*TT0GM2-(C2-C5/T)*TT0DM2)*TT0 + RDNDR = R*(-C3*TT0GM2+(C4-C6/TT0)*TT0DM2) + + END Index: trunk/psLib/src/astronomy/dcc2s.f =================================================================== --- trunk/psLib/src/astronomy/dcc2s.f (revision 3114) +++ trunk/psLib/src/astronomy/dcc2s.f (revision 3114) @@ -0,0 +1,70 @@ + SUBROUTINE sla_DCC2S (V, A, B) +*+ +* - - - - - - +* D C C 2 S +* - - - - - - +* +* Direction cosines to spherical coordinates (double precision) +* +* Given: +* V d(3) x,y,z vector +* +* Returned: +* A,B d spherical coordinates in radians +* +* The spherical coordinates are longitude (+ve anticlockwise +* looking from the +ve latitude pole) and latitude. The +* Cartesian coordinates are right handed, with the x axis +* at zero longitude and latitude, and the z axis at the +* +ve latitude pole. +* +* If V is null, zero A and B are returned. +* At either pole, zero A is returned. +* +* P.T.Wallace Starlink July 1989 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION V(3),A,B + + DOUBLE PRECISION X,Y,Z,R + + + X = V(1) + Y = V(2) + Z = V(3) + R = SQRT(X*X+Y*Y) + + IF (R.EQ.0D0) THEN + A = 0D0 + ELSE + A = ATAN2(Y,X) + END IF + + IF (Z.EQ.0D0) THEN + B = 0D0 + ELSE + B = ATAN2(Z,R) + END IF + + END Index: trunk/psLib/src/astronomy/dcs2c.f =================================================================== --- trunk/psLib/src/astronomy/dcs2c.f (revision 3114) +++ trunk/psLib/src/astronomy/dcs2c.f (revision 3114) @@ -0,0 +1,58 @@ + SUBROUTINE sla_DCS2C (A, B, V) +*+ +* - - - - - - +* D C S 2 C +* - - - - - - +* +* Spherical coordinates to direction cosines (double precision) +* +* Given: +* A,B dp spherical coordinates in radians +* (RA,Dec), (Long,Lat) etc +* +* Returned: +* V dp(3) x,y,z unit vector +* +* The spherical coordinates are longitude (+ve anticlockwise +* looking from the +ve latitude pole) and latitude. The +* Cartesian coordinates are right handed, with the x axis +* at zero longitude and latitude, and the z axis at the +* +ve latitude pole. +* +* P.T.Wallace Starlink October 1984 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION A,B,V(3) + + DOUBLE PRECISION COSB + + + + COSB=COS(B) + + V(1)=COS(A)*COSB + V(2)=SIN(A)*COSB + V(3)=SIN(B) + + END Index: trunk/psLib/src/astronomy/drange.f =================================================================== --- trunk/psLib/src/astronomy/drange.f (revision 3114) +++ trunk/psLib/src/astronomy/drange.f (revision 3114) @@ -0,0 +1,49 @@ + DOUBLE PRECISION FUNCTION sla_DRANGE (ANGLE) +*+ +* - - - - - - - +* D R A N G E +* - - - - - - - +* +* Normalize angle into range +/- pi (double precision) +* +* Given: +* ANGLE dp the angle in radians +* +* The result (double precision) is ANGLE expressed in the range +/- pi. +* +* P.T.Wallace Starlink 23 November 1995 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION ANGLE + + DOUBLE PRECISION DPI,D2PI + PARAMETER (DPI=3.141592653589793238462643D0) + PARAMETER (D2PI=6.283185307179586476925287D0) + + + sla_DRANGE=MOD(ANGLE,D2PI) + IF (ABS(sla_DRANGE).GE.DPI) + : sla_DRANGE=sla_DRANGE-SIGN(D2PI,ANGLE) + + END Index: trunk/psLib/src/astronomy/dranrm.f =================================================================== --- trunk/psLib/src/astronomy/dranrm.f (revision 3114) +++ trunk/psLib/src/astronomy/dranrm.f (revision 3114) @@ -0,0 +1,48 @@ + DOUBLE PRECISION FUNCTION sla_DRANRM (ANGLE) +*+ +* - - - - - - - +* D R A N R M +* - - - - - - - +* +* Normalize angle into range 0-2 pi (double precision) +* +* Given: +* ANGLE dp the angle in radians +* +* The result is ANGLE expressed in the range 0-2 pi (double +* precision). +* +* P.T.Wallace Starlink 23 November 1995 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION ANGLE + + DOUBLE PRECISION D2PI + PARAMETER (D2PI=6.283185307179586476925286766559D0) + + + sla_DRANRM=MOD(ANGLE,D2PI) + IF (sla_DRANRM.LT.0D0) sla_DRANRM=sla_DRANRM+D2PI + + END Index: trunk/psLib/src/astronomy/eqeqx.f =================================================================== --- trunk/psLib/src/astronomy/eqeqx.f (revision 3114) +++ trunk/psLib/src/astronomy/eqeqx.f (revision 3114) @@ -0,0 +1,74 @@ + DOUBLE PRECISION FUNCTION sla_EQEQX (DATE) +*+ +* - - - - - - +* E Q E Q X +* - - - - - - +* +* Equation of the equinoxes (IAU 1994, double precision) +* +* Given: +* DATE dp TDB (loosely ET) as Modified Julian Date +* (JD-2400000.5) +* +* The result is the equation of the equinoxes (double precision) +* in radians: +* +* Greenwich apparent ST = GMST + sla_EQEQX +* +* References: IAU Resolution C7, Recommendation 3 (1994) +* Capitaine, N. & Gontier, A.-M., Astron. Astrophys., +* 275, 645-650 (1993) +* +* Called: sla_NUTC +* +* Patrick Wallace Starlink 23 August 1996 +* +* Copyright (C) 1996 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE + +* Turns to arc seconds and arc seconds to radians + DOUBLE PRECISION T2AS,AS2R + PARAMETER (T2AS=1296000D0, + : AS2R=0.484813681109535994D-5) + + DOUBLE PRECISION T,OM,DPSI,DEPS,EPS0 + + + +* Interval between basic epoch J2000.0 and current epoch (JC) + T=(DATE-51544.5D0)/36525D0 + +* Longitude of the mean ascending node of the lunar orbit on the +* ecliptic, measured from the mean equinox of date + OM=AS2R*(450160.280D0+(-5D0*T2AS-482890.539D0 + : +(7.455D0+0.008D0*T)*T)*T) + +* Nutation + CALL sla_NUTC(DATE,DPSI,DEPS,EPS0) + +* Equation of the equinoxes + sla_EQEQX=DPSI*COS(EPS0)+AS2R*(0.00264D0*SIN(OM)+ + : 0.000063D0*SIN(OM+OM)) + + END Index: trunk/psLib/src/astronomy/geoc.f =================================================================== --- trunk/psLib/src/astronomy/geoc.f (revision 3114) +++ trunk/psLib/src/astronomy/geoc.f (revision 3114) @@ -0,0 +1,74 @@ + SUBROUTINE sla_GEOC (P, H, R, Z) +*+ +* - - - - - +* G E O C +* - - - - - +* +* Convert geodetic position to geocentric (double precision) +* +* Given: +* P dp latitude (geodetic, radians) +* H dp height above reference spheroid (geodetic, metres) +* +* Returned: +* R dp distance from Earth axis (AU) +* Z dp distance from plane of Earth equator (AU) +* +* Notes: +* 1) Geocentric latitude can be obtained by evaluating ATAN2(Z,R). +* 2) IAU 1976 constants are used. +* +* Reference: +* Green,R.M., Spherical Astronomy, CUP 1985, p98. +* +* P.T.Wallace Starlink 4th October 1989 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION P,H,R,Z + +* Earth equatorial radius (metres) + DOUBLE PRECISION A0 + PARAMETER (A0=6378140D0) + +* Reference spheroid flattening factor and useful function + DOUBLE PRECISION F,B + PARAMETER (F=1D0/298.257D0,B=(1D0-F)**2) + +* Astronomical unit in metres + DOUBLE PRECISION AU + PARAMETER (AU=1.49597870D11) + + DOUBLE PRECISION SP,CP,C,S + + + +* Geodetic to geocentric conversion + SP=SIN(P) + CP=COS(P) + C=1D0/SQRT(CP*CP+B*SP*SP) + S=B*C + R=(A0*C+H)*CP/AU + Z=(A0*S+H)*SP/AU + + END Index: trunk/psLib/src/astronomy/gmst.f =================================================================== --- trunk/psLib/src/astronomy/gmst.f (revision 3114) +++ trunk/psLib/src/astronomy/gmst.f (revision 3114) @@ -0,0 +1,77 @@ + DOUBLE PRECISION FUNCTION sla_GMST (UT1) +*+ +* - - - - - +* G M S T +* - - - - - +* +* Conversion from universal time to sidereal time (double precision) +* +* Given: +* UT1 dp universal time (strictly UT1) expressed as +* modified Julian Date (JD-2400000.5) +* +* The result is the Greenwich mean sidereal time (double +* precision, radians). +* +* The IAU 1982 expression (see page S15 of 1984 Astronomical Almanac) +* is used, but rearranged to reduce rounding errors. This expression +* is always described as giving the GMST at 0 hours UT. In fact, it +* gives the difference between the GMST and the UT, which happens to +* equal the GMST (modulo 24 hours) at 0 hours UT each day. In this +* routine, the entire UT is used directly as the argument for the +* standard formula, and the fractional part of the UT is added +* separately. Note that the factor 1.0027379... does not appear in the +* IAU 1982 expression explicitly but in the form of the coefficient +* 8640184.812866, which is 86400x36525x0.0027379... +* +* See also the routine sla_GMSTA, which delivers better numerical +* precision by accepting the UT date and time as separate arguments. +* +* Called: sla_DRANRM +* +* P.T.Wallace Starlink 14 October 2001 +* +* Copyright (C) 2001 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION UT1 + + DOUBLE PRECISION sla_DRANRM + + DOUBLE PRECISION D2PI,S2R + PARAMETER (D2PI=6.283185307179586476925286766559D0, + : S2R=7.272205216643039903848711535369D-5) + + DOUBLE PRECISION TU + + + +* Julian centuries from fundamental epoch J2000 to this UT + TU=(UT1-51544.5D0)/36525D0 + +* GMST at this UT + sla_GMST=sla_DRANRM(MOD(UT1,1D0)*D2PI+ + : (24110.54841D0+ + : (8640184.812866D0+ + : (0.093104D0-6.2D-6*TU)*TU)*TU)*S2R) + + END Index: trunk/psLib/src/astronomy/nutc.f =================================================================== --- trunk/psLib/src/astronomy/nutc.f (revision 3114) +++ trunk/psLib/src/astronomy/nutc.f (revision 3114) @@ -0,0 +1,830 @@ + SUBROUTINE sla_NUTC (DATE, DPSI, DEPS, EPS0) +*+ +* - - - - - +* N U T C +* - - - - - +* +* Nutation: longitude & obliquity components and mean obliquity, +* using the Shirai & Fukushima (2001) theory. +* +* Given: +* DATE d TDB (loosely ET) as Modified Julian Date +* (JD-2400000.5) +* Returned: +* DPSI,DEPS d nutation in longitude,obliquity +* EPS0 d mean obliquity +* +* Notes: +* +* 1 The routine predicts forced nutation (but not free core nutation) +* plus corrections to the IAU 1976 precession model. +* +* 2 Earth attitude predictions made by combining the present nutation +* model with IAU 1976 precession are accurate to 1 mas (with respect +* to the ICRF) for a few decades around 2000. +* +* 3 The sla_NUTC80 routine is the equivalent of the present routine +* but using the IAU 1980 nutation theory. The older theory is less +* accurate, leading to errors as large as 350 mas over the interval +* 1900-2100, mainly because of the error in the IAU 1976 precession. +* +* References: +* +* Shirai, T. & Fukushima, T., Astron.J. 121, 3270-3283 (2001). +* +* Fukushima, T., 1991, Astron.Astrophys. 244, L11 (1991). +* +* Simon, J. L., Bretagnon, P., Chapront, J., Chapront-Touze, M., +* Francou, G. & Laskar, J., Astron.Astrophys. 282, 663 (1994). +* +* P.T.Wallace Starlink 7 October 2001 +* +* Copyright (C) 2001 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE,DPSI,DEPS,EPS0 + +* Degrees to radians + DOUBLE PRECISION DD2R + PARAMETER (DD2R=1.745329251994329576923691D-2) + +* Arc seconds to radians + DOUBLE PRECISION DAS2R + PARAMETER (DAS2R=4.848136811095359935899141D-6) + +* Arc seconds in a full circle + DOUBLE PRECISION TURNAS + PARAMETER (TURNAS=1296000D0) + +* Reference epoch (J2000), MJD + DOUBLE PRECISION DJM0 + PARAMETER (DJM0=51544.5D0 ) + +* Days per Julian century + DOUBLE PRECISION DJC + PARAMETER (DJC=36525D0) + + INTEGER I,J + DOUBLE PRECISION T,EL,ELP,F,D,OM,VE,MA,JU,SA,THETA,C,S,DP,DE + +* Number of terms in the nutation model + INTEGER NTERMS + PARAMETER (NTERMS=194) + +* The SF2001 forced nutation model + INTEGER NA(9,NTERMS) + DOUBLE PRECISION PSI(4,NTERMS), EPS(4,NTERMS) + +* Coefficients of fundamental angles + DATA ( ( NA(I,J), I=1,9 ), J=1,10 ) / + : 0, 0, 0, 0, -1, 0, 0, 0, 0, + : 0, 0, 2, -2, 2, 0, 0, 0, 0, + : 0, 0, 2, 0, 2, 0, 0, 0, 0, + : 0, 0, 0, 0, -2, 0, 0, 0, 0, + : 0, 1, 0, 0, 0, 0, 0, 0, 0, + : 0, 1, 2, -2, 2, 0, 0, 0, 0, + : 1, 0, 0, 0, 0, 0, 0, 0, 0, + : 0, 0, 2, 0, 1, 0, 0, 0, 0, + : 1, 0, 2, 0, 2, 0, 0, 0, 0, + : 0, -1, 2, -2, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=11,20 ) / + : 0, 0, 2, -2, 1, 0, 0, 0, 0, + : -1, 0, 2, 0, 2, 0, 0, 0, 0, + : -1, 0, 0, 2, 0, 0, 0, 0, 0, + : 1, 0, 0, 0, 1, 0, 0, 0, 0, + : 1, 0, 0, 0, -1, 0, 0, 0, 0, + : -1, 0, 2, 2, 2, 0, 0, 0, 0, + : 1, 0, 2, 0, 1, 0, 0, 0, 0, + : -2, 0, 2, 0, 1, 0, 0, 0, 0, + : 0, 0, 0, 2, 0, 0, 0, 0, 0, + : 0, 0, 2, 2, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=21,30 ) / + : 2, 0, 0, -2, 0, 0, 0, 0, 0, + : 2, 0, 2, 0, 2, 0, 0, 0, 0, + : 1, 0, 2, -2, 2, 0, 0, 0, 0, + : -1, 0, 2, 0, 1, 0, 0, 0, 0, + : 2, 0, 0, 0, 0, 0, 0, 0, 0, + : 0, 0, 2, 0, 0, 0, 0, 0, 0, + : 0, 1, 0, 0, 1, 0, 0, 0, 0, + : -1, 0, 0, 2, 1, 0, 0, 0, 0, + : 0, 2, 2, -2, 2, 0, 0, 0, 0, + : 0, 0, 2, -2, 0, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=31,40 ) / + : -1, 0, 0, 2, -1, 0, 0, 0, 0, + : 0, 1, 0, 0, -1, 0, 0, 0, 0, + : 0, 2, 0, 0, 0, 0, 0, 0, 0, + : -1, 0, 2, 2, 1, 0, 0, 0, 0, + : 1, 0, 2, 2, 2, 0, 0, 0, 0, + : 0, 1, 2, 0, 2, 0, 0, 0, 0, + : -2, 0, 2, 0, 0, 0, 0, 0, 0, + : 0, 0, 2, 2, 1, 0, 0, 0, 0, + : 0, -1, 2, 0, 2, 0, 0, 0, 0, + : 0, 0, 0, 2, 1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=41,50 ) / + : 1, 0, 2, -2, 1, 0, 0, 0, 0, + : 2, 0, 0, -2, -1, 0, 0, 0, 0, + : 2, 0, 2, -2, 2, 0, 0, 0, 0, + : 2, 0, 2, 0, 1, 0, 0, 0, 0, + : 0, 0, 0, 2, -1, 0, 0, 0, 0, + : 0, -1, 2, -2, 1, 0, 0, 0, 0, + : -1, -1, 0, 2, 0, 0, 0, 0, 0, + : 2, 0, 0, -2, 1, 0, 0, 0, 0, + : 1, 0, 0, 2, 0, 0, 0, 0, 0, + : 0, 1, 2, -2, 1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=51,60 ) / + : 1, -1, 0, 0, 0, 0, 0, 0, 0, + : -2, 0, 2, 0, 2, 0, 0, 0, 0, + : 0, -1, 0, 2, 0, 0, 0, 0, 0, + : 3, 0, 2, 0, 2, 0, 0, 0, 0, + : 0, 0, 0, 1, 0, 0, 0, 0, 0, + : 1, -1, 2, 0, 2, 0, 0, 0, 0, + : 1, 0, 0, -1, 0, 0, 0, 0, 0, + : -1, -1, 2, 2, 2, 0, 0, 0, 0, + : -1, 0, 2, 0, 0, 0, 0, 0, 0, + : 2, 0, 0, 0, -1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=61,70 ) / + : 0, -1, 2, 2, 2, 0, 0, 0, 0, + : 1, 1, 2, 0, 2, 0, 0, 0, 0, + : 2, 0, 0, 0, 1, 0, 0, 0, 0, + : 1, 1, 0, 0, 0, 0, 0, 0, 0, + : 1, 0, -2, 2, -1, 0, 0, 0, 0, + : 1, 0, 2, 0, 0, 0, 0, 0, 0, + : -1, 1, 0, 1, 0, 0, 0, 0, 0, + : 1, 0, 0, 0, 2, 0, 0, 0, 0, + : -1, 0, 1, 0, 1, 0, 0, 0, 0, + : 0, 0, 2, 1, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=71,80 ) / + : -1, 1, 0, 1, 1, 0, 0, 0, 0, + : -1, 0, 2, 4, 2, 0, 0, 0, 0, + : 0, -2, 2, -2, 1, 0, 0, 0, 0, + : 1, 0, 2, 2, 1, 0, 0, 0, 0, + : 1, 0, 0, 0, -2, 0, 0, 0, 0, + : -2, 0, 2, 2, 2, 0, 0, 0, 0, + : 1, 1, 2, -2, 2, 0, 0, 0, 0, + : -2, 0, 2, 4, 2, 0, 0, 0, 0, + : -1, 0, 4, 0, 2, 0, 0, 0, 0, + : 2, 0, 2, -2, 1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=81,90 ) / + : 1, 0, 0, -1, -1, 0, 0, 0, 0, + : 2, 0, 2, 2, 2, 0, 0, 0, 0, + : 1, 0, 0, 2, 1, 0, 0, 0, 0, + : 3, 0, 0, 0, 0, 0, 0, 0, 0, + : 0, 0, 2, -2, -1, 0, 0, 0, 0, + : 3, 0, 2, -2, 2, 0, 0, 0, 0, + : 0, 0, 4, -2, 2, 0, 0, 0, 0, + : -1, 0, 0, 4, 0, 0, 0, 0, 0, + : 0, 1, 2, 0, 1, 0, 0, 0, 0, + : 0, 0, 2, -2, 3, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=91,100 ) / + : -2, 0, 0, 4, 0, 0, 0, 0, 0, + : -1, -1, 0, 2, 1, 0, 0, 0, 0, + : -2, 0, 2, 0, -1, 0, 0, 0, 0, + : 0, 0, 2, 0, -1, 0, 0, 0, 0, + : 0, -1, 2, 0, 1, 0, 0, 0, 0, + : 0, 1, 0, 0, 2, 0, 0, 0, 0, + : 0, 0, 2, -1, 2, 0, 0, 0, 0, + : 2, 1, 0, -2, 0, 0, 0, 0, 0, + : 0, 0, 2, 4, 2, 0, 0, 0, 0, + : -1, -1, 0, 2, -1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=101,110 ) / + : -1, 1, 0, 2, 0, 0, 0, 0, 0, + : 1, -1, 0, 0, 1, 0, 0, 0, 0, + : 0, -1, 2, -2, 0, 0, 0, 0, 0, + : 0, 1, 0, 0, -2, 0, 0, 0, 0, + : 1, -1, 2, 2, 2, 0, 0, 0, 0, + : 1, 0, 0, 2, -1, 0, 0, 0, 0, + : -1, 1, 2, 2, 2, 0, 0, 0, 0, + : 3, 0, 2, 0, 1, 0, 0, 0, 0, + : 0, 1, 2, 2, 2, 0, 0, 0, 0, + : 1, 0, 2, -2, 0, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=111,120 ) / + : -1, 0, -2, 4, -1, 0, 0, 0, 0, + : -1, -1, 2, 2, 1, 0, 0, 0, 0, + : 0, -1, 2, 2, 1, 0, 0, 0, 0, + : 2, -1, 2, 0, 2, 0, 0, 0, 0, + : 0, 0, 0, 2, 2, 0, 0, 0, 0, + : 1, -1, 2, 0, 1, 0, 0, 0, 0, + : -1, 1, 2, 0, 2, 0, 0, 0, 0, + : 0, 1, 0, 2, 0, 0, 0, 0, 0, + : 0, 1, 2, -2, 0, 0, 0, 0, 0, + : 0, 3, 2, -2, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=121,130 ) / + : 0, 0, 0, 1, 1, 0, 0, 0, 0, + : -1, 0, 2, 2, 0, 0, 0, 0, 0, + : 2, 1, 2, 0, 2, 0, 0, 0, 0, + : 1, 1, 0, 0, 1, 0, 0, 0, 0, + : 2, 0, 0, 2, 0, 0, 0, 0, 0, + : 1, 1, 2, 0, 1, 0, 0, 0, 0, + : -1, 0, 0, 2, 2, 0, 0, 0, 0, + : 1, 0, -2, 2, 0, 0, 0, 0, 0, + : 0, -1, 0, 2, -1, 0, 0, 0, 0, + : -1, 0, 1, 0, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=131,140 ) / + : 0, 1, 0, 1, 0, 0, 0, 0, 0, + : 1, 0, -2, 2, -2, 0, 0, 0, 0, + : 0, 0, 0, 1, -1, 0, 0, 0, 0, + : 1, -1, 0, 0, -1, 0, 0, 0, 0, + : 0, 0, 0, 4, 0, 0, 0, 0, 0, + : 1, -1, 0, 2, 0, 0, 0, 0, 0, + : 1, 0, 2, 1, 2, 0, 0, 0, 0, + : 1, 0, 2, -1, 2, 0, 0, 0, 0, + : -1, 0, 0, 2, -2, 0, 0, 0, 0, + : 0, 0, 2, 1, 1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=141,150 ) / + : -1, 0, 2, 0, -1, 0, 0, 0, 0, + : -1, 0, 2, 4, 1, 0, 0, 0, 0, + : 0, 0, 2, 2, 0, 0, 0, 0, 0, + : 1, 1, 2, -2, 1, 0, 0, 0, 0, + : 0, 0, 1, 0, 1, 0, 0, 0, 0, + : -1, 0, 2, -1, 1, 0, 0, 0, 0, + : -2, 0, 2, 2, 1, 0, 0, 0, 0, + : 2, -1, 0, 0, 0, 0, 0, 0, 0, + : 4, 0, 2, 0, 2, 0, 0, 0, 0, + : 2, 1, 2, -2, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=151,160 ) / + : 0, 1, 2, 1, 2, 0, 0, 0, 0, + : 1, 0, 4, -2, 2, 0, 0, 0, 0, + : 1, 1, 0, 0, -1, 0, 0, 0, 0, + : -2, 0, 2, 4, 1, 0, 0, 0, 0, + : 2, 0, 2, 0, 0, 0, 0, 0, 0, + : -1, 0, 1, 0, 0, 0, 0, 0, 0, + : 1, 0, 0, 1, 0, 0, 0, 0, 0, + : 0, 1, 0, 2, 1, 0, 0, 0, 0, + : -1, 0, 4, 0, 1, 0, 0, 0, 0, + : -1, 0, 0, 4, 1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=161,170 ) / + : 2, 0, 2, 2, 1, 0, 0, 0, 0, + : 2, 1, 0, 0, 0, 0, 0, 0, 0, + : 0, 0, 5, -5, 5, -3, 0, 0, 0, + : 0, 0, 0, 0, 0, 0, 0, 2, 0, + : 0, 0, 1, -1, 1, 0, 0, -1, 0, + : 0, 0, -1, 1, -1, 1, 0, 0, 0, + : 0, 0, -1, 1, 0, 0, 2, 0, 0, + : 0, 0, 3, -3, 3, 0, 0, -1, 0, + : 0, 0, -8, 8, -7, 5, 0, 0, 0, + : 0, 0, -1, 1, -1, 0, 2, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=171,180 ) / + : 0, 0, -2, 2, -2, 2, 0, 0, 0, + : 0, 0, -6, 6, -6, 4, 0, 0, 0, + : 0, 0, -2, 2, -2, 0, 8, -3, 0, + : 0, 0, 6, -6, 6, 0, -8, 3, 0, + : 0, 0, 4, -4, 4, -2, 0, 0, 0, + : 0, 0, -3, 3, -3, 2, 0, 0, 0, + : 0, 0, 4, -4, 3, 0, -8, 3, 0, + : 0, 0, -4, 4, -5, 0, 8, -3, 0, + : 0, 0, 0, 0, 0, 2, 0, 0, 0, + : 0, 0, -4, 4, -4, 3, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=181,190 ) / + : 0, 1, -1, 1, -1, 0, 0, 1, 0, + : 0, 0, 0, 0, 0, 0, 0, 1, 0, + : 0, 0, 1, -1, 1, 1, 0, 0, 0, + : 0, 0, 2, -2, 2, 0, -2, 0, 0, + : 0, -1, -7, 7, -7, 5, 0, 0, 0, + : -2, 0, 2, 0, 2, 0, 0, -2, 0, + : -2, 0, 2, 0, 1, 0, 0, -3, 0, + : 0, 0, 2, -2, 2, 0, 0, -2, 0, + : 0, 0, 1, -1, 1, 0, 0, 1, 0, + : 0, 0, 0, 0, 0, 0, 0, 0, 2 / + DATA ( ( NA(I,J), I=1,9 ), J=191,NTERMS ) / + : 0, 0, 0, 0, 0, 0, 0, 0, 1, + : 2, 0, -2, 0, -2, 0, 0, 3, 0, + : 0, 0, 1, -1, 1, 0, 0, -2, 0, + : 0, 0, -7, 7, -7, 5, 0, 0, 0 / + +* Nutation series: longitude + DATA ( ( PSI(I,J), I=1,4 ), J=1,10 ) / + : 3341.5D0, 17206241.8D0, 3.1D0, 17409.5D0, + : -1716.8D0, -1317185.3D0, 1.4D0, -156.8D0, + : 285.7D0, -227667.0D0, 0.3D0, -23.5D0, + : -68.6D0, -207448.0D0, 0.0D0, -21.4D0, + : 950.3D0, 147607.9D0, -2.3D0, -355.0D0, + : -66.7D0, -51689.1D0, 0.2D0, 122.6D0, + : -108.6D0, 71117.6D0, 0.0D0, 7.0D0, + : 35.6D0, -38740.2D0, 0.1D0, -36.2D0, + : 85.4D0, -30127.6D0, 0.0D0, -3.1D0, + : 9.0D0, 21583.0D0, 0.1D0, -50.3D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=11,20 ) / + : 22.1D0, 12822.8D0, 0.0D0, 13.3D0, + : 3.4D0, 12350.8D0, 0.0D0, 1.3D0, + : -21.1D0, 15699.4D0, 0.0D0, 1.6D0, + : 4.2D0, 6313.8D0, 0.0D0, 6.2D0, + : -22.8D0, 5796.9D0, 0.0D0, 6.1D0, + : 15.7D0, -5961.1D0, 0.0D0, -0.6D0, + : 13.1D0, -5159.1D0, 0.0D0, -4.6D0, + : 1.8D0, 4592.7D0, 0.0D0, 4.5D0, + : -17.5D0, 6336.0D0, 0.0D0, 0.7D0, + : 16.3D0, -3851.1D0, 0.0D0, -0.4D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=21,30 ) / + : -2.8D0, 4771.7D0, 0.0D0, 0.5D0, + : 13.8D0, -3099.3D0, 0.0D0, -0.3D0, + : 0.2D0, 2860.3D0, 0.0D0, 0.3D0, + : 1.4D0, 2045.3D0, 0.0D0, 2.0D0, + : -8.6D0, 2922.6D0, 0.0D0, 0.3D0, + : -7.7D0, 2587.9D0, 0.0D0, 0.2D0, + : 8.8D0, -1408.1D0, 0.0D0, 3.7D0, + : 1.4D0, 1517.5D0, 0.0D0, 1.5D0, + : -1.9D0, -1579.7D0, 0.0D0, 7.7D0, + : 1.3D0, -2178.6D0, 0.0D0, -0.2D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=31,40 ) / + : -4.8D0, 1286.8D0, 0.0D0, 1.3D0, + : 6.3D0, 1267.2D0, 0.0D0, -4.0D0, + : -1.0D0, 1669.3D0, 0.0D0, -8.3D0, + : 2.4D0, -1020.0D0, 0.0D0, -0.9D0, + : 4.5D0, -766.9D0, 0.0D0, 0.0D0, + : -1.1D0, 756.5D0, 0.0D0, -1.7D0, + : -1.4D0, -1097.3D0, 0.0D0, -0.5D0, + : 2.6D0, -663.0D0, 0.0D0, -0.6D0, + : 0.8D0, -714.1D0, 0.0D0, 1.6D0, + : 0.4D0, -629.9D0, 0.0D0, -0.6D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=41,50 ) / + : 0.3D0, 580.4D0, 0.0D0, 0.6D0, + : -1.6D0, 577.3D0, 0.0D0, 0.5D0, + : -0.9D0, 644.4D0, 0.0D0, 0.0D0, + : 2.2D0, -534.0D0, 0.0D0, -0.5D0, + : -2.5D0, 493.3D0, 0.0D0, 0.5D0, + : -0.1D0, -477.3D0, 0.0D0, -2.4D0, + : -0.9D0, 735.0D0, 0.0D0, -1.7D0, + : 0.7D0, 406.2D0, 0.0D0, 0.4D0, + : -2.8D0, 656.9D0, 0.0D0, 0.0D0, + : 0.6D0, 358.0D0, 0.0D0, 2.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=51,60 ) / + : -0.7D0, 472.5D0, 0.0D0, -1.1D0, + : -0.1D0, -300.5D0, 0.0D0, 0.0D0, + : -1.2D0, 435.1D0, 0.0D0, -1.0D0, + : 1.8D0, -289.4D0, 0.0D0, 0.0D0, + : 0.6D0, -422.6D0, 0.0D0, 0.0D0, + : 0.8D0, -287.6D0, 0.0D0, 0.6D0, + : -38.6D0, -392.3D0, 0.0D0, 0.0D0, + : 0.7D0, -281.8D0, 0.0D0, 0.6D0, + : 0.6D0, -405.7D0, 0.0D0, 0.0D0, + : -1.2D0, 229.0D0, 0.0D0, 0.2D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=61,70 ) / + : 1.1D0, -264.3D0, 0.0D0, 0.5D0, + : -0.7D0, 247.9D0, 0.0D0, -0.5D0, + : -0.2D0, 218.0D0, 0.0D0, 0.2D0, + : 0.6D0, -339.0D0, 0.0D0, 0.8D0, + : -0.7D0, 198.7D0, 0.0D0, 0.2D0, + : -1.5D0, 334.0D0, 0.0D0, 0.0D0, + : 0.1D0, 334.0D0, 0.0D0, 0.0D0, + : -0.1D0, -198.1D0, 0.0D0, 0.0D0, + : -106.6D0, 0.0D0, 0.0D0, 0.0D0, + : -0.5D0, 165.8D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=71,80 ) / + : 0.0D0, 134.8D0, 0.0D0, 0.0D0, + : 0.9D0, -151.6D0, 0.0D0, 0.0D0, + : 0.0D0, -129.7D0, 0.0D0, 0.0D0, + : 0.8D0, -132.8D0, 0.0D0, -0.1D0, + : 0.5D0, -140.7D0, 0.0D0, 0.0D0, + : -0.1D0, 138.4D0, 0.0D0, 0.0D0, + : 0.0D0, 129.0D0, 0.0D0, -0.3D0, + : 0.5D0, -121.2D0, 0.0D0, 0.0D0, + : -0.3D0, 114.5D0, 0.0D0, 0.0D0, + : -0.1D0, 101.8D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=81,90 ) / + : -3.6D0, -101.9D0, 0.0D0, 0.0D0, + : 0.8D0, -109.4D0, 0.0D0, 0.0D0, + : 0.2D0, -97.0D0, 0.0D0, 0.0D0, + : -0.7D0, 157.3D0, 0.0D0, 0.0D0, + : 0.2D0, -83.3D0, 0.0D0, 0.0D0, + : -0.3D0, 93.3D0, 0.0D0, 0.0D0, + : -0.1D0, 92.1D0, 0.0D0, 0.0D0, + : -0.5D0, 133.6D0, 0.0D0, 0.0D0, + : -0.1D0, 81.5D0, 0.0D0, 0.0D0, + : 0.0D0, 123.9D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=91,100 ) / + : -0.3D0, 128.1D0, 0.0D0, 0.0D0, + : 0.1D0, 74.1D0, 0.0D0, -0.3D0, + : -0.2D0, -70.3D0, 0.0D0, 0.0D0, + : -0.4D0, 66.6D0, 0.0D0, 0.0D0, + : 0.1D0, -66.7D0, 0.0D0, 0.0D0, + : -0.7D0, 69.3D0, 0.0D0, -0.3D0, + : 0.0D0, -70.4D0, 0.0D0, 0.0D0, + : -0.1D0, 101.5D0, 0.0D0, 0.0D0, + : 0.5D0, -69.1D0, 0.0D0, 0.0D0, + : -0.2D0, 58.5D0, 0.0D0, 0.2D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=101,110 ) / + : 0.1D0, -94.9D0, 0.0D0, 0.2D0, + : 0.0D0, 52.9D0, 0.0D0, -0.2D0, + : 0.1D0, 86.7D0, 0.0D0, -0.2D0, + : -0.1D0, -59.2D0, 0.0D0, 0.2D0, + : 0.3D0, -58.8D0, 0.0D0, 0.1D0, + : -0.3D0, 49.0D0, 0.0D0, 0.0D0, + : -0.2D0, 56.9D0, 0.0D0, -0.1D0, + : 0.3D0, -50.2D0, 0.0D0, 0.0D0, + : -0.2D0, 53.4D0, 0.0D0, -0.1D0, + : 0.1D0, -76.5D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=111,120 ) / + : -0.2D0, 45.3D0, 0.0D0, 0.0D0, + : 0.1D0, -46.8D0, 0.0D0, 0.0D0, + : 0.2D0, -44.6D0, 0.0D0, 0.0D0, + : 0.2D0, -48.7D0, 0.0D0, 0.0D0, + : 0.1D0, -46.8D0, 0.0D0, 0.0D0, + : 0.1D0, -42.0D0, 0.0D0, 0.0D0, + : 0.0D0, 46.4D0, 0.0D0, -0.1D0, + : 0.2D0, -67.3D0, 0.0D0, 0.1D0, + : 0.0D0, -65.8D0, 0.0D0, 0.2D0, + : -0.1D0, -43.9D0, 0.0D0, 0.3D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=121,130 ) / + : 0.0D0, -38.9D0, 0.0D0, 0.0D0, + : -0.3D0, 63.9D0, 0.0D0, 0.0D0, + : -0.2D0, 41.2D0, 0.0D0, 0.0D0, + : 0.0D0, -36.1D0, 0.0D0, 0.2D0, + : -0.3D0, 58.5D0, 0.0D0, 0.0D0, + : -0.1D0, 36.1D0, 0.0D0, 0.0D0, + : 0.0D0, -39.7D0, 0.0D0, 0.0D0, + : 0.1D0, -57.7D0, 0.0D0, 0.0D0, + : -0.2D0, 33.4D0, 0.0D0, 0.0D0, + : 36.4D0, 0.0D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=131,140 ) / + : -0.1D0, 55.7D0, 0.0D0, -0.1D0, + : 0.1D0, -35.4D0, 0.0D0, 0.0D0, + : 0.1D0, -31.0D0, 0.0D0, 0.0D0, + : -0.1D0, 30.1D0, 0.0D0, 0.0D0, + : -0.3D0, 49.2D0, 0.0D0, 0.0D0, + : -0.2D0, 49.1D0, 0.0D0, 0.0D0, + : -0.1D0, 33.6D0, 0.0D0, 0.0D0, + : 0.1D0, -33.5D0, 0.0D0, 0.0D0, + : 0.1D0, -31.0D0, 0.0D0, 0.0D0, + : -0.1D0, 28.0D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=141,150 ) / + : 0.1D0, -25.2D0, 0.0D0, 0.0D0, + : 0.1D0, -26.2D0, 0.0D0, 0.0D0, + : -0.2D0, 41.5D0, 0.0D0, 0.0D0, + : 0.0D0, 24.5D0, 0.0D0, 0.1D0, + : -16.2D0, 0.0D0, 0.0D0, 0.0D0, + : 0.0D0, -22.3D0, 0.0D0, 0.0D0, + : 0.0D0, 23.1D0, 0.0D0, 0.0D0, + : -0.1D0, 37.5D0, 0.0D0, 0.0D0, + : 0.2D0, -25.7D0, 0.0D0, 0.0D0, + : 0.0D0, 25.2D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=151,160 ) / + : 0.1D0, -24.5D0, 0.0D0, 0.0D0, + : -0.1D0, 24.3D0, 0.0D0, 0.0D0, + : 0.1D0, -20.7D0, 0.0D0, 0.0D0, + : 0.1D0, -20.8D0, 0.0D0, 0.0D0, + : -0.2D0, 33.4D0, 0.0D0, 0.0D0, + : 32.9D0, 0.0D0, 0.0D0, 0.0D0, + : 0.1D0, -32.6D0, 0.0D0, 0.0D0, + : 0.0D0, 19.9D0, 0.0D0, 0.0D0, + : -0.1D0, 19.6D0, 0.0D0, 0.0D0, + : 0.0D0, -18.7D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=161,170 ) / + : 0.1D0, -19.0D0, 0.0D0, 0.0D0, + : 0.1D0, -28.6D0, 0.0D0, 0.0D0, + : 4.0D0, 178.8D0,-11.8D0, 0.3D0, + : 39.8D0, -107.3D0, -5.6D0, -1.0D0, + : 9.9D0, 164.0D0, -4.1D0, 0.1D0, + : -4.8D0, -135.3D0, -3.4D0, -0.1D0, + : 50.5D0, 75.0D0, 1.4D0, -1.2D0, + : -1.1D0, -53.5D0, 1.3D0, 0.0D0, + : -45.0D0, -2.4D0, -0.4D0, 6.6D0, + : -11.5D0, -61.0D0, -0.9D0, 0.4D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=171,180 ) / + : 4.4D0, -68.4D0, -3.4D0, 0.0D0, + : 7.7D0, -47.1D0, -4.7D0, -1.0D0, + : -42.9D0, -12.6D0, -1.2D0, 4.2D0, + : -42.8D0, 12.7D0, -1.2D0, -4.2D0, + : -7.6D0, -44.1D0, 2.1D0, -0.5D0, + : -64.1D0, 1.7D0, 0.2D0, 4.5D0, + : 36.4D0, -10.4D0, 1.0D0, 3.5D0, + : 35.6D0, 10.2D0, 1.0D0, -3.5D0, + : -1.7D0, 39.5D0, 2.0D0, 0.0D0, + : 50.9D0, -8.2D0, -0.8D0, -5.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=181,190 ) / + : 0.0D0, 52.3D0, 1.2D0, 0.0D0, + : -42.9D0, -17.8D0, 0.4D0, 0.0D0, + : 2.6D0, 34.3D0, 0.8D0, 0.0D0, + : -0.8D0, -48.6D0, 2.4D0, -0.1D0, + : -4.9D0, 30.5D0, 3.7D0, 0.7D0, + : 0.0D0, -43.6D0, 2.1D0, 0.0D0, + : 0.0D0, -25.4D0, 1.2D0, 0.0D0, + : 2.0D0, 40.9D0, -2.0D0, 0.0D0, + : -2.1D0, 26.1D0, 0.6D0, 0.0D0, + : 22.6D0, -3.2D0, -0.5D0, -0.5D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=191,NTERMS ) / + : -7.6D0, 24.9D0, -0.4D0, -0.2D0, + : -6.2D0, 34.9D0, 1.7D0, 0.3D0, + : 2.0D0, 17.4D0, -0.4D0, 0.1D0, + : -3.9D0, 20.5D0, 2.4D0, 0.6D0 / + +* Nutation series: obliquity + DATA ( ( EPS(I,J), I=1,4 ), J=1,10 ) / + : 9205365.8D0, -1506.2D0, 885.7D0, -0.2D0, + : 573095.9D0, -570.2D0, -305.0D0, -0.3D0, + : 97845.5D0, 147.8D0, -48.8D0, -0.2D0, + : -89753.6D0, 28.0D0, 46.9D0, 0.0D0, + : 7406.7D0, -327.1D0, -18.2D0, 0.8D0, + : 22442.3D0, -22.3D0, -67.6D0, 0.0D0, + : -683.6D0, 46.8D0, 0.0D0, 0.0D0, + : 20070.7D0, 36.0D0, 1.6D0, 0.0D0, + : 12893.8D0, 39.5D0, -6.2D0, 0.0D0, + : -9593.2D0, 14.4D0, 30.2D0, -0.1D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=11,20 ) / + : -6899.5D0, 4.8D0, -0.6D0, 0.0D0, + : -5332.5D0, -0.1D0, 2.7D0, 0.0D0, + : -125.2D0, 10.5D0, 0.0D0, 0.0D0, + : -3323.4D0, -0.9D0, -0.3D0, 0.0D0, + : 3142.3D0, 8.9D0, 0.3D0, 0.0D0, + : 2552.5D0, 7.3D0, -1.2D0, 0.0D0, + : 2634.4D0, 8.8D0, 0.2D0, 0.0D0, + : -2424.4D0, 1.6D0, -0.4D0, 0.0D0, + : -123.3D0, 3.9D0, 0.0D0, 0.0D0, + : 1642.4D0, 7.3D0, -0.8D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=21,30 ) / + : 47.9D0, 3.2D0, 0.0D0, 0.0D0, + : 1321.2D0, 6.2D0, -0.6D0, 0.0D0, + : -1234.1D0, -0.3D0, 0.6D0, 0.0D0, + : -1076.5D0, -0.3D0, 0.0D0, 0.0D0, + : -61.6D0, 1.8D0, 0.0D0, 0.0D0, + : -55.4D0, 1.6D0, 0.0D0, 0.0D0, + : 856.9D0, -4.9D0, -2.1D0, 0.0D0, + : -800.7D0, -0.1D0, 0.0D0, 0.0D0, + : 685.1D0, -0.6D0, -3.8D0, 0.0D0, + : -16.9D0, -1.5D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=31,40 ) / + : 695.7D0, 1.8D0, 0.0D0, 0.0D0, + : 642.2D0, -2.6D0, -1.6D0, 0.0D0, + : 13.3D0, 1.1D0, -0.1D0, 0.0D0, + : 521.9D0, 1.6D0, 0.0D0, 0.0D0, + : 325.8D0, 2.0D0, -0.1D0, 0.0D0, + : -325.1D0, -0.5D0, 0.9D0, 0.0D0, + : 10.1D0, 0.3D0, 0.0D0, 0.0D0, + : 334.5D0, 1.6D0, 0.0D0, 0.0D0, + : 307.1D0, 0.4D0, -0.9D0, 0.0D0, + : 327.2D0, 0.5D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=41,50 ) / + : -304.6D0, -0.1D0, 0.0D0, 0.0D0, + : 304.0D0, 0.6D0, 0.0D0, 0.0D0, + : -276.8D0, -0.5D0, 0.1D0, 0.0D0, + : 268.9D0, 1.3D0, 0.0D0, 0.0D0, + : 271.8D0, 1.1D0, 0.0D0, 0.0D0, + : 271.5D0, -0.4D0, -0.8D0, 0.0D0, + : -5.2D0, 0.5D0, 0.0D0, 0.0D0, + : -220.5D0, 0.1D0, 0.0D0, 0.0D0, + : -20.1D0, 0.3D0, 0.0D0, 0.0D0, + : -191.0D0, 0.1D0, 0.5D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=51,60 ) / + : -4.1D0, 0.3D0, 0.0D0, 0.0D0, + : 130.6D0, -0.1D0, 0.0D0, 0.0D0, + : 3.0D0, 0.3D0, 0.0D0, 0.0D0, + : 122.9D0, 0.8D0, 0.0D0, 0.0D0, + : 3.7D0, -0.3D0, 0.0D0, 0.0D0, + : 123.1D0, 0.4D0, -0.3D0, 0.0D0, + : -52.7D0, 15.3D0, 0.0D0, 0.0D0, + : 120.7D0, 0.3D0, -0.3D0, 0.0D0, + : 4.0D0, -0.3D0, 0.0D0, 0.0D0, + : 126.5D0, 0.5D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=61,70 ) / + : 112.7D0, 0.5D0, -0.3D0, 0.0D0, + : -106.1D0, -0.3D0, 0.3D0, 0.0D0, + : -112.9D0, -0.2D0, 0.0D0, 0.0D0, + : 3.6D0, -0.2D0, 0.0D0, 0.0D0, + : 107.4D0, 0.3D0, 0.0D0, 0.0D0, + : -10.9D0, 0.2D0, 0.0D0, 0.0D0, + : -0.9D0, 0.0D0, 0.0D0, 0.0D0, + : 85.4D0, 0.0D0, 0.0D0, 0.0D0, + : 0.0D0, -88.8D0, 0.0D0, 0.0D0, + : -71.0D0, -0.2D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=71,80 ) / + : -70.3D0, 0.0D0, 0.0D0, 0.0D0, + : 64.5D0, 0.4D0, 0.0D0, 0.0D0, + : 69.8D0, 0.0D0, 0.0D0, 0.0D0, + : 66.1D0, 0.4D0, 0.0D0, 0.0D0, + : -61.0D0, -0.2D0, 0.0D0, 0.0D0, + : -59.5D0, -0.1D0, 0.0D0, 0.0D0, + : -55.6D0, 0.0D0, 0.2D0, 0.0D0, + : 51.7D0, 0.2D0, 0.0D0, 0.0D0, + : -49.0D0, -0.1D0, 0.0D0, 0.0D0, + : -52.7D0, -0.1D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=81,90 ) / + : -49.6D0, 1.4D0, 0.0D0, 0.0D0, + : 46.3D0, 0.4D0, 0.0D0, 0.0D0, + : 49.6D0, 0.1D0, 0.0D0, 0.0D0, + : -5.1D0, 0.1D0, 0.0D0, 0.0D0, + : -44.0D0, -0.1D0, 0.0D0, 0.0D0, + : -39.9D0, -0.1D0, 0.0D0, 0.0D0, + : -39.5D0, -0.1D0, 0.0D0, 0.0D0, + : -3.9D0, 0.1D0, 0.0D0, 0.0D0, + : -42.1D0, -0.1D0, 0.0D0, 0.0D0, + : -17.2D0, 0.1D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=91,100 ) / + : -2.3D0, 0.1D0, 0.0D0, 0.0D0, + : -39.2D0, 0.0D0, 0.0D0, 0.0D0, + : -38.4D0, 0.1D0, 0.0D0, 0.0D0, + : 36.8D0, 0.2D0, 0.0D0, 0.0D0, + : 34.6D0, 0.1D0, 0.0D0, 0.0D0, + : -32.7D0, 0.3D0, 0.0D0, 0.0D0, + : 30.4D0, 0.0D0, 0.0D0, 0.0D0, + : 0.4D0, 0.1D0, 0.0D0, 0.0D0, + : 29.3D0, 0.2D0, 0.0D0, 0.0D0, + : 31.6D0, 0.1D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=101,110 ) / + : 0.8D0, -0.1D0, 0.0D0, 0.0D0, + : -27.9D0, 0.0D0, 0.0D0, 0.0D0, + : 2.9D0, 0.0D0, 0.0D0, 0.0D0, + : -25.3D0, 0.0D0, 0.0D0, 0.0D0, + : 25.0D0, 0.1D0, 0.0D0, 0.0D0, + : 27.5D0, 0.1D0, 0.0D0, 0.0D0, + : -24.4D0, -0.1D0, 0.0D0, 0.0D0, + : 24.9D0, 0.2D0, 0.0D0, 0.0D0, + : -22.8D0, -0.1D0, 0.0D0, 0.0D0, + : 0.9D0, -0.1D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=111,120 ) / + : 24.4D0, 0.1D0, 0.0D0, 0.0D0, + : 23.9D0, 0.1D0, 0.0D0, 0.0D0, + : 22.5D0, 0.1D0, 0.0D0, 0.0D0, + : 20.8D0, 0.1D0, 0.0D0, 0.0D0, + : 20.1D0, 0.0D0, 0.0D0, 0.0D0, + : 21.5D0, 0.1D0, 0.0D0, 0.0D0, + : -20.0D0, 0.0D0, 0.0D0, 0.0D0, + : 1.4D0, 0.0D0, 0.0D0, 0.0D0, + : -0.2D0, -0.1D0, 0.0D0, 0.0D0, + : 19.0D0, 0.0D0, -0.1D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=121,130 ) / + : 20.5D0, 0.0D0, 0.0D0, 0.0D0, + : -2.0D0, 0.0D0, 0.0D0, 0.0D0, + : -17.6D0, -0.1D0, 0.0D0, 0.0D0, + : 19.0D0, 0.0D0, 0.0D0, 0.0D0, + : -2.4D0, 0.0D0, 0.0D0, 0.0D0, + : -18.4D0, -0.1D0, 0.0D0, 0.0D0, + : 17.1D0, 0.0D0, 0.0D0, 0.0D0, + : 0.4D0, 0.0D0, 0.0D0, 0.0D0, + : 18.4D0, 0.1D0, 0.0D0, 0.0D0, + : 0.0D0, 17.4D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=131,140 ) / + : -0.6D0, 0.0D0, 0.0D0, 0.0D0, + : -15.4D0, 0.0D0, 0.0D0, 0.0D0, + : -16.8D0, -0.1D0, 0.0D0, 0.0D0, + : 16.3D0, 0.0D0, 0.0D0, 0.0D0, + : -2.0D0, 0.0D0, 0.0D0, 0.0D0, + : -1.5D0, 0.0D0, 0.0D0, 0.0D0, + : -14.3D0, -0.1D0, 0.0D0, 0.0D0, + : 14.4D0, 0.0D0, 0.0D0, 0.0D0, + : -13.4D0, 0.0D0, 0.0D0, 0.0D0, + : -14.3D0, -0.1D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=141,150 ) / + : -13.7D0, 0.0D0, 0.0D0, 0.0D0, + : 13.1D0, 0.1D0, 0.0D0, 0.0D0, + : -1.7D0, 0.0D0, 0.0D0, 0.0D0, + : -12.8D0, 0.0D0, 0.0D0, 0.0D0, + : 0.0D0, -14.4D0, 0.0D0, 0.0D0, + : 12.4D0, 0.0D0, 0.0D0, 0.0D0, + : -12.0D0, 0.0D0, 0.0D0, 0.0D0, + : -0.8D0, 0.0D0, 0.0D0, 0.0D0, + : 10.9D0, 0.1D0, 0.0D0, 0.0D0, + : -10.8D0, 0.0D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=151,160 ) / + : 10.5D0, 0.0D0, 0.0D0, 0.0D0, + : -10.4D0, 0.0D0, 0.0D0, 0.0D0, + : -11.2D0, 0.0D0, 0.0D0, 0.0D0, + : 10.5D0, 0.1D0, 0.0D0, 0.0D0, + : -1.4D0, 0.0D0, 0.0D0, 0.0D0, + : 0.0D0, 0.1D0, 0.0D0, 0.0D0, + : 0.7D0, 0.0D0, 0.0D0, 0.0D0, + : -10.3D0, 0.0D0, 0.0D0, 0.0D0, + : -10.0D0, 0.0D0, 0.0D0, 0.0D0, + : 9.6D0, 0.0D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=161,170 ) / + : 9.4D0, 0.1D0, 0.0D0, 0.0D0, + : 0.6D0, 0.0D0, 0.0D0, 0.0D0, + : -87.7D0, 4.4D0, -0.4D0, -6.3D0, + : 46.3D0, 22.4D0, 0.5D0, -2.4D0, + : 15.6D0, -3.4D0, 0.1D0, 0.4D0, + : 5.2D0, 5.8D0, 0.2D0, -0.1D0, + : -30.1D0, 26.9D0, 0.7D0, 0.0D0, + : 23.2D0, -0.5D0, 0.0D0, 0.6D0, + : 1.0D0, 23.2D0, 3.4D0, 0.0D0, + : -12.2D0, -4.3D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=171,180 ) / + : -2.1D0, -3.7D0, -0.2D0, 0.1D0, + : -18.6D0, -3.8D0, -0.4D0, 1.8D0, + : 5.5D0, -18.7D0, -1.8D0, -0.5D0, + : -5.5D0, -18.7D0, 1.8D0, -0.5D0, + : 18.4D0, -3.6D0, 0.3D0, 0.9D0, + : -0.6D0, 1.3D0, 0.0D0, 0.0D0, + : -5.6D0, -19.5D0, 1.9D0, 0.0D0, + : 5.5D0, -19.1D0, -1.9D0, 0.0D0, + : -17.3D0, -0.8D0, 0.0D0, 0.9D0, + : -3.2D0, -8.3D0, -0.8D0, 0.3D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=181,190 ) / + : -0.1D0, 0.0D0, 0.0D0, 0.0D0, + : -5.4D0, 7.8D0, -0.3D0, 0.0D0, + : -14.8D0, 1.4D0, 0.0D0, 0.3D0, + : -3.8D0, 0.4D0, 0.0D0, -0.2D0, + : 12.6D0, 3.2D0, 0.5D0, -1.5D0, + : 0.1D0, 0.0D0, 0.0D0, 0.0D0, + : -13.6D0, 2.4D0, -0.1D0, 0.0D0, + : 0.9D0, 1.2D0, 0.0D0, 0.0D0, + : -11.9D0, -0.5D0, 0.0D0, 0.3D0, + : 0.4D0, 12.0D0, 0.3D0, -0.2D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=191,NTERMS ) / + : 8.3D0, 6.1D0, -0.1D0, 0.1D0, + : 0.0D0, 0.0D0, 0.0D0, 0.0D0, + : 0.4D0, -10.8D0, 0.3D0, 0.0D0, + : 9.6D0, 2.2D0, 0.3D0, -1.2D0 / + + + +* Interval between fundamental epoch J2000.0 and given epoch (JC). + T = (DATE-DJM0)/DJC + +* Mean anomaly of the Moon. + EL = 134.96340251D0*DD2R+ + : MOD(T*(1717915923.2178D0+ + : T*( 31.8792D0+ + : T*( 0.051635D0+ + : T*( - 0.00024470D0)))),TURNAS)*DAS2R + +* Mean anomaly of the Sun. + ELP = 357.52910918D0*DD2R+ + : MOD(T*( 129596581.0481D0+ + : T*( - 0.5532D0+ + : T*( 0.000136D0+ + : T*( - 0.00001149D0)))),TURNAS)*DAS2R + +* Mean argument of the latitude of the Moon. + F = 93.27209062D0*DD2R+ + : MOD(T*(1739527262.8478D0+ + : T*( - 12.7512D0+ + : T*( - 0.001037D0+ + : T*( 0.00000417D0)))),TURNAS)*DAS2R + +* Mean elongation of the Moon from the Sun. + D = 297.85019547D0*DD2R+ + : MOD(T*(1602961601.2090D0+ + : T*( - 6.3706D0+ + : T*( 0.006539D0+ + : T*( - 0.00003169D0)))),TURNAS)*DAS2R + +* Mean longitude of the ascending node of the Moon. + OM = 125.04455501D0*DD2R+ + : MOD(T*( - 6962890.5431D0+ + : T*( 7.4722D0+ + : T*( 0.007702D0+ + : T*( - 0.00005939D0)))),TURNAS)*DAS2R + +* Mean longitude of Venus. + VE = 181.97980085D0*DD2R+MOD(210664136.433548D0*T,TURNAS)*DAS2R + +* Mean longitude of Mars. + MA = 355.43299958D0*DD2R+MOD( 68905077.493988D0*T,TURNAS)*DAS2R + +* Mean longitude of Jupiter. + JU = 34.35151874D0*DD2R+MOD( 10925660.377991D0*T,TURNAS)*DAS2R + +* Mean longitude of Saturn. + SA = 50.07744430D0*DD2R+MOD( 4399609.855732D0*T,TURNAS)*DAS2R + +* Geodesic nutation (Fukushima 1991) in microarcsec. + DP = -153.1D0*SIN(ELP)-1.9D0*SIN(2D0*ELP) + DE = 0D0 + +* Shirai & Fukushima (2001) nutation series. + DO J=NTERMS,1,-1 + THETA = DBLE(NA(1,J))*EL+ + : DBLE(NA(2,J))*ELP+ + : DBLE(NA(3,J))*F+ + : DBLE(NA(4,J))*D+ + : DBLE(NA(5,J))*OM+ + : DBLE(NA(6,J))*VE+ + : DBLE(NA(7,J))*MA+ + : DBLE(NA(8,J))*JU+ + : DBLE(NA(9,J))*SA + C = COS(THETA) + S = SIN(THETA) + DP = DP+(PSI(1,J)+PSI(3,J)*T)*C+(PSI(2,J)+PSI(4,J)*T)*S + DE = DE+(EPS(1,J)+EPS(3,J)*T)*C+(EPS(2,J)+EPS(4,J)*T)*S + END DO + +* Change of units, and addition of the precession correction. + DPSI = (DP*1D-6-0.042888D0-0.29856D0*T)*DAS2R + DEPS = (DE*1D-6-0.005171D0-0.02408D0*T)*DAS2R + +* Mean obliquity of date (Simon et al. 1994). + EPS0 = (84381.412D0+ + : (-46.80927D0+ + : (-0.000152D0+ + : (0.0019989D0+ + : (-0.00000051D0+ + : (-0.000000025D0)*T)*T)*T)*T)*T)*DAS2R + + END Index: trunk/psLib/src/astronomy/oapqk.f =================================================================== --- trunk/psLib/src/astronomy/oapqk.f (revision 3114) +++ trunk/psLib/src/astronomy/oapqk.f (revision 3114) @@ -0,0 +1,250 @@ + SUBROUTINE sla_OAPQK (TYPE, OB1, OB2, AOPRMS, RAP, DAP) +*+ +* - - - - - - +* O A P Q K +* - - - - - - +* +* Quick observed to apparent place +* +* Given: +* TYPE c*(*) type of coordinates - 'R', 'H' or 'A' (see below) +* OB1 d observed Az, HA or RA (radians; Az is N=0,E=90) +* OB2 d observed ZD or Dec (radians) +* AOPRMS d(14) star-independent apparent-to-observed parameters: +* +* (1) geodetic latitude (radians) +* (2,3) sine and cosine of geodetic latitude +* (4) magnitude of diurnal aberration vector +* (5) height (HM) +* (6) ambient temperature (T) +* (7) pressure (P) +* (8) relative humidity (RH) +* (9) wavelength (WL) +* (10) lapse rate (TLR) +* (11,12) refraction constants A and B (radians) +* (13) longitude + eqn of equinoxes + sidereal DUT (radians) +* (14) local apparent sidereal time (radians) +* +* Returned: +* RAP d geocentric apparent right ascension +* DAP d geocentric apparent declination +* +* Notes: +* +* 1) Only the first character of the TYPE argument is significant. +* 'R' or 'r' indicates that OBS1 and OBS2 are the observed Right +* Ascension and Declination; 'H' or 'h' indicates that they are +* Hour Angle (West +ve) and Declination; anything else ('A' or +* 'a' is recommended) indicates that OBS1 and OBS2 are Azimuth +* (North zero, East is 90 deg) and zenith distance. (Zenith +* distance is used rather than elevation in order to reflect the +* fact that no allowance is made for depression of the horizon.) +* +* 2) The accuracy of the result is limited by the corrections for +* refraction. Providing the meteorological parameters are +* known accurately and there are no gross local effects, the +* predicted apparent RA,Dec should be within about 0.1 arcsec +* for a zenith distance of less than 70 degrees. Even at a +* topocentric zenith distance of 90 degrees, the accuracy in +* elevation should be better than 1 arcmin; useful results +* are available for a further 3 degrees, beyond which the +* sla_REFRO routine returns a fixed value of the refraction. +* The complementary routines sla_AOP (or sla_AOPQK) and sla_OAP +* (or sla_OAPQK) are self-consistent to better than 1 micro- +* arcsecond all over the celestial sphere. +* +* 3) It is advisable to take great care with units, as even +* unlikely values of the input parameters are accepted and +* processed in accordance with the models used. +* +* 5) "Observed" Az,El means the position that would be seen by a +* perfect theodolite located at the observer. This is +* related to the observed HA,Dec via the standard rotation, using +* the geodetic latitude (corrected for polar motion), while the +* observed HA and RA are related simply through the local +* apparent ST. "Observed" RA,Dec or HA,Dec thus means the +* position that would be seen by a perfect equatorial located +* at the observer and with its polar axis aligned to the +* Earth's axis of rotation (n.b. not to the refracted pole). +* By removing from the observed place the effects of +* atmospheric refraction and diurnal aberration, the +* geocentric apparent RA,Dec is obtained. +* +* 5) Frequently, mean rather than apparent RA,Dec will be required, +* in which case further transformations will be necessary. The +* sla_AMP etc routines will convert the apparent RA,Dec produced +* by the present routine into an "FK5" (J2000) mean place, by +* allowing for the Sun's gravitational lens effect, annual +* aberration, nutation and precession. Should "FK4" (1950) +* coordinates be needed, the routines sla_FK524 etc will also +* need to be applied. +* +* 6) To convert to apparent RA,Dec the coordinates read from a +* real telescope, corrections would have to be applied for +* encoder zero points, gear and encoder errors, tube flexure, +* the position of the rotator axis and the pointing axis +* relative to it, non-perpendicularity between the mounting +* axes, and finally for the tilt of the azimuth or polar axis +* of the mounting (with appropriate corrections for mount +* flexures). Some telescopes would, of course, exhibit other +* properties which would need to be accounted for at the +* appropriate point in the sequence. +* +* 7) The star-independent apparent-to-observed-place parameters +* in AOPRMS may be computed by means of the sla_AOPPA routine. +* If nothing has changed significantly except the time, the +* sla_AOPPAT routine may be used to perform the requisite +* partial recomputation of AOPRMS. +* +* 8) The azimuths etc used by the present routine are with respect +* to the celestial pole. Corrections from the terrestrial pole +* can be computed using sla_POLMO. +* +* Called: sla_DCS2C, sla_DCC2S, sla_REFRO, sla_DRANRM +* +* P.T.Wallace Starlink 23 June 1997 +* +* Copyright (C) 1996 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + CHARACTER*(*) TYPE + DOUBLE PRECISION OB1,OB2,AOPRMS(14),RAP,DAP + +* Breakpoint for fast/slow refraction algorithm: +* ZD greater than arctan(4), (see sla_REFCO routine) +* or vector Z less than cosine(arctan(Z)) = 1/sqrt(17) + DOUBLE PRECISION ZBREAK + PARAMETER (ZBREAK=0.242535625D0) + + CHARACTER C + DOUBLE PRECISION C1,C2,SPHI,CPHI,ST,CE,XAEO,YAEO,ZAEO,V(3), + : XMHDO,YMHDO,ZMHDO,AZ,SZ,ZDO,TZ,DREF,ZDT, + : XAET,YAET,ZAET,XMHDA,YMHDA,ZMHDA,DIURAB,F,HMA + + DOUBLE PRECISION sla_DRANRM + + + +* Coordinate type + C = TYPE(1:1) + +* Coordinates + C1 = OB1 + C2 = OB2 + +* Sin, cos of latitude + SPHI = AOPRMS(2) + CPHI = AOPRMS(3) + +* Local apparent sidereal time + ST = AOPRMS(14) + +* Standardise coordinate type + IF (C.EQ.'R'.OR.C.EQ.'r') THEN + C = 'R' + ELSE IF (C.EQ.'H'.OR.C.EQ.'h') THEN + C = 'H' + ELSE + C = 'A' + END IF + +* If Az,ZD convert to Cartesian (S=0,E=90) + IF (C.EQ.'A') THEN + CE = SIN(C2) + XAEO = -COS(C1)*CE + YAEO = SIN(C1)*CE + ZAEO = COS(C2) + ELSE + +* If RA,Dec convert to HA,Dec + IF (C.EQ.'R') THEN + C1 = ST-C1 + END IF + +* To Cartesian -HA,Dec + CALL sla_DCS2C(-C1,C2,V) + XMHDO = V(1) + YMHDO = V(2) + ZMHDO = V(3) + +* To Cartesian Az,El (S=0,E=90) + XAEO = SPHI*XMHDO-CPHI*ZMHDO + YAEO = YMHDO + ZAEO = CPHI*XMHDO+SPHI*ZMHDO + END IF + +* Azimuth (S=0,E=90) + IF (XAEO.NE.0D0.OR.YAEO.NE.0D0) THEN + AZ = ATAN2(YAEO,XAEO) + ELSE + AZ = 0D0 + END IF + +* Sine of observed ZD, and observed ZD + SZ = SQRT(XAEO*XAEO+YAEO*YAEO) + ZDO = ATAN2(SZ,ZAEO) + +* +* Refraction +* ---------- + +* Large zenith distance? + IF (ZAEO.GE.ZBREAK) THEN + +* Fast algorithm using two constant model + TZ = SZ/ZAEO + DREF = AOPRMS(11)*TZ+AOPRMS(12)*TZ*TZ*TZ + + ELSE + +* Rigorous algorithm for large ZD + CALL sla_REFRO(ZDO,AOPRMS(5),AOPRMS(6),AOPRMS(7),AOPRMS(8), + : AOPRMS(9),AOPRMS(1),AOPRMS(10),1D-8,DREF) + END IF + + ZDT = ZDO+DREF + +* To Cartesian Az,ZD + CE = SIN(ZDT) + XAET = COS(AZ)*CE + YAET = SIN(AZ)*CE + ZAET = COS(ZDT) + +* Cartesian Az,ZD to Cartesian -HA,Dec + XMHDA = SPHI*XAET+CPHI*ZAET + YMHDA = YAET + ZMHDA = -CPHI*XAET+SPHI*ZAET + +* Diurnal aberration + DIURAB = -AOPRMS(4) + F = (1D0-DIURAB*YMHDA) + V(1) = F*XMHDA + V(2) = F*(YMHDA+DIURAB) + V(3) = F*ZMHDA + +* To spherical -HA,Dec + CALL sla_DCC2S(V,HMA,DAP) + +* Right Ascension + RAP = sla_DRANRM(ST+HMA) + + END Index: trunk/psLib/src/astronomy/psAstronomyErrors.h =================================================================== --- trunk/psLib/src/astronomy/psAstronomyErrors.h (revision 3109) +++ trunk/psLib/src/astronomy/psAstronomyErrors.h (revision 3114) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-07 19:09:18 $ + * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-02-03 00:45:06 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -26,6 +26,4 @@ * DO NOT EDIT THE LINES BETWEEN //~Start and //~End! ANY CHANGES WILL BE OVERWRITTEN. */ - -#define PS_ERRORNAME_DOMAIN "psLib.astronomy." //~Start #define PS_ERRORTEXT_$1 "$2" Index: trunk/psLib/src/astronomy/psCoord.c =================================================================== --- trunk/psLib/src/astronomy/psCoord.c (revision 3109) +++ trunk/psLib/src/astronomy/psCoord.c (revision 3114) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-02 21:10:37 $ +* @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-03 00:45:06 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -599,6 +599,6 @@ XXX: How to compound errors? *****************************************************************************/ -psSphere* psSphereGetOffset(const psSphere* restrict position1, - const psSphere* restrict position2, +psSphere* psSphereGetOffset(const psSphere* position1, + const psSphere* position2, psSphereOffsetMode mode, psSphereOffsetUnit unit) @@ -683,6 +683,6 @@ *****************************************************************************/ -psSphere* psSphereSetOffset(const psSphere* restrict position, - const psSphere* restrict offset, +psSphere* psSphereSetOffset(const psSphere* position, + const psSphere* offset, psSphereOffsetMode mode, psSphereOffsetUnit unit) Index: trunk/psLib/src/astronomy/psCoord.h =================================================================== --- trunk/psLib/src/astronomy/psCoord.h (revision 3109) +++ trunk/psLib/src/astronomy/psCoord.h (revision 3114) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-01-26 20:24:16 $ +* @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-03 00:45:06 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -346,6 +346,6 @@ */ psSphere* psSphereGetOffset( - const psSphere* restrict position1, - const psSphere* restrict position2, + const psSphere* position1, + const psSphere* position2, psSphereOffsetMode mode, psSphereOffsetUnit unit @@ -364,6 +364,6 @@ */ psSphere* psSphereSetOffset( - const psSphere* restrict position, - const psSphere* restrict offset, + const psSphere* position, + const psSphere* offset, psSphereOffsetMode mode, psSphereOffsetUnit unit Index: trunk/psLib/src/astronomy/psMetadata.h =================================================================== --- trunk/psLib/src/astronomy/psMetadata.h (revision 3109) +++ trunk/psLib/src/astronomy/psMetadata.h (revision 3114) @@ -11,6 +11,6 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-01-17 20:58:20 $ +* @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-03 00:45:06 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -108,4 +108,5 @@ ); +#ifndef SWIG /** Create a metadata item with va_list. * @@ -130,4 +131,5 @@ va_list list ///< Arguments for name formatting and metadata item data. ); +#endif /** Create a metadata collection. Index: trunk/psLib/src/astronomy/psMetadataIO.c =================================================================== --- trunk/psLib/src/astronomy/psMetadataIO.c (revision 3109) +++ trunk/psLib/src/astronomy/psMetadataIO.c (revision 3114) @@ -9,6 +9,6 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-01-13 20:48:59 $ +* @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-03 00:45:06 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -760,6 +760,6 @@ psF64 tempDbl = 0.0; psS32 tempInt = 0.0; - psElemType pType; - psMetadataType mdType; + psElemType pType = PS_TYPE_PTR; + psMetadataType mdType = PS_META_PRIMITIVE; char *fileName = NULL; char *strName = NULL; @@ -909,6 +909,6 @@ psS32 status = 0; psU32 lineNumber = 0; - psElemType pType; - psMetadataType mdType; + psElemType pType = PS_TYPE_PTR; + psMetadataType mdType = PS_META_PRIMITIVE; char *strName = NULL; char *strType = NULL; Index: trunk/psLib/src/astronomy/psPhotometry.h =================================================================== --- trunk/psLib/src/astronomy/psPhotometry.h (revision 3109) +++ trunk/psLib/src/astronomy/psPhotometry.h (revision 3114) @@ -10,6 +10,6 @@ * @author George Gusciora, MHPCC * -* @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-10-27 00:57:30 $ +* @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-03 00:45:06 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -43,9 +43,9 @@ typedef struct { - const psS32 ID; ///< ID number for this photometric system - const char *name; ///< Name of photometric system - const char *camera; ///< Camera for photometric system - const char *filter; ///< Filter used for photometric system - const char *detector; ///< Detector used for photometric system + const psS32 ID; ///< ID number for this photometric system + const char *name; ///< Name of photometric system + const char *camera; ///< Camera for photometric system + const char *filter; ///< Filter used for photometric system + const char *detector; ///< Detector used for photometric system } psPhotSystem; @@ -59,13 +59,13 @@ typedef struct { - psPhotSystem src; ///< Source photometric system - psPhotSystem dst; ///< Destination photometric system - psPhotSystem pP; ///< Primary color reference - psPhotSystem pM; ///< Primary color reference - psPhotSystem sP; ///< Secondary color reference - psPhotSystem sM; ///< Secondary color reference - float pA; ///< Color offset for references - float sA; ///< Color offset for references - psPolynomial3D transform; ///< Transformation from source to destination + const psPhotSystem src; ///< Source photometric system + const psPhotSystem dst; ///< Destination photometric system + const psPhotSystem pP; ///< Primary color reference + const psPhotSystem pM; ///< Primary color reference + const psPhotSystem sP; ///< Secondary color reference + const psPhotSystem sM; ///< Secondary color reference + float pA; ///< Color offset for references + float sA; ///< Color offset for references + psPolynomial3D transform; ///< Transformation from source to destination } psPhotTransform; Index: trunk/psLib/src/astronomy/psTime.c =================================================================== --- trunk/psLib/src/astronomy/psTime.c (revision 3109) +++ trunk/psLib/src/astronomy/psTime.c (revision 3114) @@ -10,6 +10,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-17 20:47:08 $ + * @version $Revision: 1.50 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-02-03 00:45:06 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -35,8 +35,5 @@ #include "psAstronomyErrors.h" -#ifndef TIME_CONFIG_FILE -#define TIME_CONFIG_FILE "../../config/psTime.config" -#pragma warning TIME_CONFIG_FILE was not defined in the makefile. -#endif +#include "config.h" #define MAX_STRING_LENGTH 256 @@ -151,5 +148,5 @@ // Check time metadata. Function call reports errors. if(timeMetadata == NULL) { - if(!p_psTimeInit(TIME_CONFIG_FILE)) + if(!p_psTimeInit(CONFIG_FILE)) return 0.0; } @@ -708,5 +705,5 @@ // Check time metadata if(timeMetadata == NULL) { - if(!p_psTimeInit(TIME_CONFIG_FILE)) { + if(!p_psTimeInit(CONFIG_FILE)) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_FILE_NOT_FOUND, "psTime.config"); return 0.0; Index: trunk/psLib/src/astronomy/refco.f =================================================================== --- trunk/psLib/src/astronomy/refco.f (revision 3114) +++ trunk/psLib/src/astronomy/refco.f (revision 3114) @@ -0,0 +1,87 @@ + SUBROUTINE sla_REFCO (HM, TDK, PMB, RH, WL, PHI, TLR, EPS, + : REFA, REFB) +*+ +* - - - - - - +* R E F C O +* - - - - - - +* +* Determine the constants A and B in the atmospheric refraction +* model dZ = A tan Z + B tan**3 Z. +* +* Z is the "observed" zenith distance (i.e. affected by refraction) +* and dZ is what to add to Z to give the "topocentric" (i.e. in vacuo) +* zenith distance. +* +* Given: +* HM d height of the observer above sea level (metre) +* TDK d ambient temperature at the observer (deg K) +* PMB d pressure at the observer (millibar) +* RH d relative humidity at the observer (range 0-1) +* WL d effective wavelength of the source (micrometre) +* PHI d latitude of the observer (radian, astronomical) +* TLR d temperature lapse rate in the troposphere (degK/metre) +* EPS d precision required to terminate iteration (radian) +* +* Returned: +* REFA d tan Z coefficient (radian) +* REFB d tan**3 Z coefficient (radian) +* +* Called: sla_REFRO +* +* Notes: +* +* 1 Typical values for the TLR and EPS arguments might be 0.0065D0 and +* 1D-10 respectively. +* +* 2 The radio refraction is chosen by specifying WL > 100 micrometres. +* +* 3 The routine is a slower but more accurate alternative to the +* sla_REFCOQ routine. The constants it produces give perfect +* agreement with sla_REFRO at zenith distances arctan(1) (45 deg) +* and arctan(4) (about 76 deg). It achieves 0.5 arcsec accuracy +* for ZD < 80 deg, 0.01 arcsec accuracy for ZD < 60 deg, and +* 0.001 arcsec accuracy for ZD < 45 deg. +* +* P.T.Wallace Starlink 3 June 1997 +* +* Copyright (C) 1997 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION HM,TDK,PMB,RH,WL,PHI,TLR,EPS,REFA,REFB + + DOUBLE PRECISION ATN1,ATN4,R1,R2 + +* Sample zenith distances: arctan(1) and arctan(4) + PARAMETER (ATN1=0.7853981633974483D0, + : ATN4=1.325817663668033D0) + + + +* Determine refraction for the two sample zenith distances + CALL sla_REFRO(ATN1,HM,TDK,PMB,RH,WL,PHI,TLR,EPS,R1) + CALL sla_REFRO(ATN4,HM,TDK,PMB,RH,WL,PHI,TLR,EPS,R2) + +* Solve for refraction constants + REFA = (64D0*R1-R2)/60D0 + REFB = (R2-4D0*R1)/60D0 + + END Index: trunk/psLib/src/astronomy/refro.f =================================================================== --- trunk/psLib/src/astronomy/refro.f (revision 3114) +++ trunk/psLib/src/astronomy/refro.f (revision 3114) @@ -0,0 +1,401 @@ + SUBROUTINE sla_REFRO (ZOBS, HM, TDK, PMB, RH, WL, PHI, TLR, + : EPS, REF) +*+ +* - - - - - - +* R E F R O +* - - - - - - +* +* Atmospheric refraction for radio and optical/IR wavelengths. +* +* Given: +* ZOBS d observed zenith distance of the source (radian) +* HM d height of the observer above sea level (metre) +* TDK d ambient temperature at the observer (deg K) +* PMB d pressure at the observer (millibar) +* RH d relative humidity at the observer (range 0-1) +* WL d effective wavelength of the source (micrometre) +* PHI d latitude of the observer (radian, astronomical) +* TLR d temperature lapse rate in the troposphere (K/metre) +* EPS d precision required to terminate iteration (radian) +* +* Returned: +* REF d refraction: in vacuo ZD minus observed ZD (radian) +* +* Notes: +* +* 1 A suggested value for the TLR argument is 0.0065D0. The +* refraction is significantly affected by TLR, and if studies +* of the local atmosphere have been carried out a better TLR +* value may be available. The sign of the supplied TLR value +* is ignored. +* +* 2 A suggested value for the EPS argument is 1D-8. The result is +* usually at least two orders of magnitude more computationally +* precise than the supplied EPS value. +* +* 3 The routine computes the refraction for zenith distances up +* to and a little beyond 90 deg using the method of Hohenkerk +* and Sinclair (NAO Technical Notes 59 and 63, subsequently adopted +* in the Explanatory Supplement, 1992 edition - see section 3.281). +* +* 4 The code is a development of the optical/IR refraction subroutine +* AREF of C.Hohenkerk (HMNAO, September 1984), with extensions to +* support the radio case. Apart from merely cosmetic changes, the +* following modifications to the original HMNAO optical/IR refraction +* code have been made: +* +* . The angle arguments have been changed to radians. +* +* . Any value of ZOBS is allowed (see note 6, below). +* +* . Other argument values have been limited to safe values. +* +* . Murray's values for the gas constants have been used +* (Vectorial Astrometry, Adam Hilger, 1983). +* +* . The numerical integration phase has been rearranged for +* extra clarity. +* +* . A better model for Ps(T) has been adopted (taken from +* Gill, Atmosphere-Ocean Dynamics, Academic Press, 1982). +* +* . More accurate expressions for Pwo have been adopted +* (again from Gill 1982). +* +* . Provision for radio wavelengths has been added using +* expressions devised by A.T.Sinclair, RGO (private +* communication 1989). The refractivity model currently +* used is from J.M.Rueger, "Refractive Index Formulae for +* Electronic Distance Measurement with Radio and Millimetre +* Waves", in Unisurv Report S-68 (2002), School of Surveying +* and Spatial Information Systems, University of New South +* Wales, Sydney, Australia. +* +* . Various small changes have been made to gain speed. +* +* None of the changes significantly affects the optical/IR results +* with respect to the algorithm given in the 1992 Explanatory +* Supplement. For example, at 70 deg zenith distance the present +* routine agrees with the ES algorithm to better than 0.05 arcsec +* for any reasonable combination of parameters. However, the +* improved water-vapour expressions do make a significant difference +* in the radio band, at 70 deg zenith distance reaching almost +* 4 arcsec for a hot, humid, low-altitude site during a period of +* low pressure. +* +* 5 The radio refraction is chosen by specifying WL > 100 micrometres. +* Because the algorithm takes no account of the ionosphere, the +* accuracy deteriorates at low frequencies, below about 30 MHz. +* +* 6 Before use, the value of ZOBS is expressed in the range +/- pi. +* If this ranged ZOBS is -ve, the result REF is computed from its +* absolute value before being made -ve to match. In addition, if +* it has an absolute value greater than 93 deg, a fixed REF value +* equal to the result for ZOBS = 93 deg is returned, appropriately +* signed. +* +* 7 As in the original Hohenkerk and Sinclair algorithm, fixed values +* of the water vapour polytrope exponent, the height of the +* tropopause, and the height at which refraction is negligible are +* used. +* +* 8 The radio refraction has been tested against work done by +* Iain Coulson, JACH, (private communication 1995) for the +* James Clerk Maxwell Telescope, Mauna Kea. For typical conditions, +* agreement at the 0.1 arcsec level is achieved for moderate ZD, +* worsening to perhaps 0.5-1.0 arcsec at ZD 80 deg. At hot and +* humid sea-level sites the accuracy will not be as good. +* +* 9 It should be noted that the relative humidity RH is formally +* defined in terms of "mixing ratio" rather than pressures or +* densities as is often stated. It is the mass of water per unit +* mass of dry air divided by that for saturated air at the same +* temperature and pressure (see Gill 1982). +* +* 10 The algorithm is designed for observers in the troposphere. The +* supplied temperature, pressure and lapse rate are assumed to be +* for a point in the troposphere and are used to define a model +* atmosphere with the tropopause at 11km altitude and a constant +* temperature above that. However, in practice, the refraction +* values returned for stratospheric observers, at altitudes up to +* 25km, are quite usable. +* +* Called: sla_DRANGE, sla__ATMT, sla__ATMS +* +* P.T.Wallace Starlink 28 May 2002 +* +* Copyright (C) 2002 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION ZOBS,HM,TDK,PMB,RH,WL,PHI,TLR,EPS,REF + +* +* Fixed parameters +* + DOUBLE PRECISION D93,GCR,DMD,DMW,S,DELTA,HT,HS + INTEGER ISMAX +* 93 degrees in radians + PARAMETER (D93=1.623156204D0) +* Universal gas constant + PARAMETER (GCR=8314.32D0) +* Molecular weight of dry air + PARAMETER (DMD=28.9644D0) +* Molecular weight of water vapour + PARAMETER (DMW=18.0152D0) +* Mean Earth radius (metre) + PARAMETER (S=6378120D0) +* Exponent of temperature dependence of water vapour pressure + PARAMETER (DELTA=18.36D0) +* Height of tropopause (metre) + PARAMETER (HT=11000D0) +* Upper limit for refractive effects (metre) + PARAMETER (HS=80000D0) +* Numerical integration: maximum number of strips. + PARAMETER (ISMAX=16384) + + INTEGER IS,K,N,I,J + LOGICAL OPTIC,LOOP + DOUBLE PRECISION ZOBS1,ZOBS2,HMOK,TDKOK,PMBOK,RHOK,WLOK,ALPHA, + : TOL,WLSQ,GB,A,GAMAL,GAMMA,GAMM2,DELM2, + : TDC,PSAT,PWO,W, + : C1,C2,C3,C4,C5,C6,R0,TEMPO,DN0,RDNDR0,SK0,F0, + : RT,TT,DNT,RDNDRT,SINE,ZT,FT,DNTS,RDNDRP,ZTS,FTS, + : RS,DNS,RDNDRS,ZS,FS,REFOLD,Z0,ZRANGE,FB,FF,FO,FE, + : H,R,SZ,RG,DR,TG,DN,RDNDR,T,F,REFP,REFT + + DOUBLE PRECISION sla_DRANGE + +* The refraction integrand + DOUBLE PRECISION REFI + REFI(DN,RDNDR) = RDNDR/(DN+RDNDR) + + + +* Transform ZOBS into the normal range. + ZOBS1 = sla_DRANGE(ZOBS) + ZOBS2 = MIN(ABS(ZOBS1),D93) + +* Keep other arguments within safe bounds. + HMOK = MIN(MAX(HM,-1D3),HS) + TDKOK = MIN(MAX(TDK,100D0),500D0) + PMBOK = MIN(MAX(PMB,0D0),10000D0) + RHOK = MIN(MAX(RH,0D0),1D0) + WLOK = MAX(WL,0.1D0) + ALPHA = MIN(MAX(ABS(TLR),0.001D0),0.01D0) + +* Tolerance for iteration. + TOL = MIN(MAX(ABS(EPS),1D-12),0.1D0)/2D0 + +* Decide whether optical/IR or radio case - switch at 100 microns. + OPTIC = WLOK.LE.100D0 + +* Set up model atmosphere parameters defined at the observer. + WLSQ = WLOK*WLOK + GB = 9.784D0*(1D0-0.0026D0*COS(PHI+PHI)-0.00000028D0*HMOK) + IF (OPTIC) THEN + A = (287.604D0+(1.6288D0+0.0136D0/WLSQ)/WLSQ) + : *273.15D-6/1013.25D0 + ELSE + A = 77.6890D-6 + END IF + GAMAL = (GB*DMD)/GCR + GAMMA = GAMAL/ALPHA + GAMM2 = GAMMA-2D0 + DELM2 = DELTA-2D0 + TDC = TDKOK-273.15D0 + PSAT = 10D0**((0.7859D0+0.03477D0*TDC)/(1D0+0.00412D0*TDC))* + : (1D0+PMBOK*(4.5D-6+6D-10*TDC*TDC)) + IF (PMBOK.GT.0D0) THEN + PWO = RHOK*PSAT/(1D0-(1D0-RHOK)*PSAT/PMBOK) + ELSE + PWO = 0D0 + END IF + W = PWO*(1D0-DMW/DMD)*GAMMA/(DELTA-GAMMA) + C1 = A*(PMBOK+W)/TDKOK + IF (OPTIC) THEN + C2 = (A*W+11.2684D-6*PWO)/TDKOK + ELSE + C2 = (A*W+6.3938D-6*PWO)/TDKOK + END IF + C3 = (GAMMA-1D0)*ALPHA*C1/TDKOK + C4 = (DELTA-1D0)*ALPHA*C2/TDKOK + IF (OPTIC) THEN + C5 = 0D0 + C6 = 0D0 + ELSE + C5 = 375463D-6*PWO/TDKOK + C6 = C5*DELM2*ALPHA/(TDKOK*TDKOK) + END IF + +* Conditions at the observer. + R0 = S+HMOK + CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6, + : R0,TEMPO,DN0,RDNDR0) + SK0 = DN0*R0*SIN(ZOBS2) + F0 = REFI(DN0,RDNDR0) + +* Conditions in the troposphere at the tropopause. + RT = S+MAX(HT,HMOK) + CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6, + : RT,TT,DNT,RDNDRT) + SINE = SK0/(RT*DNT) + ZT = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0))) + FT = REFI(DNT,RDNDRT) + +* Conditions in the stratosphere at the tropopause. + CALL sla__ATMS(RT,TT,DNT,GAMAL,RT,DNTS,RDNDRP) + SINE = SK0/(RT*DNTS) + ZTS = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0))) + FTS = REFI(DNTS,RDNDRP) + +* Conditions at the stratosphere limit. + RS = S+HS + CALL sla__ATMS(RT,TT,DNT,GAMAL,RS,DNS,RDNDRS) + SINE = SK0/(RS*DNS) + ZS = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0))) + FS = REFI(DNS,RDNDRS) + +* +* Integrate the refraction integral in two parts; first in the +* troposphere (K=1), then in the stratosphere (K=2). +* + DO K = 1,2 + +* Initialize previous refraction to ensure at least two iterations. + REFOLD = 1D0 + +* Start off with 8 strips. + IS = 8 + +* Start Z, Z range, and start and end values. + IF (K.EQ.1) THEN + Z0 = ZOBS2 + ZRANGE = ZT-Z0 + FB = F0 + FF = FT + ELSE + Z0 = ZTS + ZRANGE = ZS-Z0 + FB = FTS + FF = FS + END IF + +* Sums of odd and even values. + FO = 0D0 + FE = 0D0 + +* First time through the loop we have to do every point. + N = 1 + +* Start of iteration loop (terminates at specified precision). + LOOP = .TRUE. + DO WHILE (LOOP) + +* Strip width. + H = ZRANGE/DBLE(IS) + +* Initialize distance from Earth centre for quadrature pass. + IF (K.EQ.1) THEN + R = R0 + ELSE + R = RT + END IF + +* One pass (no need to compute evens after first time). + DO I=1,IS-1,N + +* Sine of observed zenith distance. + SZ = SIN(Z0+H*DBLE(I)) + +* Find R (to the nearest metre, maximum four iterations). + IF (SZ.GT.1D-20) THEN + W = SK0/SZ + RG = R + DR = 1D6 + J = 0 + DO WHILE (ABS(DR).GT.1D0.AND.J.LT.4) + J=J+1 + IF (K.EQ.1) THEN + CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2, + : C1,C2,C3,C4,C5,C6,RG,TG,DN,RDNDR) + ELSE + CALL sla__ATMS(RT,TT,DNT,GAMAL,RG,DN,RDNDR) + END IF + DR = (RG*DN-W)/(DN+RDNDR) + RG = RG-DR + END DO + R = RG + END IF + +* Find the refractive index and integrand at R. + IF (K.EQ.1) THEN + CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2, + : C1,C2,C3,C4,C5,C6,R,T,DN,RDNDR) + ELSE + CALL sla__ATMS(RT,TT,DNT,GAMAL,R,DN,RDNDR) + END IF + F = REFI(DN,RDNDR) + +* Accumulate odd and (first time only) even values. + IF (N.EQ.1.AND.MOD(I,2).EQ.0) THEN + FE = FE+F + ELSE + FO = FO+F + END IF + END DO + +* Evaluate the integrand using Simpson's Rule. + REFP = H*(FB+4D0*FO+2D0*FE+FF)/3D0 + +* Has the required precision been achieved (or can't be)? + IF (ABS(REFP-REFOLD).GT.TOL.AND.IS.LT.ISMAX) THEN + +* No: prepare for next iteration. + +* Save current value for convergence test. + REFOLD = REFP + +* Double the number of strips. + IS = IS+IS + +* Sum of all current values = sum of next pass's even values. + FE = FE+FO + +* Prepare for new odd values. + FO = 0D0 + +* Skip even values next time. + N = 2 + ELSE + +* Yes: save troposphere component and terminate the loop. + IF (K.EQ.1) REFT = REFP + LOOP = .FALSE. + END IF + END DO + END DO + +* Result. + REF = REFT+REFP + IF (ZOBS1.LT.0D0) REF = -REF + + END Index: trunk/psLib/src/astronomy/refz.f =================================================================== --- trunk/psLib/src/astronomy/refz.f (revision 3114) +++ trunk/psLib/src/astronomy/refz.f (revision 3114) @@ -0,0 +1,156 @@ + SUBROUTINE sla_REFZ (ZU, REFA, REFB, ZR) +*+ +* - - - - - +* R E F Z +* - - - - - +* +* Adjust an unrefracted zenith distance to include the effect of +* atmospheric refraction, using the simple A tan Z + B tan**3 Z +* model (plus special handling for large ZDs). +* +* Given: +* ZU dp unrefracted zenith distance of the source (radian) +* REFA dp tan Z coefficient (radian) +* REFB dp tan**3 Z coefficient (radian) +* +* Returned: +* ZR dp refracted zenith distance (radian) +* +* Notes: +* +* 1 This routine applies the adjustment for refraction in the +* opposite sense to the usual one - it takes an unrefracted +* (in vacuo) position and produces an observed (refracted) +* position, whereas the A tan Z + B tan**3 Z model strictly +* applies to the case where an observed position is to have the +* refraction removed. The unrefracted to refracted case is +* harder, and requires an inverted form of the text-book +* refraction models; the formula used here is based on the +* Newton-Raphson method. For the utmost numerical consistency +* with the refracted to unrefracted model, two iterations are +* carried out, achieving agreement at the 1D-11 arcseconds level +* for a ZD of 80 degrees. The inherent accuracy of the model +* is, of course, far worse than this - see the documentation for +* sla_REFCO for more information. +* +* 2 At ZD 83 degrees, the rapidly-worsening A tan Z + B tan**3 Z +* model is abandoned and an empirical formula takes over. Over a +* wide range of observer heights and corresponding temperatures and +* pressures, the following levels of accuracy (arcsec) are +* typically achieved, relative to numerical integration through a +* model atmosphere: +* +* ZR error +* +* 80 0.4 +* 81 0.8 +* 82 1.5 +* 83 3.2 +* 84 4.9 +* 85 5.8 +* 86 6.1 +* 87 7.1 +* 88 10 +* 89 20 +* 90 40 +* 91 100 } relevant only to +* 92 200 } high-elevation sites +* +* The high-ZD model is scaled to match the normal model at the +* transition point; there is no glitch. +* +* 3 Beyond 93 deg zenith distance, the refraction is held at its +* 93 deg value. +* +* 4 See also the routine sla_REFV, which performs the adjustment in +* Cartesian Az/El coordinates, and with the emphasis on speed +* rather than numerical accuracy. +* +* P.T.Wallace Starlink 19 September 1995 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 59 Temple Place, Suite 330, +* Boston, MA 02111-1307 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION ZU,REFA,REFB,ZR + +* Radians to degrees + DOUBLE PRECISION R2D + PARAMETER (R2D=57.29577951308232D0) + +* Largest usable ZD (deg) + DOUBLE PRECISION D93 + PARAMETER (D93=93D0) + +* Coefficients for high ZD model (used beyond ZD 83 deg) + DOUBLE PRECISION C1,C2,C3,C4,C5 + PARAMETER (C1=+0.55445D0, + : C2=-0.01133D0, + : C3=+0.00202D0, + : C4=+0.28385D0, + : C5=+0.02390D0) + +* ZD at which one model hands over to the other (radians) + DOUBLE PRECISION Z83 + PARAMETER (Z83=83D0/R2D) + +* High-ZD-model prediction (deg) for that point + DOUBLE PRECISION REF83 + PARAMETER (REF83=(C1+C2*7D0+C3*49D0)/(1D0+C4*7D0+C5*49D0)) + + DOUBLE PRECISION ZU1,ZL,S,C,T,TSQ,TCU,REF,E,E2 + + + +* Perform calculations for ZU or 83 deg, whichever is smaller + ZU1 = MIN(ZU,Z83) + +* Functions of ZD + ZL = ZU1 + S = SIN(ZL) + C = COS(ZL) + T = S/C + TSQ = T*T + TCU = T*TSQ + +* Refracted ZD (mathematically to better than 1 mas at 70 deg) + ZL = ZL-(REFA*T+REFB*TCU)/(1D0+(REFA+3D0*REFB*TSQ)/(C*C)) + +* Further iteration + S = SIN(ZL) + C = COS(ZL) + T = S/C + TSQ = T*T + TCU = T*TSQ + REF = ZU1-ZL+ + : (ZL-ZU1+REFA*T+REFB*TCU)/(1D0+(REFA+3D0*REFB*TSQ)/(C*C)) + +* Special handling for large ZU + IF (ZU.GT.ZU1) THEN + E = 90D0-MIN(D93,ZU*R2D) + E2 = E*E + REF = (REF/REF83)*(C1+C2*E+C3*E2)/(1D0+C4*E+C5*E2) + END IF + +* Return refracted ZD + ZR = ZU-REF + + END