code: def linearize(wcsim, wcsima, wcsref, imcrpix, f, shift, hx=1.0, hy=1.0):
      
      x0 = imcrpix[0]
      y0 = imcrpix[1]
      p = np.asarray([[x0, y0],
                      [x0 - hx, y0],
                      [x0 - hx * 0.5, y0],
                      [x0 + hx * 0.5, y0],
                      [x0 + hx, y0],
                      [x0, y0 - hy],
                      [x0, y0 - hy * 0.5],
                      [x0, y0 + hy * 0.5],
                      [x0, y0 + hy]],
                     dtype=np.float64)
      # convert image coordinates to reference image coordinates:
      p = wcsref.wcs_world2pix(wcsim.wcs_pix2world(p, 1), 1).astype(ndfloat128)
      # apply linear fit transformation:
      p = np.dot(f, (p - shift).T).T
      # convert back to image coordinate system:
      p = wcsima.wcs_world2pix(
          wcsref.wcs_pix2world(p.astype(np.float64), 1), 1).astype(ndfloat128)
  
      # derivative with regard to x:
      u1 = ((p[1] - p[4]) + 8 * (p[3] - p[2])) / (6*hx)
      # derivative with regard to y:
      u2 = ((p[5] - p[8]) + 8 * (p[7] - p[6])) / (6*hy)
  
      return (np.asarray([u1, u2]).T, p[0])
docs: linearization using 5-point formula for first order derivative