c---------------------------------------------------------------------
c
c Jay Anderson (jay@astron.berkeley.edu)
c
c Here is a FORTRAN routine that implemenents the distortion solution
c I dervived from the inner Omega Cen calibration field, about 80 F555W
c observations of 100s.  The paper that describes the derivation and 
c use of this solution should appear in the Jan 2003 version of the PASP.
c
c
c This routine takes in an (xin,yin) position as measured on a chip
c then returns the distortion-corrected position (xout,yout); if the
c input chip is 0, no correction is performed.  k=1 corresponds to
c pc1, k=2 to wf2, etc.
c
      subroutine wfpc2gcK(xin,yin,xout,yout,k)
      implicit none
      real xin,  yin
      real xout, yout
      integer k
      real x_wfpc2gcK
      real y_wfpc2gcK

      real xinh, yinh

      xinh = xin ! doing it this way allows the calling
      yinh = yin !    routine to use the same variable for xin/xout

      xout = x_wfpc2gcK(xinh,yinh,k) 
      yout = y_wfpc2gcK(xinh,yinh,k) 

      return
      end

c--------------------------------------------------------

      real function x_wfpc2gcK(xin,yin,k)
      implicit none
      real xin,  yin
      integer k
      real x, y, d

      real pwfpc2gc(2,4,10)
      data pwfpc2gc /                                           ! ORDER
     .  0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, ! const
     .  0.000, 0.418, 0.000, 0.051, 0.000,-0.028, 0.000, 0.070, ! x
     .  0.000,-0.016, 0.000,-0.015, 0.000,-0.036, 0.000, 0.059, ! y
     . -0.525,-0.280,-0.624,-0.038,-0.349,-0.027,-0.489,-0.050, ! xx
     . -0.268,-0.292,-0.411,-0.568,-0.353,-0.423,-0.391,-0.485, ! xy
     . -0.249,-0.470,-0.092,-0.444, 0.009,-0.373,-0.066,-0.406, ! yy
     . -1.902,-0.011,-1.762, 0.003,-1.791, 0.004,-1.821,-0.015, ! xxx
     .  0.024,-1.907, 0.016,-1.832, 0.006,-1.848, 0.022,-1.890, ! xxy
     . -1.890, 0.022,-1.825, 0.011,-1.841, 0.006,-1.875, 0.022, ! xyy
     . -0.004,-1.923, 0.010,-1.730, 0.021,-1.788,-0.006,-1.821/ ! yyy
c
c        APC1  BPC1    AWF2   BWF2   AWF3   BWF3   AWF4   BWF4  
c

      x_wfpc2gcK = xin
      if (k.lt.1.or.k.gt.4) return
      if (xin.lt.000) return
      if (yin.lt.000) return
      if (xin.gt.800) return
      if (yin.gt.800) return

      x = (xin-425.0)/375.0
      y = (yin-425.0)/375.0

c
c     for x use the A coeffs...
c 
      d = pwfpc2gc(1,k,01)*1.0 +
     .    pwfpc2gc(1,k,02)*x  +
     .    pwfpc2gc(1,k,03)*y  +
     .    pwfpc2gc(1,k,04)*x**2 +
     .    pwfpc2gc(1,k,05)*x*y  +
     .    pwfpc2gc(1,k,06)*y**2 +
     .    pwfpc2gc(1,k,07)*x**3   +
     .    pwfpc2gc(1,k,08)*x**2*y +
     .    pwfpc2gc(1,k,09)*y**2*x +
     .    pwfpc2gc(1,k,10)*y**3   

      x_wfpc2gcK = xin + d 

      return
      end 

c--------------------------------------------------------

      real function y_wfpc2gcK(xin,yin,k)
      implicit none
      real xin,  yin
      integer k
      real x, y, d

      real pwfpc2gc(2,4,10)
      data pwfpc2gc /                                           ! ORDER
     .  0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, ! const
     .  0.000, 0.418, 0.000, 0.051, 0.000,-0.028, 0.000, 0.070, ! x
     .  0.000,-0.016, 0.000,-0.015, 0.000,-0.036, 0.000, 0.059, ! y
     . -0.525,-0.280,-0.624,-0.038,-0.349,-0.027,-0.489,-0.050, ! xx
     . -0.268,-0.292,-0.411,-0.568,-0.353,-0.423,-0.391,-0.485, ! xy
     . -0.249,-0.470,-0.092,-0.444, 0.009,-0.373,-0.066,-0.406, ! yy
     . -1.902,-0.011,-1.762, 0.003,-1.791, 0.004,-1.821,-0.015, ! xxx
     .  0.024,-1.907, 0.016,-1.832, 0.006,-1.848, 0.022,-1.890, ! xxy
     . -1.890, 0.022,-1.825, 0.011,-1.841, 0.006,-1.875, 0.022, ! xyy
     . -0.004,-1.923, 0.010,-1.730, 0.021,-1.788,-0.006,-1.821/ ! yyy
c
c         APC1  BPC1   AWF2   BWF2   AWF3   BWF3   AWF4   BWF4  
c

      y_wfpc2gck = yin
      if (k.lt.1.or.k.gt.4) return
      if (xin.lt.000) return
      if (yin.lt.000) return
      if (xin.gt.800) return
      if (yin.gt.800) return

      x = (xin-425.0)/375.0
      y = (yin-425.0)/375.0

c
c     for y use the B coeffs...
c 
      d = pwfpc2gc(2,k,01)*1.0  +
     .    pwfpc2gc(2,k,02)*x    +
     .    pwfpc2gc(2,k,03)*y    +
     .    pwfpc2gc(2,k,04)*x**2 +
     .    pwfpc2gc(2,k,05)*x*y  +
     .    pwfpc2gc(2,k,06)*y**2 +
     .    pwfpc2gc(2,k,07)*x**3   +
     .    pwfpc2gc(2,k,08)*x**2*y +
     .    pwfpc2gc(2,k,09)*y**2*x +
     .    pwfpc2gc(2,k,10)*y**3

      y_wfpc2gck = yin + d 

      return
      end