C 
C Utility routines for Drizzle and Blot
C
C This version for drizzle 1.3 and blot 0.6, March 1998
C The module GOLDH was transferred from drizzle as it is also
C now needed by blot.
C
C WCS routines added (xy2rd & rd2xy), September 1998
C
C Modified GTCOEF routine to also handle STIS and NICMOS,
C    January 1999
C
C Added additional code for reversing the cubic transformation added
C  June 1999
C
      SUBROUTINE SETIM(A,NX,NY,V)
C
C Set a 2D array to a specified value 
C
      IMPLICIT NONE

      INTEGER NX,NY
      REAL A(NX,NY),V

      INTEGER I,J

      DO J=1,NY
         DO I=1,NX
            A(I,J)=V
         ENDDO
      ENDDO

      RETURN
      END

      SUBROUTINE MULC(A,NX,NY,V)
C
C Multiply a 2D array to a specified value
C
      IMPLICIT NONE

      INTEGER NX,NY
      REAL A(NX,NY),V

      INTEGER I,J

      DO J=1,NY
         DO I=1,NX
            A(I,J)=A(I,J)*V
         ENDDO
      ENDDO

      RETURN
      END

      SUBROUTINE COPYIM(IN,OUT,NX,NY)
C
C Simply copy an array into another
C
      IMPLICIT NONE

      INTEGER NX,NY
      REAL IN(NX,NY),OUT(NX,NY)

      INTEGER I,J

      DO J=1,NY
         DO I=1,NX
            OUT(I,J)=IN(I,J)
         ENDDO
      ENDDO

      RETURN
      END

      REAL FUNCTION EVALCU(X,Y,CO)
C
C Evaluate a cubic geometric distortion of the form given in the 
C WFPC2 instrument handbook
C
      IMPLICIT NONE

      REAL X,Y,CO(10)

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

      RETURN
      END

      SUBROUTINE GETCO(LUN,LAM,XCO,YCO,ISTAT)
C
C Read geometric distortion coefficients in free format
C from a text file which has already been opened and return them.
C
C If there is a problem reading the numbers then the status flag will be
C returned as 1. If all is well it will be returned
C as 0.
C
C This routine caused problems with SUN SPARC SC4 compilers and was
C modified to avoid them. March 1997
C
      IMPLICIT NONE

      INTEGER LUN,ISTAT,I,J
      CHARACTER*80 LINE
      REAL LAM,MGF2,N
      REAL A,B,C
      REAL XCO(10),YCO(10)
      
C Skip blank lines and those beginning with # until we find a label
      DO WHILE(.TRUE.)
         CALL UFGLIN(LUN,LINE,ISTAT)
         IF(ISTAT.NE.0) GO TO 99

         IF(LINE.NE.' ' .AND. LINE(1:1).NE.'#') THEN

C First "Trauger" style coefficients - 60 of them
            IF(LINE(1:7).EQ.'trauger') THEN

C Now we need to calculate the coefficients which are a
C function of wavelength
               N=MGF2(LAM)

C Now we loop through extracting the coefficients, 3 per line,
C and calculating the wavelength dependence
               J=1
               DO WHILE(J.LE.20)
                  CALL UFGLIN(LUN,LINE,ISTAT)
                  IF(ISTAT.NE.0) GO TO 99

                  IF(LINE.NE.' ' .AND. LINE(1:1).NE.'#') THEN
                     READ(LINE,*,IOSTAT=ISTAT) A,B,C          
                     IF(ISTAT.NE.0) GO TO 99

                     IF(J.LE.10) THEN
                        XCO(J)=A+B*(N-1.5)+C*(N-1.5)**2
                     ELSE
                        YCO(J-10)=A+B*(N-1.5)+C*(N-1.5)**2
                     ENDIF

                     J=J+1
                  ENDIF
               ENDDO

               ISTAT=0
               GO TO 99

C Next "cubic" style - 20 coefficients, read in rows
C of five
            ELSE IF(LINE(1:5).EQ.'cubic') THEN
               J=1
               DO WHILE(J.LT.20)
                CALL UFGLIN(LUN,LINE,ISTAT)
                IF(LINE.NE.' ' .AND. LINE(1:1).NE.'#') THEN
                  IF(J.LE.10) THEN
                   READ(LINE,*,IOSTAT=ISTAT) (XCO(I),I=J,J+4)
                   IF(ISTAT.NE.0) GO TO 99
                  ELSE
                   READ(LINE,*,IOSTAT=ISTAT) (YCO(I-10),I=J,J+4)
                   IF(ISTAT.NE.0) GO TO 99
                  ENDIF
                  J=J+5
                ENDIF
               ENDDO

               ISTAT=0
               GO TO 99
            ELSE
               CALL UMSPUT('! Unknown coefficient set type',
     :                     1,0,ISTAT)
               ISTAT=1
               GO TO 99
            ENDIF
         ENDIF
      ENDDO

99    CONTINUE

      RETURN
      END

      REAL FUNCTION MGF2(LAM)
C
C Calculate the refractive index of MgF2 for a given
C wavelength (in nm) using the formula given by Trauger (1995)
C
      IMPLICIT NONE

      REAL LAM,SIG

      SIG=1.0E7/LAM

      MGF2=SQRT(1.0 + 2.590355E10/(5.312993E10-SIG*SIG) 
     :              + 4.4543708E9/(11.17083E9-SIG*SIG) 
     :              + 4.0838897E5/(1.766361E5-SIG*SIG))

      RETURN
      END

      SUBROUTINE GETWCS(ID,WCS,CTYPE1,CTYPE2,ISTAT)
C
C Get the WCS header items (8 values) and CTYPE strings
C
C The image is assumed to be available with descriptor ID
C and the status flag is returned as non-zero if a parameter
C is not found
C
      IMPLICIT NONE

      INTEGER ID,ISTAT,IS
      CHARACTER*8 CTYPE1,CTYPE2
      DOUBLE PRECISION WCS(8)

      ISTAT=0

      CALL UHDGST(ID,'CTYPE1',CTYPE1,ISTAT)
      IF(ISTAT.NE.0) CALL UMSPUT(
     :   '! Warning, unable to read CTYPE1 from header',1,0,IS)
      IF(ISTAT.NE.0) RETURN

      CALL UHDGST(ID,'CTYPE2',CTYPE2,ISTAT)
      IF(ISTAT.NE.0) CALL UMSPUT(
     :   '! Warning, unable to read CTYPE2 from header',1,0,IS)
      IF(ISTAT.NE.0) RETURN

      CALL UHDGSD(ID,'CRPIX1',WCS(1),ISTAT)
      IF(ISTAT.NE.0) CALL UMSPUT(
     :   '! Warning, unable to read CRPIX1 from header',1,0,IS)
      IF(ISTAT.NE.0) RETURN  

      CALL UHDGSD(ID,'CRVAL1',WCS(2),ISTAT)
      IF(ISTAT.NE.0) CALL UMSPUT(
     :   '! Warning, unable to read CRVAL1 from header',1,0,IS)
      IF(ISTAT.NE.0) RETURN  

      CALL UHDGSD(ID,'CRPIX2',WCS(3),ISTAT)
      IF(ISTAT.NE.0) CALL UMSPUT(
     :   '! Warning, unable to read CRPIX2 from header',1,0,IS)
      IF(ISTAT.NE.0) RETURN  

      CALL UHDGSD(ID,'CRVAL2',WCS(4),ISTAT)
      IF(ISTAT.NE.0) CALL UMSPUT(
     :   '! Warning, unable to read CRVAL2 from header',1,0,IS)
      IF(ISTAT.NE.0) RETURN  

      CALL UHDGSD(ID,'CD1_1',WCS(5),ISTAT)
      IF(ISTAT.NE.0) CALL UMSPUT(
     :   '! Warning, unable to read CD1_1 from header',1,0,IS)
      IF(ISTAT.NE.0) RETURN  

      CALL UHDGSD(ID,'CD2_1',WCS(6),ISTAT)
      IF(ISTAT.NE.0) CALL UMSPUT(
     :   '! Warning, unable to read CD2_1 from header',1,0,IS)
      IF(ISTAT.NE.0) RETURN  

      CALL UHDGSD(ID,'CD1_2',WCS(7),ISTAT)
      IF(ISTAT.NE.0) CALL UMSPUT(
     :   '! Warning, unable to read CD1_2 from header',1,0,IS)
      IF(ISTAT.NE.0) RETURN  

      CALL UHDGSD(ID,'CD2_2',WCS(8),ISTAT)
      IF(ISTAT.NE.0) CALL UMSPUT(
     :   '! Warning, unable to read CD2_2 from header',1,0,IS)

      RETURN
      END

      SUBROUTINE UWCS(ID,WCS,CTYPE1,CTYPE2)
C
C Update the WCS header items using the 8 values
C supplied in WCS and the values of CTYPE1 and CTYPE2.
C
C The image is assumed to be available with descriptor ID
C and to already have all the parameters
C
      IMPLICIT NONE

      INTEGER ID,ISTAT
      CHARACTER*8 CTYPE1,CTYPE2
      DOUBLE PRECISION WCS(8)

      CALL UHDPST(ID,'CTYPE1',CTYPE1,ISTAT)
      CALL UHDPST(ID,'CTYPE2',CTYPE2,ISTAT)
      CALL UHDPSD(ID,'CRPIX1',WCS(1),ISTAT)
      CALL UHDPSD(ID,'CRVAL1',WCS(2),ISTAT)
      CALL UHDPSD(ID,'CRPIX2',WCS(3),ISTAT)
      CALL UHDPSD(ID,'CRVAL2',WCS(4),ISTAT)
      CALL UHDPSD(ID,'CD1_1',WCS(5),ISTAT)
      CALL UHDPSD(ID,'CD2_1',WCS(6),ISTAT)
      CALL UHDPSD(ID,'CD1_2',WCS(7),ISTAT)
      CALL UHDPSD(ID,'CD2_2',WCS(8),ISTAT)

      RETURN
      END

      SUBROUTINE PUTFIL(DAT,COU,ONX,ONY,FILVAL)
C
C Put in the missing filling value where the weight is
C zero
C
      IMPLICIT NONE

      INTEGER ONX,ONY
      REAL DAT(ONX,ONY),COU(ONX,ONY),FILVAL

      INTEGER I,J

      DO J=1,ONY
         DO I=1,ONX
            IF(COU(I,J).EQ.0.0) DAT(I,J)=FILVAL
         ENDDO
      ENDDO
   
      RETURN
      END

      SUBROUTINE GTCOEF(ID,COEFFS,LAMBDA,ISTAT)
C
C This routine is invoked when it has been requested that
C the coefficients file name for geometric distortion, and
C the associated wavelength is to be extracted from the
C header. 
C
C It is currently handles WFPC2, STIS and NICMOS and assumes
C coefficients in a directory drizzle$coeffs.
C
C Extended for STIS and NICMOS, January 1999
C
      IMPLICIT NONE

      INTEGER ID,ISTAT,IDET
      REAL LAMBDA,PLAM
      CHARACTER*(*) COEFFS
      CHARACTER*8 INST,DETECT
      CHARACTER*80 CHARS

C First try to get the instrument name
      CALL UHDGST(ID,'INSTRUME',INST,ISTAT)
      IF(ISTAT.NE.0) RETURN

C Check instrument can be handled
      IF(INST.EQ.'WFPC2') THEN

C Get the other keywords we need
         CALL UHDGSI(ID,'DETECTOR',IDET,ISTAT)
         IF(ISTAT.NE.0) RETURN

         IF(IDET.LT.1 .OR. IDET.GT.4) THEN
            CALL UMSPUT('! Invalid detector number in header',
     :               1,0,ISTAT)
            ISTAT=1
            RETURN
         ENDIF

         CALL UHDGSR(ID,'PHOTPLAM',PLAM,ISTAT)
         IF(ISTAT.NE.0) RETURN

C We have to convert the wavelength from Angstroms to nm
         LAMBDA=PLAM/10.0

C Now we build the full "Trauger style" name for the coefficients
C file
         IF(IDET.EQ.1) THEN
            COEFFS='drizzle$coeffs/pc1-trauger'
         ELSE IF(IDET.EQ.2) THEN
            COEFFS='drizzle$coeffs/wf2-trauger'
         ELSE IF(IDET.EQ.3) THEN
            COEFFS='drizzle$coeffs/wf3-trauger'
         ELSE IF(IDET.EQ.4) THEN
            COEFFS='drizzle$coeffs/wf4-trauger'
         ENDIF

         WRITE(CHARS,
     : '(''-Instrument: '',A8,'' Chip: '''//
     : ',I1,'' Wavelength: '',F7.2,'' (nm)'')')
     :   INST,IDET,LAMBDA
         CALL UMSPUT(CHARS,1,0,ISTAT)
      ELSE IF(INST.EQ.'NICMOS') THEN

C Get the other keywords we need
         CALL UHDGSI(ID,'CAMERA',IDET,ISTAT)
         IF(ISTAT.NE.0) RETURN

         IF(IDET.LT.1 .OR. IDET.GT.3) THEN
            CALL UMSPUT('! Invalid camera number in header',
     :               1,0,ISTAT)
            ISTAT=1
            RETURN
         ELSE
            WRITE(COEFFS,'(''drizzle$coeffs/nic-'',I1)') IDET
         ENDIF

         WRITE(CHARS,
     : '(''-Instrument: NICMOS  Camera: '',i1)') IDET
         CALL UMSPUT(CHARS,1,0,ISTAT)
      ELSE IF(INST.EQ.'STIS') THEN

C Get the other keywords we need
         CALL UHDGST(ID,'DETECTOR',DETECT,ISTAT)
         IF(DETECT.EQ.'CCD') THEN
            COEFFS='drizzle$coeffs/stis-ccd'
         ELSE IF(DETECT.EQ.'NUV-MAMA') THEN
            COEFFS='drizzle$coeffs/stis-nuv-mama'
         ELSE IF(DETECT.EQ.'FUV-MAMA') THEN
            COEFFS='drizzle$coeffs/stis-fuv-mama'
         ELSE
            CALL UMSPUT('! Invalid STIS detector name in header',
     :               1,0,ISTAT)
            ISTAT=1
            RETURN
         ENDIF

         WRITE(CHARS,
     : '(''-Instrument: STIS    Detector: '',A8)') DETECT
         CALL UMSPUT(CHARS,1,0,ISTAT)
      ELSE
         CALL UMSPUT(
     :  '! We can currently only handle WFPC2, NICMOS or STIS headers',
     :               1,0,ISTAT)
         ISTAT=1
         RETURN
      ENDIF

      RETURN
      END

      SUBROUTINE GOLDH(ID,IC,LASTPF,LASTSC,
     :                 LASTUN,LASTOX,LASTOY,ISTAT)
C
C This routine finds out details of the last image
C to have been drizzled onto the image with identifier
C ID and returns then for comparison with the new
C ones.
C
C This version was moved from the source code of V1.2 of
C drizzle.
C
      IMPLICIT NONE

      INTEGER IC,ISTAT,ID
      CHARACTER*8 LASTUN
      CHARACTER*4 STEM
      CHARACTER*80 BUFFER
      REAL LASTPF,LASTSC,LASTOX,LASTOY

C Before doing anything check for old-style headers and
C warn the user
      CALL UHDGST(ID,'D01VER',BUFFER,ISTAT)
      IF(ISTAT.EQ.0) THEN
         CALL UMSPUT(
     : '! Warning - there is an old-style Drizzle header present,',
     : 1,0,ISTAT)
         CALL UMSPUT(
     : '! please consult the header for full processing history',
     : 1,0,ISTAT)
      ENDIF

C First find out if there are old values
C This is very ugly
      STEM='D000'

      DO IC=1,999
         WRITE(STEM(2:4),'(I3)') IC
         IF(STEM(2:2).EQ.' ') STEM(2:2)='0'
         IF(STEM(3:3).EQ.' ') STEM(3:3)='0'
         CALL UHDGST(ID,STEM//'VER',BUFFER,ISTAT)

C Check for a status code indicating that this header
C item doesn't exist
         IF(ISTAT.EQ.40) GO TO 111
      ENDDO

111   CONTINUE
      IC=IC-1

      WRITE(STEM(2:4),'(I3)') IC
      IF(STEM(2:2).EQ.' ') STEM(2:2)='0'
      IF(STEM(3:3).EQ.' ') STEM(3:3)='0'

      CALL UHDGSR(ID,STEM//'PIXF',LASTPF,ISTAT)
      IF(ISTAT.NE.0) RETURN

      CALL UHDGSR(ID,STEM//'SCAL',LASTSC,ISTAT)
      IF(ISTAT.NE.0) RETURN

      CALL UHDGST(ID,STEM//'OUUN',LASTUN,ISTAT)
      IF(ISTAT.NE.0) RETURN

      CALL UHDGSR(ID,STEM//'OUXC',LASTOX,ISTAT)
      IF(ISTAT.NE.0) RETURN

      CALL UHDGSR(ID,STEM//'OUYC',LASTOY,ISTAT)
      IF(ISTAT.NE.0) RETURN

      ISTAT=0

      RETURN
      END

      SUBROUTINE INMAT(P,Q,R,S)
C
C Invert a 2 by 2 double precision matrix
C in place. It assumes that the matrix is not
C singular.
C
      IMPLICIT NONE

      DOUBLE PRECISION P,Q,R,S
      DOUBLE PRECISION A,B,C,D,DET

      DET=P*S-Q*R

      A=S/DET
      B=-Q/DET
      C=-R/DET
      D=P/DET

      P=A
      Q=B
      R=C
      S=D

      RETURN
      END

      SUBROUTINE XY2RD(X,Y,R,D,WCS)
C
C XY2RD - convert a pixel position to an equatorial (RA,DEC)
C         position assuming a TAN projection and WCS with the
C         standard 8 elements.
C
C It is based on the xy2rd task in STSDAS
C
C Output angles are in degrees
C
C Richard Hook, STScI, 11th August 1998
C
      IMPLICIT NONE

      DOUBLE PRECISION X,Y,R,D,WCS(8)
      DOUBLE PRECISION XI,ETA,RA0,DEC0

      DOUBLE PRECISION PIBY
      PARAMETER (PIBY=3.141592653589793/180.0)

C First convert pixel coordinates to tangent plane in radians
      XI=(WCS(5)*(X-WCS(1))+WCS(7)*(Y-WCS(3)))*PIBY
      ETA=(WCS(6)*(X-WCS(1))+WCS(8)*(Y-WCS(3)))*PIBY

C And convert the reference point on the sky to radians
      RA0=WCS(2)*PIBY
      DEC0=WCS(4)*PIBY

C Now go to equatorial from tangent plane
      R=ATAN2(XI, COS(DEC0)-ETA*SIN(DEC0)) + RA0
      D=ATAN2(ETA*COS(DEC0)+SIN(DEC0),
     :        SQRT((COS(DEC0)-ETA*SIN(DEC0))**2 + XI**2))

C Convert back to degrees and check the range
      R=R/PIBY
      D=D/PIBY
      IF(R.LT.0.0) R=R+360.0

      RETURN
      END

      SUBROUTINE RD2XY(R,D,X,Y,WCS,ISTAT)
C
C RD2XY - convert an equatorial (RA,DEC) position to a pixel position
C         assuming a TAN projection and WCS with the
C         standard 8 elements.
C
C It is based on the rd2xy task in STSDAS
C
C Input angles are in degrees
C
C Richard Hook, STScI, 13th August 1998
C
      IMPLICIT NONE

      INTEGER ISTAT
      DOUBLE PRECISION X,Y,R,D,WCS(8),RA,DEC
      DOUBLE PRECISION XI,ETA,RA0,DEC0
      DOUBLE PRECISION DET,CDINV(2,2),BOTTOM
      DOUBLE PRECISION PIBY
      PARAMETER (PIBY=3.141592653589793/180.0)

C First invert the CD matrix
      DET=WCS(5)*WCS(8)-WCS(7)*WCS(6)

      IF(DET.EQ.0.0) THEN
         ISTAT=1
         RETURN
      ENDIF

      CDINV(1,1)=WCS(8)/DET
      CDINV(1,2)=-WCS(7)/DET
      CDINV(2,1)=-WCS(6)/DET
      CDINV(2,2)=WCS(5)/DET

C Translate from RA,Dec to X,Y
      RA0=WCS(2)*PIBY
      DEC0=WCS(4)*PIBY

      RA=R*PIBY
      DEC=D*PIBY

      BOTTOM=SIN(DEC)*SIN(DEC0)+COS(DEC)*COS(DEC0)*COS(RA-RA0)
      IF(BOTTOM.EQ.0) THEN
         ISTAT=1
         RETURN
      ENDIF

C Calculate tangent plane position and convert to degrees
      XI=COS(DEC)*SIN(RA-RA0)/BOTTOM/PIBY
      ETA=(SIN(DEC)*COS(DEC0)-COS(DEC)*SIN(DEC0)*COS(RA-RA0))/
     :     BOTTOM/PIBY

C Convert back to pixels using the inverse of the CD matrix
      X=CDINV(1,1)*XI+CDINV(1,2)*ETA+WCS(1)
      Y=CDINV(2,1)*XI+CDINV(2,2)*ETA+WCS(3)

      ISTAT=0
      RETURN
      END

      SUBROUTINE WCSLIN(WCSIN,WCSOUT,NONLIN,
     :            XCEN,YCEN,XCO,YCO,XC,YC,XS,YS,XT,YT)
C
C WCSLIN - derive best linear transformation coefficients to map
C          pixel positions from one world coordinate system to another.
C
C Supplied:
C
C  WCSIN - double precision 8 element array. The WCS of the input image.
C          This includes CRPIX1/2, CRVAL1/2 and a CD matrix.
C
C  WCSOUT - double precision 8 element array. The WCS of the output image.
C          This includes CRPIX1/2, CRVAL1/2 and a CD matrix.
C
C  NONLIN - logical, whether there are non-linear distortion coefficients
C
C  XCEN,YCEN - reals, the centre of the chip, for reference
C
C  XCO,YCO - double precision arrays, the non-linear distortion coefficients
C
C Returned:
C
C XC,YC,XS,YS,XT,YT - reals, linear transformation coefficients
C
C   Richard Hook, 18th September 1998
C
      IMPLICIT NONE

C Input parameters
      DOUBLE PRECISION WCSIN(8),WCSOUT(8)
      REAL XCO(*),YCO(*),XCEN,YCEN
      LOGICAL NONLIN

C Output parameters
      REAL XC,YC,XS,YS,XT,YT

C Local variables
      DOUBLE PRECISION R,D
      DOUBLE PRECISION XIN(4),YIN(4),XOUT(4),YOUT(4)
      REAL EVALCU,XX,YY
      INTEGER I,ISTAT

C Set up a central square for reference
      XIN(1)=DBLE(XCEN)
      XIN(2)=DBLE(XCEN)
      XIN(3)=DBLE(XCEN+1.0)
      XIN(4)=DBLE(XCEN+1.0)
      YIN(1)=DBLE(YCEN) 
      YIN(2)=DBLE(YCEN+1.0)
      YIN(3)=DBLE(YCEN+1.0)
      YIN(4)=DBLE(YCEN) 

C Transform these points into the WCS and then back out again
C using the target WCS
      DO I=1,4
         CALL XY2RD(XIN(I),YIN(I),R,D,WCSIN)
         CALL RD2XY(R,D,XOUT(I),YOUT(I),WCSOUT,ISTAT)
      ENDDO

C Now we apply the geometric distortion to the input points so that
C the linear transformation which we derive is appropriate after
C the distortion is corrected
      DO I=1,4
         XX=SNGL(XIN(I))
         YY=SNGL(YIN(I))
         IF(NONLIN) THEN
            XIN(I)=EVALCU(XX-XCEN,YY-YCEN,XCO)
            YIN(I)=EVALCU(XX-XCEN,YY-YCEN,YCO)
         ELSE
            XIN(I)=XX-XCEN
            YIN(I)=YY-YCEN
         ENDIF
      ENDDO

C Now we have the inputs and outputs and can derive the linear
C transform between them

C First we get the linear shift from the reference pixel
      XT=SNGL(XOUT(1))-XIN(1)
      YT=SNGL(YOUT(1))-YIN(1)

C And now the scales (ie, the CD matrix)
      XC=SNGL(0.5*((XOUT(4)-XOUT(1))/(XIN(4)-XIN(1)) +
     :        (XOUT(3)-XOUT(2))/(XIN(3)-XIN(2))))
      YS=SNGL(-0.5*((XOUT(2)-XOUT(1))/(YIN(2)-YIN(1)) +
     :        (XOUT(3)-XOUT(4))/(YIN(3)-YIN(4))))
      XS=SNGL(0.5*((YOUT(4)-YOUT(1))/(XIN(4)-XIN(1)) +
     :        (YOUT(3)-YOUT(2))/(XIN(3)-XIN(2))))
      YC=SNGL(0.5*((YOUT(2)-YOUT(1))/(YIN(2)-YIN(1)) +
     :        (YOUT(3)-YOUT(4))/(YIN(3)-YIN(4))))

      RETURN
      END

      SUBROUTINE SETWCS(OUTSCL,ORIENT,CRPIX1,CRVAL1,
     :                  CRPIX2,CRVAL2,WCS)
C
C SETWCS - convert a scale and orientation, along with
C          reference pixel and position on sky, to a WCS
C          including a CD matrix.
C
C Note: OUTSCL is in arcsecs/pix, ORIENT and CRVAL1/2
C       are in degrees.
C
C Richard Hook, 22nd September 1998
C
      IMPLICIT NONE

C Supplied:
      REAL OUTSCL,CRPIX1,CRPIX2,ORIENT
      DOUBLE PRECISION CRVAL1,CRVAL2

C Returned:
      DOUBLE PRECISION WCS(8)

C Local variables
      DOUBLE PRECISION PIBY,ROR,DSC
      PARAMETER (PIBY=3.14159265358979/180.0)

C Convert to radians for orientation
      ROR=DBLE(ORIENT)*PIBY

C and to degrees for scale
      DSC=DBLE(OUTSCL)/3600.0

C Set reference points
      WCS(1)=DBLE(CRPIX1)
      WCS(3)=DBLE(CRPIX2)
      WCS(2)=CRVAL1
      WCS(4)=CRVAL2

C and CD matrix (no skew, equal X,Y scales)
      WCS(5)=-DSC*DCOS(ROR)
      WCS(6)=DSC*DSIN(ROR)
      WCS(7)=WCS(6) 
      WCS(8)=-WCS(5)       

      RETURN
      END

      SUBROUTINE COPY1D(IN,OUT,N)
C
C Copy a 1d real array from one place to another.
C
      IMPLICIT NONE

      INTEGER N
      REAL IN(*),OUT(*)

      INTEGER I

      DO I=1,N
         OUT(I)=IN(I)
      ENDDO
 
      RETURN
      END

      SUBROUTINE LCORN(DNX,DNY,ONX,ONY,
     :           XSH,YSH,ROT,SCALE,ALIGN,ROTFIR,
     :           USEWCS,WCSIN,WCSOUT,
     :           XCO,YCO,XMIN,XMAX,YMIN,YMAX,ISTAT)
 
C
C LCORN - transform the range of input pixel coordinates covered
C         by an image onto the output and find the extreme values.
C
      IMPLICIT NONE

      INTEGER DNX,DNY,ONX,ONY
      REAL XSH,YSH,ROT,SCALE
      LOGICAL ROTFIR,USEWCS
      CHARACTER*8 ALIGN
      DOUBLE PRECISION WCSIN(8),WCSOUT(8)
      REAL XCO(10),YCO(10)
      INTEGER XMIN,XMAX,YMIN,YMAX,ISTAT
      REAL XIN(4),YIN(4)
      REAL XOUT(4),YOUT(4),XMA,YMA,XMI,YMI,XR,YR

      XIN(1)=1.0
      XIN(2)=1.0
      XIN(3)=REAL(DNX)
      XIN(4)=REAL(DNX)
      YIN(1)=1.0
      YIN(2)=REAL(DNY)
      YIN(3)=REAL(DNY)
      YIN(4)=1.0
      
C Transform onto output coordinate system
      CALL DRIVAL(XIN,YIN,4,DNX,DNY,ONX,ONY,
     :            XSH,YSH,ROT,SCALE,ALIGN,ROTFIR,
     :            USEWCS,WCSIN,WCSOUT,
     :            XCO,YCO,XOUT,YOUT)

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 5% margin of error
      XMAX=NINT(XMA+0.05*XR)
      XMIN=NINT(XMI-0.05*XR)
      YMAX=NINT(YMA+0.05*YR)
      YMIN=NINT(YMI-0.05*YR)
 
      ISTAT=0
      RETURN
      END

      SUBROUTINE DRIVAL(XIN,YIN,N,DNX,DNY,ONX,ONY,
     :            XSH,YSH,ROT,SCALE,ALIGN,ROTFIR,
     :            USEWCS,WCSIN,WCSOUT,
     :            XCO,YCO,XOUT,YOUT)

C
C DRIVAL - apply the standard Drizzle transformation
C          from input to output pixel coordinates.
C          This may optional include a polynomial
C          distortion solution or be specified using
C          input and output WCS. It duplicates the
C          logic of the DOBOX routine within Drizzle.
C
C  Richard Hook, ST-ECF, 5th October 1998
C
      IMPLICIT NONE

      INTEGER N,DNX,DNY,ONX,ONY,I
      REAL XIN(N),YIN(N),XOUT(N),YOUT(N)
      REAL XCEN,YCEN,XF,YF,XCORN,YCORN
      REAL XSH,YSH,ROT,SCALE,XCO(*),YCO(*)
      REAL XS,XC,YS,YC
      REAL EVALCU,SINTH,COSTH,XOFF,YOFF,XT,YT,XP,YP
      LOGICAL ROTFIR,USEWCS,NONLIN
      CHARACTER*8 ALIGN
      DOUBLE PRECISION WCSIN(8),WCSOUT(8)

      IF(ALIGN.EQ.'corner') THEN
         XCEN=FLOAT(DNX/2)+0.5
         YCEN=FLOAT(DNY/2)+0.5
      ELSE
         XCEN=FLOAT(DNX/2)+1.0
         YCEN=FLOAT(DNY/2)+1.0
      ENDIF
 
C If all the non-linear coefficients are zero we can
C ignore that part of things
      NONLIN=.FALSE.
      DO I=1,10
         IF(XCO(I).NE.0.0 .OR. YCO(I).NE.0.0) NONLIN=.TRUE.
      ENDDO
 
C Calculate some numbers to simplify things later
      SINTH=SIN(ROT)
      COSTH=COS(ROT)
 
C The shifts are in units of INPUT pixels.
      XOFF=XSH/SCALE
      YOFF=YSH/SCALE
      XF=1.0/SCALE
      YF=1.0/SCALE
 
C Some more economies
      XS=XF*SINTH
      XC=XF*COSTH
      YS=YF*SINTH
      YC=YF*COSTH
 
      IF(ALIGN.EQ.'corner') THEN
         XP=REAL(ONX/2)+0.5
         YP=REAL(ONY/2)+0.5
      ELSE
         XP=REAL(ONX/2)+1.0
         YP=REAL(ONY/2)+1.0
      ENDIF
 
      XT=XOFF+XP
      YT=YOFF+YP

C Here is some new code to support the WCS option
C first we work out a linear transformation from input to
C out
      IF(USEWCS) THEN
         CALL WCSLIN(WCSIN,WCSOUT,NONLIN,
     :         XCEN,YCEN,XCO,YCO,XC,YC,XS,YS,XT,YT)
         ROTFIR=.TRUE.
      ENDIF

      DO I=1,N
         IF(NONLIN) THEN
            XCORN=EVALCU(XIN(I)-XCEN,YIN(I)-YCEN,XCO)
            YCORN=EVALCU(XIN(I)-XCEN,YIN(I)-YCEN,YCO)
         ELSE
            XCORN=XIN(I)-XCEN
            YCORN=YIN(I)-YCEN
         ENDIF

C Apply the linear transform
C There are two ways this can be done - shift then
C rotate or rotate then shift
         IF(ROTFIR) THEN
            XOUT(I)=XC*XCORN-YS*YCORN+XT
            YOUT(I)=XS*XCORN+YC*YCORN+YT
         ELSE
            XOUT(I)=XC*(XCORN+XSH)-YS*(YCORN+YSH)+XP
            YOUT(I)=XS*(XCORN+XSH)+YC*(YCORN+YSH)+YP
         ENDIF
      ENDDO

      RETURN
      END
      SUBROUTINE DERICU(X,Y,CO,DX,DY)
C
C Evaluate the derivatives of a cubic polynomial distortion with
C respect to X and Y. This is mainly for inverting the distortion
C using Newton-Raphson.
C
C Richard Hook, STScI, May 1999
C
      IMPLICIT NONE

      REAL X,Y,CO(*),DX,DY

      DX=CO(2)+
     :   2.0*CO(4)*X+
     :   CO(5)*Y+
     :   3.0*CO(7)*X*X+
     :   2.0*CO(8)*X*Y+
     :   CO(9)*Y*Y
 
      DY=CO(3)+
     :   CO(5)*X+
     :   2.0*CO(6)*Y+
     :   CO(8)*X*X+
     :   2.0*CO(9)*X*Y+
     :   3.0*CO(10)*Y*Y

      RETURN
      END

      SUBROUTINE INVECU(XOUT,YOUT,XCO,YCO,ERR,XIN,YIN)
C
C Invert a cubic distortion in 2d using Newton-Raphson
C
C Supplied:
C
C Xout,Yout - real - the output position from the geometric
C                     distortion.
C
C Xco,Yco - real arrays - the cubic distortion coefficients
C
C Err - real - the accuracy required of the inversion (in X and Y)
C
C Returned:
C
C Xin, Yin - real - the position which, when distorted, ends up
C                   at Xout,Yout to accuracy Err.
C
C This method is interative and relatively slow.
C
C Richard Hook, STScI, May 1999
C
      IMPLICIT NONE

      REAL XOUT,YOUT,XCO(*),YCO(*),XIN,YIN,ERR
      REAL DXX,DXY,DYX,DYY,XO,YO,X,Y,D
      REAL EVALCU
 
C First guess
      X=2.0*XOUT-EVALCU(XOUT,YOUT,XCO)
      Y=2.0*YOUT-EVALCU(XOUT,YOUT,YCO)

      DO WHILE(.TRUE.)

C Distort current guesses
         XO=EVALCU(X,Y,XCO)
         YO=EVALCU(X,Y,YCO)

C Check for error criterion
         IF(ABS(XOUT-XO).LT.ERR .AND.
     :   ABS(YOUT-YO).LT.ERR) GO TO 999

C Calculate derivatives - there are four of these
         CALL DERICU(X,Y,XCO,DXX,DXY)
         CALL DERICU(X,Y,YCO,DYX,DYY)

C Work out the determinant
         D=DXX*DYY-DYX*DXY

C Improve guess with Newton-Raphson
         X=X+((XOUT-XO)*DYY-(YOUT-YO)*DXY)/D
         Y=Y+((YOUT-YO)*DXX-(XOUT-XO)*DYX)/D
      ENDDO

 999  CONTINUE
      XIN=X
      YIN=Y

      RETURN
      END