C; Copyright (C) 1995, 1996 C; Associated Universities, Inc. Washington DC, USA. C; C; This program is free software; you can redistribute it and/or C; modify it under the terms of the GNU General Public License as C; published by the Free Software Foundation; either version 2 of C; the License, or (at your option) any later version. C; C; This program is distributed in the hope that it will be useful, C; but WITHOUT ANY WARRANTY; without even the implied warranty of C; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C; GNU General Public License for more details. C; C; You should have received a copy of the GNU General Public C; License along with this program; if not, write to the Free C; Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C; MA 02139, USA. C; C; Correspondence concerning AIPS should be addressed as follows: C; Internet email: aipsmail@nrao.edu. C; Postal address: AIPS Project Office C; National Radio Astronomy Observatory C; 520 Edgemont Road C; Charlottesville, VA 22903-2475 USA C----------------------------------------------------------------------- C C---------------------- Equinox conversion routines from AIPS ---------- C SUBROUTINE CRDSET (CRD1, CRD2, CRDPRM, IERR) C----------------------------------------------------------------------- C! Define coefficients for transforming between coordinate systems. C# Math Coordinates C----------------------------------------------------------------------- C----------------------------------------------------------------------- C CRDSET defines the coefficients for transforming between spherical C coordinate systems. C Inputs C CRD1 C*40 Coordinate system to transform from. Options C currently supported are listed below. C CRD2 C*40 Coordinate system to transform to. C Outputs: C CRDPRM D(11) Parameters for the coordinate transformation. C Euler angles for the transformation. C 1) (PHI0) Longitude of the ascending node in the C old system, in degrees. The ascending node C is the point of intersection of the equators C of the two systems such that the equator of C the new system crosses from south to north as C viewed in the old system. C 2) (THETA) The angle between the poles of the C two systems. THETA is positive for a C positive rotation about the ascending node. C 3) (PHI) Longitude of the ascending node in the C new system, in degrees. C 4) (CTHETA) = cos(THETA). C 5) (STHETA) = sin(THETA). C The elliptic terms of aberration which were C applied to pre-IAU1976 catalogue positions. C 6) C term C 7) D term C 8) C*TAN(eccentricity) C 9-11) similarly for the output system C IERR I Error code, 0: success C 1: unrecognized coordinate system. C Algorithm: C The Euler angles for the transformation from a variety of C coordinate systems to Equatorial B1950.0 are stored or computed. C To obtain the Euler angles to go from system X to system Y, those C from X to B1950.0 are combined with the inverse of those from Y C to B1950.0. C References: C Smith, C.A., et al, 1989. Astron. J., 97, 265. C Yallop, B.D., et al, 1989. Astron. J., 97, 274. C Notes: C 1) A failing of the IAU1976 system was not to give simple and C unambiguous names to the systems of constants and C coordinates. For example, the equatorial coordinate system C loosely described as "J2000.0" usually means "the position at C epoch J2000.0 on the mean equator and equinox of J2000.0 in C the IAU1976 system of coordinates". However, the J prefix C only indicates the new convention for computing epochs, there C is no formal reason to disassociate it from the old Bessel- C Newcomb system. C C Nevertheless, one is left to assume that any coordinate C system associated with a Julian epoch refers to the IAU1976 C system, and that Besselian epochs indicate Bessel-Newcomb C precession. An epoch prefix of 'b' is interpreted as C Bessel-Newcomb without E-terms. C C 2) References here to "equatorial B1950.0" or "ecliptic J2000.0" C etc. refer strictly to space-fixed coordinate systems. These C are defined by the mean equator (or ecliptic) and equinox C specified by the Bessel-Newcomb precession formulae for C Besselian epochs, or the IAU1976 system for Julian epochs. C C In transforming star catalogue positions it should be noted C that coordinates in the "FK4" system include the E-terms of C aberration (less than 1 arcsec), whereas the "FK5" system C excludes them. The exact transformation between these two C systems is therefore not a simple rotation. C C 3) The angles for transforming from J2000.0 to B1950.0 were C computed from the matrix coefficients Mij in the SLALIB C routine FK524 according to the following C PHI0 = ATAN2(-M31,M32) C THETA = ATAN2(STH,M33) C PHI = ATAN2(M13,-M23) C where STH is an average C STH = (M31/SIN(PHI0) - M13/SIN(PHI))/2D0 C C 4) Formulae for the precession angles are conventionally given C for precession from epoch E1 to epoch E2, referenced to the C basic epoch E0 (J2000.0 or B1900.0). This is avoided here, C since it results in precession coefficients which are not C independent of one another. Instead, only precession from C the basic epoch E0 to arbitrary epoch E is used. Precession C between any two arbitrary epochs is then given by the inverse C of precession from E0 to E1 followed by precession from E0 to C E2. C C 5) If the Euler angles for a rotation are (PHI0, THETA, PHI) C then for the inverse rotation they are (PHI, -THETA, PHI0). C Author: Mark Calabretta, Australia Telescope. C Origin; 1988/Oct/06. Code last modified; 1990/Aug/15 C----------------------------------------------------------------------- CHARACTER CRD1*40, CRD2*40 DOUBLE PRECISION CRDPRM(11) INTEGER IERR C # coordinate systems supported INTEGER NCRD PARAMETER (NCRD = 7) C INTEGER CTYP(2), FOUND, I, J, K, NDX1, NDX2 DOUBLE PRECISION CRDEUL(3,NCRD), D2R, ECCENT, EPS, EPOCH, * EUL(3,2), KABERR, OBLQTY, PERIGE, PHI, PHI0, PI, T, THETA, W1, * W2, W3 CHARACTER CRD*40, CRDNAM(NCRD)*40 C PARAMETER (PI = 3.141592653589793238462643D0) PARAMETER (D2R = PI/180D0) C Euler angles defining the C transformation from each C coordinate system to equatorial C b1950.0 (PHI0, THETA, PHI). DATA (CRDNAM(J), (CRDEUL(I,J), I=1,3), J=1,NCRD) * /'EQUATORIAL b1950.0', * 0.0D0, 0.0D0, 0.0D0, * 'EQUATORIAL b1900.0', * 89.6799436805556D0, +0.278397013888889D0, 90.3201112500000D0, * 'EQUATORIAL J2000.0', * 90.3202082924616D0, -0.278405865666317D0, 89.6795393829316D0, * 'GALACTIC', * 33.0D0, -62.6D0, -77.75D0, * 'OLD GALACTIC (OHLSSON)', * -0.05634D0, -62.27390D0, -79.38616D0, * 'OLD GALACTIC (VAN TULDER)', * -0.06146D0, -62.77299D0, -78.38808D0, * 'SUPERGALACTIC (REVISED)', * +116.86867D0, +74.35592D0, +193.18941D0/ C----------------------------------------------------------------------- C Anticipate and return C immediately on error IERR = 1 DO 10 J = 1,11 CRDPRM(J) = 0D0 10 CONTINUE CRDPRM(4) = 1D0 C Check for a unity transformation CALL MATCH (CRD1, CRD2, FOUND, NDX1, NDX2) IF ((NDX1.EQ.1) .AND. (NDX2.EQ.1)) THEN IERR = 0 GO TO 999 END IF C Identify coordinates DO 30 K = 1,2 IF (K.EQ.1) CRD = CRD1 IF (K.EQ.2) CRD = CRD2 C CTYP(K) = 0 DO 20 J = 3,NCRD CALL MATCH (CRD, CRDNAM(J), FOUND, NDX1, NDX2) IF (FOUND.GT.0) THEN CTYP(K) = J EUL(1,K) = CRDEUL(1,J) EUL(2,K) = CRDEUL(2,J) EUL(3,K) = CRDEUL(3,J) GO TO 30 END IF 20 CONTINUE C Test for equatorial coordinates CALL MATCH (CRD(1:10), 'EQUATORIAL', FOUND, NDX1, NDX2) IF (FOUND.GT.0) THEN READ (CRD(13:),1020,ERR=999) EPOCH IF ((CRD(12:12).EQ.'B') .OR. (CRD(12:12).EQ.'b')) THEN C Elliptic terms of aberration IF (CRD(12:12).EQ.'B') THEN T = 1.00002135903D0*(EPOCH - 1950D0)/100D0 KABERR = 20.49552/3600D0 W1 = (0.00004193D0 + 0.000000126D0*T)*T ECCENT = 0.01673011D0 - W1 W1 = (6190.67D0 + (1.65D0 + 0.012D0*T)*T)*T PERIGE = (282.080541944D0 + W1/3600D0)*D2R W1 = (46.8495D0 + (0.00319 - 0.00181D0*T)*T)*T OBLQTY = (23.445787778D0 - W1/3600D0)*D2R J = 5 + 3*(K-1) CRDPRM(J+1) = -KABERR*ECCENT*COS(PERIGE)*COS(OBLQTY) CRDPRM(J+2) = -KABERR*ECCENT*SIN(PERIGE) CRDPRM(J+3) = CRDPRM(J+1)*TAN(OBLQTY) END IF C Bessel-Newcomb precession from C B1900.0 to epoch T = (EPOCH - 1900D0)/100D0 W1 = (2304.25D0 + (0.302D0 + 0.018D0*T)*T)*T W2 = (2004.682D0 - (0.426D0 + 0.042D0*T)*T)*T W3 = (2304.25D0 + (1.093D0 + 0.018D0*T)*T)*T PHI0 = 90D0 - W1/3600D0 THETA = W2/3600D0 PHI = 90D0 + W3/3600D0 C Epoch to B1900.0, followed by C B1900.0 to B1950.0 CALL XEULER (PHI, -THETA, PHI0, CRDEUL(1,2), CRDEUL(2,2), * CRDEUL(3,2), EUL(1,K), EUL(2,K), EUL(3,K)) C IAU1976 precession from J2000.0 C to epoch ELSE IF (CRD(12:12).EQ.'J') THEN T = (EPOCH - 2000D0)/100D0 W1 = (2306.2181D0 + (0.30188D0 + 0.017998D0*T)*T)*T W2 = (2004.3109D0 - (0.42665D0 + 0.041833D0*T)*T)*T W3 = (2306.2181D0 + (1.09468D0 + 0.018203D0*T)*T)*T PHI0 = 90D0 - W1/3600D0 THETA = W2/3600D0 PHI = 90D0 + W3/3600D0 C Epoch to J2000.0, followed by C J2000.0 to B1950.0 CALL XEULER (PHI, -THETA, PHI0, CRDEUL(1,3), CRDEUL(2,3), * CRDEUL(3,3), EUL(1,K), EUL(2,K), EUL(3,K)) C Epoch prefix not recognized ELSE GO TO 999 END IF GO TO 30 END IF C Test for ecliptic coordinates CALL MATCH (CRD(1:8), 'ECLIPTIC', FOUND, NDX1, NDX2) IF (FOUND.GT.0) THEN READ (CRD(11:),1020,ERR=999) EPOCH T = EPOCH/100D0 IF ((CRD(10:10).EQ.'B') .OR. (CRD(10:10).EQ.'b')) THEN C Bessel-Newcomb mean obliquity of C the ecliptic T = (EPOCH - 1900D0)/100D0 EPS = (84428D0 - 47D0*T)/3600D0 C Bessel-Newcomb precession from C B1900.0 to epoch W1 = (2304.25D0 + (0.302D0 + 0.018D0*T)*T)*T W2 = (2004.682D0 - (0.426D0 + 0.042D0*T)*T)*T W3 = (2304.25D0 + (1.093D0 + 0.018D0*T)*T)*T PHI0 = 90D0 - W1/3600D0 THETA = W2/3600D0 PHI = 90D0 + W3/3600D0 C Ecliptic to equatorial, followed C by epoch to B1900.0... CALL XEULER (0D0, -EPS, 0D0, PHI, -THETA, PHI0, EUL(1,K), * EUL(2,K), EUL(3,K)) C ...followed by B1900.0 to C B1950.0 CALL XEULER (EUL(1,K), EUL(2,K), EUL(3,K), CRDEUL(1,2), * CRDEUL(2,2), CRDEUL(3,2), EUL(1,K), EUL(2,K), * EUL(3,K)) ELSE IF (CRD(10:10).EQ.'J') THEN C IAU1976 mean obliquity of the C ecliptic. T = (EPOCH - 2000D0)/100D0 W1 = (46.8150D0 + (0.00059D0 - 0.001813D0*T)*T)*T EPS = (84381.448D0 - W1)/3600D0 C IAU1976 precession from J2000.0 C to epoch W1 = (2306.2181D0 + (0.30188D0 + 0.017998D0*T)*T)*T W2 = (2004.3109D0 - (0.42665D0 + 0.041833D0*T)*T)*T W3 = (2306.2181D0 + (1.09468D0 + 0.018203D0*T)*T)*T PHI0 = 90D0 - W1/3600D0 THETA = W2/3600D0 PHI = 90D0 + W3/3600D0 C Ecliptic to equatorial, followed C by epoch to J2000.0... CALL XEULER (0D0, -EPS, 0D0, PHI, -THETA, PHI0, EUL(1,K), * EUL(2,K), EUL(3,K)) C ...followed by J2000.0 to C B1950.0 CALL XEULER (EUL(1,K), EUL(2,K), EUL(3,K), CRDEUL(1,3), * CRDEUL(2,3), CRDEUL(3,3), EUL(1,K), EUL(2,K), * EUL(3,K)) C Epoch prefix not recognized ELSE GO TO 999 END IF GO TO 30 END IF C Coordinate system not recognized GO TO 999 30 CONTINUE C Determine the Euler angles C Transformation from equatorial C B1950.0 IF (CTYP(1).EQ.1) THEN PHI0 = EUL(3,2) THETA = -EUL(2,2) PHI = EUL(1,2) C Transformation to equatorial C B1950.0 ELSE IF (CTYP(2).EQ.1) THEN PHI0 = EUL(1,1) THETA = EUL(2,1) PHI = EUL(3,1) ELSE CALL XEULER (EUL(1,1), EUL(2,1), EUL(3,1), EUL(3,2), -EUL(2,2), * EUL(1,2), PHI0, THETA, PHI) END IF C CRDPRM(1) = PHI0 CRDPRM(2) = THETA CRDPRM(3) = PHI CRDPRM(4) = COS(THETA*D2R) CRDPRM(5) = SIN(THETA*D2R) C Success IERR = 0 C 999 RETURN C----------------------------------------------------------------------- 1020 FORMAT (F9.4) END SUBROUTINE CRDTRN (AZ1, POL1, CRDPRM, AZ2, POL2, ROTN) C----------------------------------------------------------------------- C! Apply Euler angle based coordinate transformation. C# Math Coordinates C----------------------------------------------------------------------- C CRDTRN applies the Euler angle based coordinate transformation. C Inputs: C AZ1 D Azimuthal angle in the old system (degrees) C POL1 D Polar angle in old system (degrees) C CRDPRM D(11) Parameters for the coordinate transformation. C Euler angles for the transformation. C 1) (PHI0) Longitude of the ascending node in the C old system, in degrees. The ascending node is C the point of intersection of the equators of the C two systems such that the equator of the new C system crosses from south to north as viewed in C the old system. C 2) (THETA) The angle between the poles of the C two systems. THETA is positive for a positive C rotation about the ascending node. C 3) (PHI) Longitude of the ascending node in the C new system, in degrees. C 4) (CTHETA) = cos(THETA). C 5) (STHETA) = sin(THETA). C The elliptic terms of aberration which were C applied to pre-IAU1976 catalogue positions. C 6) C term C 7) D term C 8) C*TAN(eccentricity) C 9-11) similarly for the output system C Outputs: C AZ2 D Azimuthal angle in the new system (degrees) C POL2 D Polar angle in new system (degrees) C ROTN D The position angle of the old pole at the C specified point, in degrees. C Called: EULROT C Algorithm: C A coordinate rotation specified by the three Euler angles. C The elliptic terms of aberration are applied as necessary, C Refer to the precursor comments to CRDSET. C Notes: C 1) Longitude at the poles in the new system is consistent with that C specified in the old system. This may be important when dealing C with map projections in which the poles are represented by finite C line segments. Such is the case for cylindrical projections for C example. C 2) Subroutine CRDSET may be called to set the transformation C parameters for any of a variety of transformations. C Author: Mark Calabretta, Australia Telescope. C Origin; 1988/Oct/05. Code last modified; 1990/Aug/15 C----------------------------------------------------------------------- DOUBLE PRECISION AZ1, AZ2, POL1, POL2, CRDPRM(11), ROTN C DOUBLE PRECISION CAZ1, CAZ2, CPOL1, CPOL2, D2R, EAZ1, EAZ2, * EPOL1, EPOL2, PI, SAZ1, SAZ2, SPOL1, SPOL2 C PARAMETER (PI = 3.141592653589793238462643D0) PARAMETER (D2R = PI/180D0) C----------------------------------------------------------------------- C Remove the E-terms (non-zero C only for FK4) CAZ1 = COS (AZ1*D2R) SAZ1 = SIN (AZ1*D2R) CPOL1 = COS (POL1*D2R) SPOL1 = SIN (POL1*D2R) IF (CPOL1.NE.0D0) THEN EAZ1 = (CRDPRM(6) * CAZ1 + CRDPRM(7) * SAZ1) / CPOL1 ELSE EAZ1 = 0D0 END IF EPOL1 = (CRDPRM(7) * CAZ1 - CRDPRM(6) * SAZ1) * SPOL1 + * CRDPRM(8) * CPOL1 C Spherical coordinate rotation CALL EULROT (AZ1-EAZ1, POL1-EPOL1, CRDPRM, AZ2, POL2, ROTN) C Apply the E-terms (non-zero only C for FK4) CAZ2 = COS (AZ2*D2R) SAZ2 = SIN (AZ2*D2R) CPOL2 = COS (POL2*D2R) SPOL2 = SIN (POL2*D2R) IF (CPOL2.NE.0D0) THEN EAZ2 = (CRDPRM(9) * CAZ2 + CRDPRM(10) * SAZ2) / CPOL2 ELSE EAZ2 = 0D0 END IF EPOL2 = (CRDPRM(10) * CAZ2 - CRDPRM(9) * SAZ2) * SPOL2 + * CRDPRM(11) * CPOL2 C AZ2 = AZ2 + EAZ2 POL2 = POL2 + EPOL2 C Normalize AZ2 = MOD (AZ2, 360D0) IF (AZ2.LT.0D0) AZ2 = AZ2 + 360D0 C 999 RETURN END SUBROUTINE MATCH (TEXT1, TEXT2, FOUND, N1, N2) C----------------------------------------------------------------------- C! Tests whether either of two strings is wholly contained in the other. C# Character Parsing C----------------------------------------------------------------------- C tests whether either of two strings is wholly contained in the other C Inputs: C TEXT1 C*(*) Input character strings. Length should not C TEXT2 C*(*) exceed 132 characters. C Outputs: C FOUND I Match results. 0: no match, C 1: case insensitive match. C 2: case sensitive match, C N1 I > 0 => the starting character position of TEXT2 C in TEXT1. C N2 I > 0 => the starting character position of TEXT1 C in TEXT2. C Author: Mark Calabretta, Australia Telescope. C Origin; 1988/Apr/28. Code last modified; 1989/Sep/05. C----------------------------------------------------------------------- INTEGER FOUND, N1, N2 CHARACTER TEXT1*(*), TEXT2*(*) C INTEGER I1, I2, J1, J2 CHARACTER STR1*132, STR2*132 C----------------------------------------------------------------------- C Find the first and last non- C blank, non-tab characters STR1 = TEXT1 STR2 = TEXT2 CALL TXTLEN (STR1, I1, I2) CALL TXTLEN (STR2, J1, J2) C Try for a case sensitive match FOUND = 2 N1 = INDEX (STR1(1:I2), STR2(J1:J2)) N2 = INDEX (STR2(1:J2), STR1(I1:I2)) C Perhaps a case insensitive match IF ((N1.LE.0) .AND. (N2.LE.0)) THEN FOUND = 1 CALL CHLTOU (132, STR1) CALL CHLTOU (132, STR2) N1 = INDEX (STR1(1:I2), STR2(J1:J2)) N2 = INDEX (STR2(1:J2), STR1(I1:I2)) END IF C No match IF ((N1.LE.0) .AND. (N2.LE.0)) FOUND = 0 C 999 RETURN END SUBROUTINE XEULER (P01, TH1, PH1, P02, TH2, PH2, PHI0, THETA, PHI) C----------------------------------------------------------------------- C! Combine the Euler angles for two consecutive rotations into one. C# Math Coordinates C----------------------------------------------------------------------- C XEULER combines the Euler angles for two consecutive rotations into C one. C Inputs: C P01 D The three Euler angles for the first C TH1 D rotation, in degrees. C PH1 D C P02 D The three Euler angles for the second C TH2 D rotation, in degrees. C PH2 D C Output: C PHI0 D The three Euler angles for the combined C THETA D rotation, in degrees. C PHI D C Called: EULROT C Algorithm: C The two rotations are applied successively to special points in C the original coordinate frame. In any rotation the north pole C will be transformed to polar angle 90 - THETA, and thus THETA, C and likewise PHI, for the combined rotation can be found easily. C Determination of PHI0 is a little trickier, and requires the C coordinates in the original system of the north pole of the new C frame. However, it is straightforward to express the inverse of C the combined rotation in terms of the inverse of the two C constituent rotations. C Notes: C 1) Note that rotations are not commutative. Therefore, the C order of the two rotations is important. C C 2) If the Euler angles for a rotation are (PHI0, THETA, PHI) C then for the inverse rotation they are (PHI, -THETA, PHI0). C C 3) If R1, R2, and R are the rotations, and r1, r2, and r are C their inverses then C R(x) = R2(R1(x)) and r(x) = r1(r2(x)). C Author: Mark Calabretta, Australia Telescope. C Origin; 1988/Oct/05, Code last modified; 1990/Jul/25 C----------------------------------------------------------------------- DOUBLE PRECISION P01, TH1, PH1, P02, TH2, PH2, PHI0, THETA, PHI C DOUBLE PRECISION AZ0, AZ1, AZ2, D2R, EULER1(5), EULER2(5), PI, * POL0, POL1, POL2, ROTN PARAMETER (PI = 3.141592653589793238462643D0) PARAMETER (D2R = PI/180D0) C----------------------------------------------------------------------- C Apply the two rotations in C succession to the north pole EULER1(1) = P01 EULER1(2) = TH1 EULER1(3) = PH1 EULER1(4) = COS (TH1*D2R) EULER1(5) = SIN (TH1*D2R) CALL EULROT (0D0, 90D0, EULER1, AZ1, POL1, ROTN) C EULER2(1) = P02 EULER2(2) = TH2 EULER2(3) = PH2 EULER2(4) = COS (TH2*D2R) EULER2(5) = SIN (TH2*D2R) CALL EULROT (AZ1, POL1, EULER2, AZ2, POL2, ROTN) C THETA is easily found THETA = 90D0 - POL2 C Now determine PHI and PHI0 C THETA = zero is special case IF (SIN(THETA*D2R).LT.1E-9) THEN CALL EULROT (0D0, 0D0, EULER1, AZ1, POL1, ROTN) CALL EULROT (AZ1, POL1, EULER2, AZ2, POL2, ROTN) PHI0 = 0D0 PHI = AZ2 C Otherwise PHI is obtained in the C same way as THETA ELSE PHI = AZ2 - 90D0 IF (PHI.LT.0D0) PHI = PHI + 360D0 C The easiest way to find PHI0 is C to apply the inverse C transformation to the north pole C of the new coordinate frame EULER2(1) = PH2 EULER2(2) = -TH2 EULER2(3) = P02 EULER2(4) = COS (TH2*D2R) EULER2(5) = -SIN (TH2*D2R) CALL EULROT (0D0, 90D0, EULER2, AZ1, POL1, ROTN) C EULER1(1) = PH1 EULER1(2) = -TH1 EULER1(3) = P01 EULER1(4) = COS (TH1*D2R) EULER1(5) = -SIN (TH1*D2R) CALL EULROT (AZ1, POL1, EULER1, AZ0, POL0, ROTN) C PHI0 = AZ0 + 90D0 IF (PHI0.GE.360D0) PHI0 = PHI0 - 360D0 END IF C 999 RETURN END SUBROUTINE EULROT (AZ1, POL1, EULER, AZ2, POL2, ROTN) C----------------------------------------------------------------------- C! Apply Euler angle based rotation of spherical coordinates. C# Math Coordinates C----------------------------------------------------------------------- C EULROT applies the Euler angle based rotation of spherical C coordinates. C Inputs C AZ1 D Azimuthal and polar angles in the old C POL1 D coordinate system, in degrees. C EULER D(5) Euler angles for the rotation. C 1) (PHI0) Longitude of the ascending node in the C old system, in degrees. The ascending node is C the point of intersection of the equators of C the two systems such that the equator of the C new system crosses from south to north as C viewed in the old system. C 2) (THETA) The angle between the poles of the two C systems. THETA is positive for a positive C rotation about the ascending node. C 3) (PHI) Longitude of the ascending node in the C new system, in degrees. C 4) (CTHETA) = cos(THETA). C 5) (STHETA) = sin(THETA). C Outputs; C AZ2 D Azimuthal and polar angles in the new coordinate C POL2 D system, in degrees. AZ2 lies between 0 and 360 C ROTN D The position angle of the old pole at the C specified point, in degrees. C Algorithm: C A coordinate rotation specified by the three Euler angles. C Notes: C 1) Longitude at the poles in the new system is consistent with C that specified in the old system. This may be important C when dealing with map projections in which the poles are C represented by finite line segments. Such is the case for C cylindrical projections for example. C Author: Mark Calabretta, Australia Telescope. C Origin; 1988/Oct/05. Code last modified; 1990/Aug/15 C----------------------------------------------------------------------- DOUBLE PRECISION AZ1, POL1, EULER(5), AZ2, POL2, ROTN C DOUBLE PRECISION AZ1P, CAZ1P, CPOL1, CTHETA, D2R, PHI0, PHI, PI, * R2D, SAZ1P, SPOL1, STHETA, THETA, X, Y PARAMETER (PI = 3.141592653589793238462643D0) PARAMETER (D2R = PI/180D0) PARAMETER (R2D = 180D0/PI) C----------------------------------------------------------------------- C The following assignments are C made for convenience PHI0 = EULER(1) THETA = EULER(2) PHI = EULER(3) CTHETA = EULER(4) STHETA = EULER(5) C Compute trigonometric functions C of the old coordinates AZ1P = AZ1 - PHI0 CAZ1P = COS (AZ1P*D2R) SAZ1P = SIN (AZ1P*D2R) CPOL1 = COS (POL1*D2R) SPOL1 = SIN (POL1*D2R) C Compute the azimuthal angle in C the new system X = CPOL1 * CAZ1P Y = CPOL1 * SAZ1P *CTHETA + SPOL1 * STHETA IF ((X.NE.0D0) .OR. (Y.NE.0D0)) THEN AZ2 = PHI + ATAN2 (Y, X) * R2D C Longitude at the poles in the C new system is consistent with C that specified in the old system ELSE AZ2 = PHI + AZ1P END IF C Compute the polar angle in the C new system POL2 = ASIN (SPOL1*CTHETA - CPOL1*STHETA*SAZ1P) * R2D C Compute the position angle of C the old pole at this point X = SAZ1P * STHETA * SPOL1 + CTHETA * CPOL1 Y = CAZ1P * STHETA ROTN = ATAN2 (Y, X) * R2D C Normalize AZ2 = MOD (AZ2, 360D0) IF (AZ2.LT.0D0) AZ2 = AZ2 + 360D0 C 999 RETURN END C C-------------- Miscellaneous routines from AIPS --------------------- C SUBROUTINE GETNUM (KB, KBPLIM, KBP, X) C----------------------------------------------------------------------- C! returns numeric field from character buffer C# Parsing FITS POPS-util C----------------------------------------------------------------------- C GETNUM converts ASCII numeric field into double precision number. C Inputs: C KB C*(*) character buffer C KBPLIM I max character position in buffer to examine C KBP I start of numeric field C Outputs: C KBP I start of next field (incl blanks) C X D numerical value: sets to magic indef = DBLANK C from DDCH.INC when overflow exponent or when C there are no numeric characters found C Common: C DERR.INC Sets ERRNUM to 27 on failures C Modified by W. D. Cotton for NVSSlist C----------------------------------------------------------------------- CHARACTER KB*(*) INTEGER KBP, KBPLIM DOUBLE PRECISION X C INTEGER EXSIGN, IEXP, IT, LEADER, IDEC, NC DOUBLE PRECISION SIGN, TEN C INCLUDE 'INCS:DERR.INC' C INCLUDE 'INCS:DDCH.INC' C----------------------------------------------------------------------- X = 0.0D0 IF (KBP.GT.KBPLIM) GO TO 980 LEADER = 0 IDEC = 0 SIGN = 1.0D0 EXSIGN = 1 TEN = 10.0D0 IEXP = 0 NC = 0 C skip leading blank,comma,+ 10 IF ((KB(KBP:KBP).NE.' ') .AND. (KB(KBP:KBP).NE.',') .AND. * (KB(KBP:KBP).NE.'+')) GO TO 15 KBP = KBP + 1 IF (KBP.GT.KBPLIM) GO TO 999 GO TO 10 15 KBP = KBP - 1 C get chars 20 KBP = KBP + 1 IF (KBP.GT.KBPLIM) GO TO 90 IF (KB(KBP:KBP).NE.'.') GO TO 25 IDEC = 1 GO TO 20 C unary minus 25 IF (KB(KBP:KBP).NE.'-') GO TO 30 IF (LEADER.NE.0) GO TO 40 SIGN = -SIGN GO TO 20 C other: leave if not num 30 CONTINUE LEADER = 1 IT = ICHAR(KB(KBP:KBP)) - ICHAR('0') IF ((IT.LT.0) .OR. (IT.GT.9)) GO TO 40 NC = NC + 1 X = X * TEN + IT IF (IDEC.NE.0) IDEC = IDEC + 1 GO TO 20 C exponential notation E or D C then (optional) blank,+,- C then integer number 40 IF ((KB(KBP:KBP).NE.'E') .AND. (KB(KBP:KBP).NE.'D') .AND. * (KB(KBP:KBP).NE.'e') .AND. (KB(KBP:KBP).NE.'d')) GO TO 90 KBP = KBP + 1 IF (KBP.GT.KBPLIM) GO TO 90 IF ((KB(KBP:KBP).NE.'-') .AND. (KB(KBP:KBP).NE.'+') .AND. * (KB(KBP:KBP).NE.' ')) GO TO 50 IF (KB(KBP:KBP).EQ.'-') EXSIGN = -1 45 KBP = KBP + 1 IF (KBP.GT.KBPLIM) GO TO 90 50 IT = ICHAR(KB(KBP:KBP)) - ICHAR('0') IF ((IT.LT.0) .OR. (IT.GT.9)) GO TO 90 IEXP = 10*IEXP + IT NC = NC + 1 GO TO 45 C got it: apply sign, decimal pl 90 IF (NC.LE.0) GO TO 980 IDEC = MAX (0, IDEC-1) IDEC = EXSIGN*IEXP - IDEC IF (IDEC.GT.18) GO TO 980 IF (IDEC.GT.-34) THEN X = X * SIGN / (TEN**(-IDEC)) ELSE X = 0.0 END IF GO TO 999 C error 980 Continue C ERRNUM = 27 C ERRLEV = ERRLEV + 1 C IF (ERRLEV.LE.10) PNAME(ERRLEV) = 'GETNUM' C X = DBLANK C 999 RETURN END INTEGER FUNCTION ITRIM (STRING) C----------------------------------------------------------------------- C! returns length of CHARACTER variable to last non-blank, non-carrage C! return C# Character Utility C----------------------------------------------------------------------- C Function to determine length of a string. I.e., it trims trailing C blanks. Use with calls like: C TRIMMED = GROSS(1:ITRIM(GROSS)) C NOTE: this does not check for NULL characters, use JTRIM for that C----------------------------------------------------------------------- CHARACTER STRING*(*), CR*1 C----------------------------------------------------------------------- C Carrage return CR = CHAR(13) ITRIM = LEN (STRING) + 1 C look backwards for non-blank 10 CONTINUE ITRIM = ITRIM - 1 IF (ITRIM.LT.1) GO TO 999 IF (STRING(ITRIM:ITRIM).EQ.' ') GO TO 10 IF (STRING(ITRIM:ITRIM).EQ.CR) GO TO 10 C 999 RETURN END SUBROUTINE DIRCOS (TYPE, RA0, DEC0, RA, DEC, L, M, IERR) C----------------------------------------------------------------------- C! determines direction cosines between ref position and test position C# Coordinates C----------------------------------------------------------------------- C This routine only supports type 2 C DIRCOS determines the direction cosines (L,M) corresponding to C the separation from (RA0,DEC0), the reference pixel, to C (RA,DEC). The direction cosine L is assumed to be positive to the C east; M is positive to the north. Several projective geometries C are supported. C Inputs: C TYPE I Type of projection: 2 => SIN, 3 => TAN, 4 => Arc, C 5 => WSRT (rel north pole), 6 Global sinusoidal, C 7 Mercator, 8 Aitoff C RA0 D Coordinate reference right ascension/longitude C DEC0 D Coordinate reference declination/latitude C RA D right ascension/longitude of point C DEC D declination/latitude of point C Outputs: C L D Cosine angle of displacement to east C M D Cosine angle of displacement to north C IERR I 0 ok, 1 out of range, 2 bad type, 3 undefined C ALL ANGLES ARE IN RADIANS. C----------------------------------------------------------------------- INTEGER IERR, TYPE DOUBLE PRECISION RA0, DEC0, RA, DEC, L, M C DOUBLE PRECISION COSS, SINS, DT, SINT DOUBLE PRECISION PI, TWOPI PARAMETER (PI = 3.141592654) C INCLUDE 'INCS:DMSG.INC' C INCLUDE 'INCS:DLOC.INC' C INCLUDE 'INCS:PSTD.INC' C DATA DEPS /1.0D-5/ C----------------------------------------------------------------------- TWOPI = 2.0D0 * PI C Use full accuracy 10 IERR = 2 IF ((TYPE.LT.2) .OR. (TYPE.GT.2)) GO TO 999 IF (ABS(DEC).GT.PI/2.0D0) GO TO 990 IF (ABS(DEC0).GT.PI/2.0D0) GO TO 990 IERR = 0 COSS = COS (DEC) SINS = SIN (DEC) DT = RA - RA0 IF (DT.GT.PI) DT = DT - TWOPI IF (DT.LT.-PI) DT = DT + TWOPI L = SIN(DT) * COSS SINT = SINS * SIN(DEC0) + COSS * COS(DEC0) * COS(DT) C SIN projection IF (SINT.LT.0.0D0) IERR = 1 M = SINS * COS(DEC0) - COSS * SIN(DEC0) * COS(DT) GO TO 999 C Undefined 990 IERR = 3 C 999 RETURN C----------------------------------------------------------------------- C 1050 FORMAT ('DIRCOS: WSRT COORDINATES CANNOT HAVE DEC0 = 0') END SUBROUTINE COORDD (N, RAD, CHM, HM, SEC) C----------------------------------------------------------------------- C! converts angles between degrees and sexagesimal format C# Utility Coordinates C----------------------------------------------------------------------- C COORDD Converts RA and DEC to and from degrees under control of N. C Inputs: C N I OPcode: 1 degrees to RA C 2 degrees to DEC C 3 RA to degrees C 4 DEC to degrees C In/out: C RAD D The angle in degrees C CHM C*1 The sign of the sexagesimal number (or a digit C for hundreds or days C HM I(2) The angle in RA/dec as C HM(1) Hours or Degrees C HM(2) Minutes of time or arc C SEC R The rest of the Ra/dec : seconds of time or arc C----------------------------------------------------------------------- CHARACTER CHM*1 INTEGER HM(2), N REAL SEC DOUBLE PRECISION RAD C CHARACTER LOCHAR(10)*1 INTEGER I DOUBLE PRECISION TEMP, CON, F15, F60 C INCLUDE 'INCS:DMSG.INC' DATA LOCHAR /'1','2','3','4','5','6','7','8','9','*'/ C----------------------------------------------------------------------- C Error in N IF ((N.LT.1) .OR. (N.GT.4)) GO TO 900 F15 = 15.0D0 F60 = 60.0D0 CON = 3600.D0 C From degrees IF (N.GT.2) GO TO 50 TEMP = RAD IF ((N.EQ.1) .AND. (TEMP.LT.0.0D0)) TEMP = TEMP + 360.0D0 IF ((N.EQ.2) .AND. (TEMP.LT.-90.0D0)) TEMP = TEMP + 360.0D0 CHM = ' ' IF (N.EQ.1) TEMP = TEMP / F15 C Check on sign IF (TEMP.GE.0.0D0) GO TO 10 TEMP = -TEMP CHM = '-' GO TO 20 C longitude leading digits 10 IF (N.EQ.1) GO TO 15 I = TEMP / 100.0D0 IF (I.LE.0) GO TO 20 IF (I.LE.9) CHM = LOCHAR(I) IF (I.GE.10) CHM = LOCHAR(10) TEMP = TEMP - 100.0D0 * I GO TO 20 C RA: days part 15 CONTINUE I = TEMP / 24.0D0 IF (I.LE.0) GO TO 20 IF (I.LE.9) CHM = LOCHAR(I) IF (I.GE.10) CHM = LOCHAR(10) TEMP = TEMP - 24.0D0 * I C Separate into h m s 20 TEMP = TEMP * CON HM(1) = TEMP / CON TEMP = TEMP - HM(1) * CON HM(2) = TEMP / F60 SEC = TEMP - HM(2) * F60 GO TO 999 C Convert to degrees 50 CONTINUE TEMP = F60 * (F60 * HM(1) + HM(2)) + SEC IF (N.EQ.3) TEMP = TEMP * F15 IF (CHM.EQ.'-') TEMP = -TEMP RAD = TEMP / CON IF ((CHM.EQ.' ') .OR. (CHM.EQ.'-')) GO TO 999 DO 55 I = 1,9 IF (CHM.EQ.LOCHAR(I)) GO TO 60 55 CONTINUE I = 0 60 IF (N.EQ.2) RAD = RAD + I * 24.0D0 IF (N.EQ.4) RAD = RAD + I * 100.0D0 GO TO 999 C Error return 900 CONTINUE C WRITE (MSGTXT,1900) N C CALL MSGWRT (6) C 999 RETURN C----------------------------------------------------------------------- C 1900 FORMAT ('COORDD: IMPROPER CONTROL PARAMETER. N= ',I4) END SUBROUTINE CHLTOU (N, STRING) C----------------------------------------------------------------------- C! converts a CHARACTER string to all upper case letters C# Character Parsing C----------------------------------------------------------------------- C CHLTOU converts any lower case characters in a CHARACTER string to C upper case. C Rewritten substantially for NVSSlist C Inputs: N I Number of characters C In/out: STRING C*(*) String to be converted. C----------------------------------------------------------------------- CHARACTER STRING*(*) INTEGER N C INTEGER I, INDX CHARACTER IT*1, ALPHAU*26, ALPHAL*26 DATA ALPHAU /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/ DATA ALPHAL /'abcdefghijklmnopqrstuvwxyz'/ C----------------------------------------------------------------------- IF (N.LE.0) GO TO 999 DO 10 I = 1,N IT = STRING(I:I) C See if it's lower case INDX = INDEX (ALPHAL, IT) C If so replace IF (INDX.GT.0) STRING(I:I) = ALPHAU(INDX:INDX) 10 CONTINUE C 999 RETURN END SUBROUTINE TXTLEN (STRING, I1, I2) C----------------------------------------------------------------------- C! Find the first and last non-blank, non-tab characters in a string. C# Character Parsing C----------------------------------------------------------------------- C finds the positions of the first and last non-blank, non-tab C characters in a string. C Inputs: C STRING C*(*) Input character string. C Outputs: C I1 I Position of the first non-blank, non-tab; one C plus the length of the string if none C I2 I Position of the last non-blank, non-tab; 0 if C none C Author: Mark Calabretta, Australia Telescope. C Origin; 1988/May/23. Code last modified; 1989/Sep/06. C----------------------------------------------------------------------- INTEGER I1, I2 CHARACTER STRING*(*) C INTEGER NCH CHARACTER TAB*1 C----------------------------------------------------------------------- TAB = CHAR(9) NCH = LEN(STRING) C Search forward for the start. DO 10 I1 = 1,NCH IF (STRING(I1:I1).EQ.' ') GO TO 10 IF (STRING(I1:I1).EQ.TAB) GO TO 10 GO TO 20 10 CONTINUE C Search backward for the end. 20 DO 30 I2 = NCH,1,-1 IF (STRING(I2:I2).EQ.' ') GO TO 30 IF (STRING(I2:I2).EQ.TAB) GO TO 30 GO TO 999 30 CONTINUE C 999 RETURN END