Using SIP Coefficients for optical distortion correction

This transformation converts pixel coordinates fitsX, fitsY from an astronomical image to celestial coordinates α, δ (right ascension and declination) using the WCS (World Coordinate System). The reference position α0, δ0 is the celestial position of the reference pixel crpix1,crpix2 stored in the header as crval1, crval2. The transformation matrix is cd1_1, cd1_2, cd2_1, cd2_2. Assuming TAN projection.

// Transformation measured coordinates fitsX,fitsY-->celestial coordinates α,δ.

// Step 1, If available apply SIP coefficients to correct for image distortions. The corrected pixel coordinates U, V are calculated using the SIP coefficients (A_i_j, B_i_j) as

u:=fitsX -crpix1
v:=fitsY -crpix2

U = u + Σ A_i_j * u^i * v^j
V = v +

Σ B_i_j * u^i * v^j

// Apply the CD matrix transformation to get the intermediate standard coordinates
ξ = (cd1_1 * U + cd1_2 * V
η

 = (cd2_1 * U + cd2_2 * V

// Step3, conversion from standard coordinates to celestial coordinates. Assume Gnomonic (Tangent Plane) Projection ( CTYPE1= 'RA---TAN')
ξ = ξ * π/180
η = η * π/180

α0 = crval1

δ0 = crval2

sincos(δ0, sinδ0,cosδ0)
delta = cosδ0 -

η * sinδ0
α = α0 + arctan2(ξ, delta)
δ = arctan((sinδ0 + η * cosδ0) / sqrt(ξ^2 + delta^2))
α = α * 180/π

δ = δ * 180/π

Alternative formulas 2:

sincos(δ0,sinδ0,cosδ0)
α = α0 + arctan2(ξ,cosδ0 - V*sinδ0);
δ = arcsin((sinδ0 + η * cosδ0) / sqrt( 1 + ξ^2 + η^2))
α = α * 180/π
δ = δ * 180/π

Alternative formulas 3:

sincos(δ0,sinδ0,cosδ0);
delta_α = arctan2(η,cosδ0 - ξ*sinδ0)
α = α0 + delta_α
δ = arctan2((η*cosδ0+cosδ0) * sin( delta_α),ξ)

This transformation converts celestial coordinates α, δ (right ascension and declination) to pixel coordinates fitsX, fitsY in an image. The reference position α0, dec0 is the celestial position of the reference pixel at crpix1,crpix2 and stored in the header as crval1, crval2. he transformation matrix is cd1_1, cd1_2, cd2_1, cd2_2.. Assuming TAN projection.

// Transformation celestial α,δ position --> fitsX,fitsY

// Step 1, convert the celestial coordinates to intermediate standard coordinates (ξ, η) using the gnomonic (tangent) projection formula.:

α0 = crval1

δ0 = crval2

sincos(δ,sinδ,cosδ)
sincos(δ0,sinδ0,cosδ0)
sincos(α-α0,sin_delta_α,cos_delta_α)

cosc =

(sinδ * sinδ0 + cosδ * cosδ0 * cosdelta_α)
ξ = (cosδ * sin_delta_α) / cosc)
η = (cosδ0 * sinδ - sinδ0 * cosδ * cos_delta_α) / cosc
ξ = ξ * 180/π
η = η * 180 / π

// Step 2, calculate the determinant of the CD matrix
detCD = cd1_1 * cd2_2 - cd1_2 * cd2_1
// Apply the CD matrix transformation to the intermediate standard coordinates and adjust with the determinant. Then add the reference pixel coordinates (crpix1, crpix2).
u = (cd2_2 * ξ - cd1_2 * η) / detCD + crpix1

v = (-cd2_1 * ξ + cd1_1 * η) / detCD + crpix2

//Step 3, If available apply SIP coefficients to correct for image distortions.. The corrected pixel coordinates are calculated using the SIP coefficients (AP_i_j, BP_i_j) as follows:
fitsX = u + Σ AP_i_j * u^i * v^j

fitsY = v + Σ BP_i_j * u^i * v^j

Note1: The summations are carried over the orders of the SIP coefficients.
Note2: The Fits origin is at [1,1]

x' = x + a₁x + a₂y + a₃x² + a₄xy + a₅y² + a₆x³ + a₇x²y + a₈xy² + a₉y³

y' = y + b₁x + b₂y + b₃x² + b₄xy + b₅y² + b₆x³ + b₇x²y + b₈xy² + b₉y³

Relation with FITS Header Keywords

The SIP coefficients are stored in the FITS header and labeled according to their order and the axis they correspond to. For example:

The A_ORDER and B_ORDER specify the order of the polynomial for the pixel-to-sky correction. The AP_ORDER and BP_ORDER keywords specify the order of the polynomial for the sky-to-pixel correction.