Skip to content

Add the RON on the the contrast curves #15

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jmilou wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Add the RON on the the contrast curves #15

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
66ef028
improvements to the contrast curves
jmilou Mar 22, 2023
7982aff
Merge branch 'master' of https://github.com/robvanholstein/IRDAP
jmilou Jul 12, 2023
678b062
saving derotation files
jmilou Mar 31, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
118 changes: 113 additions & 5 deletions irdap/irdap.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1124,6 +1124,11 @@ def check_sort_data_create_directories(frames_to_remove=[],
fh.write('%s' % [number_frames_Q, number_frames_U])
printandlog('\nWrote file ' + os.path.join(path_preprocessed_dir, name_file_root + 'number_frames_QU.txt') + '.', wrap=False)

# Write the OBJECT DIT to a .txt-file
with open(os.path.join(path_preprocessed_dir, name_file_root + 'DIT.txt'), 'w') as fh:
fh.write('%s' % object_exposure_time)
printandlog('\nWrote file ' + os.path.join(path_preprocessed_dir, name_file_root + 'DIT.txt') + '.', wrap=False)

return path_dark_files, path_flat_files, path_object_files, path_sky_files, path_center_files, \
path_object_center_files, path_flux_files, path_sky_flux_files, \
indices_to_remove_dark, indices_to_remove_flat, indices_to_remove_object, \
Expand Down Expand Up @@ -4948,10 +4953,12 @@ def compute_contrast(frame, filter_used, sigma_clip=True):
###############################################################################

def plot_contrast_extended_source(path_table_star_flux, frame_I_Q, frame_I_U, frame_Q, frame_U,
filter_used, number_frames_Q, number_frames_U, QU_QUphi='QU',
filter_used, number_frames_Q, number_frames_U,
object_exposure_time, QU_QUphi='QU',
sigma_clip=True, star_flux_in_jansky=None):
'''
Plot the contrast curve for an extended source in contrast/arcsec^2 or Jansky/arcsec^2
and save the result in a csv file

Input:
path_table_star_flux: path to CSV-file with star flux data
Expand All @@ -4977,6 +4984,10 @@ def plot_contrast_extended_source(path_table_star_flux, frame_I_Q, frame_I_U, fr
frame_Q = frame_Q * irdis_gain * number_frames_Q
frame_U = frame_U * irdis_gain * number_frames_U

# Express the ron in total number of electrons per pixel
ron_Q_e_per_px = get_ron(object_exposure_time) * np.sqrt(number_frames_Q)
ron_U_e_per_px = get_ron(object_exposure_time) * np.sqrt(number_frames_U)

# Compute mean and standard deviations as a function of separation
separation, I_Q_mean, _, fwhm = compute_contrast(frame=frame_I_Q, filter_used=filter_used, sigma_clip=sigma_clip)
_, I_U_mean, _, _ = compute_contrast(frame=frame_I_U, filter_used=filter_used, sigma_clip=sigma_clip)
Expand Down Expand Up @@ -5020,6 +5031,10 @@ def plot_contrast_extended_source(path_table_star_flux, frame_I_Q, frame_I_U, fr
Q_std = Q_std * star_flux_in_jansky / (reference_flux_Q * pixel_scale**2)
U_std = U_std * star_flux_in_jansky / (reference_flux_U * pixel_scale**2)

# Express RON in contrast (or Jansky) per arcsec^2
ron_Q = ron_Q_e_per_px * star_flux_in_jansky / (reference_flux_Q * pixel_scale**2)
ron_U = ron_U_e_per_px * star_flux_in_jansky / (reference_flux_Q * pixel_scale**2)

# Compute separation in arcseconds
separation_as = separation * pixel_scale

Expand Down Expand Up @@ -5054,6 +5069,8 @@ def plot_contrast_extended_source(path_table_star_flux, frame_I_Q, frame_I_U, fr
plt.plot(separation_as, U_std, 'r-', label=r'$1\sigma$ ' + label_U)
plt.plot(separation_as, photon_noise_I_Q, 'b:', label=r'photon noise $I_Q$')
plt.plot(separation_as, photon_noise_I_U, 'r:', label=r'photon noise $I_U$')
plt.axhline(ron_Q,0,xlim_max,color='blue',linestyle='-.',label=r'RON noise $I_Q$')
plt.axhline(ron_U,0,xlim_max,color='red',linestyle='-.',label=r'RON noise $I_U$')
plt.yscale('log')
plt.xlabel('Separation (arcsec)')
plt.ylabel(ylabel)
Expand All @@ -5066,15 +5083,30 @@ def plot_contrast_extended_source(path_table_star_flux, frame_I_Q, frame_I_U, fr
plt.close()
printandlog(os.path.join(path_pdi_figures_dir, plot_name), wrap=False)

# Save contrast curve in a csv file
csv_name = name_file_root + 'contrast_extended_source.csv'
printandlog('\nCreating csv contrast curve for the detection of extended sources.')
dico_contrast_extended_source = {'Separation [arcsec]':separation_as,\
'I_Q_mean ['+ylabel+']':I_Q_mean,\
'I_U_mean ['+ylabel+']':I_U_mean,\
'1sigma ' + label_Q:Q_std,\
'1sigma ' + label_U:U_std,\
'photon noise I_Q':photon_noise_I_Q,\
'photon noise I_U':photon_noise_I_U,\
'RON noise I_Q':ron_Q,'RON noise I_U':ron_U}
df_contrast_extended_source = pd.DataFrame(dico_contrast_extended_source)
df_contrast_extended_source.to_csv(os.path.join(path_pdi_figures_dir, csv_name),\
index=False)

###############################################################################
# plot_contrast_point_source
###############################################################################

def plot_contrast_point_source(path_table_star_flux, path_flux_left, path_flux_right, frame_I_Q,
frame_I_U, frame_Q, frame_U, filter_used, number_frames_Q,
number_frames_U, QU_QUphi='QU', sigma_clip=True):
number_frames_U,object_exposure_time, QU_QUphi='QU', sigma_clip=True):
'''
Plot the contrast curve for a point source
Plot the contrast curve for a point source and save the result in a csv file

Input:
path_table_star_flux: path to CSV-file with star flux data
Expand All @@ -5099,6 +5131,10 @@ def plot_contrast_point_source(path_table_star_flux, path_flux_left, path_flux_r
frame_I_U = frame_I_U * irdis_gain * number_frames_U
frame_Q = frame_Q * irdis_gain * number_frames_Q
frame_U = frame_U * irdis_gain * number_frames_U

# Express the ron in total number of electrons per pixel
ron_Q_e_per_px = get_ron(object_exposure_time) * np.sqrt(number_frames_Q)
ron_U_e_per_px = get_ron(object_exposure_time) * np.sqrt(number_frames_U)

# Compute mean and standard deviations as a function of separation
separation, I_Q_mean, _, fwhm = compute_contrast(frame=frame_I_Q, filter_used=filter_used, sigma_clip=sigma_clip)
Expand Down Expand Up @@ -5158,6 +5194,10 @@ def plot_contrast_point_source(path_table_star_flux, path_flux_left, path_flux_r
Q_std = Q_std / reference_flux_Q
U_std = U_std / reference_flux_U

# Express RON as contrast per element of resolution
ron_Q = ron_Q_e_per_px * (np.sqrt(np.pi)*fwhm/2) / reference_flux_Q
ron_U = ron_U_e_per_px * (np.sqrt(np.pi)*fwhm/2) / reference_flux_U

# Compute separation in arcseconds
separation_as = separation * pixel_scale

Expand Down Expand Up @@ -5192,6 +5232,8 @@ def plot_contrast_point_source(path_table_star_flux, path_flux_left, path_flux_r
plt.plot(separation_as, U_std, 'r-', label=r'$1\sigma$ ' + label_U)
plt.plot(separation_as, photon_noise_I_Q, 'b:', label=r'photon noise $I_Q$')
plt.plot(separation_as, photon_noise_I_U, 'r:', label=r'photon noise $I_U$')
plt.axhline(ron_Q,0,xlim_max,color='blue',linestyle='-.',label=r'RON noise $I_Q$')
plt.axhline(ron_U,0,xlim_max,color='red',linestyle='-.',label=r'RON noise $I_U$')
plt.yscale('log')
plt.xlabel('Separation (arcsec)')
plt.ylabel('Point-source contrast')
Expand All @@ -5204,6 +5246,21 @@ def plot_contrast_point_source(path_table_star_flux, path_flux_left, path_flux_r
plt.close()
printandlog(os.path.join(path_pdi_figures_dir, plot_name), wrap=False)

# Save contrast curve in a csv file
csv_name = name_file_root + 'contrast_point_source.csv'
printandlog('\nCreating csv contrast curve for the detection of polarized point sources.')
dico_contrast_point_source = {'Separation [arcsec]':separation_as,\
'I_Q_mean [contrast]':I_Q_mean,\
'I_U_mean [contrast]':I_U_mean,\
'1sigma ' + label_Q:Q_std,\
'1sigma ' + label_U:U_std,\
'photon noise I_Q':photon_noise_I_Q,\
'photon noise I_U':photon_noise_I_U,\
'RON noise I_Q':ron_Q,'RON noise I_U':ron_U}
df_contrast_point_source = pd.DataFrame(dico_contrast_point_source)
df_contrast_point_source.to_csv(os.path.join(path_pdi_figures_dir, csv_name),\
index=False)

###############################################################################
# apply_pdi
###############################################################################
Expand Down Expand Up @@ -5432,6 +5489,10 @@ def apply_pdi(cube_left_frames,
number_frames_U = number_frames_QU[1]
printandlog('\nRead number of frames in Q and U from ' + path_number_frames + '.')

# Read DIT from .txt-file
path_dit = os.path.join(path_preprocessed_dir, name_file_root + 'DIT.txt')
object_exposure_time = literal_eval([x for x in open(path_dit, 'r')][0])

# Compute normalized Stokes q and u, DoLP and AoLP in an annulus on the star
q_star, u_star, sigma_q_star, sigma_u_star = \
determine_star_polarization(cube_I_Q=frame_I_Q_background_subtracted,
Expand Down Expand Up @@ -5724,7 +5785,8 @@ def apply_pdi(cube_left_frames,
frame_I_U_background_subtracted,
frame_Q_phi_star_polarization_subtracted,
frame_U_phi_star_polarization_subtracted,
filter_used, number_frames_Q, number_frames_U, QU_QUphi='QUphi',
filter_used, number_frames_Q, number_frames_U,
object_exposure_time,QU_QUphi='QUphi',
sigma_clip=True, star_flux_in_jansky=None)


Expand All @@ -5734,7 +5796,8 @@ def apply_pdi(cube_left_frames,
frame_I_U_background_subtracted,
frame_Q_star_polarization_subtracted,
frame_U_star_polarization_subtracted,
filter_used, number_frames_Q, number_frames_U, QU_QUphi='QU',
filter_used, number_frames_Q, number_frames_U,
object_exposure_time,QU_QUphi='QU',
sigma_clip=True)

else:
Expand Down Expand Up @@ -5921,6 +5984,9 @@ def apply_adi(cube_left_frames, cube_right_frames,header, principal_components='
printandlog('\nWARNING, the total parallactic rotation is limited to {0:.1f} deg, which will result in severe self-subtraction.'.format(amplitude_parang_rotation))
else:
printandlog('\nThe total parallactic rotation is {0:.1f} deg.'.format(amplitude_parang_rotation))
write_fits_files(data=derotation_angle,
path=os.path.join(path_adi_dir, name_file_root + 'derotation_angles.fits'),
header=False)

# Define strings to use in file names
string_principal_components = '-'.join(map(str, principal_components))
Expand Down Expand Up @@ -6414,6 +6480,44 @@ def create_overview_headers_main(path_main_dir):
print_wrap('\nCreated overview of the headers ' + path_overview + ' and ' + \
path_overview.replace('.txt','.csv') + '.')

###############################################################################
# compute the ron
###############################################################################

def get_ron(dit):
"""
Return the ReadOut Noise of the IRDIS detector. This function was copied from
the ETC system parameter document (ETC-system-param_NEW_2014-08-14.xls)
The IRDIS Read-out noise is 10 e- rms/pix/readout for a normal double correlated
and can go down to 3 e- for non destructive sampling (as measured in lab)

Parameters
----------
dit : float
Detector Integration Time.

Returns
-------
ron in e-

File written by Julien Milli
Function status: in dev
"""
return 10.
# When using the nondest readout mode, the RON should normally be lower when
# the integration time is larger, because several readouts are done
# like on JWST. Therefore, the previous implementation of this function
# was as follow. However, this did not match the contrast floor at large separations
# so this is probably not correct and a constant 10e- RON was found to
# better match the observations.
#
# if dit<7:
# return irdis_ron
# else:
# min_dit = 0.84
# return np.max([irdis_ron/np.sqrt((np.round(dit/min_dit,decimals=0)+1)/2),2.5])


###############################################################################
# run_pipeline
###############################################################################
Expand Down Expand Up @@ -6451,6 +6555,7 @@ def run_pipeline(path_main_dir):
global pixel_scale
global msd
global irdis_gain
global irdis_ron
global path_raw_dir
global path_flat_dir
global path_bpm_dir
Expand All @@ -6462,6 +6567,7 @@ def run_pipeline(path_main_dir):
global path_pdi_no_subtr_dir
global path_pdi_subtr_dir
global path_pdi_figures_dir
global path_adi_dir
global path_adi_classical_dir
global path_adi_pca_dir
global name_file_root
Expand Down Expand Up @@ -6489,6 +6595,8 @@ def run_pipeline(path_main_dir):

# Define gain of IRDIS detector (e–/ADU) (SPHERE User Manual, 14th public release, P105 Phase 1)
irdis_gain = 1.75
# Define gain of IRDIS detector (e–) (see also the function get_ron defined below)
irdis_ron = 10

# Define paths of directories
path_raw_dir = os.path.join(path_main_dir, 'raw')
Expand Down