SUBROUTINE EIEVAL(XIN,YIN,N,
     :   CRVAL1,CRPIX1,CRVAL2,CRPIX2,
     :   CDELT1,CDELT2,MPIX,
     :   XCO,YCO,RAOFF,DECOFF,NC,
     :   RAREF,DECREF,XPOFF,YPOFF,RPP,PROJ,
     :   XOUT,YOUT,ISTAT)
C
C EIEVAL
C
C Convert X,Y positions from an input EIS frame into
C X,Y pixel positions on an output "tile" of the super-image.
C
C Supplied:
C
C XIN,YIN - double precision arrays with N elements - input frame
C           pixel positions.
C
C N - integer - the number of elements in the input arrays
C
C CRVAL1,CRPIX1,CRVAL2,CRPIX2,CDELT1,CDELT2 - double precision
C     scalars - the WCS origin and scale for the input image.
C
C MPIX - integer - a scale used by Erik Deul's pipeline. For the
C        new style (PV) coefficients this will be supplied as 0
C        and is used to determine which representation to use.
C        It can also be used to signal the use of test transformations.
C        This "fudge" was added in November 1999 for galactic coords
C        testing (MPIX=-1). In addition MPIX=-2 means that the input
C        and output grids are the same except for an integer pixel
C        shift. This is useful for fast "stacking".
C        Finally, for EIS Drizzle V0.83 the additional option MPIX=-3
C        was added for input COE images.
C
C XCO,YCO - double precision - distortion coefficients from LDAC
C           pipeline.
C
C RAOFF,DECOFF - double precisions - offsets from LDAC
C                pipeline. If the astrometric solution
C                is in the new (PV) form these values are ignored.
C 
C NC - integer - the number of coefficients in X and Y (dimensions
C                of XCO and YCO).
C
C RAREF,DECREF - double precision - the reference position of the
C                super image, ie, the reference point of the conical
C                equal area projection. 
C
C XPOFF,YPOFF - double precision - the offset, in pixels for this tile.
C
C RPP - radians/pixel - the scale of the output image.
C
C PROJ - character*3 - the projection (currently COE or TAN)
C
C Returned:
C
C XOUT,YOUT - double precision arrays, the pixel positions in the output
C             image.
C
C ISTAT - a status flag (0=good)
C
C Richard Hook, April 1997
C
C Modified to support the PV style coefficients, January 1999.
C Added the support for multiple projections using PROJ, January 1999.
C Negative MPIX "fudge" added - November 1999
C MPIX=-2 pixel-shift option added - December 1999
C MPIX=-3 COE input option added - March 2000
C
      IMPLICIT NONE

      INTEGER N
      DOUBLE PRECISION XIN(*),YIN(*)
      DOUBLE PRECISION CRVAL1,CRPIX1,CRVAL2,CRPIX2
      DOUBLE PRECISION CDELT1,CDELT2
      INTEGER MPIX
      DOUBLE PRECISION XCO(*),YCO(*),RAOFF,DECOFF
      INTEGER NC
      DOUBLE PRECISION RAREF,DECREF,XPOFF,YPOFF,RPP
      DOUBLE PRECISION XOUT(*),YOUT(*)
      CHARACTER*3 PROJ
      INTEGER ISTAT

C Local variables
      INTEGER I
      DOUBLE PRECISION S,TD,CD,X2,Y2,XD,YD,CT
      DOUBLE PRECISION EVEIP5,EVEIP9,EVEI14,EVEI20
      DOUBLE PRECISION PIBY

      PIBY=3.141592653589/180.0

C Calculate some constants to make life easier and faster
      IF(MPIX.GT.0) S=2.0/DBLE(MPIX)

      TD=TAN(CRVAL2)
      CD=COS(CRVAL2)

C First the MPIX=-2 case which is the simple one of a pixel
C offset in both X and Y
      IF(MPIX.EQ.-2) THEN

         DO I=1,N
            XOUT(I)=XIN(I)-CRPIX1-XPOFF
            YOUT(I)=YIN(I)-CRPIX2-YPOFF
         ENDDO

         GO TO 999
   
C Fudge - if MPIX=-1 we use this as a way of trying a test transformation
C         in this case galactic-cartesian
      ELSE IF(MPIX.EQ.-1) THEN

C First go from pixels to galactic coordinates
         DO I=1,N
            XOUT(I)=CRVAL1+(XIN(I)-CRPIX1)*CDELT1
            YOUT(I)=CRVAL2+(YIN(I)-CRPIX2)*CDELT2

C and then from galactic to equatorial
C This was removed as we may not have SLALIB available and it is
C unsupported anyway
CCCC            CALL SLGAEQ(XOUT(I),YOUT(I),XOUT(I),YOUT(I))
         ENDDO

C now jump to converting to output pixels (ugly)
         GO TO 33

C Case of COE input 
      ELSE IF(MPIX.EQ.-3) THEN

C Convert from input pixels directly to RA,Dec 
C note that we use CDELT2 as the scale in radians/pixel
         CALL COE2RD(XIN,YIN,N,CRVAL1,CRVAL2,CRPIX1,CRPIX2,
     :               CDELT2,XOUT,YOUT)

C Now skip to the final stage of converting to output position
         GO TO 33
      ENDIF

C Step 1 - convert from input pixels to tangent plane positions
C          in radians
C
C First handle the new-style PV coefficients
      IF(MPIX.EQ.0) THEN
         CALL EVALPV(XIN,YIN,N,CRPIX1,CRPIX2,XCO,YCO,NC,XOUT,YOUT)
      ELSE IF(MPIX.GT.0) THEN

C This section is for the old-style polynomial coefficients
        DO I=1,N

C Scale and offset
         XD=(XIN(I)-CRPIX1)*S
         YD=(YIN(I)-CRPIX2)*S

C Now the polynomial distortion
         IF(NC.EQ.2) THEN

C Linear case
C This was changed, RNH, 14th May 1998 to reverse the order
C But changed back again so that it worked correctly with linear
C LDAC coefficients as well as CD matrix
            XOUT(I)=RAOFF+XCO(1)*YD+XCO(2)*XD
            YOUT(I)=DECOFF+YCO(1)*YD+YCO(2)*XD
         ELSE IF(NC.EQ.5) THEN

C 2nd order case
            XOUT(I)=RAOFF+EVEIP5(XD,YD,XCO)
            YOUT(I)=DECOFF+EVEIP5(XD,YD,YCO)
         ELSE IF(NC.EQ.9) THEN

C 3rd order case
            XOUT(I)=RAOFF+EVEIP9(XD,YD,XCO)
            YOUT(I)=DECOFF+EVEIP9(XD,YD,YCO)

         ELSE IF(NC.EQ.14) THEN

C 4th order case
            XOUT(I)=RAOFF+EVEI14(XD,YD,XCO)
            YOUT(I)=DECOFF+EVEI14(XD,YD,YCO)

         ELSE IF(NC.EQ.20) THEN

C 5th order case
            XOUT(I)=RAOFF+EVEI20(XD,YD,XCO)
            YOUT(I)=DECOFF+EVEI20(XD,YD,YCO)
         ELSE
            ISTAT=1
            RETURN
         ENDIF
        ENDDO
      ENDIF

C Step 2 - convert from tangent plane to RA,Dec on the sky
      DO I=1,N
         X2=XOUT(I)**2/2.0
         Y2=YOUT(I)**2/2.0
   
         CT=CD*(1.0-Y2-TD*YOUT(I))
         XOUT(I)=CRVAL1+XOUT(I)/CT
         YOUT(I)=CRVAL2+YOUT(I)-TD*X2
      ENDDO

 33   CONTINUE

C Step 3 - convert from RA,Dec to pixels on the super-image
      IF(PROJ.EQ.'COE') THEN
         CALL COEPRO(XOUT,YOUT,N,RAREF,DECREF,RPP,XOUT,YOUT)
      ELSE IF(PROJ.EQ.'TAN') THEN
         CALL TANPRO(XOUT,YOUT,N,RAREF,DECREF,RPP,XOUT,YOUT)
      ELSE
         CALL UMSPUT('! Invalid projection specified',
     :               1,0,ISTAT)
         ISTAT=1
         GO TO 999
      ENDIF

C Step 4 - convert to pixels on the output tile
      DO I=1,N
         XOUT(I)=XOUT(I)-XPOFF
         YOUT(I)=YOUT(I)-YPOFF
      ENDDO

999   CONTINUE

      RETURN
      END
      SUBROUTINE EIGTCO(ID,CRVAL1,CRPIX1,CRVAL2,CRPIX2,
     :                  CDELT1,CDELT2,MPIX,FLXSCL,
     :                  XCO,YCO,RAOFF,DECOFF,NC,ISTAT)
C
C Read the coefficients of an EIS image from the LDAC
C pipeline needed for the geometric distortion.
C
C Note - all angles are converted to radians before being
C        returned.
C
C Richard Hook, ST-ECF, April 1997
C
C Revised to get sensible linear coefficients from the CD
C matrix if the LDAC keywords are missing. Richard Hook, November 1997
C
C Revised to handle the new-style PV coefficients, January 1999
C
C Added "fudge" to allow test projections, November 1999
C
C Modified to set sensible CDELT values even if they aren't in the
C header, February 2000
C
C Added MPIX=-3 case for COE input, March 2000
C
      IMPLICIT NONE

      INTEGER ID
      DOUBLE PRECISION CRVAL1,CRPIX1,CRVAL2,CRPIX2
      DOUBLE PRECISION CDELT1,CDELT2,CD11,CD12,CD21,CD22
      DOUBLE PRECISION XCO(*),YCO(*),RAOFF,DECOFF,FLXSCL
      INTEGER MPIX,NC

C Local variables
      INTEGER ISTAT,N
      DOUBLE PRECISION PIBY
      CHARACTER*8 CTYPE1,CTYPE2
      CHARACTER*7 CONAME

      PIBY=3.141592653589/180.0

C First get the basic WCS entries for the reference point
C and scales
      CALL UHDGSD(ID,'CRVAL1',CRVAL1,ISTAT)
      IF(ISTAT.NE.0) RETURN
      CRVAL1=CRVAL1*PIBY

      CALL UHDGSD(ID,'CRPIX1',CRPIX1,ISTAT)
      IF(ISTAT.NE.0) RETURN

      CALL UHDGSD(ID,'CRVAL2',CRVAL2,ISTAT)
      IF(ISTAT.NE.0) RETURN
      CRVAL2=CRVAL2*PIBY

      CALL UHDGSD(ID,'CRPIX2',CRPIX2,ISTAT)
      IF(ISTAT.NE.0) RETURN

C The CDELTs are not so important so we don't complain if they
C aren't there
      CALL UHDGSD(ID,'CDELT1',CDELT1,ISTAT)
      IF(ISTAT.NE.0) THEN
         CDELT1=0.0
      ELSE
         CDELT1=CDELT1*PIBY
      ENDIF

      CALL UHDGSD(ID,'CDELT2',CDELT2,ISTAT)
      IF(ISTAT.NE.0) THEN
         CDELT2=0.0
      ELSE
         CDELT2=CDELT2*PIBY
      ENDIF

C If there are no CDELT keywords get the CD matrix and
C use those values
      IF(CDELT1.EQ.0.0 .OR. CDELT2.EQ.0.0) THEN
         CALL UHDGSD(ID,'CD1_1',CD11,ISTAT)
         CALL UHDGSD(ID,'CD1_2',CD12,ISTAT)
         CALL UHDGSD(ID,'CD2_1',CD21,ISTAT)
         CALL UHDGSD(ID,'CD2_2',CD22,ISTAT)
         CDELT1=-SQRT(CD11**2+CD21**2)*PIBY
         CDELT2=SQRT(CD12**2+CD22**2)*PIBY
      ENDIF

C First we check to see whether we have an MPIX keyword. If we
C do we assume the "old-style" of coefficients and if not the
C new PV ones.
      CALL UHDGSI(ID,'MPIX',MPIX,ISTAT)
      IF(ISTAT.NE.0) MPIX=0

C Now for the "fudge" (added November 1999), check the CTYPE
C values
      CALL UHDGST(ID,'CTYPE1',CTYPE1,ISTAT)
      CALL UHDGST(ID,'CTYPE2',CTYPE2,ISTAT)

      IF(CTYPE1.EQ.'GLON-CAR' .AND. CTYPE2.EQ.'GLAT-CAR') THEN
         CALL UMSPUT('#+ Note - galactic cartesian coordinates',
     :               1,0,ISTAT)
         MPIX=-1
      ELSE IF(CTYPE1.EQ.'RA---COE' .AND. CTYPE2.EQ.'DEC--COE') THEN
         CALL UMSPUT('#+ Note - input image projection is COE',
     :               1,0,ISTAT)
         MPIX=-3
      ELSE 
        IF(CTYPE1.NE.'RA---TAN' .OR. CTYPE2.NE.'DEC--TAN') 
     :    CALL UMSPUT('#+ Input image WCS is '//CTYPE1//'/'//CTYPE2,
     :               1,0,ISTAT)
      ENDIF

C Now get the old-style polynomial coefficients and count them
      IF(MPIX.GT.0) THEN
        N=0
        DO WHILE(.TRUE.)
         IF(N.LT.10) THEN
            WRITE(CONAME,'(''Xpol_'',I1,'' '')') N
         ELSE
            WRITE(CONAME,'(''Xpol_'',I2)') N
         ENDIF

         CALL UHDGSD(ID,CONAME,XCO(N+1),ISTAT)
         IF(ISTAT.NE.0) GO TO 111

         IF(N.LT.10) THEN
            WRITE(CONAME,'(''Ypol_'',I1,'' '')') N
         ELSE
            WRITE(CONAME,'(''Ypol_'',I2)') N
         ENDIF

         CALL UHDGSD(ID,CONAME,YCO(N+1),ISTAT)
         IF(ISTAT.NE.0) GO TO 111

         N=N+1
        ENDDO

111     CONTINUE
        ISTAT=0
        NC=N
      ELSE
        
C PV style coefficients
        N=1
        DO WHILE(.TRUE.)
         IF(N.LT.10) THEN
            WRITE(CONAME,'(''PV'',I1,''_1  '')') N
         ELSE
            WRITE(CONAME,'(''PV'',I2,''_1 '')') N
         ENDIF

         CALL UHDGSD(ID,CONAME,XCO(N),ISTAT)
         IF(ISTAT.NE.0) GO TO 222

         IF(N.LT.10) THEN
            WRITE(CONAME,'(''PV'',I1,''_2  '')') N
         ELSE
            WRITE(CONAME,'(''PV'',I2,''_2 '')') N
         ENDIF

         CALL UHDGSD(ID,CONAME,YCO(N),ISTAT)
         IF(ISTAT.NE.0) GO TO 222

         N=N+1
        ENDDO

222     CONTINUE
        ISTAT=0
        NC=N-1
      ENDIF

C If there are no coefficients we may have a perfectly valid
C WCS which can be used for a linear transform
      IF(NC.EQ.0 .AND. (MPIX.GE.0 .OR. MPIX.LE.-1)) THEN

C Get the CD matrix
         CALL UHDGSD(ID,'CD1_1',CD11,ISTAT)
         IF(ISTAT.NE.0) THEN
            CALL UMSPUT('! Unable to find CD1_1 keyword',1,0,
     :                  ISTAT)
            ISTAT=1
            RETURN
         ENDIF

         CALL UHDGSD(ID,'CD1_2',CD12,ISTAT)
         IF(ISTAT.NE.0) THEN
            CALL UMSPUT('! Unable to find CD1_2 keyword',1,0,
     :                  ISTAT)
            CALL UMSPUT('! using default of 0.0',1,0,ISTAT)
            CD12=0.0
            RETURN
         ENDIF

         CALL UHDGSD(ID,'CD2_1',CD21,ISTAT)
         IF(ISTAT.NE.0) THEN
            CALL UMSPUT('! Unable to find CD2_1 keyword',1,0,
     :                  ISTAT)
            CALL UMSPUT('! using default of 0.0',1,0,ISTAT)
            CD21=0.0
            RETURN
         ENDIF

         CALL UHDGSD(ID,'CD2_2',CD22,ISTAT)
         IF(ISTAT.NE.0) THEN
            CALL UMSPUT('! Unable to find CD2_2 keyword',1,0,
     :                  ISTAT)
            ISTAT=1
            RETURN
         ENDIF

C Set linear coefficients (converting degrees to radians)
C Note the minus signs which were inserted for EIS drizzle V0.50
C ...and subsequently removed again.
         XCO(1)=CD12*PIBY
         XCO(2)=CD11*PIBY
         YCO(1)=CD22*PIBY
         YCO(2)=CD21*PIBY

C Fake the other numbers
         MPIX=2
         RAOFF=0.0D0
         DECOFF=0.0D0
         NC=2

      ELSE

C Now we get the RA and Dec offsets and the pixel range
       IF(MPIX.GT.0) THEN
         CALL UHDGSD(ID,'RA_OFF',RAOFF,ISTAT)
         IF(ISTAT.NE.0) RETURN

         CALL UHDGSD(ID,'DEC_OFF',DECOFF,ISTAT)
         IF(ISTAT.NE.0) RETURN
       ENDIF
      ENDIF

      CALL UHDGSD(ID,'FLXSCALE',FLXSCL,ISTAT)
      IF(ISTAT.NE.0) THEN
         FLXSCL=1.0
         ISTAT=0
         RETURN
      ENDIF

      RETURN
      END
      SUBROUTINE TANPRO(RA,DEC,N,RAREF,DECREF,RPP,X,Y)
C
C Project a position from RA,Dec expressed in radians
C to X,Y in pixels on an output master image using the 
C TAN projection.
C
C An approximate polynomial representation of the exact
C transformation is used which is accurate to better than
C 0.1 arcsecs over a 1 degree wide region.
C
C Note that this projection approximation is only suitable
C for fairly small "deep" fields.
C
C Supplied:
C
C RA - double precision array with N elements, Right Ascension
C DEC - double precision array with N elements, Declination
C N - integer, number of elements
C RAREF - double precision, reference right ascension,
C DECREF - double precision, reference declination
C RPP - double precision, scale, radians/pixel
C
C Returned:
C
C X - double precision array with N elements, X pixel offset
C Y - double precision array with N elements, Y pixel offset
C 
C Note - all angles are in RADIANS
C
C Richard Hook, January 1999, EIS project (created from COE version)
C
C Revised, 25th January 1999 - added extra terms (cubic)
C
      IMPLICIT NONE

C Inputs
      INTEGER N
      DOUBLE PRECISION RA(N),DEC(N),RAREF,DECREF,RPP

C Outputs
      DOUBLE PRECISION X(N),Y(N)

C Local variables and constants
      INTEGER I
      DOUBLE PRECISION SINDR,COSDR,SINCOS,DRA,DDEC,SIN2,COS2
      DOUBLE PRECISION DD2,DR2,BA

      SINDR=SIN(DECREF)
      SIN2=SINDR**2
      COSDR=COS(DECREF)
      COS2=COSDR**2
      SINCOS=SINDR*COSDR

C Loop over the input arrays
      DO I=1,N
         DRA=RAREF-RA(I)
         DDEC=DEC(I)-DECREF
         DD2=DDEC**2
         DR2=DRA**2
         BA=1.0-DD2/2.0-DR2/2*(COS2-DDEC*SINCOS)
         X(I)=(COSDR*(1.0-DD2/2.0)-DDEC*(1.0-DD2/6.0)*SINDR)*
     :         (1.0-DR2/6.0)*DRA/RPP/BA
         Y(I)=(DDEC+SINCOS*DR2/2.0-
     :         DDEC*DD2/6.0-DDEC*DR2/2.0*SIN2)/RPP/BA
      ENDDO

      RETURN
      END
      SUBROUTINE COEPRO(RA,DEC,N,RAREF,DECREF,RPP,X,Y)
C
C Project a position from RA,Dec expressed in radians
C to X,Y in pixels on an output master image using the 
C Conical Equal-area projection.
C
C An approximate polynomial representation of the exact
C transformation is used which is accurate to better than
C 0.01 arcsecs over a 4 degree wide region.
C
C Supplied:
C
C RA - double precision array with N elements, Right Ascension
C DEC - double precision array with N elements, Declination
C N - integer, number of elements
C RAREF - double precision, reference right ascension,
C DECREF - double precision, reference declination
C RPP - double precision, scale, radians/pixel
C
C Returned:
C
C X - double precision array with N elements, X pixel offset
C Y - double precision array with N elements, Y pixel offset
C 
C Note - all angles are in RADIANS
C
C Richard Hook, March 1997, EIS project
C RNH - change from degrees to radians, April 1997
C
      IMPLICIT NONE

C Inputs
      INTEGER N
      DOUBLE PRECISION RA(N),DEC(N),RAREF,DECREF,RPP

C Outputs
      DOUBLE PRECISION X(N),Y(N)

C Local variables and constants
      INTEGER I
      DOUBLE PRECISION RTHEA,GAMBY2,A,B,C
      DOUBLE PRECISION PHI,AN,RAP,V,TTHEA
      DOUBLE PRECISION CO1,CO2,CO3,CO4,ANSQ

      GAMBY2=SIN(DECREF)
      TTHEA=TAN(DECREF)
      RTHEA=1.0/(TTHEA)
      A=-2.0*TTHEA
      B=TTHEA**2
      C=TTHEA/3.0

      CO1=A/2.0
      CO2=-1.0/8.0*A**2+B/2.0
      CO3=-1.0/4.0*A*B+1.0/16.0*A**3+C/2
      CO4=-1.0/8.0*B**2-1.0/4.0*A*C+3.0/16.0*B*A**2-5.0/128.0*A**4

C Loop over the input arrays
      DO I=1,N
         PHI=RAREF-RA(I)
         AN=PHI*GAMBY2
         V=DEC(I)-DECREF
         RAP=RTHEA*(1.0+V*(CO1+V*(CO2+V*(CO3+V*CO4))))
         ANSQ=AN*AN
         X(I)=RAP*AN*(1.0-ANSQ/6.0)/RPP
         Y(I)=(RTHEA-RAP*(1.0-ANSQ/2.0))/RPP
      ENDDO

      RETURN
      END
      SUBROUTINE CORNER(NX,NY,
     :   CRVAL1,CRPIX1,CRVAL2,CRPIX2,
     :   CDELT1,CDELT2,MPIX,
     :   XCO,YCO,RAOFF,DECOFF,NC,
     :   RAREF,DECREF,XPOFF,YPOFF,RPP,PROJ,
     :   XMIN,XMAX,YMIN,YMAX,ISTAT)
C
C CORNER
C 
C Return an estimate of the section of an output image
C which will be covered by the transformation of an input
C image using the COE projection.
C
C The result is conservative and has a 5% margin of error
C added in both directions.
C
C The results are returned as pixel positions (integers)
C on the selected image tile.
C
C Richard Hook, ST-ECF, April 1997
C Added PROJ, January 1999
C Modified boundary margin to be 20 pixels rather than 5%
C  March 2000
C
      IMPLICIT NONE

      INTEGER NX,NY
      DOUBLE PRECISION CRVAL1,CRPIX1,CRVAL2,CRPIX2
      DOUBLE PRECISION CDELT1,CDELT2
      INTEGER MPIX
      DOUBLE PRECISION XCO(*),YCO(*),RAOFF,DECOFF
      INTEGER NC
      DOUBLE PRECISION RAREF,DECREF,XPOFF,YPOFF,RPP
      INTEGER XMIN,XMAX,YMIN,YMAX
      CHARACTER*3 PROJ
      INTEGER ISTAT

C Local variables and constants
      DOUBLE PRECISION XIN(4),YIN(4),XOUT(4),YOUT(4)
      DOUBLE PRECISION XR,YR,XMA,XMI,YMA,YMI

C First fill out the corner positions of the input image
      XIN(1)=1.0
      XIN(2)=DBLE(NX)
      XIN(3)=DBLE(NX)
      XIN(4)=1.0
      YIN(1)=1.0
      YIN(2)=1.0
      YIN(3)=DBLE(NY)
      YIN(4)=DBLE(NY)

C Find where these go on the output
      CALL EIEVAL(XIN,YIN,4,
     :   CRVAL1,CRPIX1,CRVAL2,CRPIX2,
     :   CDELT1,CDELT2,MPIX,
     :   XCO,YCO,RAOFF,DECOFF,NC,
     :   RAREF,DECREF,XPOFF,YPOFF,RPP,PROJ,
     :   XOUT,YOUT,ISTAT)
      IF(ISTAT.NE.0) RETURN

C Calculate the extreme values
      XMA=MAX(XOUT(1),XOUT(2),XOUT(3),XOUT(4))
      YMA=MAX(YOUT(1),YOUT(2),YOUT(3),YOUT(4))
      XMI=MIN(XOUT(1),XOUT(2),XOUT(3),XOUT(4))
      YMI=MIN(YOUT(1),YOUT(2),YOUT(3),YOUT(4))

      XR=XMA-XMI
      YR=YMA-YMI 

C Calculate limits allowing for a 20 pixel margin of error
C Note - this was formerly 5% which seems too generous
      XMAX=NINT(XMA+20)
      XMIN=NINT(XMI-20)
      YMAX=NINT(YMA+20)
      YMIN=NINT(YMI-20)

      ISTAT=0
      RETURN
      END
      DOUBLE PRECISION FUNCTION EVEIP5(X,Y,CO)
C
C Evaluate a polynomial in the form used by the LDAC
C pipeline. This version for 5 coefficients only.
C
C Richard Hook, ST-ECF, April 1997
C
      IMPLICIT NONE

      DOUBLE PRECISION X,Y,CO(*)

      EVEIP5=Y*(CO(1)+Y*CO(2))+X*(CO(3)+Y*CO(4)+X*CO(5))

      RETURN
      END
      DOUBLE PRECISION FUNCTION EVEIP9(X,Y,CO)
C
C Evaluate a polynomial in the form used by the LDAC
C pipeline. This version for 9 coefficients only.
C
C Richard Hook, ST-ECF, April 1997
C
      IMPLICIT NONE

      DOUBLE PRECISION X,Y,CO(*)

      EVEIP9=Y*(CO(1)+Y*(CO(2)+Y*CO(3)+X*CO(6)))
     :      +X*(CO(4)+Y*CO(5)+X*(CO(7)+Y*CO(8)+X*CO(9)))

      RETURN
      END

      DOUBLE PRECISION FUNCTION EVEI14(X,Y,CO)
C
C Evaluate a polynomial in the form used by the LDAC
C pipeline. This version for 14 coefficients only.
C
C Richard Hook, ST-ECF, April 1997
C
      IMPLICIT NONE

      DOUBLE PRECISION X,Y,CO(*)

      EVEI14=Y*(CO(1)+Y*(CO(2)+Y*(CO(3)+X*CO(8)+Y*CO(4))+
     :       X*CO(7)))
     :      +X*(CO(5)+Y*CO(6)+X*(CO(9)+Y*(CO(10)+Y*CO(11))+
     :       X*(CO(12)+Y*CO(13)+X*CO(14))))

      RETURN
      END

      DOUBLE PRECISION FUNCTION EVEI20(X,Y,CO)
C
C Evaluate a polynomial in the form used by the LDAC
C pipeline. This version for 20 coefficients only.
C
C Richard Hook, ST-ECF, July 1997
C
      IMPLICIT NONE

      DOUBLE PRECISION X,Y,CO(*)

      EVEI20=Y*(CO(1)+Y*(CO(2)+X*CO(8)+
     :      Y*(CO(3)+X*CO(9)+Y*(CO(4)+X*CO(10)+Y*CO(5)))))+
     :      X*(CO(6)+Y*CO(7)+X*(CO(11)+
     :      Y*(CO(12)+Y*(CO(13)+Y*CO(14)))+
     :      X*(CO(15)+Y*(CO(16)+Y*CO(17))+X*(CO(18)+
     :      Y*CO(19)+X*CO(20)))))

      RETURN
      END

      SUBROUTINE EVALPV(XIN,YIN,N,XOR,YOR,XCO,YCO,NC,XOUT,YOUT)
C
C Evaluate a polynomial stored in the PV style for an
C array of X,Y points.
C
C Supplied:
C
C XIN,YIN - double precision arrays of input X,Y pixel positions
C N - integer, the number of items in these arrays
C XOR,YOR - double precision, the origin (reference pixel position)
C XCO,YCO - double precision arrays with NC elements, the coefficients
C NC - the number of polynomial coefficients (up to 21 for quintic)
C
C Returned:
C
C XOUT,YOUT - double precision arrays of the output tangent plane
C             positions in radians.
C
C Richard Hook, 14th January 1999
C
      IMPLICIT NONE

      INTEGER N,NC
      DOUBLE PRECISION XIN(N),YIN(N),XOR,YOR,XCO(NC),YCO(NC)
      DOUBLE PRECISION XOUT(N),YOUT(N),X,Y

      INTEGER I

C Separate loops for different numbers of coefficients, for efficiency
      IF(NC.EQ.3) THEN
         DO I=1,N
            X=XIN(I)-XOR
            Y=YIN(I)-YOR
            XOUT(I)=XCO(1)+XCO(2)*X+XCO(3)*Y
            YOUT(I)=YCO(1)+YCO(2)*X+YCO(3)*Y
         ENDDO
      ELSE IF(NC.EQ.6) THEN
         DO I=1,N
            X=XIN(I)-XOR
            Y=YIN(I)-YOR
            XOUT(I)=XCO(1)+XCO(2)*X+XCO(3)*Y+
     :              XCO(4)*X*X+XCO(5)*X*Y+XCO(6)*Y*Y
            YOUT(I)=YCO(1)+YCO(2)*X+YCO(3)*Y+
     :              YCO(4)*X*X+YCO(5)*X*Y+YCO(6)*Y*Y
         ENDDO
      ELSE IF(NC.EQ.10) THEN
         DO I=1,N
            X=XIN(I)-XOR
            Y=YIN(I)-YOR
            XOUT(I)=XCO(1)+XCO(2)*X+XCO(3)*Y+
     :              XCO(4)*X*X+XCO(5)*X*Y+XCO(6)*Y*Y+
     :              XCO(7)*X*X*X+XCO(8)*X*X*Y+
     :              XCO(9)*X*Y*Y+XCO(10)*Y*Y*Y
            YOUT(I)=YCO(1)+YCO(2)*X+YCO(3)*Y+
     :              YCO(4)*X*X+YCO(5)*X*Y+YCO(6)*Y*Y+
     :              YCO(7)*X*X*X+YCO(8)*X*X*Y+
     :              YCO(9)*X*Y*Y+YCO(10)*Y*Y*Y
         ENDDO
      ELSE IF(NC.EQ.15) THEN
         DO I=1,N
            X=XIN(I)-XOR
            Y=YIN(I)-YOR
            XOUT(I)=XCO(1)+XCO(2)*X+XCO(3)*Y+
     :              XCO(4)*X*X+XCO(5)*X*Y+XCO(6)*Y*Y+
     :              XCO(7)*X*X*X+XCO(8)*X*X*Y+
     :              XCO(9)*X*Y*Y+XCO(10)*Y*Y*Y+
     :              XCO(11)*X**4+XCO(12)*X**3*Y+
     :              XCO(13)*X**2*Y**2+XCO(14)*X*Y**3+XCO(15)*Y**4
            YOUT(I)=YCO(1)+YCO(2)*X+YCO(3)*Y+
     :              YCO(4)*X*X+YCO(5)*X*Y+YCO(6)*Y*Y+
     :              YCO(7)*X*X*X+YCO(8)*X*X*Y+
     :              YCO(9)*X*Y*Y+YCO(10)*Y*Y*Y+
     :              YCO(11)*X**4+YCO(12)*X**3*Y+
     :              YCO(13)*X**2*Y**2+YCO(14)*X*Y**3+YCO(15)*Y**4
         ENDDO
      ELSE IF(NC.EQ.21) THEN
         DO I=1,N
            X=XIN(I)-XOR
            Y=YIN(I)-YOR
            XOUT(I)=XCO(1)+XCO(2)*X+XCO(3)*Y+
     :              XCO(4)*X*X+XCO(5)*X*Y+XCO(6)*Y*Y+
     :              XCO(7)*X*X*X+XCO(8)*X*X*Y+
     :              XCO(9)*X*Y*Y+XCO(10)*Y*Y*Y+
     :              XCO(11)*X**4+XCO(12)*X**3*Y+
     :              XCO(13)*X*X*Y*Y+XCO(14)*X*Y**3+XCO(15)*Y**4+
     :              XCO(16)*X**5+XCO(17)*X**4*Y+
     :              XCO(18)*X**3*Y*Y+XCO(19)*X*X*Y**3+
     :              XCO(20)*X*Y**4+XCO(21)*Y**5
            YOUT(I)=YCO(1)+YCO(2)*X+YCO(3)*Y+
     :              YCO(4)*X*X+YCO(5)*X*Y+YCO(6)*Y*Y+
     :              YCO(7)*X*X*X+YCO(8)*X*X*Y+
     :              YCO(9)*X*Y*Y+YCO(10)*Y*Y*Y+
     :              YCO(11)*X**4+YCO(12)*X**3*Y+
     :              YCO(13)*X**2*Y**2+YCO(14)*X*Y**3+YCO(15)*Y**4+
     :              YCO(16)*X**5+YCO(17)*X**4*Y+
     :              YCO(18)*X**3*Y*Y+YCO(19)*X*X*Y**3+
     :              YCO(20)*X*Y**4+YCO(21)*Y**5
         ENDDO
      ENDIF

      RETURN
      END

      SUBROUTINE COE2RD(X,Y,N,RAREF,DECREF,CRPIX1,CRPIX2,
     :                  RPP,RA,DEC)
C
C Convert a pixel position on an image in conic equal-area
C projection to an RA,Dec.
C
C Supplied:
C
C X,Y - double precision - the pixel position (array)
C
C N - integer - the number of values to process
C
C RAREF,DECREF - double precision - the reference point (radians)
C
C CRPIX1,CRPIX2 - double precision - the pixel position corresponding 
C                 to the reference point
C
C RPP - double precision - radians/pixel in image (scale), assumed
C       to be the same in X and Y
C
C Returned:
C
C RA,DEC - double precision - the RA,Dec in radians
C
C Richard Hook, ST-ECF, April 1997
C
      IMPLICIT NONE

      DOUBLE PRECISION X(*),Y(*),RAREF,DECREF
      DOUBLE PRECISION CRPIX1,CRPIX2,RPP
      DOUBLE PRECISION RA(*),DEC(*)
      INTEGER N

C Local variables 
      DOUBLE PRECISION TD,Y0,GAMBY,A,XOFF,YOFF,R2
      INTEGER I

      TD=TAN(DECREF)
      Y0=1.0/TD
      GAMBY=SIN(DECREF)

C Check for the south - where we have to get things
C into the right quadrant  
      IF(DECREF.LT.0.0) THEN
       DO I=1,N
         XOFF=RPP*(X(I)-CRPIX1)
         YOFF=Y0-RPP*(Y(I)-CRPIX2)
         A=ATAN2(-XOFF,-YOFF)
         RA(I)=RAREF-A/GAMBY
         R2=XOFF**2+YOFF**2
         DEC(I)=ASIN(1.0/(GAMBY*2.0)*(1.0+GAMBY**2*(1.0-R2)))
       ENDDO
      ELSE
       DO I=1,N
         XOFF=RPP*(X(I)-CRPIX1)
         YOFF=Y0-RPP*(Y(I)-CRPIX2)
         A=ATAN2(XOFF,YOFF)
         RA(I)=RAREF-A/GAMBY
         R2=XOFF**2+YOFF**2
         DEC(I)=ASIN(1.0/(GAMBY*2.0)*(1.0+GAMBY**2*(1.0-R2)))
       ENDDO
      ENDIF

      RETURN
      END

      SUBROUTINE MERCON(A,B,MAXIM,MAXEN,NEA,NEB,MT,ISTAT)
C
C Merge two context tables.
C
C Supplied:
C
C A(MAXIM,MAXEN) - INTEGER 2d array - the context table which is to 
C                  be merged with
C
C B(MAXIM,MAXEN) - INTEGER 2d array - a pre-existing table.
C
C NEA,NEB - INTEGER the size of these tables
C           (note - NEB is updated)
C
C Returned:
C
C MT - a match table showing the new context numbers in B which
C      have been assigned to the old ones in A.
C
C History:
C
C Experimental first try, Richard Hook, ST-ECF, March 5th 1999
C
      IMPLICIT NONE

      INTEGER MAXEN,MAXIM,MAXMAT
      INTEGER NEA,NEB
      PARAMETER (MAXMAT=100000)
      INTEGER A(MAXIM,MAXEN),B(MAXIM,MAXEN),MT(*)

      INTEGER I,J,II,ISTAT
      LOGICAL SAME 

C Work through the input table
      DO I=1,NEA

C Check for empty contexts and ignore them
       IF(A(2,I).NE.0) THEN

C Initialise to no match
         SAME=.FALSE.

C Compare to each entry in the existing table
         DO J=1,NEB

C Is the number of images the same?
            IF(A(1,I).EQ.B(1,J)) THEN

C Are the rest the same?
               SAME=.TRUE.
               DO II=3,A(1,I)+2
                  IF(A(II,I).NE.B(II,J)) THEN
                     SAME=.FALSE.
                     GO TO 11
                  ENDIF
               ENDDO
            ELSE
               SAME=.FALSE.
            ENDIF

 11         CONTINUE

C If we have a match we can do the merging and continue
            IF(SAME) THEN

C The next line added up the numbers of pixels of a given
C context, we often don't want to do this so it is commented
C out
               B(2,J)=B(2,J)+A(2,I)
CCC               B(2,J)=A(2,I)
               MT(I)=J
               GO TO 22
            ENDIF

         ENDDO

 22      CONTINUE

C If there is no match then we need to append a new entry
         IF(.NOT.SAME) THEN
            NEB=NEB+1
            IF(NEB.GT.MAXEN) THEN
               ISTAT=1
               RETURN
            ENDIF

            DO II=1,A(1,I)+2
               B(II,NEB)=A(II,I)
            ENDDO
            MT(I)=NEB
         ENDIF
       ENDIF
      ENDDO

      ISTAT=0
      RETURN
      END

      SUBROUTINE CONUP(IMAGE,NX,NY,UT,NUP)
C
C Update the values in an array according to an
C update table. This is normally intended for updating a
C context map after the context table has changed.
C
C Richard Hook, March 1999
C 
      IMPLICIT NONE

      INTEGER NX,NY
      INTEGER*2 IMAGE(NX,NY)
      INTEGER NUP
      INTEGER UT(NUP)

      INTEGER SWAP(65536),I,J

C Populate swap table
      DO I=1,65536
         SWAP(I)=I-32769
      ENDDO

C Add correction values
      DO I=1,NUP
         SWAP(I+32769)=UT(I)
      ENDDO

C Loop over the data and change
      DO J=1,NY
         DO I=1,NX
            IMAGE(I,J)=SWAP(IMAGE(I,J)+32769)
         ENDDO
      ENDDO

      RETURN
      END

      CHARACTER*24 FUNCTION ICTIME(T)
C
C Convert the time from an integer (seconds since the start
C of 1980) to the current date.
C
C This was added in January 2000 to support Linux which doesn't
C have "native" TIME and CTIME functions.
C
C Richard Hook, January 11th 2000
C
      IMPLICIT NONE

      INTEGER T
      INTEGER*2 OUTSTR(24)

C Get the time as an integer from the IRAF VOS
      CALL CNVTIE(T,OUTSTR,24)

C Convert to a character string
      CALL F77PAK(OUTSTR,ICTIME,24)

      RETURN
      END

      INTEGER FUNCTION IITIME()
C
C Obtain the time as an integer, seconds since the start
C of 1980, from the IRAF VOS
C
C This was added in January 2000 to support Linux which doesn't
C have "native" TIME and CTIME functions.
C
C Richard Hook, January 11th 2000
C     
      IMPLICIT NONE

      INTEGER T,CLKTIE
      IITIME=CLKTIE(T)

      RETURN
      END