function TPV_eval, xy ;+ ; NAME: ; TPV_EVAL ; PURPOSE: ; Compute distorted coordinates given TPV (Tangent + PV_ polynomial) ; coefficients. ; EXPLANATION: ; See http://fits.gsfc.nasa.gov/registry/tpvwcs.html for the TPV convention ; ; This distortion convention is used by the SCAMP software ; ( http://www.astromatic.net/software/scamp ) though SCAMP does not ; include the '-TPV' in the CTYPE keyword. ; ; The coefficients are passed via common block. This is because this ; routine is called by the intrinisc BROYDEN() function in AD2XY, and ; common blocks are the only way to pass parameters to the user supplied ; function in BROYDEN(). ; CALLING SEQUENCE: ; res = TPV_EVAL(xy) ; INPUTS: ; xy - 2 elements vector giving the undistorted X,Y position ; OUTPUTS: ; res - 2 element vector giving the distorted position ; COMMON BLOCKS: ; common broyden_coeff,pv1,ycoeff ; ; pv1, YCOEFF are both vectors giving the TPV coefficients ; REVISION HISTORY: ; Written W. Landsman Dec 2013 ; Correct several typos for 4th power terms M. Sullivan Mar 2014 ; Corrected 4th order terms once again and ; added 5th,6th and 7th order terms Arjun Dey Sep 3, 2015 ;- compile_opt idl2,hidden common broyden_coeff,pv1,pv2 Npv1 = N_elements(pv1) NPv2 = N_elements(pv2) if N_elements(xy) EQ 2 then begin x = xy[0] y = xy[1] endif else begin x = reform(xy[*,0]) y = reform(xy[*,1]) endelse x2 = x*x y2 = y*y xp = pv1[0] + pv1[1]*x + pv1[2]*y if Npv1 GT 3 && (pv1[3] NE 0.0) then xp += pv1[3]*sqrt(x2 + y2) if Npv1 GT 4 && (pv1[4] NE 0.0) then xp += pv1[4]*x2 if Npv1 GT 5 && (pv1[5] NE 0.0) then xp += pv1[5]*x*y if Npv1 GT 6 && (pv1[6] NE 0.0) then xp += pv1[6]*y2 if Npv1 GT 7 then begin if pv1[7] NE 0.0 then xp += pv1[7]*x^3 if Npv1 GT 8 && (pv1[8] NE 0.0) then xp += pv1[8]*x2*y if Npv1 GT 9 && (pv1[9] NE 0.0) then xp += pv1[9]*x*y2 if Npv1 GT 10 && (pv1[10] NE 0.0) then xp += pv1[10]*y2*y if Npv1 GT 11 && (pv1[11] NE 0.0) then xp += pv1[11]*sqrt(x2+y2)^3 if Npv1 GT 12 then begin if (pv1[12] NE 0.0) then xp += pv1[12]*x2*x2 if Npv1 GT 13 && (pv1[13] NE 0.0) then xp += pv1[13]*x2*x*y if Npv1 GT 14 && (pv1[14] NE 0.0) then xp += pv1[14]*x2*y2 if Npv1 GT 15 && (pv1[15] NE 0.0) then xp += pv1[15]*x*y2*y if Npv1 GT 16 && (pv1[16] NE 0.0) then xp += pv1[16]*y2*y2 if Npv1 GT 17 then begin x4 = x2*x2 y4 = y2*y2 if (pv1[17] NE 0.0) then xp += pv1[17]*x4*x if Npv1 GT 18 && (pv1[18] NE 0.0) then xp += pv1[18]*x4*y if Npv1 GT 19 && (pv1[19] NE 0.0) then xp += pv1[19]*x2*x*y2 if Npv1 GT 20 && (pv1[20] NE 0.0) then xp += pv1[20]*x2*y2*y if Npv1 GT 21 && (pv1[21] NE 0.0) then xp += pv1[21]*x*y4 if Npv1 GT 22 && (pv1[22] NE 0.0) then xp += pv1[22]*y4*y if Npv1 GT 23 && (pv1[23] NE 0.0) then xp += pv1[23]*sqrt(x2+y2)^5 if Npv1 GT 24 then begin if (pv1[24] NE 0.0) then xp += pv2[24]*x4*x2 if Npv1 GT 25 && (pv1[25] NE 0.0) then xp += pv1[25]*x4*x*y if Npv1 GT 26 && (pv1[26] NE 0.0) then xp += pv1[26]*x4*y2 if Npv1 GT 27 && (pv1[27] NE 0.0) then xp += pv1[27]*x2*x*y2*y if Npv1 GT 28 && (pv1[28] NE 0.0) then xp += pv1[28]*x2*y4 if Npv1 GT 29 && (pv1[29] NE 0.0) then xp += pv1[29]*x*y4*y if Npv1 GT 30 && (pv1[30] NE 0.0) then xp += pv1[30]*y4*y2 if Npv1 GT 31 then begin if (pv1[31] NE 0.0) then xp += pv1[31]*x4*x2*x if Npv1 GT 32 && (pv1[32] NE 0.0) then xp += pv1[32]*x4*x2*y if Npv1 GT 33 && (pv1[33] NE 0.0) then xp += pv1[33]*x4*x*y2 if Npv1 GT 34 && (pv1[34] NE 0.0) then xp += pv1[34]*x4*y2*y if Npv1 GT 35 && (pv1[35] NE 0.0) then xp += pv1[35]*x2*x*y4 if Npv1 GT 36 && (pv1[36] NE 0.0) then xp += pv1[36]*x2*y*y4 if Npv1 GT 37 && (pv1[37] NE 0.0) then xp += pv1[37]*x*y2*y4 if Npv1 GT 38 && (pv1[38] NE 0.0) then xp += pv1[38]*y4*y2*y if Npv1 GT 39 && (pv1[39] NE 0.0) then xp += pv1[39]*sqrt(x2+y2)^7 if Npv1 GT 40 then print,'PV1 TERMS > 40 (ORDER > 7) NOT EVALUATED' endif endif endif endif endif yp = pv2[0] + pv2[1]*y + pv2[2]*x if Npv2 GT 3 && (pv2[3] NE 0.0) then yp += pv2[3]*sqrt(x2 + y2) if NPv2 GT 4 && (pv2[4] NE 0.0) then yp += pv2[4]*y2 if NPv2 GT 5 && (pv2[5] NE 0.0) then yp += pv2[5]*x*y if NPv2 GT 6 && (pv2[6] NE 0.0) then yp += pv2[6]*x2 if NPv2 GT 7 then begin if pv2[7] NE 0.0 then yp += pv2[7]*y^3 if NPv2 GT 8 && (pv2[8] NE 0.0) then yp += pv2[8]*y2*x if NPv2 GT 9 && (pv2[9] NE 0.0) then yp += pv2[9]*y*x2 if NPv2 GT 10 && (pv2[10] NE 0.0) then yp += pv2[10]*x2*x if NPv2 GT 11 && (pv2[11] NE 0.0) then yp += pv2[11]*sqrt(x2+y2)^3 if NPv2 GT 12 then begin if (pv2[12] NE 0.0) then yp += pv2[12]*y2*y2 if NPv2 GT 13 && (pv2[13] NE 0.0) then yp += pv2[13]*y2*y*x if NPv2 GT 14 && (pv2[14] NE 0.0) then yp += pv2[14]*y2*x2 if NPv2 GT 15 && (pv2[15] NE 0.0) then yp += pv2[15]*y*x2*x if NPv2 GT 16 && (pv2[16] NE 0.0) then yp += pv2[16]*x2*x2 if Npv2 GT 17 then begin x4 = x2*x2 y4 = y2*y2 if (pv2[17] NE 0.0) then yp += pv2[17]*y4*y if Npv2 GT 18 && (pv2[18] NE 0.0) then yp += pv2[18]*y4*x if Npv2 GT 19 && (pv2[19] NE 0.0) then yp += pv2[19]*y2*y*x2 if Npv2 GT 20 && (pv2[20] NE 0.0) then yp += pv2[20]*y2*x2*x if Npv2 GT 21 && (pv2[21] NE 0.0) then yp += pv2[21]*y*x4 if Npv2 GT 22 && (pv2[22] NE 0.0) then yp += pv2[22]*x4*x if Npv2 GT 23 && (pv2[23] NE 0.0) then yp += pv2[23]*sqrt(x2+y2)^5 if Npv2 GT 24 then begin if (pv2[24] NE 0.0) then yp += pv2[24]*y4*y2 if Npv2 GT 25 && (pv2[25] NE 0.0) then yp += pv2[25]*y4*y*x if Npv2 GT 26 && (pv2[26] NE 0.0) then yp += pv2[26]*y4*x2 if Npv2 GT 27 && (pv2[27] NE 0.0) then yp += pv2[27]*y2*y*x2*x if Npv2 GT 28 && (pv2[28] NE 0.0) then yp += pv2[28]*y2*x4 if Npv2 GT 29 && (pv2[29] NE 0.0) then yp += pv2[29]*y*x4*x if Npv2 GT 30 && (pv2[30] NE 0.0) then yp += pv2[30]*x4*x2 if Npv2 GT 31 then begin if (pv2[31] NE 0.0) then yp += pv2[31]*y4*y2*y if Npv2 GT 32 && (pv2[32] NE 0.0) then yp += pv2[32]*y4*y2*x if Npv2 GT 33 && (pv2[33] NE 0.0) then yp += pv2[33]*y4*y*x2 if Npv2 GT 34 && (pv2[34] NE 0.0) then yp += pv2[34]*y4*x2*x if Npv2 GT 35 && (pv2[35] NE 0.0) then yp += pv2[35]*y2*y*x4 if Npv2 GT 36 && (pv2[36] NE 0.0) then yp += pv2[36]*y2*x*x4 if Npv2 GT 37 && (pv2[37] NE 0.0) then yp += pv2[37]*y*x2*x4 if Npv2 GT 38 && (pv2[38] NE 0.0) then yp += pv2[38]*x4*x2*x if Npv2 GT 39 && (pv2[39] NE 0.0) then yp += pv2[39]*sqrt(x2+y2)^7 if Npv2 GT 40 then print,'PV2 TERMS > 40 (ORDER > 7) NOT EVALUATED' endif endif endif endif endif return, [[xp],[yp]] end