Index: /trunk/psLib/src/astronomy/Makefile.am =================================================================== --- /trunk/psLib/src/astronomy/Makefile.am (revision 4142) +++ /trunk/psLib/src/astronomy/Makefile.am (revision 4143) @@ -15,7 +15,5 @@ psTime.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 + psAstrometry.c BUILT_SOURCES = psAstronomyErrors.h @@ -31,5 +29,4 @@ psCoord.h \ psAstrometry.h \ - psPhotometry.h \ - slalib.h + psPhotometry.h Index: unk/psLib/src/astronomy/airmas.f =================================================================== --- /trunk/psLib/src/astronomy/airmas.f (revision 4142) +++ (revision ) @@ -1,75 +1,0 @@ - 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: unk/psLib/src/astronomy/aoppa.f =================================================================== --- /trunk/psLib/src/astronomy/aoppa.f (revision 4142) +++ (revision ) @@ -1,193 +1,0 @@ - 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: unk/psLib/src/astronomy/aoppat.f =================================================================== --- /trunk/psLib/src/astronomy/aoppat.f (revision 4142) +++ (revision ) @@ -1,62 +1,0 @@ - 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: unk/psLib/src/astronomy/aopqk.f =================================================================== --- /trunk/psLib/src/astronomy/aopqk.f (revision 4142) +++ (revision ) @@ -1,259 +1,0 @@ - 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: unk/psLib/src/astronomy/atms.f =================================================================== --- /trunk/psLib/src/astronomy/atms.f (revision 4142) +++ (revision ) @@ -1,57 +1,0 @@ - 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: unk/psLib/src/astronomy/atmt.f =================================================================== --- /trunk/psLib/src/astronomy/atmt.f (revision 4142) +++ (revision ) @@ -1,71 +1,0 @@ - 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: unk/psLib/src/astronomy/dcc2s.f =================================================================== --- /trunk/psLib/src/astronomy/dcc2s.f (revision 4142) +++ (revision ) @@ -1,70 +1,0 @@ - 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: unk/psLib/src/astronomy/dcs2c.f =================================================================== --- /trunk/psLib/src/astronomy/dcs2c.f (revision 4142) +++ (revision ) @@ -1,58 +1,0 @@ - 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: unk/psLib/src/astronomy/drange.f =================================================================== --- /trunk/psLib/src/astronomy/drange.f (revision 4142) +++ (revision ) @@ -1,49 +1,0 @@ - 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: unk/psLib/src/astronomy/dranrm.f =================================================================== --- /trunk/psLib/src/astronomy/dranrm.f (revision 4142) +++ (revision ) @@ -1,48 +1,0 @@ - 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: unk/psLib/src/astronomy/eqeqx.f =================================================================== --- /trunk/psLib/src/astronomy/eqeqx.f (revision 4142) +++ (revision ) @@ -1,74 +1,0 @@ - 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: unk/psLib/src/astronomy/geoc.f =================================================================== --- /trunk/psLib/src/astronomy/geoc.f (revision 4142) +++ (revision ) @@ -1,74 +1,0 @@ - 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: unk/psLib/src/astronomy/gmst.f =================================================================== --- /trunk/psLib/src/astronomy/gmst.f (revision 4142) +++ (revision ) @@ -1,77 +1,0 @@ - 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: unk/psLib/src/astronomy/nutc.f =================================================================== --- /trunk/psLib/src/astronomy/nutc.f (revision 4142) +++ (revision ) @@ -1,830 +1,0 @@ - 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: unk/psLib/src/astronomy/oapqk.f =================================================================== --- /trunk/psLib/src/astronomy/oapqk.f (revision 4142) +++ (revision ) @@ -1,250 +1,0 @@ - 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/psAstrometry.c =================================================================== --- /trunk/psLib/src/astronomy/psAstrometry.c (revision 4142) +++ /trunk/psLib/src/astronomy/psAstrometry.c (revision 4143) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.67 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-05-31 21:47:46 $ + * @version $Revision: 1.68 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-06-08 00:26:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -29,6 +29,4 @@ #include "psTrace.h" #include "psLogMsg.h" - -#include "slalib.h" /***************************************************************************** @@ -233,5 +231,5 @@ *(float *)&exp->positionAngle = 0.0f; // XXX: need input, see Bug #207 *(float *)&exp->parallacticAngle = 0.0f; // XXX: need input, see Bug #207 - *(float *)&exp->airmass = slaAirmas(zenithDistance); + *(float *)&exp->airmass = 0.0f; // XXX: needs calculation! = slaAirmas(zenithDistance); *(float *)&exp->parallacticFactor = 0.0f; exp->cameraName = NULL; @@ -408,25 +406,22 @@ PS_ASSERT_PTR_NON_NULL(exp, NULL); - double date = psTimeToMJD(exp->time); - double dut = psTimeGetUT1Delta(exp->time,PS_IERS_A); - double elongm = exp->observatory->longitude; - double phim = exp->observatory->latitude; - double hm = exp->observatory->height; - psSphere* polarMotion = p_psTimeGetPoleCoords(exp->time); - double xp = polarMotion->r; - double yp = polarMotion->d; + + psGrommit* grommit = (psGrommit* ) psAlloc(sizeof(psGrommit)); + + *(double*)&grommit->latitude = exp->observatory->latitude; + *(double*)&grommit->longitude = exp->observatory->longitude; + *(double*)&grommit->height = exp->observatory->height; + *(double*)&grommit->abberationMag = 0.0; // XXX: need to figure out what to set here. + *(double*)&grommit->temperature = exp->temperature; + *(double*)&grommit->pressure = exp->pressure; + *(double*)&grommit->humidity = exp->humidity; + *(double*)&grommit->wavelength = exp->wavelength; + *(double*)&grommit->lapseRate = exp->observatory->tlr; + *(double*)&grommit->refractA = polarMotion->r; // XXX: need to figure out what to set here too. + *(double*)&grommit->refractB = polarMotion->d; // XXX: need to figure out what to set here too. + *(double*)&grommit->siderealTime = psTimeToMJD(exp->time); // XXX: this is probably not correct + psFree(polarMotion); - - double tdk = exp->temperature; - double pmb = exp->pressure; - double rh = exp->humidity; - double wl = exp->wavelength; - double tlr = exp->observatory->tlr; - - psGrommit* grommit = (psGrommit* ) psAlloc(sizeof(psGrommit)); - - slaAoppa(date, dut, elongm, phim, hm, xp, yp, - tdk, pmb, rh, wl, tlr, (double*)grommit); return (grommit); @@ -579,13 +574,12 @@ PS_ASSERT_PTR_NON_NULL(grommit, NULL); - double AOB = 0.0; - double ZOB = 0.0; - double HOB = 0.0; if (outSphere == NULL) { outSphere = (psSphere* ) psAlloc(sizeof(psSphere)); } - slaAopqk(tpCoord->x, tpCoord->y, (double*)grommit, - &AOB, &ZOB, &HOB, &outSphere->r, &outSphere->d); + // XXX: this was done by a SLALIB call -- needs to be reimplemented + psWarning("Warning! psCoordTPToSky functionality is no longer implemented"); + /* slaAopqk(tpCoord->x, tpCoord->y, (double*)grommit, + &AOB, &ZOB, &HOB, &outSphere->r, &outSphere->d); */ return (outSphere); @@ -693,5 +687,5 @@ PS_ASSERT_PTR_NON_NULL(grommit, NULL); - char* type = "RA"; + // char* type = "RA"; if (tpCoord == NULL) { @@ -699,5 +693,7 @@ } - slaOapqk(type, in->r, in->d, (double*)grommit, &tpCoord->x, &tpCoord->y); + // XXX: this was done by a SLALIB call -- needs to be reimplemented + psWarning("Warning! psCoordSkyToTP functionality is no longer implemented"); + /* slaOapqk(type, in->r, in->d, (double*)grommit, &tpCoord->x, &tpCoord->y); */ return(tpCoord); Index: /trunk/psLib/src/astronomy/psAstrometry.h =================================================================== --- /trunk/psLib/src/astronomy/psAstrometry.h (revision 4142) +++ /trunk/psLib/src/astronomy/psAstrometry.h (revision 4143) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-03-31 23:01:46 $ +* @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-06-08 00:26:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -44,18 +44,15 @@ typedef struct { - const double latitude; ///< geodetic latitude (radians) - const double sinLat; ///< sine of geodetic latitude - const double cosLat; ///< cosine of geodetic latitude - const double abberationMag; ///< magnitude of diurnal aberration vector - const double height; ///< height (HM) - const double temperature; ///< ambient temperature (TDK) - const double pressure; ///< pressure (PMB) - const double humidity; ///< relative humidity (RH) - const double wavelength; ///< wavelength (WL) - const double lapseRate; ///< lapse rate (TLR) - const double refractA; ///< refraction constant A (radians) - const double refractB; ///< refraction constant B (radians) - const double longitudeOffset; ///< longitude + ... (radians) - const double siderealTime; ///< local apparent sidereal time (radians) + const double latitude; ///< geodetic latitude (radians) + const double longitude; ///< longitude + ... (radians) + const double height; ///< height (HM) + const double abberationMag; ///< magnitude of diurnal aberration vector + const double temperature; ///< ambient temperature (TDK) + const double pressure; ///< pressure (PMB) + const double humidity; ///< relative humidity (RH) + const double wavelength; ///< wavelength (WL) + const double lapseRate; ///< lapse rate (TLR) + const double refractA, refractB; ///< refraction constants A and B (radians) + const double siderealTime; ///< local apparent sidereal time (radians) } psGrommit; Index: unk/psLib/src/astronomy/refco.f =================================================================== --- /trunk/psLib/src/astronomy/refco.f (revision 4142) +++ (revision ) @@ -1,87 +1,0 @@ - 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: unk/psLib/src/astronomy/refro.f =================================================================== --- /trunk/psLib/src/astronomy/refro.f (revision 4142) +++ (revision ) @@ -1,401 +1,0 @@ - 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: unk/psLib/src/astronomy/refz.f =================================================================== --- /trunk/psLib/src/astronomy/refz.f (revision 4142) +++ (revision ) @@ -1,156 +1,0 @@ - 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 Index: unk/psLib/src/astronomy/slalib.h =================================================================== --- /trunk/psLib/src/astronomy/slalib.h (revision 4142) +++ (revision ) @@ -1,129 +1,0 @@ -/** @file slalib.h -* -* @brief Simple wrapper of the required fortran SLALIB functions -* -* @ingroup AstroImage -* -* @author Robert DeSonia, MHPCC -* -* @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-17 19:26:23 $ -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -*/ - -#ifndef SLALIB_H -#define SLALIB_H - -#include "config.h" - -#define fslaAoppa F77_FUNC_(sla_aoppa,SLA_AOPPA) -extern void fslaAoppa( - const double* date, - const double* dut, - const double* elongm, - const double* phim, - const double* him, - const double* xp, - const double* yp, - const double* tdk, - const double* pmb, - const double* rh, - const double* wl, - const double* tlr, - double* AOPRMS - ); - -#define fslaAopqk F77_FUNC_(sla_aopqk,SLA_AOPQK) -extern void fslaAopqk( - const double* RAP, - const double* DAP, - double* AOPRMS, - double* AOB, - double* ZOB, - double* HOB, - double* DOB, - double* ROB - ); - -#define fslaOapqk F77_FUNC_(sla_oapqk,SLA_OAPQK) -extern void fslaOapqk( - char* TYPE, - const double* OB1, - const double* OB2, - double* AOPRMS, - double* RAP, - double* DAP - ); - -#define fslaAirmas F77_FUNC_(sla_airmas,SLA_AIRMAS) -extern double fslaAirmas( - const double* zenithDistance - ); - -/* - void slaAoppa(date, dut, elongm, phim, hm, xp, yp, tdk, pmb, rh, wl, tlr, aoprms) - - PARAMETERS - - double date - double dut - double elongm - double phim - double hm - double xp - double yp - double tdk - double pmb - double rh - double wl - double tlr - double *aoprms -*/ -#define slaAoppa(date,dut,elongm,phim,him,xp,yp,tdk,pmb,rh,wl,tlr,AOPRMS) \ -fslaAoppa(&date,&dut,&elongm,&phim,&him,&xp,&yp,&tdk,&pmb,&rh,&wl,&tlr,AOPRMS) - -/* - void slaAopqk(rap, dap, aoprms, aob, zob, hob, dob, rob) - - PARAMETERS - - double rap - double dap - double *aoprms - double *aob - double *zob - double *hob - double *dob - double *rob -*/ - -#define slaAopqk(RAP,DAP,AOPRMS,AOB,ZOB,HOB,DOB,ROB) \ -fslaAopqk(&RAP,&DAP,AOPRMS,AOB,ZOB,HOB,DOB,ROB) - -/* - void slaOapqk(type, ob1, ob2, aoprms, rap, dap) - - PARAMETERS - - char *type - double ob1 - double ob2 - double *aoprms - double *rap - double *dap -*/ - -#define slaOapqk(TYPE, OB1, OB2, AOPRMS, RAP, DAP) \ -fslaOapqk(TYPE, &OB1, &OB2, AOPRMS, RAP, DAP) - -/* - double slaAirmas(zd) - - PARAMETERS - - double zd -*/ -#define slaAirmas(ZD) fslaAirmas(&ZD) - -#endif Index: /trunk/psLib/test/astronomy/tst_psAstrometry01.c =================================================================== --- /trunk/psLib/test/astronomy/tst_psAstrometry01.c (revision 4142) +++ /trunk/psLib/test/astronomy/tst_psAstrometry01.c (revision 4143) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-05-20 00:44:12 $ +* @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-06-08 00:26:01 $ * * XXX: Must test @@ -381,8 +381,7 @@ printf("grom->latitude is %.2f\n", grom->latitude); - printf("grom->sinLat is %.2f\n", grom->sinLat); - printf("grom->cosLat is %.2f\n", grom->cosLat); + printf("grom->longitude is %.2f\n", grom->longitude); + printf("grom->height is %.2f\n", grom->height); printf("grom->abberationMag is %.2f\n", grom->abberationMag); - printf("grom->height is %.2f\n", grom->height); printf("grom->temperature is %.2f\n", grom->temperature); printf("grom->pressure is %.2f\n", grom->pressure); @@ -392,5 +391,5 @@ printf("grom->refractA is %.2f\n", grom->refractA); printf("grom->refractB is %.2f\n", grom->refractB); - printf("grom->longitudeOffset is %.2f\n", grom->longitudeOffset); + // printf("grom->siderealTime is %.2f\n", grom->siderealTime);