SUBROUTINE TRAXY  
C++
C
C TRAXY.F V0.5 Convert X,Y pixel position to drizzle output position
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 tranback.f drutil.f x_tranback.x <dir>libiraf77.a
C
C History:
C
C First quick version, 17th April 1997
C
C V0.2 modified to have the centre position as parameters, 4th Sept 97
C
C V0.3 added a rotation around the centre, 19th August 1998, STScI
C      also specified "align" and chip size.
C
C V0.31 possibility to have coeffs="" for simple linear case, 25th August 98
C
C V0.4 with full set of drizzle parameter for compatibility with
C      TRANBACK and release of STECF IRAF package. 17th December 1999
C
C V0.5 modified to use new calling sequences for some DRUTIL routines
C      and hence fix bugs, November 2001
C
C--
      IMPLICIT NONE

C Local variables
      INTEGER ISTAT,LUN,NXIN,NYIN,NXOUT,NYOUT,I
      INTEGER COTY,CONUM,COMAX
      PARAMETER (COMAX=100)

      REAL XCO(10),YCO(10),LAM,SCALE,XSH,YSH
      REAL XIN,YIN,YOUT,XOUT,XCEN,YCEN,ROT
      REAL SINTH,COSTH,XOFF,YOFF,XF,YF,XS,YC,XC,YS,XP,YP
      REAL XT,YT,EVALCU,EX,EY
      CHARACTER*80 COEFFS
      CHARACTER*132 CHARS
      CHARACTER*8 ALIGN,SHFTFR,SHFTUN
      CHARACTER*45 VERS

      REAL PIBY
      PARAMETER (PIBY = 3.141592653/180.0)

      LOGICAL VERBOSE
      LOGICAL NOCO
      LOGICAL ROTFIR
      LOGICAL NONLIN

C-- Start of executable code
      VERS='TRAXY Version 0.5 (8th November 2001)'

C First announce the version
      CALL UMSPUT('+ '//VERS,1,0,ISTAT)

      VERBOSE=.TRUE.

C Get the name of the file containing the coefficients
      CALL UCLGST('coeffs',COEFFS,ISTAT)
      IF(COEFFS.EQ.' ') THEN
         NOCO=.TRUE.
      ELSE
         NOCO=.FALSE.
      ENDIF

C Get the X,Y pixel position to be transformed 
      CALL UCLGSR('xin',XIN,ISTAT)
      CALL UCLGSR('yin',YIN,ISTAT)

C Open the coefficients file
      IF(.NOT.NOCO) THEN
         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)

C Get the wavelength (nm). This is only used for "trauger"
C style coefficients sets, 
            CALL UCLGSR('lambda',LAM,ISTAT)

C Read in the coefficients from the file
            CALL GETCO(LUN,LAM,COTY,COMAX,CONUM,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
      ENDIF

C Get the various linear transformation parameters
      CALL UCLGSR('scale',SCALE,ISTAT)
      CALL UCLGSR('xsh',XSH,ISTAT)
      CALL UCLGSR('ysh',YSH,ISTAT)
      CALL UCLGSR('rot',ROT,ISTAT)

C Get the origin of the image coordinates - in other words
C whether we treat the center or the corner of the central
C pixel as the reference point
      CALL UCLGST('align',ALIGN,ISTAT)

C Get parameter specifying which frames shifts are applied
      CALL UCLGST('shft_fr',SHFTFR,ISTAT)
      IF(SHFTFR.EQ.'input') THEN
         ROTFIR=.FALSE.
      ELSE
         ROTFIR=.TRUE.
      ENDIF

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

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

C Also get the chip size - in and out
      CALL UCLGSI('nxin',NXIN,ISTAT)
      CALL UCLGSI('nyin',NYIN,ISTAT)
      CALL UCLGSI('nxout',NXOUT,ISTAT)
      CALL UCLGSI('nyout',NYOUT,ISTAT)

C Calculate the centre position
      IF(ALIGN.EQ.'corner') THEN
         XCEN=FLOAT(NXIN/2)+0.5
         YCEN=FLOAT(NYIN/2)+0.5
      ELSE
         XCEN=FLOAT(NXIN/2)+1.0
         YCEN=FLOAT(NYIN/2)+1.0
      ENDIF

C Convert to radians
      ROT=ROT*PIBY

C Setup offsets, scales etc in the same way as Drizzle
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(NXOUT/2)+0.5
         YP=REAL(NYOUT/2)+0.5
      ELSE
         XP=REAL(NXOUT/2)+1.0
         YP=REAL(NYOUT/2)+1.0
      ENDIF

      XT=XOFF+XP
      YT=YOFF+YP

C Transform the position
C Note that this is exactly the same as that in Drizzle
      IF(NOCO) THEN
        EX=XIN-XCEN
	EY=YIN-YCEN
      ELSE
        EX=EVALCU(XIN-XCEN,YIN-YCEN,XCO)
	EY=EVALCU(XIN-XCEN,YIN-YCEN,YCO)
      ENDIF

C Apply the linear part of the transformation
      IF(ROTFIR) THEN
         XOUT=XC*EX-YS*EY+XT
	 YOUT=XS*EX+YC*EY+YT
      ELSE
         XOUT=XC*(EX+XSH)-YS*(EY+YSH)+XP
         YOUT=XS*(EX+XSH)+YC*(EY+YSH)+YP
      ENDIF

C Write out the result
      WRITE(CHARS,
     :  '('' Xin,Yin: '',2F12.5,'' Xout,Yout: '',2F12.5)')
     : XIN,YIN,XOUT,YOUT
      CALL UMSPUT(CHARS,1,0,ISTAT)

C End of main TRAXY module
 99   CONTINUE
      RETURN
      END