Tibeuleu/FOC_Reduction
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

add transmittance correction to Stokes fluxes computation

Browse Source
This commit is contained in:
Thibault Barnouin
2024-06-18 16:37:30 +02:00
parent 97479f90fb
commit 7ab87b50d1
2 changed files with 49 additions and 36 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
package/FOC_reduction.py
View File

@@ -35,24 +35,27 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
# Background estimation
error_sub_type = 'freedman-diaconis' # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51))
subtract_error = 0.50
subtract_error = 0.01
display_bkg = True
# Data binning
rebin = True
pxsize = 0.10
px_scale = 'arcsec' # pixel, arcsec or full
pxsize = 2
px_scale = 'px' # pixel, arcsec or full
rebin_operation = 'sum' # sum or average
# Alignement
align_center = 'center' # If None will not align the images
display_align = False
display_align = True
display_data = False
# Transmittance correction
transmitcorr = True
# Smoothing
smoothing_function = 'combine' # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine
smoothing_FWHM = 0.150 # If None, no smoothing is done
smoothing_scale = 'arcsec' # pixel or arcsec
smoothing_FWHM = None # If None, no smoothing is done
smoothing_scale = 'px' # pixel or arcsec
# Rotation
rotate_data = False # rotation to North convention can give erroneous results
@@ -120,10 +123,11 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
data_array, error_array, headers, background = proj_red.get_error(data_array, headers, error_array, data_mask=data_mask, sub_type=error_sub_type, subtract_error=subtract_error, display=display_bkg, savename="_".join([figname, "errors"]), plots_folder=plots_folder, return_background=True)
# Align and rescale images with oversampling.
data_array, error_array, headers, data_mask = proj_red.align_data(
data_array, headers, error_array=error_array, background=background, upsample_factor=10, ref_center=align_center, return_shifts=False)
data_array, error_array, headers, data_mask, shifts, error_shifts = proj_red.align_data(
data_array, headers, error_array=error_array, background=background, upsample_factor=10, ref_center=align_center, return_shifts=True)
if display_align:
print("Image shifts: {} \nShifts uncertainty: {}".format(shifts, error_shifts))
proj_plots.plot_obs(data_array, headers, savename="_".join([figname, str(align_center)]), plots_folder=plots_folder, norm=LogNorm(
vmin=data_array[data_array > 0.].min()*headers[0]['photflam'], vmax=data_array[data_array > 0.].max()*headers[0]['photflam']))
@@ -154,7 +158,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
# see Jedrzejewski, R.; Nota, A.; Hack, W. J., A Comparison Between FOC and WFPC2
# Bibcode : 1995chst.conf...10J
I_stokes, Q_stokes, U_stokes, Stokes_cov = proj_red.compute_Stokes(
data_array, error_array, data_mask, headers, FWHM=smoothing_FWHM, scale=smoothing_scale, smoothing=smoothing_function, transmitcorr=False)
data_array, error_array, data_mask, headers, FWHM=smoothing_FWHM, scale=smoothing_scale, smoothing=smoothing_function, transmitcorr=transmitcorr)
I_bkg, Q_bkg, U_bkg, S_cov_bkg = proj_red.compute_Stokes(background, background_error, np.array(True).reshape(
1, 1), headers, FWHM=None, scale=smoothing_scale, smoothing=smoothing_function, transmitcorr=False)

63
package/lib/reduction.py
View File

@@ -1026,7 +1026,7 @@ def polarizer_avg(data_array, error_array, data_mask, headers, FWHM=None, scale=
return polarizer_array, polarizer_cov, pol_headers
def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale='pixel', smoothing='combine', transmitcorr=False):
def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale='pixel', smoothing='combine', transmitcorr=True):
"""
Compute the Stokes parameters I, Q and U for a given data_set
----------
@@ -1060,9 +1060,9 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale
Defaults to 'combine'. Won't be used if FWHM is None.
transmitcorr : bool, optional
Weither the images should be transmittance corrected for each filter
along the line of sight. Latest calibrated data products (.c0f) does
along the line of sight. Latest calibrated data products (.c1f) does
not require such correction.
Defaults to False.
Defaults to True.
----------
Returns:
I_stokes : numpy.ndarray
@@ -1132,8 +1132,8 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale
pol_eff[(i+2) % 3]*np.cos(2.*globals()["theta"][(i+2) % 3]))*2./transmit[i]
# Normalization parameter for Stokes parameters computation
A = (coeff_stokes[0, :]*transmit/2.).sum()
coeff_stokes = coeff_stokes/A
N = (coeff_stokes[0, :]*transmit/2.).sum()
coeff_stokes = coeff_stokes/N
I_stokes = np.zeros(pol_array[0].shape)
Q_stokes = np.zeros(pol_array[0].shape)
U_stokes = np.zeros(pol_array[0].shape)
@@ -1177,29 +1177,40 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale
s_Q2_stat = np.sum([coeff_stokes[1, i]**2*sigma_flux[i]**2 for i in range(len(sigma_flux))], axis=0)
s_U2_stat = np.sum([coeff_stokes[2, i]**2*sigma_flux[i]**2 for i in range(len(sigma_flux))], axis=0)
pol_flux_corr = np.array([pf*2./t for (pf, t) in zip(pol_flux, transmit)])
coeff_stokes_corr = np.array([cs*t/2. for (cs, t) in zip(coeff_stokes.T, transmit)]).T
# Compute the derivative of each Stokes parameter with respect to the polarizer orientation
dI_dtheta1 = 2.*pol_eff[0]/A*(pol_eff[2]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0])*(pol_flux[1]-I_stokes) -
pol_eff[1]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1])*(pol_flux[2]-I_stokes))
dI_dtheta2 = 2.*pol_eff[1]/A*(pol_eff[0]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1])*(pol_flux[2]-I_stokes) -
pol_eff[2]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2])*(pol_flux[0]-I_stokes))
dI_dtheta3 = 2.*pol_eff[2]/A*(pol_eff[1]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2])*(pol_flux[0]-I_stokes) -
pol_eff[0]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0])*(pol_flux[1]-I_stokes))
dI_dtheta1 = 2.*pol_eff[0]/N*(pol_eff[2]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0])*(pol_flux_corr[1]-I_stokes) -
pol_eff[1]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1])*(pol_flux_corr[2]-I_stokes) +
coeff_stokes_corr[0, 0]*(np.sin(2.*globals()["theta"][0])*Q_stokes-np.cos(2*globals()["theta"][0])*U_stokes))
dI_dtheta2 = 2.*pol_eff[1]/N*(pol_eff[0]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1])*(pol_flux_corr[2]-I_stokes) -
pol_eff[2]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2])*(pol_flux_corr[0]-I_stokes) +
coeff_stokes_corr[0, 1]*(np.sin(2.*globals()["theta"][1])*Q_stokes-np.cos(2*globals()["theta"][1])*U_stokes))
dI_dtheta3 = 2.*pol_eff[2]/N*(pol_eff[1]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2])*(pol_flux_corr[0]-I_stokes) -
pol_eff[0]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0])*(pol_flux_corr[1]-I_stokes) +
coeff_stokes_corr[0, 2]*(np.sin(2.*globals()["theta"][2])*Q_stokes-np.cos(2*globals()["theta"][2])*U_stokes))
dI_dtheta = np.array([dI_dtheta1, dI_dtheta2, dI_dtheta3])
dQ_dtheta1 = 2.*pol_eff[0]/A*(np.cos(2.*globals()["theta"][0])*(pol_flux[1]-pol_flux[2]) - (pol_eff[2]*np.cos(-2.*globals()
["theta"][2]+2.*globals()["theta"][0]) - pol_eff[1]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1]))*Q_stokes)
dQ_dtheta2 = 2.*pol_eff[1]/A*(np.cos(2.*globals()["theta"][1])*(pol_flux[2]-pol_flux[0]) - (pol_eff[0]*np.cos(-2.*globals()
["theta"][0]+2.*globals()["theta"][1]) - pol_eff[2]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2]))*Q_stokes)
dQ_dtheta3 = 2.*pol_eff[2]/A*(np.cos(2.*globals()["theta"][2])*(pol_flux[0]-pol_flux[1]) - (pol_eff[1]*np.cos(-2.*globals()
["theta"][1]+2.*globals()["theta"][2]) - pol_eff[0]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0]))*Q_stokes)
dQ_dtheta1 = 2.*pol_eff[0]/N*(np.cos(2.*globals()["theta"][0])*(pol_flux_corr[1]-pol_flux_corr[2]) - (pol_eff[2]*np.cos(-2.*globals()
["theta"][2]+2.*globals()["theta"][0]) - pol_eff[1]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1]))*Q_stokes +
coeff_stokes_corr[1, 0]*(np.sin(2.*globals()["theta"][0])*Q_stokes-np.cos(2*globals()["theta"][0])*U_stokes))
dQ_dtheta2 = 2.*pol_eff[1]/N*(np.cos(2.*globals()["theta"][1])*(pol_flux_corr[2]-pol_flux_corr[0]) - (pol_eff[0]*np.cos(-2.*globals()
["theta"][0]+2.*globals()["theta"][1]) - pol_eff[2]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2]))*Q_stokes +
coeff_stokes_corr[1, 1]*(np.sin(2.*globals()["theta"][1])*Q_stokes-np.cos(2*globals()["theta"][1])*U_stokes))
dQ_dtheta3 = 2.*pol_eff[2]/N*(np.cos(2.*globals()["theta"][2])*(pol_flux_corr[0]-pol_flux_corr[1]) - (pol_eff[1]*np.cos(-2.*globals()
["theta"][1]+2.*globals()["theta"][2]) - pol_eff[0]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0]))*Q_stokes +
coeff_stokes_corr[1, 2]*(np.sin(2.*globals()["theta"][2])*Q_stokes-np.cos(2*globals()["theta"][2])*U_stokes))
dQ_dtheta = np.array([dQ_dtheta1, dQ_dtheta2, dQ_dtheta3])
dU_dtheta1 = 2.*pol_eff[0]/A*(np.sin(2.*globals()["theta"][0])*(pol_flux[1]-pol_flux[2]) - (pol_eff[2]*np.cos(-2.*globals()
["theta"][2]+2.*globals()["theta"][0]) - pol_eff[1]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1]))*U_stokes)
dU_dtheta2 = 2.*pol_eff[1]/A*(np.sin(2.*globals()["theta"][1])*(pol_flux[2]-pol_flux[0]) - (pol_eff[0]*np.cos(-2.*globals()
["theta"][0]+2.*globals()["theta"][1]) - pol_eff[2]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2]))*U_stokes)
dU_dtheta3 = 2.*pol_eff[2]/A*(np.sin(2.*globals()["theta"][2])*(pol_flux[0]-pol_flux[1]) - (pol_eff[1]*np.cos(-2.*globals()
["theta"][1]+2.*globals()["theta"][2]) - pol_eff[0]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0]))*U_stokes)
dU_dtheta1 = 2.*pol_eff[0]/N*(np.sin(2.*globals()["theta"][0])*(pol_flux_corr[1]-pol_flux_corr[2]) - (pol_eff[2]*np.cos(-2.*globals()
["theta"][2]+2.*globals()["theta"][0]) - pol_eff[1]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1]))*U_stokes +
coeff_stokes_corr[2, 0]*(np.sin(2.*globals()["theta"][0])*Q_stokes-np.cos(2*globals()["theta"][0])*U_stokes))
dU_dtheta2 = 2.*pol_eff[1]/N*(np.sin(2.*globals()["theta"][1])*(pol_flux_corr[2]-pol_flux_corr[0]) - (pol_eff[0]*np.cos(-2.*globals()
["theta"][0]+2.*globals()["theta"][1]) - pol_eff[2]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2]))*U_stokes +
coeff_stokes_corr[2, 1]*(np.sin(2.*globals()["theta"][1])*Q_stokes-np.cos(2*globals()["theta"][1])*U_stokes))
dU_dtheta3 = 2.*pol_eff[2]/N*(np.sin(2.*globals()["theta"][2])*(pol_flux_corr[0]-pol_flux_corr[1]) - (pol_eff[1]*np.cos(-2.*globals()
["theta"][1]+2.*globals()["theta"][2]) - pol_eff[0]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0]))*U_stokes +
coeff_stokes_corr[2, 2]*(np.sin(2.*globals()["theta"][2])*Q_stokes-np.cos(2*globals()["theta"][2])*U_stokes))
dU_dtheta = np.array([dU_dtheta1, dU_dtheta2, dU_dtheta3])
# Compute the uncertainty associated with the polarizers' orientation (see Kishimoto 1999)
@@ -1228,12 +1239,10 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale
QU_diluted_err = np.sqrt(np.sum(Stokes_cov[1, 2][mask]**2))
P_diluted = np.sqrt(Q_diluted**2+U_diluted**2)/I_diluted
P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted **
2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err)
P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted**2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err)
PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted, Q_diluted))
PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err **
2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err)
PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err**2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err)
for header in headers:
header['P_int'] = (P_diluted, 'Integrated polarization degree')