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