SUBROUTINE TRANBK 
C++
C
C TRANBACK.F V0.3 Convert X,Y positions in drizzle output frames
C                 back to input positions in the initial frames
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 Where <dir> is the location of the STSDAS f77/vos library.
C
C History:
C
C First quick version (V0.1), 7th May 1999 (Richard Hook, STScI)
C Bug fixes and tests for release of V0.2 in the IRAF stecf package
C 17th December 1999, Richard Hook
C
C Bugfix for different calling sequences, Richard Hook, 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,XX,YY
      REAL SINTH,COSTH,XF,YF,XS,YC,XC,YS,XP,YP
      REAL XCORN,YCORN,ERR
      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='TRANBACK Version 0.2 (17th December 1999)'

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

C Acceptable position error for iterative search
      ERR=0.001
      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('xout',XOUT,ISTAT)
      CALL UCLGSR('yout',YOUT,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)

C Also get the chip size - in and out
      CALL UCLGSI('nxout',NXOUT,ISTAT)
      CALL UCLGSI('nyout',NYOUT,ISTAT)
      CALL UCLGSI('nxin',NXIN,ISTAT)
      CALL UCLGSI('nyin',NYIN,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.
      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

C Transform the position
C Note that the order is the opposite of that in Drizzle

C Offset from centre
      IF(ROTFIR) THEN
	 IF(SHFTUN.EQ.'input') THEN
	    XSH=XSH/SCALE
	    YSH=YSH/SCALE
         ENDIF

         XX=XOUT-XP-XSH
         YY=YOUT-YP-YSH
         XCORN=(XX*YC+YY*YS)/(XC*YC+XS*YS)
         YCORN=(YY*XC-XX*XS)/(XC*YC+XS*YS)
      ELSE
	 IF(SHFTUN.EQ.'output') THEN
	   XSH=XSH*SCALE
	   YSH=YSH*SCALE
	 ENDIF

         XX=XOUT-XP
         YY=YOUT-YP
         XCORN=(XX*YC+YY*YS)/(XC*YC+XS*YS)-XSH
         YCORN=(YY*XC-XX*XS)/(XC*YC+XS*YS)-YSH
      ENDIF

C Now we have to do the geometric transform, backwards
      IF(.NOT.NOCO) THEN
        CALL INVECU(XCORN,YCORN,XCO,YCO,ERR,XX,YY) 
      ELSE
	XX=XCORN
	YY=YCORN
      ENDIF

C Offset from the centre of the input image
      XIN=XX+XCEN
      YIN=YY+YCEN

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 TRANBACK module
 99   CONTINUE
      RETURN
      END