;+
; NAME: PIXEL_PHASE_CORRECT_GAUSS
;
;
;
; PURPOSE: Correct for pixel phase the aperture photometry flux(es)
; OBSERVED_FLUX of point sources observed with Warm IRAC,
; given a moment centroid at pixel position(s) X,Y. (The
; coordinate system here is defined with the center of a
; pixel having integer values and the edges having
; half-integer values). The model of pixel phase response is
; a "double-gaussian", summation of gaussians in the
; orthogonal pixel phase directions:
;
; RESPONSE(xphase,yphase) = F0 + deltaF_x*exp(-WRAP(xphase-X0,0.5)^2/(2*sigma_x^2)) +
; deltaF_y*exp(-WRAP(yphase-Y0,0.5)^2/(2*sigma_y^2))
;
; Here, xphase and yphase are the observed pixel phase
; (X-ROUND(X) and Y-ROUND(Y)), deltaF_x and
; deltaF_y are the peak to baseline offsets of the x and y
; gaussians, X0 and Y0 are the central x and y pixel phases
; of the gaussians, sigma_x and sigma_y are the gaussian
; sigmas, and F0_pp is the baseline relative flux (relative
; flux at infinity). The function WRAP(z,0.5) wraps the pixel
; phase offset so that it is periodic at +/- 0.5 pixels
; (z=+0.5 becomes z=-0.5, continuing to 0, and vice
; versa). The center of the wrap (z=0) is the phase function
; peak. The response function is normalized such that its
; integral over a pixel is unity. In other words, for random
; pointing, the average correction is 1.
;
;
; CATEGORY: Photometry, corrections
;
;
; CALLING SEQUENCE: IDL> corrected_flux = PIXEL_PHASE_CORRECT_GAUSS(observed_flux,x,y,channel [,APER]
; [,PARAMS=PARAMS]
; [,/SUBARRAY] [,/HELP])
;
;
;
; INPUTS: OBSERVED_FLUX - A scalar or vector containing observed
; aperture flux measurements in Jy or
; proportional units.
; X, Y - Scalars or vectors with the same number of
; elements as OBSERVED_FLUX, giving the X
; and Y pixel centroid location(s) of the
; source for each flux measurement. For the
; correction to be valid, centroids should
; be obtained using the moment method.
; Alternately, X and/or Y can be the X and Y
; pixel phases of the measurements (the
; pixel phase correction does not currently depend on
; absolute position on the array, only pixel
; phase). The coordinate system is defined
; with the center of a pixel having a whole
; number pixel value (pixel phase 0.0) and
; the edge having half-integer values (pixel
; phase +/- 0.5).
; CHANNEL -The IRAC channel number for the
; observations. If a scalar value, then all
; observations will be assumed to have the
; same channel. If a vector the it should
; have the same number of elements as
; OBSERVED_FLUX (otherwise only the first
; N_ELEMENTS(channel) elements of
; OBSERVED_FLUX will be corrected).
; Allowable values are 1 and 2. Can be an
; integer, floating point or string
; variable.
;
;
; OPTIONAL INPUTS: APER -An optional scalar string describing
; the aperture used to perform the
; centroiding and photometry. This entry
; consists of three integers, separated
; by underscores (_):
; (1) the radius of the aperture, in
; pixels;
; (2) the inner radius of the sky
; annulus, in pixels; and
; (3) the outer radius of the sky
; annulus.
;
; Currently available values are:
;
; APER RADIUS ANNULUS
; =============================
; 2_2_6 2 2 - 6
; 2_12_20 12 - 20
; 3_3_7 3 3 - 7
; 3_12_20 12 - 20
; 5_5_10 5 5 - 10
; 4_4_8 4 4 - 8
; 4_12_20 4 12 - 20
; 5_12_20 12 - 20
; 8_12_20 8 12 - 20
; 9_12_20 9 12 - 20
; 10_12_20 10 12 - 20 ***DEFAULT***
; 11_12_20 11 12 - 20
; 12_12_20 12 12 - 20
; 15_15_25 15 15 - 25
;
; If APER is not passed, or is not passed with one
; of the values above, it will default to '10_12_20'.
;
;
; KEYWORD PARAMETERS: PARAMS - An optional structure variable that allows you
; to override the default pixel phase
; response parameters. The structure is
; initialized/defined using the
; following syntax (or similar):
; IDL> params = {x0:fltarr(2),y0:fltarr(2),sigma_x:fltarr(2),sigma_y:fltarr(2),$
; deltaF_x:fltarr(2),deltaF_y:fltarr(2),F0:fltarr(2),theta:fltarr(2)}
;
; Each tag of PARAMS specifies a
; 2-element vector, whose 2 elements
; give the pixel phase response
; function parameters for Channels 1
; and 2:
; PARAMS.X0: Peak pixel phase of the X gaussian
; PARAMS.Y0: Peak pixel phase of the Y gaussian
; PARAMS.SIGMA_X: Sigma width of the X gaussian
; PARAMS.SIGMA_Y: Sigma width of the Y gaussian
; PARAMS.deltaF_X: Peak-to-baseline offset of
; the X gaussian
; PARAMS.deltaF_Y: Peak-to-baseline offset of
; the Y gaussian
; PARAMS.F0: Baseline response at infinity
; PARAMS.THETA: Angle (in degrees) that the X
; Gaussian makes with the X pixel
; axis. (Y Gaussian is orthogonal
; to X Gaussian). In the default
; model, the angle is assumed to be
; 0.
; If you set PARAMS, then APER will be ignored
;
; /SUBARRAY - use the parameter set for subarray mode
;
; /HELP - print this header
;
;
; OUTPUTS: This routine returns the corrected flux value(s) in the same
; units as OBSERVED_FLUX. The result will have the same
; number of elements and type as OBSERVED_FLUX.
;
;
;
; OPTIONAL OUTPUTS: None
;
;
;
; COMMON BLOCKS: None
;
;
;
; SIDE EFFECTS: The response function is normalized such
; that its integral over a pixel is unity. In other words, for
; random pointing, the average result would be no correction.
;
;
;
; RESTRICTIONS: (1) Inputs OBSERVED_FLUX, X, and Y should have
; the same number of elements. If CHANNEL is a vector,
; then it should also have the same number of elements
; as OBSERVED_FLUX, X, and Y.
;
;
; PROCEDURE:
;
;
;
; EXAMPLE:
;
;
;
; MODIFICATION HISTORY: 2/24/10: Created. Parameters based on Focus
; Check data. J. Ingalls (SSC)
; 4/20/10: Introduced PARAMS structure to
; override default parameters. Changed
; default parameters to reflect Full
; Array measurements, with more robust
; normalization. J. Ingalls (SSC)
; 5/03/10: Introduced APER parameter to accomodate
; variations of pixel phase response with
; aperture. Added /HELP keyword.
; J. Ingalls (SSC)
; 5/18/10: Introduced /SUBARRAY parameter to
; correct pixel phase response in
; subarray FOV.
; J. Ingalls (SSC)
; 10/27/11: Added 4-pixel radius aperture (with 4,8 and 12,20 background),
; Recomputed parameters for all apertures using
; updated centroiding technique. Now all apertures have
; valid results.
; J. Ingalls (SSC)
;
; Copyright (C) <2011>
; This program is free software: you can redistribute it and/or modify
; it under the terms of the GNU General Public License as published by
; the Free Software Foundation, either version 3 of the License, or
; any later version.
;
; This program is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
; GNU General Public License for more details.
;
; You should have received a copy of the GNU General Public License
; along with this program. If not, see .
;
;-
FUNCTION pixel_phase_double_gauss,xphase,yphase,x0,y0,sigma_x,sigma_y,deltaF_x,deltaF_y,F0,theta
dx = xphase-x0
dy = yphase-y0
IF Theta EQ 0 THEN BEGIN
dx_prime = dx
dy_prime = dy
ENDIF ELSE BEGIN
dx_prime = dx * cos(theta) + dy * sin(theta)
dy_prime = -dx * sin(theta) + dy * cos(theta)
ENDELSE
iadj_x = where(abs(dx_prime) GT 0.5,nadj_x)
IF nadj_x NE 0 THEN dx_prime[iadj_x] = 1.0 - abs(dx_prime[iadj_x])
iadj_y = where(abs(dy_prime) GT 0.5,nadj_y)
IF nadj_y NE 0 THEN dy_prime[iadj_y] = 1.0 - abs(dy_prime[iadj_y])
relative_flux = deltaF_x * exp( -dx_prime^2/(2*sigma_x^2) ) + deltaF_y * exp( -dy_prime^2/(2*sigma_y^2) ) + F0
RETURN,relative_flux
END
FUNCTION pixel_phase_correct_gauss,observed_flux,x,y,channel,aper,params=params,subarray=subarray,help=help
;;
;; Correct the flux(es) OBSERVED_FLUX of measured objects with moment centroid at
;; pixel position(s) X,Y. Use the pixel phase response double
;; gaussian fit. OBSERVED_FLUX,X,Y, and CHANNEL are arrays or scalars with the
;; same number of elements
;;
IF KEYWORD_SET(HELP) THEN doc_library,'pixel_phase_correct_gauss'
AP_ARRAY = ['15_15_25','12_12_20','11_12_20','10_12_20','9_12_20','8_12_20','5_12_20','5_5_10','4_12_20','4_4_8','3_12_20','3_3_7','2_12_20','2_2_6']
IF N_ELEMENTS(APER) EQ 0 THEN aper = '10_12_20'
nap = 0
WHILE nap EQ 0 DO BEGIN
iap = where(AP_ARRAY EQ APER,nap)
IF NAP EQ 0 THEN aper = '10_12_20'
ENDWHILE
nflux = N_ELEMENTS(observed_flux)
nx = N_ELEMENTS(x)
ny = N_ELEMENTS(y)
nch = N_ELEMENTS(channel)
IF nch EQ 1 THEN BEGIN
IF Nflux GT 1 THEN channel = FIX(MAKE_ARRAY(SIZE=size(observed_flux),VALUE=channel)) ELSE channel = FIX(channel)
nch = N_ELEMENTS(channel)
ENDIF
IF nflux NE nx OR nflux NE ny OR nflux NE nch THEN BEGIN
print,'PIXEL_PHASE_CORRECT_GAUSS: Input array sizes don''t match. Exiting.'
RETALL
ENDIF
IF N_ELEMENTS(PARAMS) EQ 0 then BEGIN
;;; Parameters constrained such that average over pixel is unity.
IF KEYWORD_SET(SUBARRAY) EQ 0 THEN BEGIN
;;; From run of 10/25/11, adding 4-pixel radius apertures and using new centroiding function
X0_CH1 = [-0.075,0.192,0.192,0.192,0.191,0.192,0.191,0.191,0.194,0.193,0.190,0.186,0.208,0.210]
Y0_CH1 = [0.051,0.029,0.029,0.030,0.029,0.030,0.027,0.027,0.023,0.023,0.022,0.021,0.036,0.034]
SIG_X_CH1 = [0.990,0.175,0.187,0.187,0.189,0.189,0.189,0.190,0.190,0.190,0.189,0.191,0.198,0.197]
SIG_Y_CH1 = [0.191,0.169,0.164,0.163,0.162,0.161,0.163,0.163,0.163,0.164,0.164,0.164,0.163,0.163]
DELF_X_CH1 = [0.0396,0.0279,0.0286,0.0288,0.0290,0.0291,0.0299,0.0298,0.0309,0.0308,0.0327,0.0330,0.0351,0.0350]
DELF_Y_CH1 = [0.0433,0.0451,0.0448,0.0453,0.0453,0.0461,0.0475,0.0479,0.0487,0.0492,0.0544,0.0540,0.0633,0.0634]
F0_CH1 = [0.941,0.969,0.968,0.968,0.968,0.968,0.967,0.966,0.966,0.965,0.962,0.962,0.957,0.957]
X0_CH2 = [-0.044,0.089,0.089,0.089,0.088,0.089,0.088,0.091,0.083,0.087,0.070,0.071,0.047,0.051]
Y0_CH2 = [-0.179,0.100,0.099,0.100,0.098,0.098,0.095,0.091,0.096,0.095,0.098,0.099,0.095,0.104]
SIG_X_CH2 = [0.014,0.218,0.217,0.215,0.212,0.210,0.209,0.209,0.206,0.207,0.203,0.201,0.100,0.104]
SIG_Y_CH2 = [0.018,0.227,0.227,0.228,0.224,0.221,0.214,0.209,0.215,0.210,0.216,0.218,0.177,0.191]
DELF_X_CH2 = [0.1457,0.0172,0.0173,0.0172,0.0170,0.0171,0.0178,0.0176,0.0173,0.0177,0.0175,0.0179,0.0078,0.0083]
DELF_Y_CH2 = [0.0312,0.0188,0.0189,0.0189,0.0187,0.0185,0.0187,0.0183,0.0188,0.0191,0.0186,0.0192,0.0260,0.0270]
F0_CH2 = [0.994,0.980,0.980,0.980,0.981,0.981,0.981,0.981,0.981,0.981,0.981,0.981,0.987,0.985]
ENDIF ELSE BEGIN
;;; From run of 10/25/11, adding 4-pixel radius apertures
X0_CH1 = [0.125,0.125,0.129,0.129,0.130,0.130,0.130,0.129,0.131,0.130,0.128,0.126,0.130,0.131]
Y0_CH1 = [0.039,0.036,0.033,0.033,0.033,0.033,0.032,0.032,0.030,0.032,0.027,0.027,0.039,0.037]
SIG_X_CH1 = [0.194,0.187,0.184,0.184,0.184,0.184,0.183,0.185,0.183,0.184,0.183,0.182,0.177,0.176]
SIG_Y_CH1 = [0.192,0.194,0.191,0.192,0.192,0.193,0.194,0.196,0.194,0.197,0.195,0.197,0.191,0.191]
DELF_X_CH1 = [0.0358,0.0362,0.0359,0.0360,0.0363,0.0366,0.0384,0.0388,0.0396,0.0394,0.0416,0.0419,0.0425,0.0430]
DELF_Y_CH1 = [0.0434,0.0438,0.0434,0.0437,0.0438,0.0447,0.0468,0.0472,0.0476,0.0488,0.0520,0.0521,0.0596,0.0596]
F0_CH1 = [0.962,0.962,0.963,0.963,0.962,0.962,0.960,0.959,0.959,0.958,0.956,0.956,0.953,0.953]
X0_CH2 = [0.071,0.074,0.070,0.067,0.065,0.064,0.052,0.051,0.048,0.049,0.035,0.037,0.009,0.011]
Y0_CH2 = [0.106,0.105,0.107,0.113,0.111,0.110,0.099,0.098,0.089,0.088,0.091,0.095,0.131,0.133]
SIG_X_CH2 = [0.175,0.175,0.172,0.164,0.167,0.168,0.173,0.177,0.187,0.188,0.196,0.198,0.171,0.174]
SIG_Y_CH2 = [0.265,0.291,0.284,0.267,0.263,0.257,0.212,0.207,0.192,0.180,0.179,0.174,0.177,0.179]
DELF_X_CH2 = [0.0167,0.0182,0.0184,0.0185,0.0182,0.0182,0.0185,0.0187,0.0193,0.0193,0.0215,0.0220,0.0205,0.0214]
DELF_Y_CH2 = [0.0169,0.0201,0.0195,0.0189,0.0192,0.0187,0.0167,0.0167,0.0162,0.0161,0.0162,0.0169,0.0223,0.0236]
F0_CH2 = [0.982,0.979,0.979,0.981,0.980,0.981,0.983,0.983,0.983,0.984,0.982,0.982,0.981,0.980]
ENDELSE
;;; CH1 CH2
x0 = [x0_ch1[iap[0]],x0_ch2[iap[0]]]
y0 = [y0_ch1[iap[0]],y0_ch2[iap[0]]]
sigma_x = [sig_x_ch1[iap[0]],sig_x_ch2[iap[0]]]
sigma_y = [sig_y_ch1[iap[0]],sig_y_ch2[iap[0]]]
deltaF_x = [delf_x_ch1[iap[0]],delf_x_ch2[iap[0]]]
deltaF_y = [delf_y_ch1[iap[0]],delf_y_ch2[iap[0]]]
F0 = [F0_ch1[iap[0]],F0_ch2[iap[0]]]
Theta = [0.,0.]
ENDIF ELSE BEGIN
x0 = params.x0
y0 = params.y0
sigma_x = params.sigma_x
sigma_y = params.sigma_y
deltaF_x = params.deltaF_x
deltaF_y = params.deltaF_y
F0 = params.F0
Theta = params.Theta
ENDELSE
chan = STRCOMPRESS(channel,/REMOVE_ALL)
ch = FIX(channel)
;;; Determine pixel X and Y phase
xphase = x - ROUND(x)
yphase = y - ROUND(y)
NPTS = N_ELEMENTS(xphase)
IF NFLUX GT 1 THEN response = make_array(size=SIZE(observed_flux)) ELSE response = fltarr(1)
ich1 = WHERE(ch EQ 1,nch1)
ich2 = WHERE(ch EQ 2,nch2)
IF nch1 NE 0 THEN response[ich1] = pixel_phase_double_gauss(xphase[ich1],yphase[ich1],x0[0],y0[0],$
sigma_x[0],sigma_y[0],$
deltaF_x[0],deltaF_y[0],F0[0],theta[0]/!radeg)
IF nch2 NE 0 THEN response[ich2] = pixel_phase_double_gauss(xphase[ich2],yphase[ich2],x0[1],y0[1],$
sigma_x[1],sigma_y[1],$
deltaF_x[1],deltaF_y[1],F0[1],theta[1]/!radeg)
RETURN,observed_flux/response
END