docs: def sky2xypv(ra, dec, crpix1, crpix2, crval1, crval2, dc, pv, nord, maxiter=300):
      """
      Transforms from celestial coordinates to pixel coordinates to taking
      non-linear distortion into account with the World Coordinate System
      FITS keywords as used in MegaPipe.
  
      For the inverse operation see xy2sky.
  
      Reference material:
      http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/megapipe/docs/CD_PV_keywords.pdf
  
      Args:
        ra: float
          Right ascension
        dec: float
          Declination
        crpix1: float
          Tangent point x, pixels
        crpix2: float
          Tangent point y, pixels
        crval1: float
          Tangent point RA, degrees
        crval2:
          Tangent point Dec, degrees
        dc: 2d array
          Expresses the scale, the rotation and any possible skew of the image
          with respect to the sky.  It is the inverse of the cd matrix in xy2sky.
        pv: 2d array
        nord: int
          order of the fit
  
      Returns:
        x, y: float
          Pixel coordinates
      """
      if math.fabs(ra - crval1) > 100:
          if crval1 < 180:
              ra -= 360
          else:
              ra += 360
  
      ra /= PI180
      dec /= PI180
  
      tdec = math.tan(dec)
      ra0 = crval1 / PI180
      dec0 = crval2 / PI180
      ctan = math.tan(dec0)
      ccos = math.cos(dec0)
  
      traoff = math.tan(ra - ra0)
      craoff = math.cos(ra - ra0)
      etar = (1 - ctan * craoff / tdec) / (ctan + craoff / tdec)
      xir = traoff * ccos * (1 - etar * ctan)
      xi = xir * PI180
      eta = etar * PI180
  
      if nord < 0:
          # The simple solution
          x = xi
          y = eta
      else:
          # Reverse by Newton's method
          tolerance = 0.001 / 3600
  
          # Initial guess
          x = xi
          y = eta
          iteration = 0
  
          converged = False
          while not converged:
              assert nord >= 0
  
              # Estimates
              f = pv[0][0]
              g = pv[1][0]
  
              # Derivatives
              fx = 0
              fy = 0
              gx = 0
              gy = 0
  
              if nord >= 1:
                  r = math.sqrt(x ** 2 + y ** 2)
                  f += pv[0][1] * x + pv[0][2] * y + pv[0][3] * r
                  g += pv[1][1] * y + pv[1][2] * x + pv[1][3] * r
                  fx += pv[0][1] + pv[0][3] * x / r
                  fy += pv[0][2] + pv[0][3] * y / r
                  gx += pv[1][2] + pv[1][3] * x / r
                  gy += pv[1][1] + pv[1][3] * y / r
  
              if nord >= 2:
                  x2 = x ** 2
                  xy = x * y
                  y2 = y ** 2
  
                  f += pv[0][4] * x2 + pv[0][5] * xy + pv[0][6] * y2
                  g += pv[1][4] * y2 + pv[1][5] * xy + pv[1][6] * x2
                  fx += pv[0][4] * 2 * x + pv[0][5] * y
                  fy += pv[0][5] * x + pv[0][6] * 2 * y
                  gx += pv[1][5] * y + pv[1][6] * 2 * x
                  gy += pv[1][4] * 2 * y + pv[1][5] * x
  
              if nord >= 3:
                  x3 = x ** 3
                  x2y = x2 * y
                  xy2 = x * y2
                  y3 = y ** 3
  
                  f += pv[0][7] * x3 + pv[0][8] * x2y + pv[0][9] * xy2 + pv[0][10] * y3
                  g += pv[1][7] * y3 + pv[1][8] * xy2 + pv[1][9] * x2y + pv[1][10] * x3
                  fx += pv[0][7] * 3 * x2 + pv[0][8] * 2 * xy + pv[0][9] * y2
                  fy += pv[0][8] * x2 + pv[0][9] * 2 * xy + pv[0][10] * 3 * y2
                  gx += pv[0][8] * y2 + pv[1][9] * 2 * xy + pv[1][10] * 3 * x2
                  gy += pv[1][7] * 3 * y2 + pv[0][8] * 2 * xy + pv[1][9] * x2
  
              f -= xi
              g -= eta
              dx = (-f * gy + g * fy) / (fx * gy - fy * gx)
              dy = (-g * fx + f * gx) / (fx * gy - fy * gx)
              x += dx
              y += dy
  
              if math.fabs(dx) < tolerance and math.fabs(dy) < tolerance:
                  converged = True
  
              iteration += 1
  
              if iteration > maxiter:
                  break
  
      xp = dc[0][0] * x + dc[0][1] * y
      yp = dc[1][0] * x + dc[1][1] * y
  
      x = xp + crpix1
      y = yp + crpix2
  
      return x, y