notebook.community

Edit and run 





In [1]:



    


# Makes print and division act like Python 3
from __future__ import print_function, division

# Import the usual libraries
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

# Enable inline plotting
%matplotlib inline

from IPython.display import display, Latex, clear_output
from matplotlib.backends.backend_pdf import PdfPages



    











In [2]:



    


import pynrc
from pynrc import nrc_utils
from pynrc.nrc_utils import S, source_spectrum
# from pynrc.obs_nircam import model_to_hdulist, obs_hci
# from pynrc.obs_nircam import plot_contrasts, plot_contrasts_mjup, planet_mags, plot_planet_patches

pynrc.setup_logging('WARNING', verbose=False)



    


















    






[     pynrc:INFO]   jwst_backgrounds is not installed and will not be used for bg estimates.

















In [3]:



    


# Observation Definitions
from pynrc.nb_funcs import make_key, model_info, obs_wfe, obs_optimize

# Functions to run a series of operations
from pynrc.nb_funcs import do_opt, do_contrast, do_gen_hdus, do_sat_levels

# Plotting routines
from pynrc.nb_funcs import do_plot_contrasts, plot_images, plot_images_swlw

Define Sources and their Reference PSF Stars





In [4]:



    


# Various Bandpasses
bp_v = S.ObsBandpass('v')
bp_k = pynrc.bp_2mass('k')
bp_w1 = pynrc.bp_wise('w1')
bp_w2 = pynrc.bp_wise('w2')



    











In [5]:



    


#                source,      dist, age, sptype, Teff, [Fe/H], log_g, mag, band, fov
args_sources = [('HD10647',  17.34, 1400, 'F9V', 5954, +0.00, 4.7, 4.34, bp_k, 13),
                ('HD107146', 27.47,  200, 'G2V', 5850, +0.00, 4.5, 5.54, bp_k, 13),
                ('HD181327', 48.21,   12, 'F6V', 6449, +0.29, 4.4, 5.91, bp_k,  7),
                ('HD61005',  36.49,  100, 'G8V', 5500, +0.00, 4.5, 6.45, bp_k,  7),
                ('HD32297', 132.79,   30, 'A7V', 7800, -0.76, 3.8, 7.59, bp_k,  5),
               ]

ref_sources  = [('iotHor',   'F8V',  6080, +0.15, 4.5, 4.14, bp_k),
                ('HD111398', 'G5V',  5689, +0.07, 4.5, 5.53, bp_k),
                ('HR7297',   'F7V',  6500, -0.10, 4.2, 5.10, bp_k),
                ('HD56161',  'G5IV', 5337, +0.00, 4.3, 4.91, bp_k),
                ('HD31411',  'A0V',  9500, +0.00, 4.0, 6.42, bp_k)]



    











In [6]:



    


# Directory housing VOTables 
# http://vizier.u-strasbg.fr/vizier/sed/
votdir = 'votables/'

# Directory to save plots and figures
outdir  = 'DebrisDisks/'



    











In [7]:



    


# List of filters
args_filter = [('F182M', 'MASK335R', 'CIRCLYOT'),
               ('F210M', 'MASK335R', 'CIRCLYOT'),
               ('F250M', 'MASK335R', 'CIRCLYOT'),
               ('F300M', 'MASK335R', 'CIRCLYOT'),
               ('F335M', 'MASK335R', 'CIRCLYOT'),
               ('F444W', 'MASK335R', 'CIRCLYOT')]

subsize = 320

#args_filter = [('F335M', 'MASK335R', 'CIRCLYOT'),
#               ('F444W', 'MASK335R', 'CIRCLYOT')]

filt_keys = []
for filt,mask,pupil in args_filter:
    filt_keys.append(make_key(filt, mask=mask, pupil=pupil))

HD 10647





In [8]:



    


# Fit spectrum to SED photometry
i=0
name_sci, dist_sci, age, spt_sci, Teff_sci, feh_sci, logg_sci, mag_sci, bp_sci, fov = args_sources[i]
vot = votdir + name_sci.replace(' ' ,'') + '.vot'

args = (name_sci, spt_sci, mag_sci, bp_sci, vot)
kwargs = {'Teff':Teff_sci, 'metallicity':feh_sci, 'log_g':logg_sci}
src = source_spectrum(*args, **kwargs)

src.fit_SED(use_err=False, robust=True, wlim=[1,4])

# Final source spectrum
sp_sci = src.sp_model



    


















    






Warning, 1 of 43 bins contained negative fluxes; they have been set to zero.
[1.0135502]

















In [9]:



    


# Do the same for the reference source
name_ref, spt_ref, Teff_ref, feh_ref, logg_ref, mag_ref, bp_ref = ref_sources[i]
vot = votdir + name_ref.replace(' ' ,'') + '.vot'

args = (name_ref, spt_ref, mag_ref, bp_ref, vot)
kwargs = {'Teff':Teff_ref, 'metallicity':feh_ref, 'log_g':logg_ref}
ref = nrc_utils.source_spectrum(*args, **kwargs)

ref.fit_SED(use_err=True, robust=True)

# Final reference spectrum
sp_ref = ref.sp_model



    


















    






Warning, 1 of 31 bins contained negative fluxes; they have been set to zero.
[1.01410814]

















In [10]:



    


# Plot spectra 

fig, axes = plt.subplots(1,2, figsize=(14,4.5))
src.plot_SED(ax=axes[0], xr=[0.5,30])
ref.plot_SED(ax=axes[1], xr=[0.5,30])

axes[0].set_title('Science Specta -- {} ({})'.format(src.name, spt_sci))
axes[1].set_title('Reference Specta -- {} ({})'.format(ref.name, spt_ref))

#for ax in axes:
#    ax.set_xscale('linear')
#    ax.xaxis.set_minor_locator(AutoMinorLocator())

fig.tight_layout()
fig.subplots_adjust(top=0.85, bottom=0.1 , left=0.05, right=0.97)

fig.savefig(outdir+'{}_SEDs.pdf'.format(name_sci.replace(' ','')))



    


















    
























In [11]:



    


# Plot the two spectra
fig, ax = plt.subplots(1,1, figsize=(8,5))

xr = [2.5,5.5]

bp = pynrc.read_filter(*args_filter[-1])
for sp in [sp_sci, sp_ref]:
    w = sp.wave / 1e4
    o = S.Observation(sp, bp, binset=bp.wave)
    sp.convert('Jy')
    f = sp.flux / o.effstim('Jy')

    ind = (w>=xr[0]) & (w<=xr[1])
    ax.plot(w[ind], f[ind], lw=1, label=sp.name)
    ax.set_ylabel('Flux (Jy) normalized over bandpass')
    sp.convert('flam')

ax.set_xlim(xr)
ax.set_ylim([0,ax.get_ylim()[1]])
ax.set_xlabel(r'Wavelength ($\mu m$)')
ax.set_title('{} Spectra'.format(sp_sci.name))

# Overplot Filter Bandpass
ax2 = ax.twinx()
ax2.plot(bp.wave/1e4, bp.throughput, color='C2', label=bp.name+' Bandpass')
ax2.set_ylim([0,0.8])
ax2.set_xlim(xr)
ax2.set_ylabel('Bandpass Throughput')

ax.legend(loc='upper left')
ax2.legend(loc='upper right')

fig.tight_layout()
fig.savefig(outdir+'{}_2SEDs.pdf'.format(name_sci.replace(' ','')))



    


















    
























In [12]:



    


# Create a dictionary that holds the obs_coronagraphy class for each filter
wfe_drift = 0
obs_dict = obs_wfe(wfe_drift, args_filter, sp_sci, dist_sci, sp_ref=sp_ref, args_disk='auto', 
                   wind_mode='WINDOW', subsize=subsize, verbose=False)



    


















    






F182M_MASK335R_CIRCLYOT
F210M_MASK335R_CIRCLYOT
F250M_MASK335R_CIRCLYOT
F300M_MASK335R_CIRCLYOT
F335M_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [13]:



    


# Optimize readout parameters
tacq = 4200
# patterns=['BRIGHT2']
# do_opt(tacq, patterns=patterns, ng_min=5, ng_max=10, tacq_frac=0.1, well_levels=[2], even_nints=True)



    











In [14]:



    


# Update multiaccum info
for key in filt_keys:
    obs = obs_dict[key]

    read_mode='MEDIUM8'
    ng, nint_sci, nint_ref = (8,21,20)
    # Double up for SW, because fewer number of paired fitlers
    if obs.bandpass.avgwave()/1e4 < 2.5: 
        nint_sci *= 2
        nint_ref *= 2
    
    obs.update_detectors(read_mode=read_mode, ngroup=ng, nint=nint_sci)
    obs.nrc_ref.update_detectors(read_mode=read_mode, ngroup=ng, nint=nint_ref)

#     print(key)
#     print(obs.multiaccum_times)
#     _ = obs.sensitivity(nsig=5, units='vegamag', verbose=True)
#     print('')



    











In [15]:



    


# Max Saturation Values
dmax = []
for k in filt_keys:
    print(k)
    obs = obs_dict[k]
    
    imsci = obs.gen_slope_image(exclude_noise=True, quick_PSF=True)
    imref = obs.gen_slope_image(exclude_noise=True, quick_PSF=True, do_ref=True)
    
    ng1 = 2
    ng2 = obs.multiaccum.ngroup
    sci_sat1 = obs.saturation_levels(ngroup=ng1, image=imsci)
    sci_sat2 = obs.saturation_levels(ngroup=ng2, image=imsci)
    ref_sat1 = obs.saturation_levels(ngroup=ng1, image=imref, do_ref=True)
    ref_sat2 = obs.saturation_levels(ngroup=ng2, image=imref, do_ref=True)
    
    print('Max Well NG={}: {:.2f} {:.2f};   Max Well NG={}: {:.2f} {:.2f}'\
          .format(ng1,sci_sat1.max(),ref_sat1.max(),ng2,sci_sat2.max(),ref_sat2.max()))



    


















    






F182M_MASK335R_CIRCLYOT
Max Well NG=2: 0.54 0.65;   Max Well NG=8: 2.34 2.82
F210M_MASK335R_CIRCLYOT
Max Well NG=2: 0.60 0.73;   Max Well NG=8: 2.62 3.15
F250M_MASK335R_CIRCLYOT
Max Well NG=2: 0.32 0.38;   Max Well NG=8: 1.38 1.67
F300M_MASK335R_CIRCLYOT
Max Well NG=2: 0.23 0.27;   Max Well NG=8: 0.99 1.19
F335M_MASK335R_CIRCLYOT
Max Well NG=2: 0.18 0.21;   Max Well NG=8: 0.77 0.92
F444W_MASK335R_CIRCLYOT
Max Well NG=2: 0.12 0.15;   Max Well NG=8: 0.52 0.63

















In [16]:



    


# Disk Images
wfe_drift = 5
wfe_roll_drift = 2
hdu_dict = do_gen_hdus(obs_dict, filt_keys, wfe_drift, wfe_roll_drift, 
                       PA1=-5, PA2=+5, opt_diff=False, use_cmask=True,
                       verbose=True)



    


















    






F182M_MASK335R_CIRCLYOT
F210M_MASK335R_CIRCLYOT
F250M_MASK335R_CIRCLYOT
F300M_MASK335R_CIRCLYOT
F335M_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [17]:



    


# Save FITS files
for k in filt_keys:
    hdul = hdu_dict[k]
    hdul.writeto(outdir+'{}_{}.fits'.format(name_sci,k), overwrite=True)



    











In [18]:



    


save_fig = True
if len(filt_keys)==2:
    plot_images(obs_dict, hdu_dict, filt_keys, wfe_drift, fov=fov, save_fig=save_fig, outdir=outdir)
else:
    plot_images_swlw(obs_dict, hdu_dict, filt_keys, wfe_drift, fov=fov, save_fig=save_fig, outdir=outdir)



    


















    
























In [19]:



    


# Determine contrast curves for various WFE drift values
wfe_list = [0,2,5,10]
nsig = 5
roll = 10

# (Roll1 - Ref) + (Roll2 - Ref)
curves_dict = do_contrast(obs_dict, wfe_list, filt_keys[-1:], nsig=nsig, roll_angle=roll, exclude_disk=False)
curves_F444W = curves_dict[filt_keys[-1]]

# Roll1 - Roll2
curves_dict = do_contrast(obs_dict, wfe_list, filt_keys[-1:], nsig=nsig, roll_angle=roll, exclude_disk=False,
                          no_ref=True)
curves_F444W2 = curves_dict[filt_keys[-1]]



    


















    






F444W_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [20]:



    


sat_rad = 0
obs = obs_dict[filt_keys[-1]]
do_plot_contrasts(curves_F444W, curves_F444W2, nsig, wfe_list, obs, age, 
                  save_fig=True, outdir=outdir, yr2=[1,100])

TODO: EVERYTHING BELOW HERE NEEDS UPDATING (5/5/2019)

HD 107146





In [21]:



    


# Fit spectrum to SED photometry
from pynrc.nrc_utils import source_spectrum

i=1
name_sci, dist_sci, age, spt_sci, Teff_sci, feh_sci, logg_sci, mag_sci, bp_sci, fov = args_sources[i]
vot = votdir + name_sci.replace(' ' ,'') + '.vot'

args = (name_sci, spt_sci, mag_sci, bp_sci, vot)
kwargs = {'Teff':Teff_sci, 'metallicity':feh_sci, 'log_g':logg_sci}
src = source_spectrum(*args, **kwargs)

src.fit_SED(use_err=False, robust=True, wlim=[1,4])

# Final source spectrum
sp_sci = src.sp_model



    


















    






[1.00911793]

















In [22]:



    


# Do the same for the reference source
name_ref, spt_ref, Teff_ref, feh_ref, logg_ref, mag_ref, bp_ref = ref_sources[i]
vot = votdir + name_ref.replace(' ' ,'') + '.vot'

args = (name_ref, spt_ref, mag_ref, bp_ref, vot)
kwargs = {'Teff':Teff_ref, 'metallicity':feh_ref, 'log_g':logg_ref}
ref = nrc_utils.source_spectrum(*args, **kwargs)

ref.fit_SED(use_err=True, robust=True)

# Final reference spectrum
sp_ref = ref.sp_model



    


















    






[1.04390456]

















In [23]:



    


# Plot spectra 

fig, axes = plt.subplots(1,2, figsize=(14,4.5))
src.plot_SED(ax=axes[0], xr=[0.5,30])
ref.plot_SED(ax=axes[1], xr=[0.5,30])

axes[0].set_title('Science Specta -- {} ({})'.format(src.name, spt_sci))
axes[1].set_title('Reference Specta -- {} ({})'.format(ref.name, spt_ref))

#for ax in axes:
#    ax.set_xscale('linear')
#    ax.xaxis.set_minor_locator(AutoMinorLocator())

fig.tight_layout()
fig.subplots_adjust(top=0.85, bottom=0.1 , left=0.05, right=0.97)

fig.savefig(outdir+'{}_SEDs.pdf'.format(name_sci.replace(' ','')))



    


















    
























In [24]:



    


# Plot the two spectra
fig, ax = plt.subplots(1,1, figsize=(8,5))

xr = [2.5,5.5]

for sp in [sp_sci, sp_ref]:
    w = sp.wave / 1e4
    ind = (w>=xr[0]) & (w<=xr[1])
    ind2 = (w>=3.9) & (w<=4.1)
    sp.convert('Jy')
    f = sp.flux / np.mean(sp.flux[ind2])
    ax.plot(w[ind], f[ind], lw=1, label=sp.name)
    ax.set_ylabel('Flux (Jy) normalized at 4 $\mu m$')
    sp.convert('flam')

ax.set_xlim(xr)
ax.set_ylim([0,ax.get_ylim()[1]])
ax.set_xlabel(r'Wavelength ($\mu m$)')
ax.set_title('{} Spectra'.format(sp_sci.name))

# Overplot Filter Bandpass
bp = pynrc.read_filter(*args_filter[-1])
ax2 = ax.twinx()
ax2.plot(bp.wave/1e4, bp.throughput, color='C2', label=bp.name+' Bandpass')
ax2.set_ylim([0,0.8])
ax2.set_xlim(xr)
ax2.set_ylabel('Bandpass Throughput')

ax.legend(loc='upper left')
ax2.legend(loc='upper right')

fig.tight_layout()
fig.savefig(outdir+'{}_2SEDs.pdf'.format(name_sci.replace(' ','')))



    


















    
























In [25]:



    


# Create a dictionary that holds the obs_coronagraphy class for each filter
wfe_drift = 0
obs_dict = obs_wfe(wfe_drift, args_filter, sp_sci, dist_sci, sp_ref=sp_ref, args_disk='auto', 
                   wind_mode='WINDOW', subsize=subsize, verbose=False)



    


















    






F182M_MASK335R_CIRCLYOT
F210M_MASK335R_CIRCLYOT
F250M_MASK335R_CIRCLYOT
F300M_MASK335R_CIRCLYOT
F335M_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [26]:



    


# Optimize readout parameters
tacq = 4200
#do_opt(tacq, patterns='MEDIUM8', ng_min=8, ng_max=8, tacq_frac=0.1, well_levels=[2], even_nints=True)



    











In [27]:



    


for key in filt_keys:
    obs = obs_dict[key]

    read_mode='MEDIUM8'
    ng, nint_sci, nint_ref = (10,17,15)
    if obs.bandpass.avgwave()/1e4 < 2.5: 
        nint_sci *= 2
        nint_ref *= 2
    
    obs.update_detectors(read_mode=read_mode, ngroup=ng, nint=nint_sci)
    obs.nrc_ref.update_detectors(read_mode=read_mode, ngroup=ng, nint=nint_ref)

    #print(key)
    #print(obs.multiaccum_times)
    #_ = obs.sensitivity(nsig=5, units='vegamag', verbose=True)
    #print('')



    











In [28]:



    


# Saturation Levels
# Only want F444W
#obs = obs_dict[filt_keys[-1]]
#sat_rad = do_sat_levels(obs, plot=True)

# Max Saturation Values
dmax = []
for k in filt_keys:
    print(k)
    obs = obs_dict[k]
    
    imsci = obs.gen_slope_image(exclude_noise=True, quick_PSF=True)
    imref = obs.gen_slope_image(exclude_noise=True, quick_PSF=True, do_ref=True)
    
    ng1 = 2
    ng2 = obs.multiaccum.ngroup
    sci_sat1 = obs.saturation_levels(ngroup=ng1, image=imsci)
    sci_sat2 = obs.saturation_levels(ngroup=ng2, image=imsci)
    ref_sat1 = obs.saturation_levels(ngroup=ng1, image=imref, do_ref=True)
    ref_sat2 = obs.saturation_levels(ngroup=ng2, image=imref, do_ref=True)
    
    print('Max Sat NG={}: {:.2f} {:.2f};   Max Sat NG={}: {:.2f} {:.2f}'\
          .format(ng1,sci_sat1.max(),ref_sat1.max(),ng2,sci_sat2.max(),ref_sat2.max()))



    


















    






F182M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.18 0.19;   Max Sat NG=10: 0.97 1.01
F210M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.20 0.21;   Max Sat NG=10: 1.09 1.13
F250M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.11 0.11;   Max Sat NG=10: 0.57 0.59
F300M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.08 0.08;   Max Sat NG=10: 0.41 0.43
F335M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.06 0.06;   Max Sat NG=10: 0.32 0.33
F444W_MASK335R_CIRCLYOT
Max Sat NG=2: 0.04 0.04;   Max Sat NG=10: 0.22 0.22

















In [29]:



    


# Disk Images
wfe_drift = 5
wfe_roll_drift = 2
hdu_dict = do_gen_hdus(obs_dict, filt_keys, wfe_drift, wfe_roll_drift, 
                       PA1=-5, PA2=+5, opt_diff=False, use_cmask=True,
                       verbose=True)



    


















    






F182M_MASK335R_CIRCLYOT
F210M_MASK335R_CIRCLYOT
F250M_MASK335R_CIRCLYOT
F300M_MASK335R_CIRCLYOT
F335M_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [30]:



    


# Save FITS files
for k in filt_keys:
    hdul = hdu_dict[k]
    hdul.writeto(outdir+'{}_{}.fits'.format(name_sci,k), overwrite=True)



    











In [31]:



    


save_fig = True
if len(filt_keys)==2:
    plot_images(obs_dict, hdu_dict, filt_keys, wfe_drift, fov=fov, save_fig=save_fig, outdir=outdir)
else:
    plot_images_swlw(obs_dict, hdu_dict, filt_keys, wfe_drift, fov=fov, save_fig=save_fig, outdir=outdir)



    


















    
























In [32]:



    


# Determine contrast curves for various WFE drift values
wfe_list = [0,2,5,10]
nsig = 5
roll = 10

# (Roll1 - Ref) + (Roll2 - Ref)
curves_dict = do_contrast(obs_dict, wfe_list, filt_keys[-1:], nsig=nsig, roll_angle=roll, exclude_disk=False)
curves_F444W = curves_dict[filt_keys[-1]]

# Roll1 - Roll2
curves_dict = do_contrast(obs_dict, wfe_list, filt_keys[-1:], nsig=nsig, roll_angle=roll, exclude_disk=False,
                          no_ref=True)
curves_F444W2 = curves_dict[filt_keys[-1]]



    


















    






F444W_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [33]:



    


sat_rad = 0
obs = obs_dict[filt_keys[-1]]
#do_plot_contrasts(curves_F444W, curves_F444W2, obs, save_fig=True)
do_plot_contrasts(curves_F444W, curves_F444W2, nsig, wfe_list, obs, age, 
                  save_fig=True, outdir=outdir, yr2=[0.1,100])

HD 181327





In [34]:



    


# Fit spectrum to SED photometry
from pynrc.nrc_utils import source_spectrum

i=2
name_sci, dist_sci, age, spt_sci, Teff_sci, feh_sci, logg_sci, mag_sci, bp_sci, fov = args_sources[i]
vot = votdir + name_sci.replace(' ' ,'') + '.vot'

args = (name_sci, spt_sci, mag_sci, bp_sci, vot)
kwargs = {'Teff':Teff_sci, 'metallicity':feh_sci, 'log_g':logg_sci}
src = source_spectrum(*args, **kwargs)

src.fit_SED(use_err=False, robust=True, wlim=[1,4])

# Final source spectrum
sp_sci = src.sp_model



    


















    






[0.97349543]

















In [35]:



    


# Do the same for the reference source
name_ref, spt_ref, Teff_ref, feh_ref, logg_ref, mag_ref, bp_ref = ref_sources[i]
vot = votdir + name_ref.replace(' ' ,'') + '.vot'

args = (name_ref, spt_ref, mag_ref, bp_ref, vot)
kwargs = {'Teff':Teff_ref, 'metallicity':feh_ref, 'log_g':logg_ref}
ref = nrc_utils.source_spectrum(*args, **kwargs)

ref.fit_SED(use_err=True, robust=True)

# Final reference spectrum
sp_ref = ref.sp_model



    


















    






[0.92729568]

















In [36]:



    


# Plot spectra 

fig, axes = plt.subplots(1,2, figsize=(14,4.5))
src.plot_SED(ax=axes[0], xr=[0.5,30])
ref.plot_SED(ax=axes[1], xr=[0.5,30])

axes[0].set_title('Science Specta -- {} ({})'.format(src.name, spt_sci))
axes[1].set_title('Reference Specta -- {} ({})'.format(ref.name, spt_ref))

#for ax in axes:
#    ax.set_xscale('linear')
#    ax.xaxis.set_minor_locator(AutoMinorLocator())

fig.tight_layout()
fig.subplots_adjust(top=0.85, bottom=0.1 , left=0.05, right=0.97)

fig.savefig(outdir+'{}_SEDs.pdf'.format(name_sci.replace(' ','')))



    


















    
























In [37]:



    


# Plot the two spectra
fig, ax = plt.subplots(1,1, figsize=(8,5))

xr = [2.5,5.5]

for sp in [sp_sci, sp_ref]:
    w = sp.wave / 1e4
    ind = (w>=xr[0]) & (w<=xr[1])
    ind2 = (w>=3.9) & (w<=4.1)
    sp.convert('Jy')
    f = sp.flux / np.mean(sp.flux[ind2])
    ax.plot(w[ind], f[ind], lw=1, label=sp.name)
    ax.set_ylabel('Flux (Jy) normalized at 4 $\mu m$')
    sp.convert('flam')

ax.set_xlim(xr)
ax.set_ylim([0,ax.get_ylim()[1]])
ax.set_xlabel(r'Wavelength ($\mu m$)')
ax.set_title('{} Spectra'.format(sp_sci.name))

# Overplot Filter Bandpass
bp = pynrc.read_filter(*args_filter[-1])
ax2 = ax.twinx()
ax2.plot(bp.wave/1e4, bp.throughput, color='C2', label=bp.name+' Bandpass')
ax2.set_ylim([0,0.8])
ax2.set_xlim(xr)
ax2.set_ylabel('Bandpass Throughput')

ax.legend(loc='upper left')
ax2.legend(loc='upper right')

fig.tight_layout()
fig.savefig(outdir+'{}_2SEDs.pdf'.format(name_sci.replace(' ','')))



    


















    
























In [38]:



    


# Create a dictionary that holds the obs_coronagraphy class for each filter
wfe_drift = 0
obs_dict = obs_wfe(wfe_drift, args_filter, sp_sci, dist_sci, sp_ref=sp_ref, args_disk='auto', 
                   wind_mode='WINDOW', subsize=subsize, verbose=False)



    


















    






F182M_MASK335R_CIRCLYOT
F210M_MASK335R_CIRCLYOT
F250M_MASK335R_CIRCLYOT
F300M_MASK335R_CIRCLYOT
F335M_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [39]:



    


# Optimize readout parameters
tacq = 4200
#do_opt(tacq, patterns='MEDIUM8', ng_min=8, ng_max=8, tacq_frac=0.1, well_levels=[2], even_nints=True)



    











In [40]:



    


for key in filt_keys:
    obs = obs_dict[key]

    read_mode='MEDIUM8'
    ng, nint_sci, nint_ref = (10,11,15)
    if obs.bandpass.avgwave()/1e4 < 2.5: 
        nint_sci *= 2
        nint_ref *= 2
    
    obs.update_detectors(read_mode=read_mode, ngroup=ng, nint=nint_sci)
    obs.nrc_ref.update_detectors(read_mode=read_mode, ngroup=ng, nint=nint_ref)

    #print(key)
    #print(obs.multiaccum_times)
    #_ = obs.sensitivity(nsig=5, units='vegamag', verbose=True)
    #print('')



    











In [41]:



    


# Saturation Levels
# Only want F444W
#obs = obs_dict[filt_keys[-1]]
#sat_rad = do_sat_levels(obs, plot=True)

# Max Saturation Values
dmax = []
for k in filt_keys:
    print(k)
    obs = obs_dict[k]
    
    imsci = obs.gen_slope_image(exclude_noise=True, quick_PSF=True)
    imref = obs.gen_slope_image(exclude_noise=True, quick_PSF=True, do_ref=True)
    
    ng1 = 2
    ng2 = obs.multiaccum.ngroup
    sci_sat1 = obs.saturation_levels(ngroup=ng1, image=imsci)
    sci_sat2 = obs.saturation_levels(ngroup=ng2, image=imsci)
    ref_sat1 = obs.saturation_levels(ngroup=ng1, image=imref, do_ref=True)
    ref_sat2 = obs.saturation_levels(ngroup=ng2, image=imref, do_ref=True)
    
    print('Max Sat NG={}: {:.2f} {:.2f};   Max Sat NG={}: {:.2f} {:.2f}'\
          .format(ng1,sci_sat1.max(),ref_sat1.max(),ng2,sci_sat2.max(),ref_sat2.max()))



    


















    






F182M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.12 0.25;   Max Sat NG=10: 0.66 1.34
F210M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.14 0.27;   Max Sat NG=10: 0.74 1.50
F250M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.07 0.15;   Max Sat NG=10: 0.39 0.79
F300M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.05 0.10;   Max Sat NG=10: 0.28 0.56
F335M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.04 0.08;   Max Sat NG=10: 0.22 0.43
F444W_MASK335R_CIRCLYOT
Max Sat NG=2: 0.03 0.06;   Max Sat NG=10: 0.15 0.30

















In [42]:



    


# Disk Images
wfe_drift = 5
wfe_roll_drift = 2
hdu_dict = do_gen_hdus(obs_dict, filt_keys, wfe_drift, wfe_roll_drift, 
                       PA1=-5, PA2=+5, opt_diff=False, use_cmask=True,
                       verbose=True)



    


















    






F182M_MASK335R_CIRCLYOT
F210M_MASK335R_CIRCLYOT
F250M_MASK335R_CIRCLYOT
F300M_MASK335R_CIRCLYOT
F335M_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [43]:



    


# Save FITS files
for k in filt_keys:
    hdul = hdu_dict[k]
    hdul.writeto(outdir+'{}_{}.fits'.format(name_sci,k), overwrite=True)



    











In [44]:



    


save_fig = True
if len(filt_keys)==2:
    plot_images(obs_dict, hdu_dict, filt_keys, wfe_drift, fov=fov, save_fig=save_fig, outdir=outdir)
else:
    plot_images_swlw(obs_dict, hdu_dict, filt_keys, wfe_drift, fov=fov, save_fig=save_fig, outdir=outdir)



    


















    
























In [45]:



    


# Determine contrast curves for various WFE drift values
wfe_list = [0,2,5,10]
nsig = 5
roll = 10

# (Roll1 - Ref) + (Roll2 - Ref)
curves_dict = do_contrast(obs_dict, wfe_list, filt_keys[-1:], nsig=nsig, roll_angle=roll, exclude_disk=False)
curves_F444W = curves_dict[filt_keys[-1]]

# Roll1 - Roll2
curves_dict = do_contrast(obs_dict, wfe_list, filt_keys[-1:], nsig=nsig, roll_angle=roll, exclude_disk=False,
                          no_ref=True)
curves_F444W2 = curves_dict[filt_keys[-1]]



    


















    






F444W_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [46]:



    


sat_rad = 0
obs = obs_dict[filt_keys[-1]]
#do_plot_contrasts(curves_F444W, curves_F444W2, obs, save_fig=True, xr=[0,5], xr2=[0,5])
do_plot_contrasts(curves_F444W, curves_F444W2, nsig, wfe_list, obs, age, 
                  save_fig=True, outdir=outdir, xr=[0,5], xr2=[0,5], yr2=[0.01,100])

HD 61005





In [47]:



    


# Fit spectrum to SED photometry
from pynrc.nrc_utils import source_spectrum

i=3
name_sci, dist_sci, age, spt_sci, Teff_sci, feh_sci, logg_sci, mag_sci, bp_sci, fov = args_sources[i]
vot = votdir + name_sci.replace(' ' ,'') + '.vot'

args = (name_sci, spt_sci, mag_sci, bp_sci, vot)
kwargs = {'Teff':Teff_sci, 'metallicity':feh_sci, 'log_g':logg_sci}
src = source_spectrum(*args, **kwargs)

src.fit_SED(use_err=False, robust=True, wlim=[1,4])

# Final source spectrum
sp_sci = src.sp_model



    


















    






[0.96681562]

















In [48]:



    


# Do the same for the reference source
name_ref, spt_ref, Teff_ref, feh_ref, logg_ref, mag_ref, bp_ref = ref_sources[i]
vot = votdir + name_ref.replace(' ' ,'') + '.vot'

args = (name_ref, spt_ref, mag_ref, bp_ref, vot)
kwargs = {'Teff':Teff_ref, 'metallicity':feh_ref, 'log_g':logg_ref}
ref = nrc_utils.source_spectrum(*args, **kwargs)

ref.fit_SED(use_err=True, robust=True)

# Final reference spectrum
sp_ref = ref.sp_model



    


















    






[0.9238596]

















In [49]:



    


# Plot spectra 

fig, axes = plt.subplots(1,2, figsize=(14,4.5))
src.plot_SED(ax=axes[0], xr=[0.5,30])
ref.plot_SED(ax=axes[1], xr=[0.5,30])

axes[0].set_title('Science Specta -- {} ({})'.format(src.name, spt_sci))
axes[1].set_title('Reference Specta -- {} ({})'.format(ref.name, spt_ref))

#for ax in axes:
#    ax.set_xscale('linear')
#    ax.xaxis.set_minor_locator(AutoMinorLocator())

fig.tight_layout()
fig.subplots_adjust(top=0.85, bottom=0.1 , left=0.05, right=0.97)

fig.savefig(outdir+'{}_SEDs.pdf'.format(name_sci.replace(' ','')))



    


















    
























In [50]:



    


# Plot the two spectra
fig, ax = plt.subplots(1,1, figsize=(8,5))

xr = [2.5,5.5]

for sp in [sp_sci, sp_ref]:
    w = sp.wave / 1e4
    ind = (w>=xr[0]) & (w<=xr[1])
    ind2 = (w>=3.9) & (w<=4.1)
    sp.convert('Jy')
    f = sp.flux / np.mean(sp.flux[ind2])
    ax.plot(w[ind], f[ind], lw=1, label=sp.name)
    ax.set_ylabel('Flux (Jy) normalized at 4 $\mu m$')
    sp.convert('flam')

ax.set_xlim(xr)
ax.set_ylim([0,ax.get_ylim()[1]])
ax.set_xlabel(r'Wavelength ($\mu m$)')
ax.set_title('{} Spectra'.format(sp_sci.name))

# Overplot Filter Bandpass
bp = pynrc.read_filter(*args_filter[-1])
ax2 = ax.twinx()
ax2.plot(bp.wave/1e4, bp.throughput, color='C2', label=bp.name+' Bandpass')
ax2.set_ylim([0,0.8])
ax2.set_xlim(xr)
ax2.set_ylabel('Bandpass Throughput')

ax.legend(loc='upper left')
ax2.legend(loc='upper right')

fig.tight_layout()
fig.savefig(outdir+'{}_2SEDs.pdf'.format(name_sci.replace(' ','')))



    


















    
























In [51]:



    


# Create a dictionary that holds the obs_coronagraphy class for each filter
wfe_drift = 0
obs_dict = obs_wfe(wfe_drift, args_filter, sp_sci, dist_sci, sp_ref=sp_ref, args_disk='auto', 
                   wind_mode='WINDOW', subsize=subsize, verbose=False)



    


















    






F182M_MASK335R_CIRCLYOT
F210M_MASK335R_CIRCLYOT
F250M_MASK335R_CIRCLYOT
F300M_MASK335R_CIRCLYOT
F335M_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [52]:



    


# Optimize readout parameters
tacq = 4200
#do_opt(tacq, patterns='MEDIUM8', ng_min=8, ng_max=8, tacq_frac=0.1, well_levels=[2], even_nints=True)



    











In [53]:



    


for key in filt_keys:
    obs = obs_dict[key]

    read_mode='MEDIUM8'
    ng, nint_sci, nint_ref = (10,14,15)
    if obs.bandpass.avgwave()/1e4 < 2.5: 
        nint_sci *= 2
        nint_ref *= 2
    
    obs.update_detectors(read_mode=read_mode, ngroup=ng, nint=nint_sci)
    obs.nrc_ref.update_detectors(read_mode=read_mode, ngroup=ng, nint=nint_ref)

    #print(key)
    #print(obs.multiaccum_times)
    #_ = obs.sensitivity(nsig=5, units='vegamag', verbose=True)
    #print('')



    











In [54]:



    


# Saturation Levels
# Only want F444W
#obs = obs_dict[filt_keys[-1]]
#sat_rad = do_sat_levels(obs, plot=True)

# Max Saturation Values
dmax = []
for k in filt_keys:
    print(k)
    obs = obs_dict[k]
    
    imsci = obs.gen_slope_image(exclude_noise=True, quick_PSF=True)
    imref = obs.gen_slope_image(exclude_noise=True, quick_PSF=True, do_ref=True)
    
    ng1 = 2
    ng2 = obs.multiaccum.ngroup
    sci_sat1 = obs.saturation_levels(ngroup=ng1, image=imsci)
    sci_sat2 = obs.saturation_levels(ngroup=ng2, image=imsci)
    ref_sat1 = obs.saturation_levels(ngroup=ng1, image=imref, do_ref=True)
    ref_sat2 = obs.saturation_levels(ngroup=ng2, image=imref, do_ref=True)
    
    print('Max Sat NG={}: {:.2f} {:.2f};   Max Sat NG={}: {:.2f} {:.2f}'\
          .format(ng1,sci_sat1.max(),ref_sat1.max(),ng2,sci_sat2.max(),ref_sat2.max()))



    


















    






F182M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.07 0.29;   Max Sat NG=10: 0.40 1.58
F210M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.08 0.33;   Max Sat NG=10: 0.45 1.78
F250M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.04 0.17;   Max Sat NG=10: 0.24 0.93
F300M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.03 0.12;   Max Sat NG=10: 0.17 0.67
F335M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.02 0.10;   Max Sat NG=10: 0.13 0.52
F444W_MASK335R_CIRCLYOT
Max Sat NG=2: 0.02 0.06;   Max Sat NG=10: 0.09 0.35

















In [55]:



    


# Disk Images
wfe_drift = 5
wfe_roll_drift = 2
hdu_dict = do_gen_hdus(obs_dict, filt_keys, wfe_drift, wfe_roll_drift, 
                       PA1=-5, PA2=+5, opt_diff=False, use_cmask=True,
                       verbose=True)



    


















    






F182M_MASK335R_CIRCLYOT
F210M_MASK335R_CIRCLYOT
F250M_MASK335R_CIRCLYOT
F300M_MASK335R_CIRCLYOT
F335M_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [56]:



    


# Save FITS files
for k in filt_keys:
    hdul = hdu_dict[k]
    hdul.writeto(outdir+'{}_{}.fits'.format(name_sci,k), overwrite=True)



    











In [57]:



    


save_fig = True
if len(filt_keys)==2:
    plot_images(obs_dict, hdu_dict, filt_keys, wfe_drift, fov=fov, save_fig=save_fig, outdir=outdir)
else:
    plot_images_swlw(obs_dict, hdu_dict, filt_keys, wfe_drift, fov=fov, save_fig=save_fig, outdir=outdir)



    


















    
























In [58]:



    


# Determine contrast curves for various WFE drift values
wfe_list = [0,2,5,10]
nsig = 5
roll = 10

# (Roll1 - Ref) + (Roll2 - Ref)
curves_dict = do_contrast(obs_dict, wfe_list, filt_keys[-1:], nsig=nsig, roll_angle=roll, exclude_disk=False)
curves_F444W = curves_dict[filt_keys[-1]]

# Roll1 - Roll2
curves_dict = do_contrast(obs_dict, wfe_list, filt_keys[-1:], nsig=nsig, roll_angle=roll, exclude_disk=False,
                          no_ref=True)
curves_F444W2 = curves_dict[filt_keys[-1]]



    


















    






F444W_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [59]:



    


sat_rad = 0
obs = obs_dict[filt_keys[-1]]
#do_plot_contrasts(curves_F444W, curves_F444W2, obs, save_fig=True, xr=[0,5], xr2=[0,5])
do_plot_contrasts(curves_F444W, curves_F444W2, nsig, wfe_list, obs, age, 
                  save_fig=True, outdir=outdir, xr=[0,5], xr2=[0,5], yr2=[0.1,100])

HD 32297





In [60]:



    


# Fit spectrum to SED photometry
from pynrc.nrc_utils import source_spectrum

i=4
name_sci, dist_sci, age, spt_sci, Teff_sci, feh_sci, logg_sci, mag_sci, bp_sci, fov = args_sources[i]
vot = votdir + name_sci.replace(' ' ,'') + '.vot'

args = (name_sci, spt_sci, mag_sci, bp_sci, vot)
kwargs = {'Teff':Teff_sci, 'metallicity':feh_sci, 'log_g':logg_sci}
src = source_spectrum(*args, **kwargs)

src.fit_SED(use_err=False, robust=True, wlim=[1,4])

# Final source spectrum
sp_sci = src.sp_model



    


















    






[0.99412904]

















In [61]:



    


# Do the same for the reference source
name_ref, spt_ref, Teff_ref, feh_ref, logg_ref, mag_ref, bp_ref = ref_sources[i]
vot = votdir + name_ref.replace(' ' ,'') + '.vot'

args = (name_ref, spt_ref, mag_ref, bp_ref, vot)
kwargs = {'Teff':Teff_ref, 'metallicity':feh_ref, 'log_g':logg_ref}
ref = nrc_utils.source_spectrum(*args, **kwargs)

ref.fit_SED(use_err=True, robust=True)

# Final reference spectrum
sp_ref = ref.sp_model



    


















    






[0.93365695]

















In [62]:



    


# Plot spectra 

fig, axes = plt.subplots(1,2, figsize=(14,4.5))
src.plot_SED(ax=axes[0], xr=[0.5,30])
ref.plot_SED(ax=axes[1], xr=[0.5,30])

axes[0].set_title('Science Specta -- {} ({})'.format(src.name, spt_sci))
axes[1].set_title('Reference Specta -- {} ({})'.format(ref.name, spt_ref))

#for ax in axes:
#    ax.set_xscale('linear')
#    ax.xaxis.set_minor_locator(AutoMinorLocator())

fig.tight_layout()
fig.subplots_adjust(top=0.85, bottom=0.1 , left=0.05, right=0.97)

fig.savefig(outdir+'{}_SEDs.pdf'.format(name_sci.replace(' ','')))



    


















    
























In [63]:



    


# Plot the two spectra
fig, ax = plt.subplots(1,1, figsize=(8,5))

xr = [2.5,5.5]

for sp in [sp_sci, sp_ref]:
    w = sp.wave / 1e4
    ind = (w>=xr[0]) & (w<=xr[1])
    ind2 = (w>=3.9) & (w<=4.1)
    sp.convert('Jy')
    f = sp.flux / np.mean(sp.flux[ind2])
    ax.plot(w[ind], f[ind], lw=1, label=sp.name)
    ax.set_ylabel('Flux (Jy) normalized at 4 $\mu m$')
    sp.convert('flam')

ax.set_xlim(xr)
ax.set_ylim([0,ax.get_ylim()[1]])
ax.set_xlabel(r'Wavelength ($\mu m$)')
ax.set_title('{} Spectra'.format(sp_sci.name))

# Overplot Filter Bandpass
bp = pynrc.read_filter(*args_filter[-1])
ax2 = ax.twinx()
ax2.plot(bp.wave/1e4, bp.throughput, color='C2', label=bp.name+' Bandpass')
ax2.set_ylim([0,0.8])
ax2.set_xlim(xr)
ax2.set_ylabel('Bandpass Throughput')

ax.legend(loc='upper left')
ax2.legend(loc='upper right')

fig.tight_layout()
fig.savefig(outdir+'{}_2SEDs.pdf'.format(name_sci.replace(' ','')))



    


















    
























In [64]:



    


# Create a dictionary that holds the obs_coronagraphy class for each filter
wfe_drift = 0
obs_dict = obs_wfe(wfe_drift, args_filter, sp_sci, dist_sci, sp_ref=sp_ref, args_disk='auto', 
                   wind_mode='WINDOW', subsize=subsize, verbose=False)



    


















    






F182M_MASK335R_CIRCLYOT
F210M_MASK335R_CIRCLYOT
F250M_MASK335R_CIRCLYOT
F300M_MASK335R_CIRCLYOT
F335M_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [65]:



    


# Optimize readout parameters
tacq = 4200
#do_opt(tacq, patterns='MEDIUM8', ng_min=8, ng_max=8, tacq_frac=0.1, well_levels=[2], even_nints=True)



    











In [66]:



    


for key in filt_keys:
    obs = obs_dict[key]

    read_mode='MEDIUM8'
    ng, nint_sci, nint_ref = (10,15,15)
    if obs.bandpass.avgwave()/1e4 < 2.5:
        nint_sci *= 2
        nint_ref *= 2
    
    obs.update_detectors(read_mode=read_mode, ngroup=ng, nint=nint_sci)
    obs.nrc_ref.update_detectors(read_mode=read_mode, ngroup=ng, nint=nint_ref)

    #print(key)
    #print(obs.multiaccum_times)
    #_ = obs.sensitivity(nsig=5, units='vegamag', verbose=True)
    #print('')



    











In [67]:



    


# Saturation Levels
# Only want F444W
#obs = obs_dict[filt_keys[-1]]
#sat_rad = do_sat_levels(obs, plot=True)

# Max Saturation Values
dmax = []
for k in filt_keys:
    print(k)
    obs = obs_dict[k]
    
    imsci = obs.gen_slope_image(exclude_noise=True, quick_PSF=True)
    imref = obs.gen_slope_image(exclude_noise=True, quick_PSF=True, do_ref=True)
    
    ng1 = 2
    ng2 = obs.multiaccum.ngroup
    sci_sat1 = obs.saturation_levels(ngroup=ng1, image=imsci)
    sci_sat2 = obs.saturation_levels(ngroup=ng2, image=imsci)
    ref_sat1 = obs.saturation_levels(ngroup=ng1, image=imref, do_ref=True)
    ref_sat2 = obs.saturation_levels(ngroup=ng2, image=imref, do_ref=True)
    
    print('Max Sat NG={}: {:.2f} {:.2f};   Max Sat NG={}: {:.2f} {:.2f}'\
          .format(ng1,sci_sat1.max(),ref_sat1.max(),ng2,sci_sat2.max(),ref_sat2.max()))



    


















    






F182M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.03 0.07;   Max Sat NG=10: 0.15 0.41
F210M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.03 0.08;   Max Sat NG=10: 0.16 0.45
F250M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.02 0.04;   Max Sat NG=10: 0.09 0.24
F300M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.01 0.03;   Max Sat NG=10: 0.06 0.17
F335M_MASK335R_CIRCLYOT
Max Sat NG=2: 0.01 0.02;   Max Sat NG=10: 0.05 0.13
F444W_MASK335R_CIRCLYOT
Max Sat NG=2: 0.01 0.02;   Max Sat NG=10: 0.03 0.09

















In [68]:



    


# Disk Images
wfe_drift = 5
wfe_roll_drift = 2
hdu_dict = do_gen_hdus(obs_dict, filt_keys, wfe_drift, wfe_roll_drift, 
                       PA1=-5, PA2=+5, opt_diff=False, use_cmask=True,
                       verbose=True)



    


















    






F182M_MASK335R_CIRCLYOT
F210M_MASK335R_CIRCLYOT
F250M_MASK335R_CIRCLYOT
F300M_MASK335R_CIRCLYOT
F335M_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [69]:



    


# Save FITS files
for k in filt_keys:
    hdul = hdu_dict[k]
    hdul.writeto(outdir+'{}_{}.fits'.format(name_sci,k), overwrite=True)



    











In [70]:



    


save_fig = True
if len(filt_keys)==2:
    plot_images(obs_dict, hdu_dict, filt_keys, wfe_drift, fov=fov, save_fig=save_fig, outdir=outdir)
else:
    plot_images_swlw(obs_dict, hdu_dict, filt_keys, wfe_drift, fov=fov, save_fig=save_fig, outdir=outdir)



    


















    
























In [71]:



    


# Determine contrast curves for various WFE drift values
wfe_list = [0,2,5,10]
nsig = 5
roll = 10

# (Roll1 - Ref) + (Roll2 - Ref)
curves_dict = do_contrast(obs_dict, wfe_list, filt_keys[-1:], nsig=nsig, roll_angle=roll, exclude_disk=False)
curves_F444W = curves_dict[filt_keys[-1]]

# Roll1 - Roll2
curves_dict = do_contrast(obs_dict, wfe_list, filt_keys[-1:], nsig=nsig, roll_angle=roll, exclude_disk=False,
                          no_ref=True)
curves_F444W2 = curves_dict[filt_keys[-1]]



    


















    






F444W_MASK335R_CIRCLYOT
F444W_MASK335R_CIRCLYOT

















In [72]:



    


sat_rad = 0
obs = obs_dict[filt_keys[-1]]
do_plot_contrasts(curves_F444W, curves_F444W2, nsig, wfe_list, obs, age, 
                  save_fig=True, outdir=outdir, xr=[0,5], xr2=[0,5], yr2=[0.3,100])



    


















    
























In [73]:



    


obs.multiaccum_times



    


















    
Out[73]:








{'t_frame': 1.06904,
 't_group': 10.6904,
 't_int': 104.76592,
 't_exp': 1571.4887999999999,
 't_acq': 1587.8232,
 't_int_tot1': 105.85488,
 't_int_tot2': 105.85488}

















In [ ]:

Content source: JarronL/pynrc

Similar notebooks:

notebook.community | gallery | about