SUBROUTINE RESWCS C C Reset the header of an image to be TAN projection with the reference C pixel at the centre. C C The input can be any image with a COE or TAN WCS C C It is assumed that a linear WCS in TAN projection with the C reference pixel at the centre is adequate for the output. C C The information is read from the header and written back there. C C Richard Hook, ST-ECF, for the EIS project, July 1999 C C Removed an inappropriate NC check - Richard Hook, May 2000 C IMPLICIT NONE C Local variables CHARACTER*80 IMAGE INTEGER NX,NY,ISTAT,ID,N INTEGER DATTYP,NDIMS,DIMS(7) DOUBLE PRECISION CRVAL1,CRPIX1,CRVAL2,CRPIX2 DOUBLE PRECISION RAREF,DECREF,RPP DOUBLE PRECISION RAOFF,DECOFF,FLXSCL,WCS(8) DOUBLE PRECISION CP1,CP2 DOUBLE PRECISION CDELT1,CDELT2,CD11,CD12,CD21,CD22 DOUBLE PRECISION XCO(100),YCO(100) DOUBLE PRECISION CD,CX,CY CHARACTER*3 PROJ CHARACTER*8 CTYPE1,CTYPE2 INTEGER NC,MPIX LOGICAL IMOPEN C Constant DOUBLE PRECISION PIBY PARAMETER (PIBY=3.141592653589793/180.0) C Get the name of the input image CALL UCLGST('image',IMAGE,ISTAT) C Open the input image IMOPEN=.FALSE. CALL UIMOPN(IMAGE,2,ID,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Unable to open input image', : 1,0,ISTAT) GO TO 999 ELSE IMOPEN=.TRUE. ENDIF C Get the size and shape C Get the array size and shape CALL UIMGID(ID,DATTYP,NDIMS,DIMS,ISTAT) IF(NDIMS.NE.2) THEN CALL UMSPUT('! Input image is not two dimensional', : 1,0,ISTAT) GO TO 999 ELSE NX=DIMS(1) NY=DIMS(2) ENDIF C Save the WCS CALL SAVWCS(ID,ISTAT) C Read the astrometric header CALL EIGTCO(ID,RAREF,CRPIX1,DECREF,CRPIX2, : CDELT1,CDELT2,MPIX,FLXSCL, : XCO,YCO,RAOFF,DECOFF,NC,ISTAT) IF(ISTAT.NE.0) GO TO 777 C Get the projection CALL UHDGST(ID,'CTYPE1',CTYPE1,ISTAT) IF(ISTAT.NE.0) GO TO 777 CALL UHDGST(ID,'CTYPE2',CTYPE2,ISTAT) 777 CONTINUE IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to get all necessary header values', : 1,0,ISTAT) GO TO 999 ENDIF C Check we have a projection we can handle IF(CTYPE1.EQ.'RA---TAN' .AND. CTYPE2.EQ.'DEC--TAN') THEN PROJ='tan' ELSE IF(CTYPE1.EQ.'RA---COE' .AND. CTYPE2.EQ.'DEC--COE') THEN PROJ='coe' ELSE CALL UMSPUT('! Input image projection must be TAN or COE', : 1,0,ISTAT) GO TO 999 ENDIF C Calculate new CRVAL and CRPIX keywords CP1=FLOAT(NX)/2.0 CP2=FLOAT(NY)/2.0 N=1 C Now handle the COE and TAN cases separately. C In both cases exact expressions are used, not approximations. IF(PROJ.EQ.'tan') THEN C Setup the WCS 8 element array for XY2RD (angles in degrees) WCS(1)=CRPIX1 WCS(2)=RAREF/PIBY WCS(3)=CRPIX2 WCS(4)=DECREF/PIBY WCS(5)=XCO(2)/PIBY WCS(6)=YCO(2)/PIBY WCS(7)=XCO(1)/PIBY WCS(8)=YCO(1)/PIBY C xy2rd gives degrees so we have to convert back to radians CALL XY2RD(CP1,CP2,CRVAL1,CRVAL2,WCS) CD=COS(CRVAL2*PIBY) CALL XY2RD(CP1+1.0D0,CP2,CX,CY,WCS) CD11=(CX-CRVAL1)*CD*PIBY CD21=(CY-CRVAL2)*PIBY CALL XY2RD(CP1,CP2+1.0D0,CX,CY,WCS) CD12=(CX-CRVAL1)*CD*PIBY CD22=(CY-CRVAL2)*PIBY CRVAL1=CRVAL1*PIBY CRVAL2=CRVAL2*PIBY C COE case - use packaged routine ELSE C Deduce the scale RPP=ABS(XCO(2)) C New origin on sky CALL COE2RD(CP1,CP2,N,RAREF,DECREF,CRPIX1,CRPIX2, : RPP,CRVAL1,CRVAL2) CD=COS(CRVAL2) CALL COE2RD(CP1+1.0D0,CP2,N,RAREF,DECREF,CRPIX1,CRPIX2, : RPP,CX,CY) CD11=(CX-CRVAL1)*CD CD21=CY-CRVAL2 CALL COE2RD(CP1,CP2+1.0D0,N,RAREF,DECREF,CRPIX1,CRPIX2, : RPP,CX,CY) CD12=(CX-CRVAL1)*CD CD22=CY-CRVAL2 ENDIF C Fill the WCS array (converting to degrees from radians) WCS(1)=CP1 WCS(2)=CRVAL1/PIBY WCS(3)=CP2 WCS(4)=CRVAL2/PIBY WCS(5)=CD11/PIBY WCS(6)=CD21/PIBY WCS(7)=CD12/PIBY WCS(8)=CD22/PIBY CTYPE1='RA---TAN' CTYPE2='DEC--TAN' C Write the new WCS values and C the TAN projection information to the header CALL UPWCS(ID,WCS,CTYPE1,CTYPE2,0.0) 999 CONTINUE C Close the image IF(IMOPEN) CALL UIMCLO(ID,ISTAT) RETURN END SUBROUTINE UPWCS(ID,WCS,CTYPE1,CTYPE2,DECRED) C C Update the WCS header items using the 8 values C supplied in WCS and the values of CTYPE1 and CTYPE2. C C This is the EIS version which doesn't update the CDELT C C The image is assumed to be available with descriptor ID C and to already have all the parameters C C Modified so that the CD matrix is added - not modified - C Richard Hook, ST-ECF, July 1997 C C Modified so that full WCS is supplied and the items are C deleted before being added, Dec 97 C C Specific version for RESWCS C IMPLICIT NONE INTEGER ID,ISTAT CHARACTER*8 CTYPE1,CTYPE2 DOUBLE PRECISION WCS(8),DECRED C Write a full WCS for the EIS COE projection CALL UHDDSP(ID,'CTYPE1',ISTAT) CALL UHDAST(ID,'CTYPE1',CTYPE1, : 'WCS projection type',0,ISTAT) CALL UHDDSP(ID,'CTYPE2',ISTAT) CALL UHDAST(ID,'CTYPE2',CTYPE2, : 'WCS projection type',0,ISTAT) CALL UHDDSP(ID,'CRPIX1',ISTAT) CALL UHDASD(ID,'CRPIX1',WCS(1), : 'Reference pixel (X)',0,ISTAT) CALL UHDDSP(ID,'CRVAL1',ISTAT) CALL UHDASD(ID,'CRVAL1',WCS(2), : 'Reference point on sky (RA,degrees)',0,ISTAT) CALL UHDDSP(ID,'CRPIX2',ISTAT) CALL UHDASD(ID,'CRPIX2',WCS(3), : 'Reference pixel (Y)',0,ISTAT) CALL UHDDSP(ID,'CRVAL2',ISTAT) CALL UHDASD(ID,'CRVAL2',WCS(4), : 'Reference point on sky (dec,degrees)',0,ISTAT) CALL UHDDSP(ID,'CD1_1',ISTAT) CALL UHDASD(ID,'CD1_1',WCS(5), : 'CD1_1 matrix element',0,ISTAT) CALL UHDDSP(ID,'CD2_1',ISTAT) CALL UHDASD(ID,'CD2_1',WCS(6), : 'CD2_1 matrix element',0,ISTAT) CALL UHDDSP(ID,'CD1_2',ISTAT) CALL UHDASD(ID,'CD1_2',WCS(7), : 'CD1_2 matrix element',0,ISTAT) CALL UHDDSP(ID,'CD2_2',ISTAT) CALL UHDASD(ID,'CD2_2',WCS(8), : 'CD2_2 matrix element',0,ISTAT) IF(CTYPE1.EQ.'RA---COE') THEN CALL UHDDSP(ID,'PROJP1',ISTAT) CALL UHDASD(ID,'PROJP1',DECRED, : 'COE projection reference declination (degrees)',0,ISTAT) CALL UHDDSP(ID,'PROJP2',ISTAT) CALL UHDASD(ID,'PROJP2',0.0D0, : 'COE projection - separation of reference decs', : 0,ISTAT) ENDIF RETURN END SUBROUTINE SAVWCS(ID,ISTAT) C C Save the WCS before altering it. The 10 parameters C CTYPE1/2, CRVAL1/2, CRPIX1/2, CD1/2_1/2 are copied C to the same names with 'O' at the start. C C It is assumed that the header has been opened. C IMPLICIT NONE INTEGER ID,ISTAT DOUBLE PRECISION CD11,CD12,CD21,CD22 DOUBLE PRECISION CRVAL1,CRVAL2,CRPIX1,CRPIX2 CHARACTER*8 CTYPE1,CTYPE2 CALL UHDGSD(ID,'CRVAL1',CRVAL1,ISTAT) IF(ISTAT.EQ.0) : CALL UHDASD(ID,'OCRVAL1',CRVAL1,'Original CRVAL1',0,ISTAT) CALL UHDGSD(ID,'CRVAL2',CRVAL2,ISTAT) IF(ISTAT.EQ.0) : CALL UHDASD(ID,'OCRVAL2',CRVAL2,'Original CRVAL2',0,ISTAT) CALL UHDGSD(ID,'CRPIX1',CRPIX1,ISTAT) IF(ISTAT.EQ.0) : CALL UHDASD(ID,'OCRPIX1',CRPIX1,'Original CRPIX1',0,ISTAT) CALL UHDGSD(ID,'CRPIX2',CRPIX2,ISTAT) IF(ISTAT.EQ.0) : CALL UHDASD(ID,'OCRPIX2',CRPIX2,'Original CRPIX2',0,ISTAT) CALL UHDGSD(ID,'CD1_1',CD11,ISTAT) IF(ISTAT.EQ.0) : CALL UHDASD(ID,'OCD1_1',CD11,'Original CD1_1',0,ISTAT) CALL UHDGSD(ID,'CD1_2',CD12,ISTAT) IF(ISTAT.EQ.0) : CALL UHDASD(ID,'OCD1_2',CD12,'Original CD1_2',0,ISTAT) CALL UHDGSD(ID,'CD2_1',CD21,ISTAT) IF(ISTAT.EQ.0) : CALL UHDASD(ID,'OCD2_1',CD21,'Original CD2_1',0,ISTAT) CALL UHDGSD(ID,'CD2_2',CD22,ISTAT) IF(ISTAT.EQ.0) : CALL UHDASD(ID,'OCD2_2',CD22,'Original CD2_2',0,ISTAT) CALL UHDGST(ID,'CTYPE1',CTYPE1,ISTAT) IF(ISTAT.EQ.0) : CALL UHDAST(ID,'OCTYPE1',CTYPE1,'Original CTYPE1',0,ISTAT) CALL UHDGST(ID,'CTYPE2',CTYPE2,ISTAT) IF(ISTAT.EQ.0) : CALL UHDAST(ID,'OCTYPE2',CTYPE2,'Original CTYPE2',0,ISTAT) RETURN END