C Position class: name = 'POSITION' C----------------------------------------------------------------------- C! Object Oriented FITSAIPS Fortran "position" class library C# Map-util Utility Object-Oriented C----------------------------------------------------------------------- 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 this software should be addressed as follows: C; Internet email: bcotton@nrao.edu. C; Postal address: Bill Cotton C; National Radio Astronomy Observatory C; 520 Edgemont Road C; Charlottesville, VA 22903-2475 USA C----------------------------------------------------------------------- C C C Public functions: C PSNCVT (init, in, xyzi, out, xyzo, ierr) C Convert pixel position in one image to pixel position in C another. C PSNPIX (init, in, x, y, xpix, ypix, ierr) C Determine the pixel location for a given x,y coordinate. C PSNVAL (init, in, xpix, ypix, x, y, z, ierr) C Determine thelocation for a given pixel. C C Private functions: C HEDGET (in, ierr) C Read file header geometry into common C SETLOC (depth, swapok) C Fill geometry information into position common C AXSTRN (axtyp, axval, axptr, nch, str) C Encodes axis type and value in a string C XYPIX (x, y, xpix, ypix, ierr) C Determines the pixel location corresponding to a specified C coordinate value. After setup by SETLOC. C XYVAL (Xpix, ypix, x, y, z, ierr) C Determines the coordinate value (X,Y,Z) corresponding to the C pixel location (XPIX,YPIX). After setup by SETLOC. C COORDT (iti, ito, longi, lati, epoki, epoko, longo, lato, ierr) C Translates between types of coordinates. C DIRCOS (type, ra0, dec0, ra, dec, l, m, ierr) C Determines the direction cosines (L,M) corresponding to C the separation from (RA0,DEC0), the reference pixel, to C (RA,DEC). C DIRDEC (type, dy, ra, ra0, dec0, rota, dx, dec, ierr) C Will find the latitude and longitude-like pixel position C of a point given the longitude and the latitude-pixel of the C point. C DIRRA (type, dx, dec, ra0, dec0, rota, dy, ra, ierr) C Will find the longitude and the latitude pixel position of a C point given the latitude and the longitude pixel of the C point. C DIRDEC (type, dy, ra, ra0, dec0, rota, dx, dec, ierr) C Will find the latitude and longitude-like pixel position C of a point given the longitude and the latitude-pixel of the C point. C NEWPOS (type, ra0, dec0, l, m, raout, decout, ierr) C Determines the coordinates (RAOUT,DECOUT) corresponding to a C displacement (L,M) given by the direction cosines from C coordinates (RA0,DEC0). C GALPOL (epoch, ragp, decgp, loncp) C Returns the coordinates of the galactic center for the epoch C of observation. C B2JPOS (dir, rain, decin, raout, decout) C Converts between B1950 positions and J2000 positions. C----------------------------------------------------------------------- LOCAL INCLUDE 'POSITION.INC' C Position info INTEGER ITI, ITO REAL EPOK, CDELT1(7), CDELT2(7), CRPIX1(7), CRPIX2(7) DOUBLE PRECISION RA0, DEC0 COMMON /PSNCOM/ RA0, DEC0, * EPOK, CDELT1, CDELT2, CRPIX1, CRPIX2, * ITI, ITO LOCAL END LOCAL INCLUDE 'DLOC.INC' C Position conversion common C Dimensionality limit INTEGER NUMLOC PARAMETER (NUMLOC = 3) DOUBLE PRECISION RPVAL(4,NUMLOC), COND2R, AXDENU(NUMLOC), * GEOMD1(NUMLOC), GEOMD2(NUMLOC), GEOMD3(NUMLOC), GEOMD4(NUMLOC) CHARACTER CTYP(4,NUMLOC)*20, CPREF(2,NUMLOC)*5, * SAXLAB(2,NUMLOC)*20 REAL RPLOC(4,NUMLOC), AXINC(4,NUMLOC), ROT(NUMLOC), * REPOCH(NUMLOC) INTEGER ZDEPTH(5,NUMLOC), ZAXIS(NUMLOC), AXTYP(NUMLOC), * CORTYP(NUMLOC), LABTYP(NUMLOC), SGNROT(NUMLOC), * AXFUNC(7,NUMLOC), KLOCL(NUMLOC), KLOCM(NUMLOC), KLOCF(NUMLOC), * KLOCS(NUMLOC), KLOCA(NUMLOC), KLOCB(NUMLOC), NCHLAB(2,NUMLOC), * LOCNUM COMMON /PSNLOC/ CTYP, CPREF, SAXLAB COMMON /PSNLOI/ RPVAL, COND2R, AXDENU, GEOMD1, GEOMD2, GEOMD3, * GEOMD4, RPLOC, AXINC, ROT, ZDEPTH, * ZAXIS, AXTYP, CORTYP, LABTYP, SGNROT, AXFUNC, KLOCL, KLOCM, * KLOCF, KLOCS, KLOCA, KLOCB, NCHLAB, REPOCH, LOCNUM C End DLOC. LOCAL END LOCAL INCLUDE 'DCAT.INC' C File header common INTEGER MAXDIM PARAMETER (MAXDIM = 7) INTEGER NAXIS, INAXES(MAXDIM) REAL CDELT(MAXDIM), CRPIX(MAXDIM), CROTA(MAXDIM), EPOCH, * ALTRP CHARACTER OBSDAT*8, BUNIT*8, CTYPE(MAXDIM)*8 DOUBLE PRECISION OBSRA, OBSDEC, RESTF, CRVAL(MAXDIM) COMMON /PSNCAI/ OBSRA, OBSDEC, RESTF, CRVAL, * CDELT, CRPIX, CROTA, EPOCH, ALTRP, * NAXIS, INAXES COMMON /PSNCAC/ OBSDAT, BUNIT, CTYPE C End DCAT. LOCAL END SUBROUTINE PSNCVT (INIT, IN, XYZI, OUT, XYZO, IERR) C----------------------------------------------------------------------- C Public C Convert pixel position in one image to pixel position in another. C Can convert between B1950 and J2000. C Inputs: C INIT L If true initialize position structures, Must be C .TRUE. on first call. C IN I Input FITS i/o unit C XYZI R(7) Pixel position in IN C OUT I Output FITS i/o unit C Outputs: C XYZO R(7) Pixel position in OUT. C IERR I Error code, 0=OK. C----------------------------------------------------------------------- LOGICAL INIT INTEGER IN, OUT REAL XYZI(7), XYZO(7) INTEGER IERR C INTEGER I, IDEP(5), IROUND, IEPOKI, IEPOKO LOGICAL SWAPOK REAL EPOKI, EPOKO DOUBLE PRECISION X, Y, Z, RA, DEC, RAT, DECT LOGICAL DOPRE INCLUDE 'DCAT.INC' INCLUDE 'DLOC.INC' INCLUDE 'POSITION.INC' DOUBLE PRECISION DEGRAD PARAMETER (DEGRAD=1.745329251994329576D-2) DATA SWAPOK /.FALSE./ SAVE EPOKI, EPOKO, DOPRE, IEPOKI, IEPOKO C----------------------------------------------------------------------- IERR = 0 C Get info for input image if C necessary. IF (INIT) THEN C Get catalog header CALL HEDGET(IN, IERR) IF (IERR.NE.0) GO TO 990 DO 10 I = 1,5 IDEP(I) = IROUND (XYZI(I+2)) 10 CONTINUE LOCNUM = 1 CALL SETLOC (IDEP, SWAPOK) ITI = 1 IF (LABTYP(LOCNUM).LT.44) ITI = 2 IF (LABTYP(LOCNUM).LT.22) ITI = 3 EPOKI = EPOCH IEPOKI = EPOKI + 0.5 C Get info for input image if C necessary. C Get catalog header CALL HEDGET (OUT, IERR) IF (IERR.NE.0) GO TO 990 DO 20 I = 1,5 IDEP(I) = IROUND (XYZO(I+2)) 20 CONTINUE LOCNUM = 2 CALL SETLOC (IDEP, SWAPOK) ITO = 1 IF (LABTYP(LOCNUM).LT.44) ITO = 2 IF (LABTYP(LOCNUM).LT.22) ITO = 3 IF (IERR.NE.0) GO TO 990 EPOKO = EPOCH IEPOKO = EPOKO + 0.5 C Do we need precession DOPRE = (ABS (EPOKI-EPOKO).GT.0.01) .AND. (EPOKI.GT.0.1) * .AND. (EPOKO.GT.0.1) END IF C Get position of input LOCNUM = 1 CALL XYVAL (XYZI(1), XYZI(2), X, Y, Z, IERR) IF (IERR.NE.0) GO TO 990 C Precess if necessary IF (DOPRE) THEN RA = X * DEGRAD DEC = Y * DEGRAD C B1950 to/from J2000 special IF ((IEPOKI.EQ.1950) .AND. (IEPOKO.EQ.2000)) THEN CALL B2JPOS (1, RA, DEC, RAT, DECT) C J2000 to B1950 ELSE IF ((IEPOKI.EQ.2000) .AND. (IEPOKO.EQ.1950)) THEN CALL B2JPOS (-1, RA, DEC, RAT, DECT) ELSE C General case -can't cope IERR = 2 WRITE (0,2001) IEPOKI, IEPOKO GO TO 999 END IF X = RAT / DEGRAD Y = DECT / DEGRAD END IF C Convert types if necessary IF (ITI.NE.ITO) THEN CALL COORDT (ITI, ITO, X, Y, EPOKO, X, Y, IERR) IF (IERR.NE.0) GO TO 990 END IF C Get corresponding pixel in OUT LOCNUM = 2 CALL XYPIX (X, Y, XYZO(1), XYZO(2), IERR) IF (IERR.NE.0) GO TO 990 GO TO 999 C Error 990 WRITE (0,2000) 'PSNCVT: ERROR CONVERTING PIXEL POSITIONS' C 999 RETURN 2000 FORMAT (A) 2001 FORMAT ('PSNCVT: CANNOT precess ',I4, ' to ', I4) END SUBROUTINE PSNPIX (INIT, IN, X, Y, XPIX, YPIX, IERR) C----------------------------------------------------------------------- C Public C Determine pixel position of a given x,y coordinate (e.g. RA, Dec) C Inputs: C INIT L If true initialize position structures C IN I Input FITS i/o unit C X D First axis coordinate in header units (degrees) C Y D Sedond axis coordinate in header units (degrees) C Outputs: C XPIX R First axis pixel number corresponding to (X,Y) C YPIX R Second axis pixel number corresponding to (X,Y) C IERR I Error code, 0=OK. C----------------------------------------------------------------------- LOGICAL INIT INTEGER IN, IERR DOUBLE PRECISION X, Y REAL XPIX, YPIX C INTEGER IDEP(5) LOGICAL SWAPOK INCLUDE 'DCAT.INC' INCLUDE 'DLOC.INC' INCLUDE 'POSITION.INC' DATA SWAPOK /.FALSE./ C----------------------------------------------------------------------- IERR = 0 C Get info for input image if C necessary. IF (INIT) THEN C Get catalog header CALL HEDGET (IN, IERR) IF (IERR.NE.0) GO TO 990 CALL FILL (5, 1, IDEP) LOCNUM = 1 CALL SETLOC (IDEP, SWAPOK) ITI = 1 IF (LABTYP(LOCNUM).LT.44) ITI = 2 IF (LABTYP(LOCNUM).LT.22) ITI = 3 END IF C Get pixel value LOCNUM = 1 CALL XYPIX (X, Y, XPIX, YPIX, IERR) IF (IERR.NE.0) GO TO 990 GO TO 999 C Error 990 WRITE(0,2000) 'PSNPIX: ERROR FINDING PIXEL' C 999 RETURN 2000 FORMAT (A) END SUBROUTINE PSNVAL (INIT, IN, XPIX, YPIX, X, Y, Z, IERR) C----------------------------------------------------------------------- C Public C Determine the coordinate value of a given pixel in an image. C Inputs: C INIT L If true initialize position structures C IN I Input FITS i/o unit number C XPIX R Pixel location, x-coordinate C YPIX R Pixel location, y-coordinate C Outputs: C X D X-coordinate value at pixel location C Y D Y-coordinate value at pixel location C Z D Z-coordinate value (if part of a position C pair with either X or Y) C IERR I Error code, 0=OK. C----------------------------------------------------------------------- LOGICAL INIT INTEGER IN, IERR DOUBLE PRECISION X, Y, Z REAL XPIX, YPIX C INTEGER IDEP(5) LOGICAL SWAPOK INCLUDE 'DCAT.INC' INCLUDE 'DLOC.INC' INCLUDE 'POSITION.INC' DATA SWAPOK /.FALSE./ C----------------------------------------------------------------------- IERR = 0 C Get info for input image if C necessary. IF (INIT) THEN C Get catalog header CALL HEDGET (IN, IERR) IF (IERR.NE.0) GO TO 990 CALL FILL (5, 1, IDEP) LOCNUM = 1 CALL SETLOC (IDEP, SWAPOK) ITI = 1 IF (LABTYP(LOCNUM).LT.44) ITI = 2 IF (LABTYP(LOCNUM).LT.22) ITI = 3 END IF C Get coordinate value LOCNUM = 1 CALL XYVAL (XPIX, YPIX, X, Y, Z, IERR) IF (IERR.NE.0) GO TO 990 GO TO 999 C Error 990 WRITE (0,2000) 'PSNVAL: ERROR FINDING COORDINATES ' C 999 RETURN 2000 FORMAT (A) END SUBROUTINE HEDGET (IN, IERR) C----------------------------------------------------------------------- C Read FITS header info from IN and save in common C Inputs: C IN I Input FITS fitsio unit number C Output: C IERR I Error code: 0 => ok C Output in common: C NAXIS I Number of axes in image C INAXES I(*) Array dimensionality C OBSDAT C DATE-OBS C BUNIT C BUNIT C CTYPE C(*) CTYPEn C EPOCH R EPOCH C ALTRP R ALTRPIX C CDELT R(*) CDELTn C CRPIX R(*) CRPIXn C OBSRA D OBSRA C RESTF D RESTF (Rest frequency) C OBSDEC D OBSDEC C CRVAL D(*) CRVALn C----------------------------------------------------------------------- INTEGER IN, IERR C INCLUDE 'DCAT.INC' INTEGER NFOUND, I INTEGER BITPIX, PCOUNT, GCOUNT LOGICAL SIMPLE, EXTEND CHARACTER CDATA(7)*8, COMMNT*80 C----------------------------------------------------------------------- C Don't know it's name RESTF = 0.0D0 C Read keyword values C Input size CALL FTGHPR (IN, MAXDIM, SIMPLE, BITPIX, NAXIS, INAXES, * PCOUNT, GCOUNT, EXTEND, IERR) IF (IERR.NE.0) THEN WRITE(0,1000) 'FTGHPR error' GO TO 999 END IF OBSDAT = ' ' CALL FTGKYS (IN, 'DATE-OBS', CDATA, COMMNT, IERR) IF (IERR.EQ.0) OBSDAT = CDATA(1) IF (IERR.EQ.202) IERR = 0 BUNIT = ' ' CALL FTGKYS (IN, 'BUNIT', CDATA, COMMNT, IERR) IF (IERR.EQ.0) BUNIT = CDATA(1) IF (IERR.EQ.202) IERR = 0 DO 10 I = 1,MAXDIM CTYPE(I) = ' ' 10 CONTINUE CALL FTGKNS (IN, 'CTYPE', 1, MAXDIM, CDATA, NFOUND, IERR) IF (IERR.EQ.0) THEN DO 20 I = 1,NFOUND CTYPE(I) = CDATA(I) 20 CONTINUE END IF IF (IERR.EQ.202) IERR = 0 EPOCH = 0.0 CALL FTGKYE (IN, 'EPOCH', EPOCH, COMMNT, IERR) IF (IERR.EQ.202) IERR = 0 ALTRP = 1.0 CALL FTGKYE (IN, 'ALTRPIX', ALTRP, COMMNT, IERR) IF (IERR.EQ.202) IERR = 0 CALL RFILL (7, 1.0, CDELT) CALL FTGKNE (IN, 'CDELT', 1, MAXDIM, CDELT, NFOUND, IERR) IF (IERR.EQ.202) IERR = 0 CALL RFILL (7, 1.0, CRPIX) CALL FTGKNE (IN, 'CRPIX', 1, MAXDIM, CRPIX, NFOUND, IERR) IF (IERR.EQ.202) IERR = 0 CALL RFILL (7, 1.0, CROTA) CALL FTGKNE (IN, 'CROTA', 1, MAXDIM, CROTA, NFOUND, IERR) IF (IERR.EQ.202) IERR = 0 OBSRA = 0.0D0 CALL FTGKYD (IN, 'OBSRA', OBSRA, COMMNT, IERR) IF (IERR.EQ.202) IERR = 0 OBSDEC = 0.0D0 CALL FTGKYD (IN, 'OBSDEC', OBSDEC, COMMNT, IERR) IF (IERR.EQ.202) IERR = 0 DO 60 I = 1,MAXDIM CRVAL(I) = 0.0D0 60 CONTINUE CALL FTGKND (IN, 'CRVAL', 1, MAXDIM, CRVAL, NFOUND, IERR) IF (IERR.EQ.202) IERR = 0 IF (IERR.NE.0) WRITE(0,1000) 'ERROR reading input FITS header' C 999 RETURN C----------------------------------------------------------------------- 1000 FORMAT(A) END SUBROUTINE SETLOC (DEPTH, SWAPOK) C----------------------------------------------------------------------- C SETLOC uses the catalog header to build the values of the position C common /LOCATI/ for use by position finding and axis labeling C routines (at least). C Inputs: C DEPTH I(5) Position of map plane axes 3 - 7 C SWAPOK L T => okay to swap axes if rotation near 90 C Common: C DCAT.INC catalog block (not modified) C DLOC.INC position parms - created here C----------------------------------------------------------------------- INTEGER DEPTH(5) LOGICAL SWAPOK C CHARACTER AXP(2,5)*4, STOKES*4, FRAXP(3)*4, GEOM(8)*4, TS*20, * CTEST1*4, CTEST2*4 INTEGER NANG, KL, JX, I, J, NAX, NGEOM, II, JF, JY REAL SWAP DOUBLE PRECISION DT, DX, DY, COSR, SINR, X, Y, VELITE, DOT INCLUDE 'DLOC.INC' INCLUDE 'DCAT.INC' DATA NANG /5/ DATA AXP /'ELON','ELAT', 'GLON','GLAT', 'RA--','DEC-', * 'RA ','DEC ', 'LL ','MM '/ DATA STOKES, FRAXP /'STOK', 'FREQ','VELO','FELO'/ DATA NGEOM, GEOM /8, '-SIN','-TAN','-ARC','-NCP','-GLS','-MER', * '-AIT','-STG'/ DATA VELITE /2.997924562D8 / C----------------------------------------------------------------------- C Basic parms: assume VLA normal C CALL CHECKL ('SETLOC') COND2R = 3.14159265358979D0 / 1.80D2 C REPOCH(LOCNUM) = CATR(KREPO) REPOCH(LOCNUM) = EPOCH C RPVAL(1,LOCNUM) = CATD(KDCRV) RPVAL(1,LOCNUM) = CRVAL(1) C RPVAL(2,LOCNUM) = CATD(KDCRV+1) RPVAL(2,LOCNUM) = CRVAL(2) C RPLOC(1,LOCNUM) = CATR(KRCRP) RPLOC(1,LOCNUM) = CRPIX(1) C RPLOC(2,LOCNUM) = CATR(KRCRP+1) RPLOC(2,LOCNUM) = CRPIX(2) C AXINC(1,LOCNUM) = CATR(KRCIC) AXINC(1,LOCNUM) = CDELT(1) IF (AXINC(1,LOCNUM).EQ.0.0) AXINC(1,LOCNUM) = 1.0 C AXINC(2,LOCNUM) = CATR(KRCIC+1) AXINC(2,LOCNUM) = CDELT(2) IF (AXINC(2,LOCNUM).EQ.0.0) AXINC(2,LOCNUM) = 1.0 C CALL H2CHR (8, 1, CATH(KHCTP), CTYP(1,LOCNUM)) CTYP(1,LOCNUM) = CTYPE(1) C CALL H2CHR (8, 1, CATH(KHCTP+2), CTYP(2,LOCNUM)) CTYP(2,LOCNUM) = CTYPE(2) C ROT(LOCNUM) = CATR(KRCRT+1) ROT(LOCNUM) = CROTA(2) CPREF(1,LOCNUM) = ' ' CPREF(2,LOCNUM) = ' ' CTYP(3,LOCNUM) = ' ' CTYP(4,LOCNUM) = ' ' CALL COPY (5, DEPTH, ZDEPTH(1,LOCNUM)) ZAXIS(LOCNUM) = 0 AXTYP(LOCNUM) = 1 CORTYP(LOCNUM) = 1 LABTYP(LOCNUM) = 65 SGNROT(LOCNUM) = 1 GEOMD1(LOCNUM) = 1.0D0 GEOMD2(LOCNUM) = 1.0D0 GEOMD3(LOCNUM) = 0.0D0 GEOMD4(LOCNUM) = 0.0D0 C Get pointers to axis types KLOCL(LOCNUM) = -1 KLOCM(LOCNUM) = -1 KLOCF(LOCNUM) = -1 KLOCS(LOCNUM) = -1 KLOCA(LOCNUM) = -1 KLOCB(LOCNUM) = -1 C NAX = CATBLK(KIDIM) NAX = NAXIS JX = 0 JY = 0 JF = 0 DO 20 I = 1,NAX C II = (I-1) * 2 + KHCTP C CALL H2CHR (4, 1, CATR(II), CTEST1) CTEST1 = CTYPE(I)(1:4) DO 10 J = 1,NANG IF (AXP(1,J).EQ.CTEST1) THEN JX = J KLOCL(LOCNUM) = I - 1 GO TO 20 ELSE IF (AXP(2,J).EQ.CTEST1) THEN JY = J KLOCM(LOCNUM) = I - 1 GO TO 20 END IF 10 CONTINUE DO 15 J = 1,3 IF (FRAXP(J).EQ.CTEST1) THEN JF = J KLOCF(LOCNUM) = I - 1 GO TO 20 END IF 15 CONTINUE IF (STOKES.EQ.CTEST1) KLOCS(LOCNUM) = I - 1 20 CONTINUE C Non-linear axes? CALL FILL (7, -1, AXFUNC(1,LOCNUM)) CALL FILL (NAX, 0, AXFUNC(1,LOCNUM)) C velocity IF (JF.NE.3) GO TO 25 AXFUNC(KLOCF(LOCNUM)+1,LOCNUM) = 1 AXDENU(LOCNUM) = 0.0D0 C IF (CATD(KDRST).GT.0.0D0) AXDENU(LOCNUM) = C * -CATR(KRCIC+KLOCF(LOCNUM)) / (VELITE + C * CATD(KDCRV+KLOCF(LOCNUM))) IF (RESTF.GT.0.0D0) AXDENU(LOCNUM) = * -CDELT(1+KLOCF(LOCNUM)) / (VELITE + * CRVAL(1+KLOCF(LOCNUM))) C Projective geometries 25 IF ((KLOCL(LOCNUM).GE.0) .AND. (KLOCM(LOCNUM).GE.0)) THEN IF ((JX.EQ.5) .AND. (JY.EQ.5)) THEN AXFUNC(KLOCL(LOCNUM)+1,LOCNUM) = 2 AXFUNC(KLOCM(LOCNUM)+1,LOCNUM) = 2 END IF JX = MIN (JX, 3) JY = MIN (JY, 3) C CALL H2CHR (4, 5, CATH(KHCTP+KLOCL(LOCNUM)*2), CTEST1) CTEST1 = CTYPE(1+KLOCL(LOCNUM))(5:8) IF (CTEST1.NE.' ') AXFUNC(KLOCL(LOCNUM)+1,LOCNUM) = -100 C CALL H2CHR (4, 5, CATH(KHCTP+KLOCM(LOCNUM)*2), CTEST2) CTEST2 = CTYPE(1+KLOCM(LOCNUM))(5:8) IF (CTEST2.NE.' ') AXFUNC(KLOCM(LOCNUM)+1,LOCNUM) = -101 DO 35 I = 1,NGEOM IF (GEOM(I).EQ.CTEST1) AXFUNC(KLOCL(LOCNUM)+1,LOCNUM) = * I + 1 IF (GEOM(I).EQ.CTEST2) AXFUNC(KLOCM(LOCNUM)+1,LOCNUM) = * I + 1 35 CONTINUE END IF C Rational angle pair? IF ((JX.EQ.JY) .AND. (AXFUNC(KLOCL(LOCNUM)+1,LOCNUM).EQ. * AXFUNC(KLOCM(LOCNUM)+1,LOCNUM))) GO TO 45 IF (((KLOCL(LOCNUM).GE.0) .AND. * (AXFUNC(KLOCL(LOCNUM)+1,LOCNUM).LT.-10)) .OR. * ((KLOCM(LOCNUM).GE.0) .AND. * (AXFUNC(KLOCM(LOCNUM)+1,LOCNUM).LT.-10))) * THEN WRITE (0,1041) ELSE WRITE (0,1040) END IF IF (KLOCL(LOCNUM).GE.0) AXFUNC(KLOCL(LOCNUM)+1,LOCNUM) = 0 IF (KLOCM(LOCNUM).GE.0) AXFUNC(KLOCM(LOCNUM)+1,LOCNUM) = 0 KLOCL(LOCNUM) = -1 KLOCM(LOCNUM) = -1 JX = 0 JY = 0 GO TO 55 C Advance parms for geometries C Mercator 45 IF (AXFUNC(KLOCL(LOCNUM)+1,LOCNUM).NE.7) GO TO 50 C ROT(LOCNUM) = CATR(KRCRT+KLOCM(LOCNUM)) ROT(LOCNUM) = CROTA(1+KLOCM(LOCNUM)) C DT = CATR(KRCIC+KLOCM(LOCNUM))*COS(ROT(LOCNUM)*COND2R) + C * CATR(KRCIC+KLOCL(LOCNUM))* SIN(ROT(LOCNUM)*COND2R) DT = CDELT(1+KLOCM(LOCNUM))*COS(ROT(LOCNUM)*COND2R) + * CDELT(1+KLOCL(LOCNUM))* SIN(ROT(LOCNUM)*COND2R) IF (DT.EQ.0.0D0) DT = 1.0D0 C DY = (CATD(KDCRV+KLOCM(LOCNUM))/2.0D0 + 45.0D0) * COND2R DY = (CRVAL(1+KLOCM(LOCNUM))/2.0D0 + 45.0D0) * COND2R DX = DY + DT / 2.0D0 * COND2R DY = LOG (TAN (DY)) DX = LOG (TAN (DX)) GEOMD2(LOCNUM) = DT * COND2R / (DX - DY) GEOMD3(LOCNUM) = GEOMD2(LOCNUM) * DY C GEOMD1(LOCNUM) = COS (CATD(KDCRV+KLOCM(LOCNUM))*COND2R) GEOMD1(LOCNUM) = COS (CRVAL(1+KLOCM(LOCNUM))*COND2R) IF (GEOMD1(LOCNUM).LE.0.0D0) GEOMD1(LOCNUM) = 1.0D0 C Aitoff 50 IF (AXFUNC(KLOCL(LOCNUM)+1,LOCNUM).NE.8) GO TO 55 C DOT = CATR(KRCRT+KLOCM(LOCNUM)) * COND2R DOT = CROTA(1+KLOCM(LOCNUM)) * COND2R C DT = CATR(KRCIC+KLOCM(LOCNUM))*COS(DOT) + C * CATR(KRCIC+KLOCL(LOCNUM))*SIN(DOT) DT = CDELT(1+KLOCM(LOCNUM))*COS(DOT) + * CDELT(1+KLOCL(LOCNUM))*SIN(DOT) IF (DT.EQ.0.0D0) DT = 1.0D0 DT = DT * COND2R C DY = CATD(KDCRV+KLOCM(LOCNUM)) * COND2R DY = CRVAL(1+KLOCM(LOCNUM)) * COND2R DX = SIN(DY+DT)/SQRT((1.0D0+COS(DY+DT))/2.0D0) - * SIN(DY)/SQRT((1.0D0+COS(DY))/2.0D0) IF (DX.EQ.0.0D0) DX = 1.0D0 GEOMD2(LOCNUM) = DT / DX C DT = CATR(KRCIC+KLOCL(LOCNUM))*COS(DOT) - C * CATR(KRCIC+KLOCM(LOCNUM))* SIN(DOT) DT = CDELT(1+KLOCL(LOCNUM))*COS(DOT) - * CDELT(1+KLOCM(LOCNUM))* SIN(DOT) IF (DT.EQ.0.0D0) DT = 1.0D0 DT = DT * COND2R DX = 2.0D0 * COS(DY) * SIN(DT/2.0D0) IF (DX.EQ.0.0D0) DX = 1.0D0 GEOMD1(LOCNUM) = DT * SQRT((1.0D0+COS(DY)*COS(DT/2.0D0))/2.0D0) * / DX GEOMD3(LOCNUM) = GEOMD2(LOCNUM) * SIN(DY) / * SQRT((1.0D0+COS(DY))/2.0D0) C Other parms C No angle pair 55 IF ((KLOCL(LOCNUM).GE.0) .AND. (KLOCM(LOCNUM).GE.0)) GO TO 60 ROT(LOCNUM) = 0.0 AXTYP(LOCNUM) = 0 CORTYP(LOCNUM) = 0 LABTYP(LOCNUM) = 0 IF (KLOCL(LOCNUM).GE.0) LABTYP(LOCNUM) = (2 * MIN (JX, 3) - 1) * * (10 ** KLOCL(LOCNUM)) IF (KLOCM(LOCNUM).GE.0) LABTYP(LOCNUM) = (2 * MIN (JX, 3)) * * (10 ** KLOCM(LOCNUM)) GO TO 100 C normal order C 60 ROT(LOCNUM) = CATR(KRCRT+KLOCM(LOCNUM)) 60 ROT(LOCNUM) = CROTA(1+KLOCM(LOCNUM)) IF ((KLOCL(LOCNUM).NE.0) .OR. (KLOCM(LOCNUM).NE.1)) GO TO 65 AXTYP(LOCNUM) = 1 CORTYP(LOCNUM) = 1 LABTYP(LOCNUM) = 22 * JX - 1 GO TO 100 C reverse order 65 IF ((KLOCL(LOCNUM).NE.1) .OR. (KLOCM(LOCNUM).NE.0)) GO TO 70 AXTYP(LOCNUM) = 1 CORTYP(LOCNUM) = 2 SGNROT(LOCNUM) = -1 LABTYP(LOCNUM) = 22 * JX - 10 GO TO 100 C longitude with z 70 IF (KLOCL(LOCNUM).GT.1) GO TO 75 AXTYP(LOCNUM) = 2 + KLOCL(LOCNUM) CORTYP(LOCNUM) = 3 + 2 * KLOCL(LOCNUM) LABTYP(LOCNUM) = (2 * JX - 1) * (10 ** KLOCL(LOCNUM)) ZAXIS(LOCNUM) = KLOCM(LOCNUM) + 1 KLOCA(LOCNUM) = KLOCM(LOCNUM) GO TO 100 C latitude with z 75 IF (KLOCM(LOCNUM).GT.1) GO TO 80 AXTYP(LOCNUM) = 2 + KLOCM(LOCNUM) CORTYP(LOCNUM) = 4 + 2 * KLOCM(LOCNUM) LABTYP(LOCNUM) = (2 * JX) * (10 ** KLOCM(LOCNUM)) ZAXIS(LOCNUM) = KLOCL(LOCNUM) + 1 KLOCA(LOCNUM) = KLOCL(LOCNUM) SGNROT(LOCNUM) = -1 GO TO 100 C 2 "z" axes paired 80 CONTINUE AXTYP(LOCNUM) = 4 KLOCA(LOCNUM) = KLOCL(LOCNUM) KLOCB(LOCNUM) = KLOCM(LOCNUM) CORTYP(LOCNUM) = 1 IF (KLOCL(LOCNUM).GT.KLOCM(LOCNUM)) CORTYP(LOCNUM) = 2 IF (CORTYP(LOCNUM).EQ.2) SGNROT(LOCNUM) = -1 LABTYP(LOCNUM) = 0 C Other axes to display? 100 IF (NAX.LE.2) GO TO 110 IF ((KLOCA(LOCNUM).GE.2) .AND. (KLOCB(LOCNUM).GE.2)) GO TO 110 DO 105 I = 3,NAX C IF (CATBLK(KINAX+I-1).LE.1) GO TO 105 IF (INAXES(I).LE.1) GO TO 105 IF ((KLOCA(LOCNUM).GE.2) .AND. (KLOCA(LOCNUM).NE.I-1)) * KLOCB(LOCNUM) = I - 1 IF (KLOCB(LOCNUM).GE.2) GO TO 110 KLOCA(LOCNUM) = I - 1 105 CONTINUE IF (KLOCA(LOCNUM).LT.2) KLOCA(LOCNUM) = KLOCL(LOCNUM) IF (KLOCA(LOCNUM).LT.2) KLOCA(LOCNUM) = KLOCM(LOCNUM) IF (KLOCA(LOCNUM).LT.2) KLOCA(LOCNUM) = KLOCS(LOCNUM) IF (KLOCA(LOCNUM).LT.2) KLOCA(LOCNUM) = KLOCF(LOCNUM) IF ((KLOCB(LOCNUM).LT.2) .OR. (KLOCB(LOCNUM).EQ.KLOCA(LOCNUM))) * KLOCB(LOCNUM) = KLOCL(LOCNUM) IF ((KLOCB(LOCNUM).LT.2) .OR. (KLOCB(LOCNUM).EQ.KLOCA(LOCNUM))) * KLOCB(LOCNUM) = KLOCM(LOCNUM) IF ((KLOCB(LOCNUM).LT.2) .OR. (KLOCB(LOCNUM).EQ.KLOCA(LOCNUM))) * KLOCB(LOCNUM) = KLOCS(LOCNUM) IF ((KLOCB(LOCNUM).LT.2) .OR. (KLOCB(LOCNUM).EQ.KLOCA(LOCNUM))) * KLOCB(LOCNUM) = KLOCF(LOCNUM) IF (KLOCB(LOCNUM).EQ.KLOCA(LOCNUM)) KLOCB(LOCNUM) = -1 IF (KLOCA(LOCNUM).LT.2) KLOCA(LOCNUM) = -1 C get header values 110 IF (KLOCA(LOCNUM).LT.2) GO TO 115 C RPVAL(3,LOCNUM) = CATD(KDCRV+KLOCA(LOCNUM)) RPVAL(3,LOCNUM) = CRVAL(1+KLOCA(LOCNUM)) C RPLOC(3,LOCNUM) = CATR(KRCRP+KLOCA(LOCNUM)) RPLOC(3,LOCNUM) = CRPIX(1+KLOCA(LOCNUM)) C AXINC(3,LOCNUM) = CATR(KRCIC+KLOCA(LOCNUM)) AXINC(3,LOCNUM) = CDELT(1+KLOCA(LOCNUM)) IF (AXINC(3,LOCNUM).EQ.0.0) AXINC(3,LOCNUM) = 1.0 C CALL H2CHR (8, 1, CATH(KHCTP+KLOCA(LOCNUM)), CTYP(3,LOCNUM)) CTYP(3,LOCNUM) = CTYPE (1+KLOCA(LOCNUM)) 115 IF (KLOCB(LOCNUM).LT.2) GO TO 120 C RPVAL(4,LOCNUM) = CATD(KDCRV+KLOCB(LOCNUM)) RPVAL(4,LOCNUM) = CRVAL(1+KLOCB(LOCNUM)) C RPLOC(4,LOCNUM) = CATR(KRCRP+KLOCB(LOCNUM)) RPLOC(4,LOCNUM) = CRPIX(1+KLOCB(LOCNUM)) C AXINC(4,LOCNUM) = CATR(KRCIC+KLOCB(LOCNUM)) AXINC(4,LOCNUM) = CDELT(1+KLOCB(LOCNUM)) IF (AXINC(4,LOCNUM).EQ.0.0) AXINC(4,LOCNUM) = 1.0 C CALL H2CHR (8, 1, CATH(KHCTP+KLOCB(LOCNUM)*2), CTYP(4,LOCNUM)) CTYP(4,LOCNUM) = CTYPE (1+KLOCB(LOCNUM)) C Check ROT 120 IF (ROT(LOCNUM).GT.180.0) ROT(LOCNUM) = ROT(LOCNUM) - 360.0 IF (ROT(LOCNUM).LT.-180.0) ROT(LOCNUM) = ROT(LOCNUM) + 360.0 IF (ABS(ROT(LOCNUM)).LE.90.0) GO TO 130 I = 1 J = 2 IF (AXTYP(LOCNUM).GT.1) J = 3 IF (AXTYP(LOCNUM).EQ.3) I = 2 IF (AXTYP(LOCNUM).EQ.4) J = 4 IF (AXTYP(LOCNUM).EQ.4) I = 3 AXINC(I,LOCNUM) = -AXINC(I,LOCNUM) AXINC(J,LOCNUM) = -AXINC(J,LOCNUM) ROT(LOCNUM) = ROT(LOCNUM) - SIGN(180.0, ROT(LOCNUM)) C Swap axes on ROT 90 130 IF (.NOT.SWAPOK) GO TO 200 IF ((AXTYP(LOCNUM).LT.1) .OR. (AXTYP(LOCNUM).GT.3)) GO TO 200 IF ((ABS(ROT(LOCNUM)-90.0).GT.20.0) .AND. * (ABS(ROT(LOCNUM)+90.0).GT.20.0)) GO TO 200 SWAP = SIGN (1.0, ROT(LOCNUM)) I = 1 J = 2 IF (AXTYP(LOCNUM).GT.1) J = 3 IF (AXTYP(LOCNUM).EQ.3) I = 2 IF (AXTYP(LOCNUM).EQ.4) J = 4 IF (AXTYP(LOCNUM).EQ.4) I = 3 DT = RPVAL(I,LOCNUM) RPVAL(I,LOCNUM) = RPVAL(J,LOCNUM) RPVAL(J,LOCNUM) = DT TS = CTYP(I,LOCNUM) CTYP(I,LOCNUM) = CTYP(J,LOCNUM) CTYP(J,LOCNUM) = TS ROT(LOCNUM) = ROT(LOCNUM) - 90.0 * SWAP AXINC(I,LOCNUM) = SWAP * AXINC(I,LOCNUM) AXINC(J,LOCNUM) = - SWAP * AXINC(J,LOCNUM) SGNROT(LOCNUM) = -SGNROT(LOCNUM) IF (CORTYP(LOCNUM).GT.0) CORTYP(LOCNUM) = CORTYP(LOCNUM) - * ((-1)**CORTYP(LOCNUM)) I = KLOCL(LOCNUM) KLOCL(LOCNUM) = KLOCM(LOCNUM) KLOCM(LOCNUM) = I IF (AXTYP(LOCNUM).EQ.1) LABTYP(LOCNUM) = 10 * * MOD(LABTYP(LOCNUM),10) + LABTYP(LOCNUM)/10 IF ((AXTYP(LOCNUM).LE.1) .OR. (AXTYP(LOCNUM).EQ.4)) GO TO 200 IF (KLOCL(LOCNUM).LT.2) LABTYP(LOCNUM) = (2*JX-1) * (10 ** * KLOCL(LOCNUM)) IF (KLOCM(LOCNUM).LT.2) LABTYP(LOCNUM) = (2*JX) * (10 ** * KLOCM(LOCNUM)) C Set up fixed axis values 200 SAXLAB(1,LOCNUM) = ' ' SAXLAB(2,LOCNUM) = ' ' NCHLAB(1,LOCNUM) = 0 NCHLAB(2,LOCNUM) = 0 KL = KLOCA(LOCNUM) J = 1 C two extra labels loop 210 IF (KL.LT.2) GO TO 999 C Linear axis IF (AXFUNC(KL+1,LOCNUM).GT.0) GO TO 220 C DT = CATD(KDCRV+KL) + (ZDEPTH(KL-1,LOCNUM) - CATR(KRCRP+KL)) C * * CATR(KRCIC+KL) DT = CRVAL(1+KL) + (ZDEPTH(KL-1,LOCNUM) - CRPIX(1+KL)) * * CDELT(1+KL) CALL AXSTRN (CTYP(J+2,LOCNUM), DT, KL, NCHLAB(J,LOCNUM), * SAXLAB(J,LOCNUM)) GO TO 240 C Felocity 220 IF (AXFUNC(KL+1,LOCNUM).GT.1) GO TO 230 C DT = ZDEPTH(KL-1,LOCNUM) - CATR(KRCRP+KL) DT = ZDEPTH(KL-1,LOCNUM) - CRPIX(1+KL) C DT = CATD(KDCRV+KL) + CATR(KRCIC+KL)*DT / (1.0D0 + C * AXDENU(LOCNUM)*DT) DT = CRVAL(1+KL) + CDELT(1+KL)*DT / (1.0D0 + * AXDENU(LOCNUM)*DT) CALL AXSTRN (CTYP(J+2,LOCNUM), DT, KL, NCHLAB(J,LOCNUM), * SAXLAB(J,LOCNUM)) GO TO 240 C Positions: only 2 z axes now 230 IF (AXTYP(LOCNUM).NE.4) GO TO 240 DX = (ZDEPTH(KLOCL(LOCNUM)-1,LOCNUM) - RPLOC(3,LOCNUM)) * * AXINC(3,LOCNUM) DY = (ZDEPTH(KLOCM(LOCNUM)-1,LOCNUM) - RPLOC(4,LOCNUM)) * * AXINC(4,LOCNUM) COSR = COS (ROT(LOCNUM) * COND2R) SINR = SGNROT(LOCNUM) * SIN (ROT(LOCNUM) * COND2R) DT = (DX * COSR - DY * SINR) * COND2R DY = (DY * COSR + DX * SINR) * COND2R DX = DT COSR = RPVAL(3,LOCNUM) * COND2R SINR = RPVAL(4,LOCNUM) * COND2R CALL NEWPOS (AXFUNC(KLOCL(LOCNUM)+1,LOCNUM), COSR, SINR, DX, * DY, X, Y, II) IF (II.NE.0) GO TO 999 X = X / COND2R Y = Y / COND2R CALL AXSTRN (CTYP(3,LOCNUM), X, KLOCL(LOCNUM), * NCHLAB(1,LOCNUM), SAXLAB(1,LOCNUM)) CALL AXSTRN (CTYP(4,LOCNUM), Y, KLOCM(LOCNUM), * NCHLAB(2,LOCNUM), SAXLAB(2,LOCNUM)) GO TO 999 C loop for second label 240 IF (J.EQ.2) GO TO 999 J = 2 KL = KLOCB(LOCNUM) GO TO 210 C 999 RETURN C----------------------------------------------------------------------- 1040 FORMAT ('ANGLE-PAIR AXES INCONSISTENT OR PARTLY MISSING') 1041 FORMAT ('ANGLE-PAIR AXES PROJECTION TYPE NOT RECOGNIZED') END SUBROUTINE AXSTRN (AXTYP, AXVAL, AXPTR, NCH, STR) C----------------------------------------------------------------------- C AXSTRN returns a 20-character, left-justified string encoding the C axis type and value appropriately. C Inputs: C AXTYP C*(*) Axis type -- 1st 8 real chars used C AXVAL D Value on that axis C AXPTR I 0-relative pointer of which axis N.B. C Output: C NCH I Number of characters in STR (can be rather C less than 20 - 16 preferred for pos.) C STR C*20 Axis descriptor string C Common: DCAH.INC C----------------------------------------------------------------------- CHARACTER AXTYP*(*), STR*(*) DOUBLE PRECISION AXVAL, X INTEGER AXPTR, NCH C C INTEGER CATBLK(256) INTEGER NA, NTY, IFRM, NB, IT, HM(2), NC REAL Y, SEC C REAL CATR(256) CHARACTER XSTR*20, WST*20, SXTYP(15)*4, CHSTOK(16)*4, CHM*1 INCLUDE 'DCAT.INC' C COMMON /MAPHDR/ CATBLK C EQUIVALENCE (CATBLK, CATR) DATA NTY, SXTYP /15, 'LL ','RA ','RA--','GLON','ELON', * 'MM ','DEC ','DEC-','GLAT','ELAT','STOK','FREQ','VELO', * 'FELO','TIME'/ DATA CHSTOK /'BEAM','IPOL','QPOL','UPOL','VPOL','PPOL','FPOL', * 'PANG','SPIX','OPTD','ROTM','????','RR ','LL ','RL ','LR '/ C----------------------------------------------------------------------- STR = ' ' C is it a special type? DO 10 IT = 1,NTY IF (AXTYP(1:4).EQ.SXTYP(IT)) GO TO 100 10 CONTINUE C No: pack AXTYP XSTR = AXTYP CALL XTRIM (XSTR, 8, XSTR, NA) NA = MIN (NA + 2, 11) IF (NA.GT.1) THEN XSTR(NA:) = ' ' NA = NA + 1 END IF C Select a format 15 X = ABS (AXVAL) IFRM = 1 IF (X.LT.1.0D0) IFRM = 2 IF (X.LT.1.D-3) IFRM = 3 IF (X.LT.1.D-5) IFRM = 5 IF (X.GT.1.D3) IFRM = 4 IF (X.GT.1.D5) IFRM = 5 20 IF (IFRM.EQ.1) WRITE (WST,1020,ERR=880) AXVAL IF (IFRM.EQ.2) WRITE (WST,1021,ERR=880) AXVAL IF (IFRM.EQ.3) WRITE (WST,1022,ERR=880) AXVAL IF (IFRM.EQ.4) WRITE (WST,1023,ERR=880) AXVAL IF (IFRM.NE.5) GO TO 900 WRITE (WST,1024,ERR=880) AXVAL IF (NA.GT.9) NA = 9 XSTR(NA-1:) = ' ' GO TO 900 C Special types: branch 100 NA = IT - 10 IF (NA.GT.0) GO TO (130, 140, 150, 150, 160), NA C An angle axis: fix prefix IF ((IT.EQ.1) .OR. (IT.EQ.3)) IT = 2 IF ((IT.EQ.6) .OR. (IT.EQ.8)) IT = 7 NA = 5 IF (IT.EQ.2) NA = 3 IF (IT.EQ.7) NA = 4 XSTR = SXTYP(IT) XSTR(NA:NA) = ' ' NA = NA + 1 C Check field of view X = AXVAL C Y = ABS (CATR(KRCIC+AXPTR)) * 3600.0 Y = ABS (CDELT(1+AXPTR)) * 3600.0 IF (Y.LE.0.0) Y = 121. IF (MOD(IT-1,5).GT.2) Y = Y * 20. C Do in degrees IF (Y.LE.1200.) GO TO 110 IF (IT.EQ.2) X = X / 15.0D0 IF (IT.NE.2) WRITE (WST,1020,ERR=880) X IF (IT.EQ.2) WRITE (WST,1025,ERR=880) X GO TO 900 C Do in sexagesimal 110 CONTINUE IF (IT.EQ.2) CALL COORDD (1, X, CHM, HM, SEC) IF (IT.NE.2) CALL COORDD (2, X, CHM, HM, SEC) C Y = ABS(CATR(KRCIC+AXPTR)) * 3600.0 Y = ABS(CDELT(1+AXPTR)) * 3600.0 IF ((Y.GT.0.1) .OR. (Y.LE.0.0)) THEN IF (IT.EQ.2) WRITE (WST,1110,ERR=880) CHM, HM, SEC IF (IT.EQ.7) WRITE (WST,1111,ERR=880) CHM, HM, SEC IF ((IT.NE.2) .AND. (IT.NE.7)) * WRITE (WST,1112,ERR=880) CHM, HM, SEC ELSE IF (IT.EQ.2) WRITE (WST,1115,ERR=880) CHM, HM, SEC IF (IT.NE.2) WRITE (WST,1116,ERR=880) CHM, HM, SEC END IF IF (WST(8:8).EQ.' ') WST(8:8) = '0' IF (WST(9:9).EQ.' ') WST(9:9) = '0' GO TO 900 C Stokes 130 CONTINUE NC = AXVAL + 1.5 IF ((NC.LE.0) .OR. (NC.GT.11)) NC = 12 IF ((AXVAL.LT.-0.5) .AND. (AXVAL.GT.-4.5)) NC = -AXVAL + 12.5 STR = CHSTOK(NC) NCH = 4 - 2 * (NC / 13) GO TO 999 C Frequency 140 CONTINUE NA = 1 IF (ABS(AXVAL).LE.1.E6) THEN WRITE (WST,1140,ERR=880) AXVAL ELSE IF(ABS(AXVAL).LE.1.E11) THEN X = AXVAL / 1.E6 WRITE (WST,1145,ERR=880) X ELSE WRITE (WST,1147,ERR=880) AXVAL END IF GO TO 900 C Velocity 150 CONTINUE NA = 1 IF (ABS(AXVAL).GT.1.E3) GO TO 155 WRITE (WST,1150,ERR=880) AXVAL GO TO 900 155 CONTINUE X = AXVAL / 1.E3 WRITE (WST,1155,ERR=880) X GO TO 900 C Time 160 CONTINUE NA = 1 IF (ABS(AXVAL).GT.8.64D4) GO TO 165 WRITE (WST,1160,ERR=880) AXVAL GO TO 900 165 CONTINUE X = AXVAL / 8.64D4 WRITE (WST,1165,ERR=880) X GO TO 900 C Error Handeling 880 CONTINUE C One more Try WRITE(WST,1200,ERR=890) AXVAL 890 CONTINUE C Final string packing 900 CALL XTRIM (WST, 20, WST, NB) IF (NB.GT.21-NA) NB = 21 - NA XSTR(NA:NA+NB-1) = WST(1:NB) NCH = NA + NB - 1 STR = XSTR(1:NCH) GO TO 999 C 999 RETURN C----------------------------------------------------------------------- 1020 FORMAT (F8.3) 1021 FORMAT (F8.5) 1022 FORMAT (F8.6) 1023 FORMAT (F8.0) 1024 FORMAT (1PE11.4) 1025 FORMAT (F8.5) 1110 FORMAT (A1,I2.2,I3.2,F7.3) 1111 FORMAT (A1,I2.2,I3.2,F6.2) 1112 FORMAT (A1,I2.2,I3.2,F5.1) 1115 FORMAT (A1,I2.2,I3.2,F9.5) 1116 FORMAT (A1,I2.2,I3.2,F8.4) 1140 FORMAT (F8.1,' HZ') 1145 FORMAT (F10.3,' MHZ') 1147 FORMAT (1PE11.4,' HZ') 1150 FORMAT (F7.2,' M/S') 1155 FORMAT (F8.1,' KM/S') 1160 FORMAT (F7.1,' SEC') 1165 FORMAT (F7.3,' DAYS') 1200 FORMAT (1PE11.4,' ') END SUBROUTINE XYPIX (X, Y, XPIX, YPIX, IERR) C----------------------------------------------------------------------- C XYPIX determines the pixel location corresponding to a specified C coordinate value. The pixel location is not necessarily an C integer. The position parms are provided by the common DLOC.INC C which requires a previous call to SETLOC. C Inputs: C X D X-coordinate value (header units) C Y D Y-coordinate value (header units) C Output: C XPIX R x-coordinate pixel location C YPIX R y-coordinate pixel location C IERR I 0 ok, 1 out of range, 2 bad type, 3 undefined C----------------------------------------------------------------------- REAL XPIX, YPIX DOUBLE PRECISION X, Y C DOUBLE PRECISION DX, DY, DZ, A, B, C, R, AX, AY, AZ, COSR, SINR INTEGER I, IERR INCLUDE 'DLOC.INC' C----------------------------------------------------------------------- DZ = 0.0D0 R = ROT(LOCNUM) * COND2R IERR = 0 C linear on both axes I = AXFUNC(KLOCL(LOCNUM)+1,LOCNUM) IF ((CORTYP(LOCNUM).GT.0) .AND. (AXTYP(LOCNUM).NE.4) .AND. * (I.GT.1)) GO TO 10 DX = X - RPVAL(1,LOCNUM) DY = Y - RPVAL(2,LOCNUM) IF ((AXTYP(LOCNUM).EQ.2) .OR. (AXTYP(LOCNUM).EQ.3)) DZ = * AXINC(3,LOCNUM) * (ZDEPTH(ZAXIS(LOCNUM)-2,LOCNUM) - * RPLOC(3,LOCNUM)) GO TO 50 C Partly non-linear 10 A = COND2R * RPVAL(1,LOCNUM) B = COND2R * RPVAL(2,LOCNUM) C = COND2R * RPVAL(3,LOCNUM) AX = X * COND2R AY = Y * COND2R IF ((AXTYP(LOCNUM).EQ.2) .OR. (AXTYP(LOCNUM).EQ.3)) DZ = * (ZDEPTH(ZAXIS(LOCNUM)-2,LOCNUM) - RPLOC(3,LOCNUM)) * * AXINC(3,LOCNUM) * COND2R GO TO (15, 20, 25, 30, 35, 40), CORTYP(LOCNUM) C x, y => L, M 15 CALL DIRCOS (I, A, B, AX, AY, DX, DY, IERR) GO TO 45 C x, y => M, L 20 CONTINUE CALL DIRCOS (I, B, A, AY, AX, DY, DX, IERR) GO TO 45 C x, z => L, M 25 CONTINUE CALL DIRDEC (I, DZ, AX, A, C, R, DX, AZ, IERR) DY = AY - B GO TO 45 C x, z => M, L 30 CONTINUE CALL DIRRA (I, DZ, AX, C, A, R, DX, AZ, IERR) DY = AY - B GO TO 45 C y, z => L, M 35 CONTINUE CALL DIRDEC (I, DZ, AY, B, C, R, DY, AZ, IERR) DX = AX - A GO TO 45 C y, z => M, L 40 CONTINUE CALL DIRRA (I, DZ, AY, C, B, R, DY, AZ, IERR) DX = AX - A C back to degrees 45 IF (IERR.NE.0) GO TO 999 DX = DX / COND2R DY = DY / COND2R DZ = DZ / COND2R IF (CORTYP(LOCNUM).GT.2) GO TO 60 C Correct for rotation C (DIRRA, DIRDEC already did) 50 IF ((ROT(LOCNUM).EQ.0.0) .OR. (AXTYP(LOCNUM).EQ.0)) GO TO 60 COSR = COS (R) SINR = SGNROT(LOCNUM) * SIN (R) IF (AXTYP(LOCNUM).EQ.2) DX = DX*COSR + DZ*SINR IF (AXTYP(LOCNUM).EQ.3) DY = DY*COSR + DZ*SINR IF (AXTYP(LOCNUM).EQ.1) THEN DZ = DX*COSR + DY*SINR DY = DY*COSR - DX*SINR DX = DZ END IF C Go to pixels 60 XPIX = DX / AXINC(1,LOCNUM) + RPLOC(1,LOCNUM) YPIX = DY / AXINC(2,LOCNUM) + RPLOC(2,LOCNUM) C Felocity IF ((KLOCF(LOCNUM).EQ.0) .AND. (AXFUNC(1,LOCNUM).EQ.1)) XPIX = DX * / (AXINC(1,LOCNUM) - AXDENU(LOCNUM)*DX) + RPLOC(1,LOCNUM) IF ((KLOCF(LOCNUM).EQ.1) .AND. (AXFUNC(2,LOCNUM).EQ.1)) YPIX = DY * / (AXINC(2,LOCNUM) - AXDENU(LOCNUM)*DY) + RPLOC(2,LOCNUM) C 999 RETURN END SUBROUTINE XYVAL (XPIX, YPIX, X, Y, Z, IERR) C----------------------------------------------------------------------- C XYVAL determines the coordinate value (X,Y,Z) corresponding to the C pixel location (XPIX,YPIX). The pixel values need not be integers. C The necessary map header data is passed via common DLOC.INC C requiring a previous call to SETLOC. This program is the inverse of C XYPIX. C Inputs: C XPIX R Pixel location, x-coordinate C YPIX R Pixel location, y-coordinate C Outputs: C X D X-coordinate value at pixel location C Y D Y-coordinate value at pixel location C Z D Z-coordinate value (if part of a position C pair with either X or Y) C IERR I 0 ok, 1 out of range, 2 bad type, 3 undefined C Common inputs: C DLOC.INC position parms deduced from the map header by C subroutine SETLOC C UNITS ARE AS IN THE MAPHEADER: degrees for position coords C----------------------------------------------------------------------- DOUBLE PRECISION X, Y, Z REAL XPIX, YPIX C REAL COSR, SINR DOUBLE PRECISION A, B, C, DX, DY, DZ, TEMP INTEGER IERR, I INCLUDE 'DLOC.INC' C----------------------------------------------------------------------- IERR = 0 C Offset from ref pixel DX = (XPIX-RPLOC(1,LOCNUM)) * AXINC(1,LOCNUM) DY = (YPIX-RPLOC(2,LOCNUM)) * AXINC(2,LOCNUM) IF ((AXTYP(LOCNUM).EQ.2) .OR. (AXTYP(LOCNUM).EQ.3)) DZ = * (ZDEPTH(ZAXIS(LOCNUM)-2,LOCNUM) - RPLOC(3,LOCNUM)) * * AXINC(3,LOCNUM) C Take out rotation IF ((ROT(LOCNUM).EQ.0.0) .OR. (AXTYP(LOCNUM).LE.0) .OR. * (AXTYP(LOCNUM).GE.4)) GO TO 50 COSR = COS(ROT(LOCNUM)*COND2R) SINR = SGNROT(LOCNUM) * SIN(ROT(LOCNUM)*COND2R) GO TO (10, 20, 30), AXTYP(LOCNUM) 10 CONTINUE TEMP = DX * COSR - DY * SINR DY = DY * COSR + DX * SINR DX = TEMP GO TO 50 20 CONTINUE TEMP = DX * COSR - DZ * SINR DZ = DZ * COSR + DX * SINR DX = TEMP GO TO 50 30 CONTINUE TEMP = DY * COSR - DZ * SINR DZ = DZ * COSR + DY * SINR DY = TEMP C Linear transformation 50 Z = 0.0 I = AXFUNC(KLOCL(LOCNUM)+1,LOCNUM) IF ((CORTYP(LOCNUM).GT.0) .AND. (AXTYP(LOCNUM).NE.4) .AND. * (I.GT.1)) GO TO 60 IF ((AXTYP(LOCNUM).GT.1) .AND. (AXTYP(LOCNUM).LT.4)) Z = * RPVAL(3,LOCNUM) + DZ X = RPVAL(1,LOCNUM) + DX Y = RPVAL(2,LOCNUM) + DY GO TO 100 C Partly non-linear 60 A = RPVAL(1,LOCNUM) * COND2R B = RPVAL(2,LOCNUM) * COND2R C = RPVAL(3,LOCNUM) * COND2R DX = DX * COND2R DY = DY * COND2R DZ = DZ * COND2R GO TO (65, 70, 75, 80, 85, 90), CORTYP(LOCNUM) C x, y => L,M 65 CALL NEWPOS (I, A, B, DX, DY, X, Y, IERR) GO TO 95 C x, y => M, L 70 CONTINUE CALL NEWPOS (I, B, A, DY, DX, Y, X, IERR) GO TO 95 C x, z => L, M 75 CONTINUE CALL NEWPOS (I, A, C, DX, DZ, X, Z, IERR) Y = B + DY GO TO 95 C x, z => M, L 80 CONTINUE CALL NEWPOS (I, C, A, DZ, DX, Z, X, IERR) Y = B + DY GO TO 95 C y, z => L, M 85 CONTINUE CALL NEWPOS (I, B, C, DY, DZ, Y, Z, IERR) X = A + DX GO TO 95 C y, z => M, L 90 CONTINUE CALL NEWPOS (I, C, B, DZ, DY, Z, Y, IERR) X = A + DX C correct units back 95 X = X / COND2R Y = Y / COND2R Z = Z / COND2R DX = DX / COND2R DY = DY / COND2R C Felocity 100 IF ((KLOCF(LOCNUM).EQ.0) .AND. (AXFUNC(1,LOCNUM).EQ.1)) X = * RPVAL(1,LOCNUM) + DX / (1.0D0 + AXDENU(LOCNUM) * * (XPIX-RPLOC(1,LOCNUM))) IF ((KLOCF(LOCNUM).EQ.1) .AND. (AXFUNC(2,LOCNUM).EQ.1)) Y = * RPVAL(2,LOCNUM) + DY / (1.0D0 + AXDENU(LOCNUM) * * (YPIX-RPLOC(2,LOCNUM))) C 999 RETURN END SUBROUTINE COORDT (ITI, ITO, LONGI, LATI, EPOK, LONGO, LATO, IERR) C----------------------------------------------------------------------- C COORDT translates between types of coordinates: C Inputs: ITI I Input type (1 Ra, Dec; 2 gal, 3 ecliptic) C ITO I Output type C LONGI D Input longitude C LATI D Input latitude C EPOK R Epoch of coordinates (no precession done) C 1950 or 2000 assumed in Galactic conversions!!! C Output: LONGO D Output longitude C LATI D Output latitude C IERR I error: 0 ok, 1 input error, 2 conversion fails C----------------------------------------------------------------------- INTEGER ITI, ITO, IERR REAL EPOK DOUBLE PRECISION LONGI, LATI, LONGO, LATO C define PI/180 and 2PI DOUBLE PRECISION DEGRAD, CIRC PARAMETER (DEGRAD=1.745329251994329576D-2, * CIRC = 6.2831853071795865788D0) DOUBLE PRECISION RAGP, DECGP, LONCP, DEC, RA, DT, * GLAT, GLON, ELAT, ELON, RAGP0, DECGP0, LONCP0, DEPS0(4), EPS DATA DEPS0 /23.452294D0, -0.0130125D0, -1.64D-6, 5.03D-7/ C----------------------------------------------------------------------- IERR = 1 IF ((ITI.LT.1) .OR. (ITI.GT.3)) GO TO 999 IF ((ITO.LT.1) .OR. (ITO.GT.3)) GO TO 999 IERR = 0 C null cases LONGO = LONGI LATO = LATI IF (ITI.EQ.ITO) GO TO 999 DT = 0.0D0 IF ((ITI.EQ.3) .OR. (ITO.EQ.3)) DT = (EPOK - 1900.0) / 100.0 EPS = DEPS0(1) + DT * (DEPS0(2) + DT*(DEPS0(3) + DT*DEPS0(4))) EPS = EPS * DEGRAD C Get galactic coords for Epoch CALL GALPOL(EPOK, RAGP0, DECGP0, LONCP0) RAGP = RAGP0 * DEGRAD DECGP = DECGP0 * DEGRAD LONCP = LONCP0 * DEGRAD GO TO (100, 200, 300), ITI C Ra, Dec in 100 DEC = LATI * DEGRAD RA = LONGI * DEGRAD C Galactic out 110 IF (ITO.EQ.3) GO TO 150 GLAT = ASIN (SIN(DEC)*SIN(DECGP) + COS(DEC)*COS(DECGP)* * COS(RA-RAGP)) GLON = LONCP + ATAN2 (COS(DEC)*SIN(RAGP-RA), SIN(DEC)* * COS(DECGP) - COS(DEC)*SIN(DECGP)*COS(RA-RAGP)) IF (GLON.GE.CIRC) GLON = GLON - CIRC IF (GLON.LT.0.0D0) GLON = GLON + CIRC LATO = GLAT / DEGRAD LONGO = GLON / DEGRAD GO TO 999 C Ecliptic out 150 CONTINUE ELAT = ASIN (SIN(DEC)*COS(EPS) - COS(DEC)*SIN(EPS)* * SIN(RA)) ELON = ATAN2 (SIN(DEC)*SIN(EPS) - * COS(DEC)*COS(EPS)*SIN(RA), COS(DEC)*COS(RA)) IF (ELON.GE.CIRC) ELON = ELON - CIRC IF (ELON.LT.0.0D0) ELON = ELON + CIRC LATO = ELAT / DEGRAD LONGO = ELON / DEGRAD GO TO 999 C Galactic in 200 GLAT = LATI * DEGRAD GLON = LONGI * DEGRAD DEC = ASIN (SIN(GLAT)*SIN(DECGP) + COS(GLAT)*COS(DECGP)* * COS(GLON-LONCP)) RA = RAGP + ATAN2 (COS(GLAT)*SIN(LONCP-GLON), SIN(GLAT)* * COS(DECGP) - COS(GLAT)*SIN(DECGP)*COS(GLON-LONCP)) C Ra, Dec out 210 IF (ITO.EQ.3) GO TO 150 IF (RA.GE.CIRC) RA = RA - CIRC IF (RA.LT.0.0D0) RA = RA + CIRC LONGO = RA / DEGRAD LATO = DEC / DEGRAD GO TO 999 C Ecliptic in 300 ELAT = LATI * DEGRAD ELON = LONGI * DEGRAD DEC = ASIN (SIN(ELAT)*COS(EPS) + * COS(ELAT)*SIN(EPS)*SIN(ELON)) RA = ATAN2 (COS(ELAT)*COS(EPS)*SIN(ELON) - * SIN(ELAT)*SIN(EPS), COS(ELAT)*COS(ELON)) C Ra, Dec out IF (ITO.EQ.2) GO TO 110 GO TO 210 C 999 RETURN 2000 FORMAT (A) END SUBROUTINE DIRCOS (TYPE, RA0, DEC0, RA, DEC, L, M, IERR) C----------------------------------------------------------------------- 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 DEPS, COSS, SINS, DT, DA, DD, SINT, PI, TWOPI INCLUDE 'DLOC.INC' PARAMETER (PI = 3.141592653589793238462643D0) PARAMETER (TWOPI = 2.0D0 * PI) DATA DEPS /1.0D-5/ C----------------------------------------------------------------------- C Use full accuracy 10 IERR = 2 IF ((TYPE.LT.2) .OR. (TYPE.GT.9)) 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) GO TO (20, 20, 30, 40, 50, 60, 70, 80, 90), TYPE C SIN projection 20 CONTINUE IF (SINT.LT.0.0D0) IERR = 1 M = SINS * COS(DEC0) - COSS * SIN(DEC0) * COS(DT) GO TO 999 C TAN projection 30 CONTINUE IF (SINT.LE.0.0D0) THEN IERR = 1 ELSE M = SINS * SIN(DEC0) + COSS * COS(DEC0) * COS(DT) L = L / M M = (SINS * COS(DEC0) - COSS * SIN(DEC0) * COS(DT)) / M END IF GO TO 999 C Arc projection 40 CONTINUE M = SINS * SIN(DEC0) + COSS * COS(DEC0) * COS(DT) M = MIN (1.0D0, MAX (-1.0D0, M)) M = ACOS (M) IF (M.NE.0.0D0) THEN M = M / SIN(M) ELSE M = 1.0D0 END IF L = L * M M = (SINS * COS(DEC0) - COSS * SIN(DEC0) * COS(DT)) * M GO TO 999 C WSRT projection 50 CONTINUE IF (DEC0.NE.0.0D0) GO TO 55 WRITE (0,1050) L = 0.0D0 M = 0.0D0 GO TO 990 55 M = (COS(DEC0) - COSS * COS(DT)) / SIN(DEC0) GO TO 999 C Global sinusoidal 60 CONTINUE M = DEC - DEC0 L = DT * COSS GO TO 999 C Mercator 70 CONTINUE L = GEOMD1(LOCNUM) * DT DT = DEC / 2.0D0 + PI / 4.0D0 DT = TAN (DT) IF (DT.LT.DEPS) GO TO 990 M = GEOMD2(LOCNUM) * LOG (DT) - GEOMD3(LOCNUM) GO TO 999 C Aitoff 80 CONTINUE L = 0.0D0 M = 0.0D0 DA = DT / 2.0D0 DT = SQRT ((1.0D0 + COS(DEC) * COS(DA))/2.0D0) IF (ABS(DT).LT.DEPS) GO TO 990 L = 2.0D0 * GEOMD1(LOCNUM) * COS(DEC) * SIN(DA) / DT M = GEOMD2(LOCNUM) * SIN(DEC) / DT - GEOMD3(LOCNUM) GO TO 999 C Stereographic 90 CONTINUE DD = 1.0D0 + SINS * SIN(DEC0) + COSS * COS(DEC0) * COS(DT) IF (ABS(DD).LT.DEPS) GO TO 990 DD = 2.0D0 / DD L = L * DD M = DD * (SINS * COS(DEC0) - COSS * SIN(DEC0) * COS(DT)) GO TO 999 C Undefined 990 IERR = 3 C 999 RETURN C----------------------------------------------------------------------- 1050 FORMAT ('DIRCOS: WSRT COORDINATES CANNOT HAVE DEC0 = 0') END SUBROUTINE DIRDEC (TYPE, DY, RA, RA0, DEC0, ROTA, DX, DEC, IERR) C----------------------------------------------------------------------- C DIRDEC will find the latitude and longitude-like pixel position C of a point given the longitude and the latitude-pixel of the point. C ALL ANGLES ARE IN RADIANS. C Inputs: C TYPE I 2 => SIN, 3 => TAN, 4 => Arc, 5 => WSRT geometry C 6 Global sinusoidal, 7 Mercator, 8 Aitoff C DY D the Y pixel pos rel to ref pixel C RA D right ascension/longitude C RA0 D RA/longitude of ref pixel C DEC0 D DEC/latitude of ref pixel C ROT D Rotation of coords +dec into +ra C Output: C DEC D declination/latitude C DX D RA-like pixel position rel to ref pixel C IERR I 0 ok, 1 out of range, 2 bad type, 3 undefined C----------------------------------------------------------------------- DOUBLE PRECISION DY, RA, RA0, DEC0, ROTA, DX, DEC INTEGER IERR, TYPE C INTEGER NC, NCLIM DOUBLE PRECISION COSR, SINR, DA, DB, DC, DD, COS0, SIN0, COSA, * SINA, DEPS, DE, TWOPI, DT, TRA, FG, DFDD, DSTEP LOGICAL DOIT INCLUDE 'DLOC.INC' DATA TWOPI, DEPS /6.28318530717959D0, 1.D-5/ DATA NCLIM /1000/ C----------------------------------------------------------------------- DX = 0.0D0 DEC = 0.0D0 IERR = 2 IF ((TYPE.LT.2) .OR. (TYPE.GT.9)) GO TO 999 IERR = 0 COSR = COS(ROTA) SINR = SIN(ROTA) COS0 = COS(DEC0) SIN0 = SIN(DEC0) COSA = COS(RA-RA0) SINA = SIN(RA-RA0) NC = 0 GO TO (20, 20, 30, 40, 50, 60, 70, 80, 90), TYPE C SIN projection 20 CONTINUE DA = COSR * COS0 DB = COSR * SIN0 * COSA + SINR * SINA IF ((ABS(DA).LT.DEPS) .AND. (ABS(DB).LT.DEPS)) GO TO 990 DT = DY / SQRT(DB*DB + DA*DA) IF (ABS(DT).GT.1.0D0) GO TO 990 DEC = ASIN (DT) + ATAN2 (DB, DA) DX = SINR * SIN(DEC) * COS0 - COS(DEC) * * (SINR * SIN0 * COSA - COSR * SINA) GO TO 900 C TAN projection 30 CONTINUE DA = COS0 * COSR - DY * SIN0 IF (ABS(DA).LT.DEPS) GO TO 990 DB = (SIN0 * COSA * COSR + SINA * SINR + DY * COS0 * COSA) / DA DEC = ATAN (DB) DX = (SINA*COSR + DB*COS0*SINR - SIN0*COSA*SINR) / * (DB * SIN0 + COS0 * COSA) GO TO 900 C Arc projection 40 CONTINUE IF (ABS(COSR).LT.DEPS) GO TO 990 DD = ATAN2 (SIN0, COS0 * COSA) DC = SQRT (1.D0 - COS0*COS0*SINA*SINA) DX = (DY * SINR + COS0 * SINA) / COSR NC = 0 45 DA = SQRT (DX*DX + DY*DY) NC = NC + 1 DE = DX DB = 1.0D0 IF (SIN(DA).NE.0.0D0) DB = DA / SIN(DA) DA = COS(DA) / DC DEC = DD IF (ABS(DA).LT.1.0D0) DEC = DEC + ACOS (DA) IF ((DY.LT.0.0D0) .AND. (DEC.GT.DEC0)) DEC = -DEC * + 2.0D0 * DD DX = (DY * SINR + DB * COS(DEC) * SINA) / COSR IF ((ABS(DX-DE).LT.1.0D-9) .OR. (NC.GT.NCLIM)) GO TO 900 IF (NC.GT.10) DX = (DX + DE) / 2.0D0 IF ((NC.GT.25) .AND. ((NC/2)*2.EQ.NC)) DX = (DX+DE) / 2.0D0 GO TO 45 C WSRT projection 50 CONTINUE IF (ABS(SIN0).LT.DEPS) GO TO 990 DA = COSA * COSR + SINA * SIN0 * SINR IF (ABS(DA).LT.DEPS) GO TO 990 DT = (COS0 * COSR - DY * SIN0) / DA IF (ABS(DT).GT.1.0D0) GO TO 990 DEC = ACOS (DT) IF (DEC0.LT.0.0D0) DEC = -DEC DX = COS(DEC) * SINA * COSR + SINR * (COS0 - COS(DEC)*COSA) * / SIN0 GO TO 900 C Global sinusoidal 60 CONTINUE IF (ABS(COSR).LT.DEPS) GO TO 990 TRA = RA - RA0 IF (ABS(TRA).GT.TWOPI/2.0D0) GO TO 990 NC = 0 DEC = DEC0 + DY / COSR IF (ABS(SINR).LT.DEPS) GO TO 67 DC = -100.0D0 65 NC = NC + 1 DSTEP = 0.05D0 IF (NC.LT.4) DSTEP = 0.20D0 / NC FG = DY + TRA*SINR*COS(DEC) - (DEC-DEC0)*COSR DFDD = COSR + TRA*SINR*SIN(DEC) DOIT = (ABS(DFDD).GT.DEPS) .OR. ((ABS(DFDD).GT.DABS(FG)) * .AND. (DFDD.NE.0.0D0)) IF ((.NOT.DOIT) .AND. (DC.LT.-50.0D0)) GO TO 990 IF (.NOT.DOIT) DEC = (DEC + DC) / 2.0D0 IF (DOIT) DC = DEC IF (DOIT) DEC = DEC + FG / DFDD IF (DEC-DC.GT.DSTEP) DEC = DC + DSTEP IF (DEC-DC.LT.-DSTEP) DEC = DC - DSTEP IF (NC.GT.NCLIM) GO TO 990 IF (ABS(DEC-DC).GT.1.0D-9) GO TO 65 67 DX = TRA * COS(DEC) * COSR + (DEC - DEC0) * SINR IF (ABS(DEC-DEC0).GT.TWOPI/2.0D0) GO TO 990 IF (ABS(DEC).GT.TWOPI/4.0D0) GO TO 990 GO TO 900 C Mercator 70 CONTINUE IF (ABS(COSR).LT.DEPS) GO TO 990 TRA = RA - RA0 IF (ABS(TRA).GT.TWOPI/2.0D0) GO TO 990 DA = (DY + GEOMD1(LOCNUM) * TRA * SINR) / COSR DX = GEOMD1(LOCNUM) * TRA * COSR + DA * SINR DA = (DA + GEOMD3(LOCNUM)) / GEOMD2(LOCNUM) DEC = 2.0D0 * ATAN (EXP(DA)) - TWOPI / 4.0D0 GO TO 900 C Aitoff 80 CONTINUE DE = (RA - RA0) / 2.0D0 IF (ABS(DE).GT.TWOPI/4.0D0) GO TO 990 COS0 = COS (DE) SIN0 = SIN (DE) IF (ABS(COSR).LT.DEPS) GO TO 990 DD = DEC0 IF (DY.GT.1.0D0) DD = 1.0D0 + DEC0 IF (DY.LT.-1.0D0) DD = -1.0D0 + DEC0 NC = 0 DC = -100.0 85 DB = SQRT ((1.0D0 + COS(DD) * COS0)/2.0D0) IF (ABS(DB).LT.DEPS) GO TO 990 NC = NC + 1 DSTEP = 0.05D0 IF (NC.LT.4) DSTEP = 0.20D0 / NC FG = (DY + GEOMD3(LOCNUM) * COSR) * DB + 2.0D0 * SINR * SIN0 * * COS(DD) * GEOMD1(LOCNUM) - GEOMD2(LOCNUM) * COSR * * SIN(DD) DFDD = (DY + GEOMD3(LOCNUM) * COSR) * COS0 * SIN(DD) / * (4.0D0 * DB) + 2.0D0 * GEOMD1(LOCNUM) * SINR * SIN0 * * SIN(DD) + GEOMD2(LOCNUM) * COSR * COS(DD) DOIT = (ABS(DFDD).GT.DEPS) .OR. ((ABS(DFDD).GT.DABS(FG)) * .AND. (DFDD.NE.0.0D0)) IF ((.NOT.DOIT) .AND. (DC.LT.-50.0D0)) GO TO 990 IF (.NOT.DOIT) DD = (DD + DC) / 2.0D0 IF (DOIT) DC = DD IF (DOIT) DD = DD + FG / DFDD IF (DD-DC.GT.DSTEP) DD = DC + DSTEP IF (DD-DC.LT.-DSTEP) DD = DC - DSTEP IF (NC.GT.NCLIM) GO TO 990 IF (ABS(DD-DC).GT.1.0D-9) GO TO 85 87 DEC = DD IF (ABS(DEC).GT.TWOPI/4.0D0) GO TO 990 IF (DB.LT.DEPS) GO TO 990 DX = (2.0D0 * SIN0 * COS(DD) * GEOMD1(LOCNUM) / DB + DY * SINR) * / COSR GO TO 900 C Stereographic 90 CONTINUE DA = 2.0D0 * COS0 * COSR - DY * SIN0 DB = COSA * (2.0D0 * SIN0 * COSR + DY * COS0) + * 2.0D0 * SINA * SINR DT = SQRT (DA * DA + DB * DB) IF (DT.LT.DEPS) GO TO 990 DT = DY / DT IF (ABS(DT).GT.1.0D0) GO TO 990 DEC = ATAN2 (DB, DA) + ASIN (DT) C check DC = 1.0D0 + SIN(DEC)*SIN0 + COS(DEC)*COS0*COSA FG = -1.0D10 IF (ABS(DC).GT.DEPS) FG = 2.0D0 * ((SIN(DEC)*COS0 - * COS(DEC)*SIN0*COSA)*COSR - COS(DEC)*SINA*SINR) / DC IF (ABS(FG-DY).GT.DEPS) DEC = ATAN2 (DB, DA) - ASIN (DT) * + TWOPI/2.0D0 DX = 1.0D0 + SIN(DEC) * SIN0 + COS(DEC) * COS0 * COSA IF (ABS(DX).LT.DEPS) GO TO 990 FG = -1.0D10 IF (ABS(DC).GT.DEPS) FG = 2.0D0 * ((SIN(DEC)*COS0 - * COS(DEC)*SIN0*COSA)*COSR - COS(DEC)*SINA*SINR) / DX IF (ABS(FG-DY).GT.DEPS) GO TO 990 DX = 2.0D0 * (COS(DEC) * SINA * COSR + * (SIN(DEC)*COS0 - COS(DEC)*SIN0*COSA) * SINR) / DX GO TO 900 C Check the loop 900 IF (NC.LE.NCLIM/2) GO TO 999 IF (NC.GT.NCLIM) THEN WRITE (0,1901) NC ELSE WRITE (0,1900) NC END IF IF (NC.GT.NCLIM) IERR = 3 GO TO 999 C Undefined 990 IERR = 3 C WRITE (0,1990) C 999 RETURN C----------------------------------------------------------------------- 1900 FORMAT ('DIRDEC: ITERATIONS REQUIRED = ',I5) 1901 FORMAT ('DIRDEC: FAILED TO CONVERGE AFTER',I5,' ITERATIONS') 1990 FORMAT ('DIRDEC: OPERATION UNDEFINED AT THE POLE OR BECAUSE OF', * ' ROTATION') END SUBROUTINE DIRRA (TYPE, DX, DEC, RA0, DEC0, ROTA, DY, RA, IERR) C----------------------------------------------------------------------- C DIRRA will find the longitude and the latitude pixel position of a C point given the latitude and the longitude pixel of the point. C For use with four projective geometries. C ALL ANGLES ARE IN RADIANS. C Inputs: C TYPE I 2 => SIN, 3 => TAN, 4 => Arc, 5 => WSRT geometry C 6 Global sinusoidal, 7 Mercator, 8 Aitoff C DX D X pixel value from ref pixel C DEC D declination/ latitude C RA0 D RA / longitude of ref pixel C DEC0 D DEC / latitude of ref pixel C ROTA D Rotation of coordinates C Output: C DY D Y pixel position C RA D right ascension/ longitude C IERR D 0 ok, 1 out of range, 2 bad type, 3 undefined answer C----------------------------------------------------------------------- INTEGER IERR, TYPE DOUBLE PRECISION DX, DEC, RA0, DEC0, ROTA, DY, RA C INTEGER NC, NCLIM DOUBLE PRECISION COSR, SINR, DA, DB, DC, COS0, SIN0, COSD, * SIND, TWOPI, DEPS, DT, FG, DFDA, DSTEP LOGICAL DOIT INCLUDE 'DLOC.INC' DATA TWOPI /6.28318530717959D0/ DATA NCLIM, DEPS /1000, 1.0D-5/ C----------------------------------------------------------------------- DY = 0.0D0 RA = 0.0D0 IERR = 2 IF ((TYPE.LT.2) .OR. (TYPE.GT.9)) GO TO 999 IERR = 0 COSR = COS(ROTA) SINR = SIN(ROTA) COSD = COS(DEC) SIND = SIN(DEC) COS0 = COS(DEC0) SIN0 = SIN(DEC0) NC = 0 IF (ABS(COSD).LT.DEPS) GO TO 990 GO TO (20, 20, 30, 40, 50, 60, 70, 80, 90), TYPE C SIN projection 20 CONTINUE DA = COSR DB = SINR * SIN0 DT = SQRT (DA*DA + DB*DB) IF (DT.LT.DEPS) GO TO 990 DT = (DX - SINR*SIND*COS0) / DT / COSD IF (ABS(DT).GT.1.0D0) GO TO 990 RA = ASIN (DT) + ATAN2 (DB, DA) + RA0 DY = SIND * COS0 * COSR - COSD * (COSR * SIN0 * COS(RA-RA0) * + SIN(RA-RA0) * SINR) GO TO 900 C TAN projection 30 CONTINUE DA = COSR DB = SINR * SIN0 + DX * COS0 DT = SQRT (DA*DA + DB*DB) IF (DT.LT.DEPS) GO TO 990 DT = (DX*SIN0 - COS0*SINR) / DT * SIND / COSD IF (ABS(DT).GT.1.0D0) GO TO 990 RA = ASIN (DT) + ATAN2 (DB, DA) + RA0 DT = SIND * SIN0 + COSD * COS0 * COS(RA-RA0) IF (ABS(DT).LT.DEPS) GO TO 990 DY = (SIND*COS0*COSR - COSD*SIN0*COS(RA-RA0)*COSR - * COSD*SINR*SIN(RA-RA0)) / DT GO TO 900 C Arc projection 40 CONTINUE IF (ABS(COSR).LT.DEPS) GO TO 990 IF (ABS(COS0).GT.DEPS) GO TO 42 IF (ABS(SIN0).LT.DEPS) GO TO 990 DY = ACOS (SIND / SIN0) DY = DY*DY - DX*DX IF (DY.LT.0.0D0) GO TO 990 DY = -SQRT (DY) DA = SQRT (DX*DX + DY*DY) DB = 1.0D0 IF (DA.NE.0.0D0) DB = SIN(DA) / DA GO TO 47 42 DY = 0.0D0 NC = 0 45 DA = SQRT (DX*DX + DY*DY) NC = NC + 1 DC = DY DB = 1.0D0 IF (DA.NE.0.0D0) DB = SIN(DA) / DA DA = COS(DA) DY = (SIND - SIN0*DA - DX*SINR*COS0*DB) / (COSR * COS0 * * DB) DA = ABS (DY - DC) IF ((DA.LT.1.0D-9) .OR. (NC.GT.NCLIM)) GO TO 47 IF (NC.GT.10) DY = (DY + DC) / 2.0D0 GO TO 45 47 DT = DB * (DX*COSR - DY*SINR) / COSD IF (ABS(DT).GT.1.0D0) GO TO 990 RA = RA0 + ASIN (DT) GO TO 900 C WSRT projection 50 CONTINUE IF (ABS(SIN0).LT.DEPS) GO TO 990 DA = COSR DB = SINR / SIN0 DT = (DX - COS0*DB) / SQRT (DA*DA + DB*DB) / COSD IF (ABS(DT).GT.1.0D0) GO TO 990 RA = ASIN (DT) + ATAN2 (DB, DA) + RA0 DY = (COS0 - COSD*COS(RA-RA0)) * COSR / SIN0 - * COSD * SINR * SIN(RA-RA0) GO TO 900 C Global sinusoidal 60 CONTINUE IF (ABS(COSR).LT.DEPS) GO TO 990 DY = (DEC - DEC0 - DX * SINR) / COSR RA = RA0 + (DX * COSR - DY * SINR) / COSD IF (ABS(RA-RA0).GT.TWOPI/2.0D0) GO TO 990 GO TO 900 C Mercator 70 CONTINUE IF (ABS(COSR).LT.DEPS) GO TO 990 DT = TAN (DEC/2.0D0+TWOPI/8.0D0) IF (DT.LT.DEPS) GO TO 990 DT = GEOMD2(LOCNUM) * LOG (DT) - GEOMD3(LOCNUM) DA = (DX - DT * SINR) / COSR DY = DT * COSR - DA * SINR RA = RA0 + DA / GEOMD1(LOCNUM) GO TO 900 C Aitoff 80 CONTINUE NC = 0 COSD = COS (DEC) SIND = SIN (DEC) IF ((ABS(COSR).LT.DEPS) .OR. (ABS(COSD).LT.DEPS)) GO TO 990 DA = 0.0D0 IF (DX.GT.1.0D0) DA = 1.0D0 IF (DX.LT.-1.0D0) DA = -1.0D0 DC = -100.0D0 85 DB = SQRT ((1.0D0 + COSD * COS(DA))/2.0D0) NC = NC + 1 DSTEP = 0.05D0 IF (NC.LT.5) DSTEP = 0.25D0 / NC IF (ABS(DB).LT.DEPS) GO TO 990 FG = DB * DX + SINR * (GEOMD3(LOCNUM) * DB - SIND * * GEOMD2(LOCNUM)) - (2.0D0 * GEOMD1(LOCNUM) * COSR * COSD) * * SIN(DA) DFDA = (DX + SINR * GEOMD3(LOCNUM)) * COSD * SIN(DA) / * (4.0D0 * DB) + 2.0D0 * GEOMD1(LOCNUM) * COSD * COSR * * COS(DA) DOIT = (ABS(DFDA).GT.DEPS) .OR. ((ABS(DFDA).GT.DABS(FG)) * .AND. (DFDA.NE.0.0D0)) IF ((.NOT.DOIT) .AND. (DC.LT.-50.0D0)) GO TO 990 IF (.NOT.DOIT) DA = (DA + DC) / 2.0D0 IF (DOIT) DC = DA IF (DOIT) DA = DA + FG / DFDA IF (DA-DC.GT.DSTEP) DA = DC + DSTEP IF (DA-DC.LT.-DSTEP) DA = DC - DSTEP IF (NC.GT.NCLIM) GO TO 990 IF (ABS(DA-DC).GT.1.0D-9) GO TO 85 87 RA = RA0 + 2.0D0 * DA IF (ABS(RA-RA0).GT.TWOPI/2.0D0) GO TO 990 IF (DB.LT.DEPS) GO TO 990 DY = (SIND * GEOMD2(LOCNUM) / DB - DX * SINR) / COSR GO TO 900 C Stereographic 90 CONTINUE DA = 2.0D0 * COSD * COSR DB = 2.0D0 * COSD * SIN0 * SINR + DX * COSD * COS0 DT = SQRT (DA*DA + DB*DB) IF (DT.LT.DEPS) GO TO 990 DT = (DX + DX*SIND*SIN0 - 2.0D0*SIND*COS0*SINR) / DT IF (ABS(DT).GT.1.0D0) GO TO 990 RA = ASIN (DT) + ATAN2 (DB, DA) DY = 1.0D0 + SIND*SIN0 + COSD*COS0*COS(RA) IF (ABS(DY).LE.DEPS) GO TO 990 DY = 2.0D0 * ((SIND*COS0 - COSD*SIN0*COS(RA)) * COSR - * COSD*SIN(RA) * SINR) / DY RA = RA + RA0 GO TO 900 C RA in range: 900 IF (RA-RA0.LT.-TWOPI/2.0D0) RA = RA + TWOPI IF (RA-RA0.GT.TWOPI/2.0D0) RA = RA - TWOPI IF (NC.LE.NCLIM/2) GO TO 999 IF (NC.GT.NCLIM) THEN WRITE (0,1901) NC ELSE WRITE (0,1900) NC END IF IF (NC.GT.NCLIM) IERR = 3 GO TO 999 C Undefined 990 IERR = 3 C 999 RETURN C----------------------------------------------------------------------- 1900 FORMAT ('DIRRA: ITERATIONS REQUIRED = ',I5) 1901 FORMAT ('DIRRA: FAILED TO CONVERGE AFTER',I5,' ITERATIONS') END SUBROUTINE NEWPOS (TYPE, RA0, DEC0, L, M, RAOUT, DECOUT, IERR) C----------------------------------------------------------------------- C NEWPOS determines the coordinates (RAOUT,DECOUT) corresponding to a C displacement (L,M) given by the direction cosines from coordinates C (RA0,DEC0). The direction cosine L is assumed to be positive to C the east; M is positive to the north. The routine works for C 4 kinds of projective geometries and for Celestial, Ecliptic, or C Galactic coordinate systems. C This subroutine always uses an accurate computation. C Inputs: C TYPE I 2 = SIN projection, 3 = TAN projection, 4 = arc C projection, 5 = projection to north pole oriented C plane (i.e. WSRT) C RA0 D Coordinate reference right ascension (longitude) C DEC0 D Coordinate reference declination (latitude) C L D Cosine angle of displacement to east C M D Cosine angle of displacement to north C Outputs: C RAOUT D right ascension or longitude at (L,M) C DECOUT D declination or latitude at (L,M) C IERR I Error condition: 0 = ok, 1 = L,M crazy, 2 = bad C type, 3 = answer undefined C ALL ANGLES IN THIS SUBROUTINE ARE IN RADIANS. C----------------------------------------------------------------------- INTEGER TYPE, IERR DOUBLE PRECISION RA0, DEC0, L, M, RAOUT, DECOUT C DOUBLE PRECISION SINS, COSS, TWOPI, DECT, RAT, DT, DEPS, MG, DA, * DD, DZ, COS0, SIN0 INCLUDE 'DLOC.INC' DATA TWOPI, DEPS /6.28318530717959D0, 1.D-5/ C----------------------------------------------------------------------- IERR = 2 IF ((TYPE.GE.2) .AND. (TYPE.LE.9)) GO TO 5 WRITE (0,1000) TYPE GO TO 999 5 IERR = 1 SINS = L*L + M*M IF (SINS.LE.1.0D0) GO TO 10 IF (TYPE.EQ.9) GO TO 10 IF (TYPE.EQ.6) GO TO 10 IF ((TYPE.EQ.4) .AND. (SINS.LT.TWOPI*TWOPI/4.0D0)) GO TO 10 WRITE (0,1005) L, M GO TO 999 10 IERR = 3 DECOUT = 0.D0 RAOUT = 0.D0 COS0 = COS(DEC0) SIN0 = SIN(DEC0) C Use accurate solution GO TO (20, 20, 30, 40, 50, 60, 70, 80, 90), TYPE C SIN projection 20 CONTINUE COSS = SQRT (1.0D0 - SINS) DT = SIN0 * COSS + COS0 * M IF (ABS(DT).GT.1.0D0) GO TO 999 DECT = ASIN (DT) RAT = COS0 * COSS - SIN0 * M IF ((RAT.EQ.0.D0) .AND. (L.EQ.0.0D0)) GO TO 999 RAT = ATAN2 (L, RAT) + RA0 GO TO 900 C TAN projection 30 CONTINUE DECT = COS0 - M * SIN0 IF (DECT.EQ.0.0D0) GO TO 999 RAT = RA0 + ATAN2 (L, DECT) DECT = ATAN (COS(RAT-RA0) * (M * COS0 + SIN0) / DECT) GO TO 900 C Arc projection 40 CONTINUE SINS = SQRT(SINS) COSS = COS (SINS) IF (SINS.NE.0.0D0) THEN SINS = SIN (SINS) / SINS ELSE SINS = 1.0D0 END IF DT = M * COS0 * SINS + SIN0 * COSS IF (ABS(DT).GT.1.0D0) GO TO 999 DECT = ASIN (DT) DA = COSS - DT * SIN0 DT = L * SINS * COS0 IF ((DA.EQ.0.0D0) .AND. (DT.EQ.0.0D0)) GO TO 999 RAT = RA0 + ATAN2 (DT, DA) GO TO 900 C WSRT (North pole projection) 50 CONTINUE DECT = COS0 - M * SIN0 IF (DECT.EQ.0.0D0) GO TO 999 RAT = RA0 + ATAN2 (L, DECT) DT = COS (RAT-RA0) IF (DT.EQ.0.0D0) GO TO 999 DECT = DECT / DT IF (ABS(DECT).GT.1.0D0) GO TO 999 DECT = ACOS (DECT) IF (DEC0.LT.0.0D0) DECT = -DECT GO TO 900 C Global sinusoid 60 CONTINUE DECT = DEC0 + M IF (ABS(DECT).GT.TWOPI/4.0D0) GO TO 999 COSS = COS (DECT) IF (ABS(L).GT.TWOPI*COSS/2.0D0) GO TO 999 RAT = RA0 IF (COSS.GT.DEPS) RAT = RAT + L / COSS GO TO 900 C Mercator 70 CONTINUE RAT = L / GEOMD1(LOCNUM) + RA0 IF (ABS(RAT-RA0).GT.TWOPI) GO TO 999 DT = 0.0D0 IF (GEOMD2(LOCNUM).NE.0.0D0) DT = (M + GEOMD3(LOCNUM)) / * GEOMD2(LOCNUM) DT = EXP (DT) DECT = 2.0D0 * ATAN (DT) - TWOPI / 4.0D0 GO TO 900 C Aitoff 80 CONTINUE RAT = RA0 DECT = DEC0 IF ((L.EQ.0.0D0) .AND. (M.EQ.0.0D0)) GO TO 900 DZ = 4.0D0 - L*L/(4.0D0*GEOMD1(LOCNUM)*GEOMD1(LOCNUM)) - * ((M+GEOMD3(LOCNUM))/GEOMD2(LOCNUM))**2 IF ((DZ.GT.4.0D0) .OR. (DZ.LT.2.0D0)) GO TO 999 DZ = 0.5D0 * SQRT (DZ) DD = (M+GEOMD3(LOCNUM)) * DZ / GEOMD2(LOCNUM) IF (ABS(DD).GT.1.0D0) GO TO 999 DD = ASIN (DD) IF (ABS(COS(DD)).LT.DEPS) GO TO 999 DA = L * DZ / (2.0D0 * GEOMD1(LOCNUM) * COS(DD)) IF (ABS(DA).GT.1.0D0) GO TO 999 DA = ASIN (DA) RAT = RA0 + 2.0D0 * DA DECT = DD IF (ABS(DECT).GT.TWOPI/4.0D0) GO TO 999 IF (ABS(DA).GT.TWOPI/4.0D0) GO TO 999 GO TO 900 C Stereographic 90 CONTINUE DZ = (4.0D0 - SINS) / (4.0D0 + SINS) IF (ABS(DZ).GT.1.0D0) GO TO 999 DECT = DZ * SIN0 + M * COS0 * (1.0D0+DZ) / 2.0D0 IF (ABS(DECT).GT.1.0D0) GO TO 999 DECT = ASIN (DECT) RAT = COS(DECT) IF (ABS(RAT).LT.DEPS) GO TO 999 RAT = L * (1.0D0+DZ) / (2.0D0 * RAT) IF (ABS(RAT).GT.1.0D0) GO TO 999 RAT = ASIN (RAT) MG = 1.0D0 + SIN(DECT) * SIN0 + * COS(DECT) * COS0 * COS(RAT) IF (ABS(MG).LT.DEPS) GO TO 999 MG = 2.0D0 * (SIN(DECT) * COS0 - COS(DECT) * SIN0 * * COS(RAT)) / MG IF (ABS(MG-M).GT.DEPS) RAT = TWOPI/2.0D0 - RAT RAT = RA0 + RAT GO TO 900 C Return: in range RA 900 RAOUT = RAT DECOUT = DECT IERR = 0 IF (RAOUT-RA0.GT.TWOPI/2.0D0) RAOUT = RAOUT - TWOPI IF (RAOUT-RA0.LT.-TWOPI/2.0D0) RAOUT = RAOUT + TWOPI C 999 RETURN C----------------------------------------------------------------------- 1000 FORMAT ('NEWPOS: ILLEGAL PROJECTION TYPE',I4,' INPUT') 1005 FORMAT ('NEWPOS: ANGLE IS TOO LARGE. L,M=',2E11.3) END SUBROUTINE GALPOL (EPOCH, RAGP, DECGP, LONCP) C----------------------------------------------------------------------- C GALPOL returns the coordinates of the galactic center for the epoch C of observation. C Inputs: C EPOCH R Date of observation years (ie 1950, 2000 etc) C Output: C RAGP D Ra of north galactic pole (degrees) C DECGP D Dec of north galactic pole (degrees) C LONCP D Longitude of North Celestial pole (degrees) C----------------------------------------------------------------------- REAL EPOCH DOUBLE PRECISION RAGP, DECGP, LONCP C REAL DAY PARAMETER (DAY=1./365.25) C----------------------------------------------------------------------- C If less than a day from 1950.0 IF (ABS(EPOCH-1950.).LT.DAY) THEN RAGP = 192.25D0 DECGP = 27.40D0 LONCP = 123.00D0 C these assume 2000.0 ELSE RAGP = 192.859375D0 DECGP = 27.128643D0 LONCP = 122.932018D0 END IF C 999 RETURN END SUBROUTINE B2JPOS (DIR, RAIN, DECIN, RAOUT, DECOUT) C----------------------------------------------------------------------- C Routine to convert between B1950 positions and J2000 positions. C Positions assumed to be at standard mean epoch (i.e. 1950 or 2000). C Note: Proper motions, parallax and radial velocity assumed 0. C Taken from the Supplement to the Astronomical Almanac 1984. C Inputs: C DIR I 1 => B1950 to J2000; -1 => J2000 to B1950. C RAIN D Input Right ascension (radians). C DECIN D Input declination. C Outputs: C RAOUT D Output Right ascension C DECOUT D Output declination C----------------------------------------------------------------------- INTEGER DIR DOUBLE PRECISION RAIN, DECIN, RAOUT, DECOUT C INTEGER LOOP, LOOP2 DOUBLE PRECISION M(6,6), MP(6,6), MR0(6), MR1(6), MA(3), MR(6), * MS(6), MS1(6), X, Y, Z, R, TEMP, SUM, TWOPI, XR1, M3P(3,3) DATA MA /-1.62557D-6, -.31919D-6, -.13843D-6/ DATA M /0.9999256782D0, -.0111820611D0, -.0048579477D0, * 0.00000242395018D0, -.00000002710663D0, -.00000001177656D0, * 0.0111820610D0, 0.9999374784D0, -.0000271765D0, * 0.00000002710663D0, 0.00000242397878D0, -.00000000006587D0, * 0.0048579479D0, -.0000271474D0, 0.9999881997D0, * 0.00000001177656D0, -.00000000006582D0, 0.00000242410173D0, * -.000551D0, -.238565D0, 0.435739D0, * 0.99994704D0, -.01118251D0, -.00485767D0, * 0.238514D0, -.002667D0, -.008541D0, * 0.01118251D0, 0.99995883D0, -.00002718D0, * -.435623D0, 0.012254D0, 0.002117D0, * 0.00485767D0, -.00002714D0, 1.00000956D0/ DATA MP /0.9999256795D0, 0.0111814828D0, 0.0048590039D0, * -.00000242389840D0, -.00000002710544D0, -.00000001177742D0, * -.0111814828D0, 0.9999374849D0, -.0000271771D0, * 0.00000002710544D0, -.00000242392702D0, 0.00000000006585D0, * -.0048590040D0, -.0000271557D0, 0.9999881946D0, * 0.00000001177742D0, 0.00000000006585D0, -.00000242404995D0, * -.000551D0, 0.238509D0, -.435614D0, * 0.99990432D0, 0.01118145D0, 0.00485852D0, * -.238560D0, -.002667D0, 0.012254D0, * -.01118145D0, 0.99991613D0, -.00002717D0, * 0.435730D0, -.008541D0, 0.002117D0, * -.00485852D0, -.00002716D0, 0.99996684D0/ C M3P is the inverse of M as a C 3x3 matrix, MP is the inverse C as a true 6x6 matrix. DATA M3P /0.999925678130111D0, 0.0111820610293105D0, * 0.004857947862766843D0, -0.01118206107068339D0, * 0.999937478463956D0, -0.00002714738732342485D0, * -0.004857947737353646D0, -0.00002717648788286833D0, * 0.999988199766155D0/ C----------------------------------------------------------------------- TWOPI = 8.D0 * DATAN (1.D0) C Which direction? IF (DIR.GT.0) THEN C B1950 => J2000 MR0(1) = COS (RAIN) * COS (DECIN) MR0(2) = SIN (RAIN) * COS (DECIN) MR0(3) = SIN (DECIN) C Remove abberation "E-terms" TEMP = MR0(1)*MA(1) + MR0(2)*MA(2) + MR0(3)*MA(3) MR1(1) = MR0(1) * (1.0D0 + TEMP) - MA(1) MR1(2) = MR0(2) * (1.0D0 + TEMP) - MA(2) MR1(3) = MR0(3) * (1.0D0 + TEMP) - MA(3) C Rotate DO 110 LOOP = 1,3 SUM = 0.0D0 DO 100 LOOP2 = 1,3 SUM = SUM + M(LOOP2,LOOP) * MR1(LOOP2) 100 CONTINUE MR(LOOP) = SUM 110 CONTINUE C Convert to polar coordinates. X = MR(1) Y = MR(2) Z = MR(3) R = SQRT (X*X + Y*Y + Z*Z) DECOUT = ASIN (Z / R) RAOUT = ATAN2 (Y, X) IF (RAOUT.LT.0.0) RAOUT = RAOUT + TWOPI C J2000 to B1950 ELSE MR0(1) = COS (RAIN) * COS (DECIN) MR0(2) = SIN (RAIN) * COS (DECIN) MR0(3) = SIN (DECIN) C Rotate DO 220 LOOP = 1,3 SUM = 0.0D0 DO 210 LOOP2 = 1,3 SUM = SUM + M3P(LOOP2,LOOP) * MR0(LOOP2) 210 CONTINUE MR1(LOOP) = SUM 220 CONTINUE C Include abberation "E-terms" XR1 = 1.0D0 / SQRT (MR1(1)*MR1(1) + MR1(2)*MR1(2) + * MR1(3)*MR1(3)) MS1(1) = MR1(1) * XR1 MS1(2) = MR1(2) * XR1 MS1(3) = MR1(3) * XR1 MS(1) = MS1(1) MS(2) = MS1(2) MS(3) = MS1(3) C Iterate a few times DO 250 LOOP = 1,5 TEMP = MS(1)*MA(1) + MS(2)*MA(2) + MS(3)*MA(3) MR(1) = MS1(1) + MA(1) - (MS(1)*TEMP) MR(2) = MS1(2) + MA(2) - (MS(2)*TEMP) MR(3) = MS1(3) + MA(3) - (MS(3)*TEMP) XR1 = 1.0D0 / SQRT (MR(1)*MR(1) + MR(2)*MR(2) + * MR(3)*MR(3)) MS(1) = MR(1) * XR1 MS(2) = MR(2) * XR1 MS(3) = MR(3) * XR1 250 CONTINUE C Convert to polar coordinates. X = MR(1) Y = MR(2) Z = MR(3) R = SQRT (X*X + Y*Y + Z*Z) DECOUT = ASIN (Z / R) RAOUT = ATAN2 (Y, X) IF (RAOUT.LT.0.0) RAOUT = RAOUT + TWOPI END IF C 999 RETURN END