c---------------------------------------
c
c developed by Jay at home Xmas 2009
c
      program wfc3ir_gc_2fits
      implicit none

      integer i,j
      real*8 xr,yr
      real*8 xc,yc

      real ximg(1014,1014)
      real yimg(1014,1014)

      do i = 0001, 1014
      do j = 0001, 1014
         xr = i
         yr = j
         call wfc3ir_gcX(xr,yr,xc,yc)
         ximg(i,j) = xc
         yimg(i,j) = yc
         enddo
         enddo

      call writfits_r4('wfc3ir_gcx.fits',ximg,1014,1014)
      call writfits_r4('wfc3ir_gcy.fits',yimg,1014,1014)

      return
      end

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

#define _LINUX_ .true.
#include "/home/jayander/FORTRAN/MAC_OSX/FITS/writfits_r4.f"

      subroutine wfc3ir_gcX(xr,yr,xc,yc)
      implicit none

      real*8 xr,yr
      real*8 xc,yc

      real xru,yru

      real x, y
      integer kinu

      real dxcor, dycor

      real polyx(10)
      common /IRpolyx_/polyx
      data polyx / 
     .  507.000, 567.453,   1.767, -0.065,  7.171, 
     .    0.033,  -0.002,   0.007,  0.003, -0.012/

      real polyy(10)
      common /IRpolyxyx_/polyy
      data polyy / 
     .  507.000,   0.000, 507.000,  1.805, -0.049,
     .    7.866,   0.015,  -0.022,  0.032, -0.012/

      real xcor(11,11)
      common /IRxcor_/xcor
      data xcor /
     .0.1070,0.0915,0.0738,0.0580,0.0333,0.0178,
     .       0.0032,-.0196,-.0343,-.0387,-.0326,
     .0.0815,0.0691,0.0547,0.0387,0.0144,0.0040,
     .       -.0156,-.0389,-.0477,-.0479,-.0376,
     .0.0589,0.0498,0.0386,0.0244,0.0032,-.0044,
     .       -.0260,-.0462,-.0508,-.0468,-.0323,
     .0.0343,0.0328,0.0285,0.0243,0.0081,-.0089,
     .       -.0303,-.0391,-.0442,-.0346,-.0142,
     .0.0146,0.0084,0.0065,0.0175,0.0076,-.0127,
     .       -.0350,-.0348,-.0243,-.0165,-.0041,
     .0.0039,-.0020,-.0027,0.0071,0.0017,-.0049,
     .       -.0175,-.0209,-.0111,-.0005,0.0163,
     .-.0204,-.0184,-.0137,-.0036,-.0070,-.0009,
     .       -.0016,-.0085,-.0012,0.0149,0.0396,
     .-.0370,-.0347,-.0274,-.0081,-.0016,0.0032,
     .       0.0068,0.0141,0.0235,0.0401,0.0616,
     .-.0470,-.0455,-.0390,-.0224,-.0071,0.0043,
     .       0.0078,0.0223,0.0412,0.0543,0.0686,
     .-.0491,-.0480,-.0419,-.0229,-.0067,0.0074,
     .       0.0165,0.0343,0.0514,0.0656,0.0810,
     .-.0452,-.0444,-.0387,-.0196,-.0055,0.0127,
     .       0.0324,0.0497,0.0583,0.0736,0.0901/

      real ycor(11,11)
      common /IRycor_/ycor
      data ycor /
     .0.1085,0.0819,0.0633,0.0621,0.0369,0.0438,
     .       0.0580,0.0933,0.1052,0.1388,0.1840,
     .0.0703,0.0479,0.0336,0.0389,0.0233,0.0340,
     .       0.0343,0.0532,0.0650,0.0919,0.1305,
     .0.0360,0.0178,0.0077,0.0158,0.0050,0.0172,
     .       0.0097,0.0176,0.0252,0.0455,0.0773,
     .0.0100,-.0094,-.0184,-.0146,-.0162,-.0056,
     .       -.0190,-.0171,-.0240,-.0039,0.0324,
     .0.0013,-.0239,-.0373,-.0308,-.0363,-.0289,
     .       -.0336,-.0329,-.0470,-.0272,0.0105,
     .0.0193,-.0198,-.0428,-.0377,-.0466,-.0465,
     .       -.0521,-.0507,-.0556,-.0353,0.0008,
     .0.0367,-.0045,-.0324,-.0496,-.0577,-.0528,
     .       -.0624,-.0599,-.0466,-.0249,0.0115,
     .0.0634,0.0252,-.0033,-.0346,-.0489,-.0462,
     .       -.0542,-.0497,-.0351,-.0114,0.0275,
     .0.0930,0.0612,0.0357,0.0086,-.0063,-.0081,
     .       -.0110,-.0046,0.0092,0.0289,0.0584,
     .0.1118,0.0855,0.0653,0.0459,0.0339,0.0313,
     .       0.0298,0.0350,0.0435,0.0594,0.0851,
     .0.1303,0.1093,0.0945,0.0889,0.0813,0.0766,
     .       0.0771,0.0781,0.0802,0.0923,0.1143/


      real    rx, ry
      integer ix, iy
      real    fx, fy

      xc = xr
      yc = yr

      xru = xr
      yru = yr

      x = (xru-507.00)/507.0
      y = (yru-507.00)/507.0
      xc = polyx(01) 
     .   + polyx(02)*x
     .   + polyx(03)*y
     .   + polyx(04)*x*x
     .   + polyx(05)*x*y
     .   + polyx(06)*y*y
     .   + polyx(07)*x*x*x
     .   + polyx(08)*x*x*y
     .   + polyx(09)*x*y*y
     .   + polyx(10)*y*y*y
      yc = polyy(01) 
     .   + polyy(02)*x
     .   + polyy(03)*y
     .   + polyy(04)*x*x
     .   + polyy(05)*x*y
     .   + polyy(06)*y*y
     .   + polyy(07)*x*x*x
     .   + polyy(08)*x*x*y
     .   + polyy(09)*x*y*y
     .   + polyy(10)*y*y*y


      rx = 1 + xr/100
      ry = 1 + yr/100
      ix = int(rx)
      iy = int(ry)
      if (ix.lt.01) ix = 01
      if (iy.lt.01) iy = 01
      if (ix.gt.10) ix = 10
      if (iy.gt.10) iy = 10
      fx = rx-ix
      fy = ry-iy

      dxcor = (1-fx)*(1-fy)*xcor(ix  ,iy  )
     .      + (1-fx)*( fy )*xcor(ix  ,iy+1) 
     .      + ( fx )*(1-fy)*xcor(ix+1,iy  ) 
     .      + ( fx )*( fy )*xcor(ix+1,iy+1) 

      dycor = (1-fx)*(1-fy)*ycor(ix  ,iy  )
     .      + (1-fx)*( fy )*ycor(ix  ,iy+1) 
     .      + ( fx )*(1-fy)*ycor(ix+1,iy  ) 
     .      + ( fx )*( fy )*ycor(ix+1,iy+1) 

      xc = xc - dxcor
      yc = yc - dycor
  
      return 
      end