docs: def xy2skypv(x, y, crpix1, crpix2, crval1, crval2, cd, pv, nord):
      """
      Transforms from pixel coordinates to celestial coordinates taking
      non-linear distortion into account with the World Coordinate System
      FITS keywords as used in MegaPipe.
  
      For the inverse operation see sky2xy
  
      Reference material:
      http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/megapipe/docs/CD_PV_keywords.pdf
  
      Args:
        x, y: float
          Input pixel coordinate
        crpix1: float
          Tangent point x, pixels
        crpix2: float
          Tangent point y, pixels
        crval1: float
          Tangent point RA, degrees
        crval2:
          Tangent point Dec, degrees
        cd: 2d array
          Expresses the scale, the rotation and any possible skew of the image
          with respect to the sky.
        pv: 2d array
        nord: int
          order of the fit
  
      Returns:
        ra: float
          Right ascension
        dec: float
          Declination
      """
      xp = x - crpix1
      yp = y - crpix2
  
      # IMPORTANT NOTE: 0-based indexing in Python means indexing for values
      # in cd and pv will be shifted from in the paper.
      x_deg = cd[0][0] * xp + cd[0][1] * yp
      y_deg = cd[1][0] * xp + cd[1][1] * yp
  
      if nord < 0:
          xi = x
          eta = y
      else:
          xi = pv[0][0]
          eta = pv[1][0]
  
      if nord >= 1:
          r = numpy.sqrt(x_deg ** 2 + y_deg ** 2)
          xi += pv[0][1] * x_deg + pv[0][2] * y_deg + pv[0][3] * r
          eta += pv[1][1] * y_deg + pv[1][2] * x_deg + pv[1][3] * r
  
          if nord >= 2:
              x2 = x_deg ** 2
              xy = x_deg * y_deg
              y2 = y_deg ** 2
  
              xi += pv[0][4] * x2 + pv[0][5] * xy + pv[0][6] * y2
              eta += pv[1][4] * y2 + pv[1][5] * xy + pv[1][6] * x2
  
              if nord >= 3:
                  x3 = x_deg ** 3
                  x2y = x2 * y_deg
                  xy2 = x_deg * y2
                  y3 = y_deg ** 3
  
                  xi += pv[0][7] * x3 + pv[0][8] * x2y + pv[0][9] * xy2 + pv[0][10] * y3
                  eta += pv[1][7] * y3 + pv[1][8] * xy2 + pv[1][9] * x2y + pv[1][10] * x3
  
      xir = xi / PI180
      etar = eta / PI180
  
      ra0 = crval1 / PI180
      dec0 = crval2 / PI180
  
      ctan = numpy.tan(dec0)
      ccos = numpy.cos(dec0)
      raoff = numpy.arctan2(xir / ccos, 1 - etar * ctan)
      ra = raoff + ra0
      dec = numpy.arctan(numpy.cos(raoff) / ((1 - (etar * ctan)) / (etar + ctan)))
  
      ra *= PI180
      ra = numpy.where(ra < 0 , ra + 360, ra)
      ra = numpy.where(ra > 360, ra - 360, ra)
      # f = numpy.int(ra < 0) + 360
      # print f
      # if ra < 0:
      #     ra += 360
      # if ra >= 360:
      #     ra -= 360
  
      dec *= PI180
  
      return ra * units.degree, dec * units.degree