SUBROUTINE BLOT
C++
C
C BLOT.F Version 0.6   "Reverse drizzling" 
C
C It is based on ideas of Andy Fruchter and was coded by Richard Hook. 
C
C Version 0.6 is compatible and tested with Drizzle V1.2.
C
C This task uses coefficients held in a text file to define
C geometrical distortions.
C
C This code uses the F77/VOS for IRAF i/o.
C It may be compiled and linked using a command like:
C
C xc -p stsdas (-z) blot.f drutil.f x_blot.x <dir>libiraf77.a <dir>libiminterp.a
C
C Note that the optional -z switch (which tells IRAF to create a static
C executable which doesn't use the IRAF sharable library) may be
C necessary to avoid a virtual memory use limitation of around 167MB in
C IRAF versions up to 2.10.4.
C
C History:
C
C Created first try (V0.1) from DRIZZLE.F V0.97, 9th July 1996
C Modified for new origin of coordinates etc (V0.2), 17th July 1996
C Modified to match up with DRIZZLE V1.05, 14th November 1996
C (rotfir parameter added, few bug fixes, change coordinate origin)
C Match up to Drizzle 1.08 - V0.4 of BLOT, 6th December 1996
C Correct shift units bug, V0.41, 16th December 1996
C Add ALIGN option, V0.42, 16th January 1997
C Use "drutil.f" utilities, no -z, 4th February 1997
C Prepared V0.5 for release, matches V1.097 of Drizzle, 17th March 1997
C
C First public release of Blot, 12th May 1997
C
C Added in_un/out_un switches for compatibility with V1.2 of Drizzle.
C This was for the first release of drizzle/blot in the dither package
C within STSDAS. 25th February 1998.
C
C Added automatic output image size selection, 27th February 1998
C--
      IMPLICIT NONE

C Iraf global storage (real in this program)
      REAL MEMR(1) !I
      COMMON /MEM/MEMR         !I

C Local variables
      INTEGER ISTAT,I,LUN
      INTEGER NDIMS,DNX,DNY,ONX,ONY,DDIMS(7),ODIMS(7),DATTYP
      INTEGER IDD,IDND
      INTEGER PDATA,PNDAT
      INTEGER INTYPE,IC
      REAL XCO(10),YCO(10),XSH,YSH,LAM,MISVAL
      REAL SCALE,ROT,EF,EXPIN,EXPOUT
      REAL LASTPF,LASTSC,LASTOX,LASTOY
      DOUBLE PRECISION WCS(8)
      CHARACTER*80 DATA,OUTDAT,CHARS
      CHARACTER*80 COEFFS
      CHARACTER*8 CTYPE1,CTYPE2,SHFTFR,SHFTUN,ALIGN
      CHARACTER*8 EXPKEY,INUN,OUTUN,LASTUN
      CHARACTER*80 EXPSTR
      CHARACTER*45 VERS
      CHARACTER*7 INTERP

C Logical flags
      LOGICAL VERBOSE
      LOGICAL NOCO
      LOGICAL GOTWCS
      LOGICAL ROTFIR
      LOGICAL HCOEFF
      LOGICAL INCPS
      LOGICAL OUTCPS

C-- Start of executable code

C First announce the version
      VERS='BLOT Version 0.6 (February 1998)'

      CALL UMSPUT('+ '//VERS,1,0,ISTAT)

      VERBOSE=.TRUE.

C Get reference frame for shifts - specifying whether shifts are done before
C or after rotation
      CALL UCLGST('shft_fr',SHFTFR,ISTAT)
      IF(SHFTFR.EQ.'output') THEN
         ROTFIR=.FALSE.
      ELSE
         ROTFIR=.TRUE.
      ENDIF

C Get the alignment parameter
      CALL UCLGST('align',ALIGN,ISTAT)

C Get the exposure related parameters
      CALL UCLGST('expkey',EXPKEY,ISTAT)
      CALL UCLGST('expout',EXPSTR,ISTAT)

C Get the switch value for whether to read in CPS or
C counts (latter is the default)
      CALL UCLGST('in_un',INUN,ISTAT)
      IF(INUN.EQ.'cps') THEN
         INCPS=.TRUE.
      ELSE
         INCPS=.FALSE.
      ENDIF

C Get the switch value for whether to write out in CPS or
C counts (latter is the default)
      CALL UCLGST('out_un',OUTUN,ISTAT)
      IF(OUTUN.EQ.'cps') THEN
         OUTCPS=.TRUE.
      ELSE
         OUTCPS=.FALSE.
      ENDIF

C Get the name of the file containing the coefficients
      CALL UCLGST('coeffs',COEFFS,ISTAT)

C Get the scale and rotation for the final image
      CALL UCLGSR('scale',SCALE,ISTAT)
      CALL UCLGSR('rot',ROT,ISTAT)

C Check for invalid scale
      IF(SCALE.EQ.0.0) THEN
         CALL UMSPUT('! Invalid scale',1,0,ISTAT)
         GO TO 99
      ENDIF

C Convert the rot to radians
      ROT=ROT*3.1415926/180.0

C Get the X,Y final shift
      CALL UCLGSR('xsh',XSH,ISTAT)
      CALL UCLGSR('ysh',YSH,ISTAT)

C Also get the binary flag for whether units are input or
C output pixels
      CALL UCLGST('shft_un',SHFTUN,ISTAT)

      IF(SHFTUN.EQ.'input') THEN
         XSH=XSH*SCALE
         YSH=YSH*SCALE
      ENDIF

C Note that the shifts are in INPUT pixel units

C Get the missing value flag
      CALL UCLGSR('fillval',MISVAL,ISTAT)

C Get the name of the input data image
      CALL UCLGST('data',DATA,ISTAT)

C Try to access the input data image
      CALL UIMOPN(DATA,1,IDD,ISTAT)
      IF(ISTAT.NE.0) THEN
         CALL UMSPUT('! Unable to open data image',
     :                    1,0,ISTAT)
         GO TO 99
      ELSE
         IF(VERBOSE) CALL UMSPUT(
     :  '-Opening input data image: '//DATA,
     :                   1,0,ISTAT)
      ENDIF

C Get the array size and shape
      CALL UIMGID(IDD,DATTYP,NDIMS,DDIMS,ISTAT)
      IF(NDIMS.NE.2) THEN
         CALL UMSPUT('! Input image is not two dimensional',
     :                    1,0,ISTAT)
         GO TO 99
      ELSE
         DNX=DDIMS(1)
         DNY=DDIMS(2)
      ENDIF

C We now get some of the previous drizzle information from
C the header of the current image and check for funny
C combinations
      CALL GOLDH(IDD,IC,LASTPF,LASTSC,
     :           LASTUN,LASTOX,LASTOY,ISTAT)

C Check for incompatibilities
      IF(INCPS) THEN
         IF(LASTUN.EQ.'counts' .OR.
     :                LASTUN.EQ.'COUNTS') THEN
          CALL UMSPUT('! Warning input image header gives units',
     :               1,0,ISTAT)
          CALL UMSPUT('! as counts NOT cps! Assuming cps',
     :               1,0,ISTAT)
         ENDIF
      ELSE
         IF(LASTUN.EQ.'cps' .OR.
     :                LASTUN.EQ.'CPS') THEN
          CALL UMSPUT('! Warning input image header gives units',
     :               1,0,ISTAT)
          CALL UMSPUT('! as cps NOT counts! Assuming counts',
     :               1,0,ISTAT)
         ENDIF
      ENDIF

C Interpret the coefficients file name and act appropriately
      NOCO=.FALSE.
      IF(COEFFS.EQ.' ') NOCO=.TRUE.

C Now check for the special string "header"
      IF(COEFFS.EQ.'header') THEN
         HCOEFF=.TRUE.

C and try to read these from the input header
         CALL GTCOEF(IDD,COEFFS,LAM,ISTAT)
         IF(ISTAT.NE.0) THEN
            CALL UMSPUT(
     :   '! Unable to get necessary header information',
     :      1,0,ISTAT)
            GO TO 99
         ENDIF
      ELSE
         HCOEFF=.FALSE.
      ENDIF

C if there is no coeffients file then set the coeffs to defaults 
      IF(NOCO) THEN
         DO I=1,10
            XCO(I)=0.0
            YCO(I)=0.0
         ENDDO
      ELSE

C Open the coefficients file
         CALL UFOPEN(COEFFS,1,LUN,ISTAT)
         IF(ISTAT.NE.0) THEN
            CALL UMSPUT('! Unable to open coefficients file',1,0,
     :                  ISTAT)
            GO TO 99
         ELSE
            IF(VERBOSE) CALL UMSPUT(
     :         '-Opening coefficients file: '//COEFFS,
     :                      1,0,ISTAT)
         ENDIF

C Get the wavelength (nm). This is only used for "trauger"
C style coefficients sets, if we are using the header we
C will already have this
         IF(.NOT.HCOEFF) CALL UCLGSR('lambda',LAM,ISTAT)

C Read in the coefficients from the file
         CALL GETCO(LUN,LAM,XCO,YCO,ISTAT)
         IF(ISTAT.NE.0) THEN
            CALL UMSPUT('! Unable to read coefficients',
     :                       1,0,ISTAT)
            CALL UFCLOS(LUN,ISTAT)
            GO TO 99
         ENDIF

C Close the coefficients file
         CALL UFCLOS(LUN,ISTAT)
      ENDIF

C We get the name of the output image here
      CALL UCLGST('outdata',OUTDAT,ISTAT)

C Check that a name is specified
      IF(OUTDAT.EQ.' ') THEN
         CALL UMSPUT(
     :       '! No output data array name specified',
     :                   1,0,ISTAT)
           GO TO 99
      ENDIF

C We have to create the output data image 
C so we need to find out their dimension from the CL
      CALL UCLGSI('outnx',ONX,ISTAT)
      CALL UCLGSI('outny',ONY,ISTAT)

      ODIMS(1)=ONX
      ODIMS(2)=ONY

C Now we need to get the name for the output image
      CALL UCLGST('outdata',OUTDAT,ISTAT)
      CALL UIMCRE(OUTDAT,6,2,ODIMS,IDND,ISTAT)
      IF(ISTAT.NE.0) THEN
         CALL UMSPUT(
     :     '! Unable to create output data image',
     :                 1,0,ISTAT)
         GO TO 99
      ENDIF

C Now we allocate space for all the arrays we will need, note that they
C are all real
      CALL UDMGET(DNX*DNY,6,PDATA,ISTAT)
      IF(ISTAT.NE.0) THEN
         CALL UMSPUT(
     :          '! Unable to allocate memory for data array',
     :                 1,0,ISTAT)
         GO TO 99
      ENDIF

      CALL UDMGET(ONX*ONY,6,PNDAT,ISTAT)
      IF(ISTAT.NE.0) THEN
         CALL UMSPUT(
     : '! Unable to allocate memory for output data array',
     :                 1,0,ISTAT)
         GO TO 99
      ENDIF

C Now we can actually read the data into the arrays
C Note that we read a line at a time to try to avoid
C excessive memory use
      DO I=1,DNY
         CALL UIGS2R(IDD,1,DNX,I,I,MEMR(PDATA+(I-1)*DNX),ISTAT)
         IF(ISTAT.NE.0) THEN
            CALL UMSPUT(
     :         '! Failed to read input data image',
     :                      1,0,ISTAT)
            GO TO 99
         ENDIF
      ENDDO

C Before we close the data array we copy the header to
C the output, this also transfers the exposure time information
      CALL UHDCPY(IDD,IDND,ISTAT)

C Also try to get the header items for the WCS, note that
C this is not the full WCS but just 8 header items commonly used
C for HST images plus CTYPE1 & CTYPE2 
      CALL GETWCS(IDD,WCS,CTYPE1,CTYPE2,ISTAT)
      IF(ISTAT.EQ.0) THEN
         GOTWCS=.TRUE.
      ELSE
         CALL UMSPUT(
     :    '! Warning, unable to read WCS information from header',
     :    1,0,ISTAT)
      ENDIF

C Here we set the WCS matrix to something sensible anyway
      IF(.NOT.GOTWCS) THEN
         CTYPE1='PIXELS'
         CTYPE2='PIXELS'

         DO I=1,8
            WCS(I)=0.0D0
         ENDDO
      ENDIF

C Get the exposure time from the header, if possible
      CALL UHDGSR(IDD,EXPKEY,EXPIN,ISTAT)
      IF(ISTAT.NE.0) THEN
         CALL UMSPUT(
     : '! Warning, failed to get exposure time keyword value',
     :      1,0,ISTAT)
            GO TO 99
      ENDIF

C Check for the "input" value for the output exposure,
C this is the standard behaviour
      IF(EXPSTR.EQ.'INPUT' .OR. EXPSTR.EQ.'input') THEN
         EXPOUT=EXPIN
      ELSE
         READ(EXPSTR,*,IOSTAT=ISTAT) EXPOUT
         IF(ISTAT.NE.0) THEN
            CALL UMSPUT('! Invalid output exposure value',
     :                  1,0,ISTAT)
            GO TO 99
         ENDIF
      ENDIF

C How we re-scale depends on the units, there are four cases
      IF(INCPS) THEN
         IF(OUTCPS) THEN
            EF=1.0
         ELSE
            EF=EXPOUT
         ENDIF
      ELSE
         IF(OUTCPS) THEN
            EF=1.0/EXPIN
         ELSE
            EF=EXPOUT/EXPIN
         ENDIF
      ENDIF

      CALL UHDPSR(IDND,EXPKEY,EXPOUT,ISTAT)
      WRITE(CHARS,
     :    '(''- Input exposure '',F10.3,'' Output exposure: '''//
     :    'F10.3)') EXPIN,EXPOUT
      CALL UMSPUT(CHARS,1,0,ISTAT)

C Close the data array now we have it
      CALL UIMCLO(IDD,ISTAT)
 
C The final case is where the output is empty.
C Before we start blotting we need to drain out the output arrays
C if there are no old arrays 
      CALL SETIM(MEMR(PNDAT),ONX,ONY,0.0)

C Before we start also get the interpolation type and convert
C to a numerical code
      CALL UCLGST('interpol',INTERP,ISTAT)
      IF(INTERP.EQ.'nearest') THEN
         INTYPE=1
      ELSE IF(INTERP.EQ.'linear') THEN
         INTYPE=2
      ELSE IF(INTERP.EQ.'poly3') THEN
         INTYPE=3
      ELSE IF(INTERP.EQ.'poly5') THEN
         INTYPE=4
      ELSE IF(INTERP.EQ.'spline3') THEN
         INTYPE=5
      ELSE 
         CALL UMSPUT('! Invalid interpolant specified',
     :               1,0,ISTAT)
         GO TO 99
      ENDIF  

C Now do the actual combination using interpolation
      CALL DOBLOT(MEMR(PDATA),MEMR(PNDAT),INTYPE,MISVAL,
     :         DNX,DNY,ROTFIR,EF,ALIGN,
     :         ONX,ONY,XCO,YCO,SCALE,ROT,XSH,YSH,WCS)

C Write out the new image
      IF(VERBOSE)
     :   CALL UMSPUT('-Writing output blotted image: '//
     :               OUTDAT,1,0,ISTAT)
         
      CALL UIPS2R(IDND,1,ONX,1,ONY,MEMR(PNDAT),ISTAT)

C Update the header
      CALL BUHEAD(IDND,VERS,DATA,OUTDAT,COEFFS,SHFTUN,MISVAL,
     :     EXPKEY,EXPSTR,ALIGN,INUN,OUTUN,
     :     SHFTFR,INTERP,LAM,SCALE,ROT,XSH,YSH,ONX,ONY)

      IF(GOTWCS) THEN
         CALL UWCS(IDND,WCS,CTYPE1,CTYPE2)
      ENDIF

      CALL UIMCLO(IDND,ISTAT)
      CALL UDMFRE(PNDAT,6,ISTAT)

C Close down everything and exit
      CALL UDMFRE(PDATA,6,ISTAT)

 99   CONTINUE
      END

      SUBROUTINE BUHEAD(ID,VERS,DATA,OUTDAT,COEFFS,SHFTUN,MISVAL,
     :      EXPKEY,EXPSTR,ALIGN,INUN,OUTUN,
     :      SHFTFR,INTERP,LAM,SCALE,ROT,XSH,YSH,OUTNX,OUTNY)
C
C Update the header of the output image with all the
C parameters of the current BLOT run.
C
      IMPLICIT NONE

      INTEGER ID,OUTNX,OUTNY,ISTAT
      CHARACTER*80 DATA
      CHARACTER*80 OUTDAT,COEFFS
      CHARACTER*7 INTERP
      CHARACTER*8 SHFTFR,SHFTUN,EXPKEY,ALIGN,INUN,OUTUN
      CHARACTER*80 EXPSTR
      CHARACTER*(*) VERS
      REAL LAM,SCALE,ROT
      REAL DROT,XSH,YSH,XSHI,YSHI,MISVAL

      CALL UHDAST(ID,'BLOTVER',VERS,
     :            'Blot, task version',0,ISTAT)
      
      CALL UHDAST(ID,'BLOTDATA',DATA,
     :            'Blot, input data image',0,ISTAT)

      CALL UHDAST(ID,'BLOTOUTDA',OUTDAT,
     :            'Blot, output data image',0,ISTAT)

      CALL UHDAST(ID,'BLOTCOEFF',COEFFS,
     :            'Blot, coefficients file name ',0,ISTAT)

      CALL UHDASR(ID,'BLOTLAM',LAM,
     : 'Blot, wavelength applied for transformation (nm)',
     :  0,ISTAT)

      CALL UHDASR(ID,'BLOTSCALE',SCALE,
     : 'Blot, scale (pixel size) of input image',0,ISTAT)

C Convert the rotation angle back to degrees
      DROT=ROT/3.1415926*180.0
      CALL UHDASR(ID,'BLOTROT',DROT,
     :  'Blot, rotation angle, degrees clockwise',0,ISTAT)

C Check the SCALE and units
      IF(SHFTUN.EQ.'input') THEN
         XSHI=XSH/SCALE
         YSHI=YSH/SCALE
      ELSE
         XSHI=XSH
         YSHI=YSH
      ENDIF

      CALL UHDASR(ID,'BLOTXSH',XSHI,
     :  'Blot, X shift applied',0,ISTAT)

      CALL UHDASR(ID,'BLOTYSH',YSHI,
     :  'Blot, Y shift applied',0,ISTAT)

      CALL UHDAST(ID,'BLOTSHFU',SHFTUN,
     : 'Blot, units used for shifts (output or input pixels)',
     :            0,ISTAT)

      CALL UHDAST(ID,'BLOTINUN',INUN,
     : 'Blot, units of input image (counts or cps)',
     :            0,ISTAT)

      CALL UHDAST(ID,'BLOTOUUN',OUTUN,
     : 'Blot, units for output image (counts or cps)',
     :            0,ISTAT)

      CALL UHDAST(ID,'BLOTSHFF',SHFTFR,
     : 'Blot, frame in which shifts were applied',0,ISTAT)

      CALL UHDASI(ID,'BLOTONX',OUTNX,
     :  'Blot, X size of output image',0,ISTAT)

      CALL UHDASI(ID,'BLOTONY',OUTNY,
     :  'Blot, Y size of output image',0,ISTAT)

      CALL UHDAST(ID,'BLOTINT',INTERP,
     :  'Blot, interpolation method used',0,ISTAT)

      CALL UHDASR(ID,'BLOTFILVA',MISVAL,
     :  'Blot, value assigned to invalid pixels',0,ISTAT)

      CALL UHDAST(ID,'BLOTEXKY',EXPKEY,
     :  'Blot, exposure time keyword in input header',0,ISTAT)

      CALL UHDAST(ID,'BLOTEXOU',EXPSTR,
     :  'Blot, exposure time for output image',0,ISTAT)

      RETURN
      END

      SUBROUTINE DOBLOT(DATA,NDAT,INTYPE,MISVAL,DNX,DNY,
     :     ROTFIR,EF,ALIGN,ONX,ONY,XCO,YCO,SCALE,ROT,XSH,YSH,WCS)
C
C This routine does the interpolation of the input array
C
      IMPLICIT NONE

      INTEGER DNX,DNY,ONX,ONY
      INTEGER I,J,NMISS,ISTAT
      INTEGER INTYPE
      REAL DATA(DNX,DNY),V
      REAL NDAT(ONX,ONY)
      REAL XCO(10),YCO(10)
      DOUBLE PRECISION WCS(8)
      REAL X,Y,XOUT,YOUT,XSH,YSH,XIN,YIN,XV,YV
      REAL EVALCU,SINTH,COSTH,ROT,XF,YF,XOFF,YOFF,SCALE
      REAL XS,XC,YS,YC,XT,YT,S2,XDC,YDC
      REAL XCEN,YCEN,XP,YP,MISVAL
      REAL MRIEVL,EF,A,B,C,D
      CHARACTER*80 CHARS
      CHARACTER*8 ALIGN
      LOGICAL NONLIN
      LOGICAL ROTFIR
C--

C Some initial settings - note that the reference pixel
C is offset from the "centre":
C
      NMISS=0

      IF(ALIGN.EQ.'corner') THEN
         XCEN=FLOAT(ONX/2)+0.5
         YCEN=FLOAT(ONY/2)+0.5
      ELSE
         XCEN=FLOAT(ONX/2)+1.0
         YCEN=FLOAT(ONY/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 the position of the first whole pixel
      Y=-YCEN

C Calculate some numbers to simplify things later
      SINTH=SIN(ROT)
      COSTH=COS(ROT)

C Note that XSH is in INPUT pixels and XOFF etc are in 
C OUTPUT pixels
      XOFF=XSH/SCALE
      YOFF=YSH/SCALE
      XF=1.0/SCALE
      YF=1.0/SCALE

C Some more economies
      S2=1.0/(XF*YF)
      XS=XF*SINTH
      XC=XF*COSTH
      YS=YF*SINTH
      YC=YF*COSTH

      IF(ALIGN.EQ.'corner') THEN
         XDC=REAL(DNX/2)+0.5
         YDC=REAL(DNY/2)+0.5
      ELSE
         XDC=REAL(DNX/2)+1.0
         YDC=REAL(DNY/2)+1.0
      ENDIF

      XT=XOFF+XDC
      YT=YOFF+YDC

C Outer loop over output image pixels (X,Y)
      DO J=1,ONY
       Y=Y+1.0

       X=-XCEN
       DO I=1,ONX
          X=X+1.0

C Calculate where this pixel lands up
C The angles are assumed to be in radians at this point

          IF(NONLIN) THEN
             XP=EVALCU(X,Y,XCO)
             YP=EVALCU(X,Y,YCO)
          ELSE
             XP=X
             YP=Y
          ENDIF

C Apply the linear transform
C Here we have to take care about the order of
C operations
          IF(ROTFIR) THEN
             XOUT=XC*XP-YS*YP+XT
             YOUT=XS*XP+YC*YP+YT
          ELSE
             XOUT=XC*(XP+XSH)-YS*(YP+YSH)+XDC
             YOUT=XS*(XP+XSH)+YC*(YP+YSH)+YDC
          ENDIF

C Check it is on the input image
         IF(XOUT.GE.1.0 .AND. XOUT.LE.REAL(DNX) .AND.
     :      YOUT.GE.1.0 .AND. YOUT.LE.REAL(DNY)) THEN

C Do the interpolation
C and allow for stretching because of scale change
            V=MRIEVL(XOUT,YOUT,DATA,DNX,DNY,DNX,INTYPE)

            NDAT(I,J)=V*EF/S2
         ELSE

C If there is nothing for us then set the output to missing
C value flag
            NDAT(I,J)=MISVAL

C Count cases where the pixel is off the output image
            NMISS=NMISS+1
         ENDIF
        ENDDO
      ENDDO

C Report number of points which were outside the input frame
      IF(NMISS.GT.0) THEN
       WRITE(CHARS,'(''! Warning, '',I7,'' points were outside '','//
     :                '''the output image.'')') NMISS
         CALL UMSPUT(CHARS,1,0,ISTAT)
      ENDIF

C Project the central position in the output array and make
C it the reference pixel for the output WCS
      IF(NONLIN) THEN
         XP=XCO(1)
         YP=YCO(1)
      ELSE
         XP=0.0         
         YP=0.0        
      ENDIF

C Now the centre position
      IF(ROTFIR) THEN
         XIN=DBLE(XC*XP-YS*YP+XT)
         YIN=DBLE(XS*XP+YC*YP+YT)
      ELSE
         XIN=DBLE(XC*(XP+XSH)-YS*(YP+YSH)+XDC)
         YIN=DBLE(XS*(XP+XSH)+YC*(YP+YSH)+YDC)
      ENDIF

C Now we convert the pixel position to a sky position assuming
C a small field and using the old WCS
      XV=((XIN-REAL(WCS(1)))*REAL(WCS(5))+
     :   (YIN-REAL(WCS(3)))*REAL(WCS(7)))/
     :   COS(REAL(WCS(4))*3.1415926/180.0)+REAL(WCS(2))
      YV=(XIN-REAL(WCS(1)))*REAL(WCS(6))+
     :   (YIN-REAL(WCS(3)))*REAL(WCS(8))+REAL(WCS(4))

      WCS(1)=DBLE(XCEN)
      WCS(2)=DBLE(XV)
      WCS(3)=DBLE(YCEN)
      WCS(4)=DBLE(YV)

C Now we scale and rotate the CD matrix (WCS(5-8))
      A=COSTH*WCS(5)+SINTH*WCS(7)
      B=COSTH*WCS(6)+SINTH*WCS(8)
      C=-SINTH*WCS(5)+COSTH*WCS(7)
      D=-SINTH*WCS(6)+COSTH*WCS(8)

      WCS(5)=DBLE(A*XF)
      WCS(6)=DBLE(B*XF)
      WCS(7)=DBLE(C*XF)
      WCS(8)=DBLE(D*XF)

      RETURN
      END