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
def sky2xypv ( ra , dec , crpix1 , crpix2 , crval1 , crval2 , dc , pv , nord , maxiter = 300 ) : 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 : x = xi y = eta else : tolerance = 0.001 / 3600 x = xi y = eta iteration = 0 converged = False while not converged : assert nord >= 0 f = pv [ 0 ] [ 0 ] g = pv [ 1 ] [ 0 ] 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
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