C 
C Utility routines for Drizzle and Blot
C
C This version for drizzle 1.2 and blot 0.6, 25th February 1998
C The module GOLDH was transferred from drizzle as it is also
C now needed by blot.
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. It is currently specific for WFPC2 and assumes
C Trauger coefficients in a directory drizzle$coeffs.
C
      IMPLICIT NONE

      INTEGER ID,ISTAT,IDET
      REAL LAMBDA,PLAM
      CHARACTER*(*) COEFFS
      CHARACTER*8 INST
      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.NE.'WFPC2') THEN
         CALL UMSPUT('! We can currently only handle WFPC2 headers',
     :               1,0,ISTAT)
         ISTAT=1
         RETURN
      ELSE

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)
      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