pro pharo_25mas_distortion, x, y, ha, dec, crN=crN, $
	xdp=xdp, ydp=ydp, xp=xp, yp=yp, xc=xc, yc=yc, mag=mag, $
	angle=angle

; pharo_25mas_distortion	- Produce the corrected (xdp, ydp) 
;			coordinates at a given location (x,y) on the 
;			PHARO 25mas/pix array, for given
;			hour angle (ha), and declination (dec).
;
; x, y		- Input pixel coordinates.
; ha		- Telescope hour angle (west is negative).
; dec		- Telescope declination.
; crN		- Cassegrain ring angle at which celestial North and the
;		  detector y-axis are aligned (default=334.92 degrees).
; xdp, ydp	- Output distortion-corrected coordinates.
; xp, yp	- Output coordinates before the beam tilt correction.
; xc, yc	- Origin of the (x,y) system; (default = 512,512).
; angle		- Orientation of the beam tilt (default=82.5 degrees).
; mag		- Beam tilt magnification along the angle+90 axis
;		  (default=1.0069). Equation (4.11) in Metchev (2006)

if n_params(0) lt 4 then $
    message,'Usage: CORRECT_DISTORTION, x, y, hour_angle, dec'+$
	'[, crN=crN, xdp=xdp, ydp=ydp, xp=xp, yp=yp, xc=xc,'+$
	' mag=mag, angle=angle]'
if n_elements(xc) lt 1 then xc = 512.
if n_elements(yc) lt 1 then yc = 512.
if n_elements(mag) eq 0 then mag = 1.0069	; Eq. (4.11) of Metchev (2006)
if n_elements(angle) eq 0 then angle = 83.3	; Eq. (4.10) of Metchev (2006)
if n_elements(crN) eq 0 then crN= 334.92	; the value listed in Table
						; 4.6 of Metchev (2006)

; Convert to radians
pi = 3.141592653D
beamtilt_ang = angle/180.*pi
northPA = crN/180.*pi
cassring_ang = northPA

x = x-xc
y = y-yc

; Set the xp coefficients in xp(x,y) = a0+a1*y+a2*y^2+a3*x+a5x^2
; Table 4.5 in Metchev (2006).
a = dblarr(6)
b = dblarr(3)
a1 = get_a1 (ha, dec)
a[0] = 0.01839547
a[1] = a1
a[2] = -4.510687e-06
a[3] = 0.9998296 
a[4] = 0
a[5] = -4.365845e-06

; Set the yp coefficients in yp(x,y) = b0+b1*y+b2*x
; Table 4.5 in Metchev (2006).
b[0] = 0.01932465
b[1] = 1.00001000
b2 = get_b2 (ha, dec)
b[2] = b2

; Correct for the detector distortion
; See Equations (4.1-4.2) in Metchev (2006).
xp = a[0] + a[1]*y + a[2]*y^2 + a[3]*x + a[4]*y*x + a[5]*x^2 
yp = b[0] + b[1]*y + b[2]*x

; Correct for the distortion due to the primary+secondary beams.
; Project along the (beamtilt_ang,beamtilt_ang+90) coordinate set and apply 
; the supplied magnification factor mag along the beamtilt_ang+90 axis.  
; Then de-project without de-magnifying.  The beam tilt angle relative
; to the detector also depends on the CR angle.  Tested - it seems that
; most of the distortion is unrelated to the primary+secondary, but comes
; from the positioning of the detector in the focal plane.
;
; See Equations (4.5) and (4.11) in Metchev (2006)
xr = xp*cos(beamtilt_ang)+yp*sin(beamtilt_ang)
yr = (-xp*sin(beamtilt_ang)+yp*cos(beamtilt_ang)) * mag
; See Equation (4.6) in Metchev (2006)
xdp = xr*cos(beamtilt_ang)-yr*sin(beamtilt_ang)
ydp = xr*sin(beamtilt_ang)+yr*cos(beamtilt_ang)

xdp += xc
ydp += yc

x += xc
y += yc

end


function get_a1, ha, dec

; get_a1	- Calculate the coefficient (a1) in front of y in the
;		dependence of the undistorted coordinates (xp,yp) on the 
;		measured ones (x,y): xp(x,y) = a0+a1*y+a2*y^2+a3*x+a5*x^2
;		These are coefficients alpha_0 to alpha_9 in Table 4.4 
;		of Metchev (2006).

par = dblarr(3)
a = dblarr(20)
a_err = dblarr(20)

par[0] = ha
par[1] = dec

; Based on 300 measurements at all hour angles, decs and CR = 333.5
a[0] = -0.0002496325
a[1] =  7.222127e-06
a[2] =  9.158649e-08
a[3] = -5.467939e-10
a[4] = -8.138857e-06
a[5] =  0
a[6] =  0
a[7] =  2.137211e-08
a[8] = -1.030719e-09
a[9] =  3.659832e-10
a_err[0] =  5.378355e-06
a_err[1] =  1.469128e-07
a_err[2] =  6.465733e-09
a_err[3] =  6.509733e-11
a_err[4] =  1.060632e-07
a_err[5] =  0
a_err[6] =  0
a_err[7] =  1.687558e-09
a_err[8] =  3.126475e-11
a_err[9] =  1.681164e-11

a1 = 0
k = 0
; Equation (4.3) in Metchev (2006)
for i=0,3 do for j=0,3-i do begin 
	a1 += a[k] * par[0]^i * par[1]^j
	k += 1
endfor

return, a1

end


function get_b2, ha, dec

; get_b2	- Calculate the coefficient (b2) in front of x in the
;		dependence of the undistorted coordinates (xp,yp) on the 
;		measured ones (x,y): yp(x,y) = b0+b1*y+b2*x
;		These are coefficients beta_0 to beta_9 in Table 4.4 
;		of Metchev (2006).

par = dblarr(3)
a = dblarr(20)
a_err = dblarr(20)

par[0] = ha
par[1] = dec

; Based on 300 measurements at all hour angles, decs and CR = 333.5
a[0] =  0.0001996558
a[1] = -1.054056e-05
a[2] = -7.407272e-08
a[3] =  5.37468e-10
a[4] =  7.580522e-06
a[5] =  0
a[6] =  0
a[7] =  1.065442e-08
a[8] =  8.868426e-10
a[9] = -2.833509e-10
a_err[0] = 5.278138e-06
a_err[1] = 1.457331e-07
a_err[2] = 6.395554e-09
a_err[3] = 6.431538e-11
a_err[4] = 1.042847e-07
a_err[5] = 0
a_err[6] = 0
a_err[7] = 1.655009e-09
a_err[8] = 3.064197e-11
a_err[9] = 1.65091e-11

b2 = 0
k = 0
; Equation (4.4) in Metchev (2006)
for i=0,3 do for j=0,3-i do begin 
	b2 += a[k] * par[0]^i * par[1]^j
	k += 1
endfor

return, b2

end