<jupyter_start><jupyter_text>First import all the things and setup the plot things<jupyter_code>from psisim import telescope,instrument,observation,spectrum,universe,plots,signal
import time
import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt
from importlib import reload
import speclite.filters
from scipy.interpolate import interp1d, RegularGridInterpolator as rgi
from scipy.ndimage import gaussian_filter
from scipy.ndimage.filters import median_filter
from scipy.signal import medfilt, correlate
from numpy.random import poisson, randn
import copy 

SMALL_SIZE = 32
MEDIUM_SIZE = 40
BIGGER_SIZE = 48

plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title<jupyter_output><empty_output><jupyter_text>Initial Setup Setup the path to the peripheral files and load in the filter objects from speclite<jupyter_code>path = '/scr3/dmawet/ETC/'
filters = spectrum.load_filters()<jupyter_output><empty_output><jupyter_text>Telescope Setup A *telescope* object contains the telecope parameters, the site characteristics and the current conditions.<jupyter_code>keck = telescope.Keck(path=path)
keck.airmass=1.0
keck.water_vapor=1.6
keck.seeing = keck.median_seeing<jupyter_output><empty_output><jupyter_text>Instrument Setup An _instrument_ object contains both the static instrument parameters, as well as the current observing setup<jupyter_code>#Initialize a hispec object and set the current filter
hispec = instrument.hispec(telescope=keck)
hispec.set_current_filter("TwoMASS-J")

#Get the set of wavelengths based on the current instrument setup
wavelengths = hispec.get_wavelength_range()

# Set the observing mode: Exposure time (per exposure), Number of Exposures,filter name, wavelength array
hispec.set_observing_mode(3600,1,'TwoMASS-J', wavelengths,mode="off-axis")<jupyter_output><empty_output><jupyter_text>Generate some Spectra We can now generate a host spectrum<jupyter_code>#First set the host properties for a Phoenix model. 
host_properties = {"StarLogg":4.00*u.dex(u.cm/ u.s**2),"StarTeff":10000*u.K,"StarZ":'-0.0',"StarAlpha":"0.0",
                   "StarRadialVelocity":100*u.km/u.s,
                    "StarVsini":10*u.km/u.s,
                    "StarLimbDarkening":0.8}

#Now setup the user parameters that a Phoenix model needs: (path, object_filter, magnitude_in_object_filter,
# filters_object,current_filter). 
host_user_params = (path,'TwoMASS-J',5.0,filters,hispec.current_filter)

#Generate the spectrum! (Here we apply a doppler shift and rotationally broaden)
host_spectrum = spectrum.get_stellar_spectrum(host_properties,wavelengths,hispec.current_R,model="Phoenix",
                                              user_params=host_user_params,doppler_shift=True,broaden=True,
                                              delta_wv=hispec.current_dwvs)<jupyter_output><empty_output><jupyter_text>Plot the Spectrum!<jupyter_code>plt.figure(figsize=(20,10))

plt.semilogy(wavelengths,host_spectrum)

plt.xlabel("Wavelength [{}]".format(wavelengths.unit))
plt.ylabel("Spectrum [{}]".format(host_spectrum.unit))<jupyter_output><empty_output><jupyter_text>Now we can generate a target spectrum as well<jupyter_code>obj_properties = {"StarLogg":3.25*u.dex(u.cm/ u.s**2),"StarTeff":700*u.K,"StarRadialVelocity":20*u.km/u.s,"StarVsini":10*u.km/u.s,
                  "StarLimbDarkening":0.9}

obj_user_params = (path,'TwoMASS-K',20,filters,hispec.current_filter)

obj_spectrum = spectrum.get_stellar_spectrum(obj_properties,wavelengths,hispec.current_R,model="Sonora",
                                             user_params=obj_user_params,doppler_shift=True,broaden=True,
                                             delta_wv=hispec.current_dwvs)
plt.figure(figsize=(30,10))

plt.semilogy(wavelengths,obj_spectrum)

plt.xlabel("Wavelength [{}]".format(wavelengths.unit))
plt.ylabel("Spectrum [{}]".format(obj_spectrum.unit))
plt.ylim(1e-10,1e-4)<jupyter_output><empty_output><jupyter_text>Because of the original PSISIM design considerations, we want the object spectrum in contrast units<jupyter_code>obj_spectrum /= host_spectrum
plt.semilogy(wavelengths,obj_spectrum)<jupyter_output><empty_output><jupyter_text>Simulate an observation First set some extra parameters<jupyter_code># The angular separation of the companion, in milliarcsecond
host_properties['AngSep'] = 400 *u.mas

#Get the host star magnitude in the AO filter
host_properties["StarAOmag"] = spectrum.get_model_ABmags(host_properties,[hispec.ao_filter], model='Phoenix',
                                                         verbose=False,user_params = host_user_params)
hispec.ao_mag = host_properties["StarAOmag"]

#Hispec doesn't care about the spectral type, but we need to include teh paramter
host_properties['StarSpT'] = None<jupyter_output><empty_output><jupyter_text>Actually simulate the observation<jupyter_code>obj_spec,total_noise,stellar_spec,thermal_spec,noise_components = observation.simulate_observation(keck,hispec,host_properties,
                                                                                      obj_spectrum,wavelengths,1e5,
                                                                                      inject_noise=False,verbose=True,
                                                                                       post_processing_gain = 10,
                                                                                      return_noise_components=True,
                                                                                      stellar_spec=host_spectrum,
                                                                                      apply_lsf=True,
                                                                                      integrate_delta_wv=False)<jupyter_output><empty_output><jupyter_text>Plot it!<jupyter_code>plt.figure(figsize=(30,10))

plt.semilogy(wavelengths,obj_spec)

plt.xlabel("Wavelength [{}]".format(wavelengths.unit))
plt.ylabel("Spectrum [{}]".format(obj_spec.unit))
plt.ylim(1e-4,2e2)<jupyter_output><empty_output><jupyter_text>Now let's cycle through the available filters and collect all the results<jupyter_code>reload(observation)
all_wavelengths = []
full_host_spectrum = []
full_obj_spectrum = []
full_obj_spectrum_nosky = []
full_total_noise = []
full_thermal_spec = []
full_noise_components = []
full_speckle_noise = []


for hispec_filter in hispec.filters:
    
    #Setup the instrument
    hispec.set_current_filter(hispec_filter)
    wavelengths = hispec.get_wavelength_range()
    hispec.set_observing_mode(3600,1,hispec_filter, wavelengths) 
    
    
    host_user_params = (path,'TwoMASS-J',5.0,filters,hispec.current_filter)
    host_spectrum = spectrum.get_stellar_spectrum(host_properties,wavelengths,hispec.current_R,
                                                  model="Phoenix",user_params=host_user_params,
                                                  doppler_shift=True,broaden=True,delta_wv=hispec.current_dwvs)
    
    obj_user_params = (path,'TwoMASS-K',20,filters,hispec.current_filter)
    obj_spectrum = spectrum.get_stellar_spectrum(obj_properties,wavelengths,hispec.current_R,model="Sonora",
                                                 user_params=obj_user_params,doppler_shift=True,broaden=True,
                                                 delta_wv=hispec.current_dwvs)
    
    obj_spectrum /= host_spectrum
    
    obj_spec,total_noise,stellar_spec,thermal_spec,noise_components= observation.simulate_observation(keck,hispec,
                                                                                      host_properties,
                                                                                      obj_spectrum,wavelengths,1e5,
                                                                                      inject_noise=False,verbose=True,
                                                                                      post_processing_gain = np.inf,
                                                                                      return_noise_components=True,
                                                                                      stellar_spec=host_spectrum,
                                                                                      apply_lsf=True,
                                                                                      integrate_delta_wv=False,
                                                                                      plot=False,
                                                                                      sky_on=True)
    obj_spec_no_sky,_,_,_ = observation.simulate_observation(keck,hispec,
                                                                                      host_properties,
                                                                                      obj_spectrum,wavelengths,1e5,
                                                                                      inject_noise=False,verbose=True,
                                                                                      post_processing_gain = np.inf,
                                                                                      return_noise_components=False,
                                                                                      stellar_spec=host_spectrum,
                                                                                      apply_lsf=True,
                                                                                      integrate_delta_wv=False,
                                                                                      plot=False,
                                                                                      sky_on=False)

    full_speckle_noise.append(hispec.get_speckle_noise(0.4*u.arcsecond,host_properties["StarAOmag"],hispec.current_filter,
                                             wavelengths,host_properties['StarSpT'],keck)[0])
    all_wavelengths.append(wavelengths)
    full_host_spectrum.append(stellar_spec)
    full_obj_spectrum.append(obj_spec)
    full_obj_spectrum_nosky.append(obj_spec_no_sky)
    full_total_noise.append(total_noise)
    full_thermal_spec.append(thermal_spec)
    full_noise_components.append(noise_components)

all_wavelengths = np.hstack(all_wavelengths).value*wavelengths.unit
full_host_spectrum = np.hstack(full_host_spectrum).value*stellar_spec.unit
full_obj_spectrum = np.hstack(full_obj_spectrum).value*obj_spec.unit
full_obj_spectrum_nosky = np.hstack(full_obj_spectrum_nosky).value*obj_spec_no_sky.unit
full_total_noise = np.hstack(full_total_noise).value*total_noise.unit
full_noise_components = np.hstack(full_noise_components)*obj_spec.unit
full_thermal_spec = np.hstack(full_thermal_spec).value*thermal_spec.unit
full_speckle_noise = np.hstack(full_speckle_noise)
plt.figure(figsize=(30,10))

plt.plot(all_wavelengths,full_obj_spectrum,label="Spectrum")
plt.plot(all_wavelengths,full_total_noise,label="Total Stastistical Noise Level")

plt.legend()

plt.xlabel("Wavelength [{}]".format(all_wavelengths.unit))
plt.ylabel("Spectrum [{}]".format(full_obj_spectrum.unit))
plt.title(r"HISPEC Observation of Sonora grid $T_{{eff}}$ = {}K,logg = {}, Exp. Time = {}, $N_{{exp}}$= {}".format(
obj_properties["StarTeff"].value,obj_properties["StarLogg"].value,hispec.exposure_time,hispec.n_exposures))
plt.ylim(1e-4,1e2)
plt.grid()<jupyter_output><empty_output><jupyter_text>We can calculate the signal to noise now<jupyter_code>plt.figure(figsize=(30,10))

plt.plot(all_wavelengths,full_obj_spectrum/full_total_noise)

plt.xlabel("Wavelength [{}]".format(all_wavelengths.unit))
plt.ylabel("SNR")

plt.title(r"HISPEC Observation of Sonora grid $T_{{eff}}$ = {}K,logg = {}, Exp. Time = {}, $N_{{exp}}$= {}".format(
obj_properties["StarTeff"].value,obj_properties["StarLogg"].value,hispec.exposure_time,hispec.n_exposures))

# plt.ylim(-0.01,0.4)
plt.xlim(0.9,2.4)
plt.grid()
plt.ylim(-0.1,2.)<jupyter_output><empty_output><jupyter_text>But let's not forget about systematics. We will estimate our ability to remove the host star spectrum and the telluric lines at the 1% level.<jupyter_code>cal = 0.01 #Telluric calibration accuracy of 1%

#The photon flux at the object will be the stellar flux multipled by the contrast there: 
# full_host_spectrum
host_flux_at_obj = full_host_spectrum*hispec.get_speckle_noise(0.4*u.arcsecond,host_properties["StarAOmag"],
                                                               "",all_wavelengths,None,keck)[0]

systematics = (cal*(host_flux_at_obj+full_thermal_spec))**2 #Variance of systematics

noise_plus_systematics = np.sqrt(full_total_noise**2+systematics)
plt.figure(figsize=(30,10))

plt.plot(all_wavelengths,full_obj_spectrum/noise_plus_systematics)

plt.xlabel("Wavelength [{}]".format(all_wavelengths.unit))
plt.ylabel("SNR with systematics")

plt.title(r"HISPEC Observation of Sonora grid $T_{{eff}}$ = {}K,logg = {}, Exp. Time = {}, $N_{{exp}}$= {}".format(
obj_properties["StarTeff"].value,obj_properties["StarLogg"].value,hispec.exposure_time,hispec.n_exposures))

# plt.ylim(-0.01,0.4)
plt.xlim(0.9,2.4)
plt.grid()
plt.ylim(-0.1,2.)<jupyter_output><empty_output><jupyter_text>Calculate the CCF-SNR  We'll do this with in the K-band (the last one from when we cycled through the filters above) Since we're doing this in k-band we'll want to recalculate the noise with the systematics just for this band<jupyter_code>cal = 0.01 #Telluric calibration accuracy of 1%

#The photon flux at the object will be the stellar flux multipled by the contrast there: 
# full_host_spectrum
host_flux_at_obj = stellar_spec*hispec.get_speckle_noise(0.4*u.arcsecond,host_properties["StarAOmag"],
                                                               "",wavelengths,None,keck)[0]

systematics = (cal*(host_flux_at_obj+thermal_spec))**2 #Variance of systematics

noise_plus_systematics = np.sqrt(total_noise**2+systematics)

plt.plot(np.sqrt(systematics),label="Systematic Noise")
plt.plot(total_noise,label="Random Noise")
plt.legend(fontsize=10)<jupyter_output><empty_output><jupyter_text>Calculate the sky transmission lines convolved with the line spread function<jupyter_code>#Get the wavelength spacing
dwvs = np.abs(wavelengths - np.roll(wavelengths, 1))
dwvs[0] = dwvs[1]
dwv_mean = np.mean(dwvs)
lsf_fwhm = (hispec.lsf_width/dwv_mean).decompose() #Get the lsf_fwhm in units of current wavelength spacing
lsf_sigma = lsf_fwhm/(2*np.sqrt(2*np.log(2))) #Convert to sigma

#Calculate the 
sky_transmission_lsf = gaussian_filter(keck.get_atmospheric_transmission(wavelengths),
                                                                           lsf_sigma.value)<jupyter_output><empty_output><jupyter_text>Now get the ccf snr<jupyter_code>ccf_snr = signal.compute_ccf_snr_matchedfilter(obj_spec, obj_spec_no_sky, noise_plus_systematics, 
                                            sky_transmission_lsf,
                                               systematics_residuals = cal, 
                                               norm_cutoff=0.8)
print(ccf_snr)<jupyter_output><empty_output><jupyter_text>We can also calculate the ccf SNF for all the bands  Some other things we can look at to compare with Dimitri's original notebook "HISPEC MODHIS ETC HDC v3"  The thermal spectrum of the sky, telecsope and instrument.<jupyter_code>plt.figure(figsize=(30,10))

# plt.semilogy(all_wavelengths,full_host_spectrum)
plt.semilogy(all_wavelengths,full_thermal_spec)

plt.xlabel("Wavelength [{}]".format(all_wavelengths.unit))
plt.ylabel('Thermal photon flux [{}]'.format(str((full_thermal_spec).unit)))
# plt.ylim(1e,1e7)
# plt.ylim(1e4,1e8)
plt.xlim(0.9,2.4)<jupyter_output><empty_output><jupyter_text>What about the relative contributions of the different noise components<jupyter_code>plt.figure(figsize=(30,10))

plt.plot(all_wavelengths,full_noise_components[0],label="Speckle Noise")
plt.plot(all_wavelengths,full_noise_components[1],label="Read Noise")
# plt.plot(all_wavelengths,full_noise_components[2],label="Dark Current Noise")
plt.plot(all_wavelengths,full_noise_components[2],label="Photon Noise")

plt.plot(all_wavelengths,full_total_noise,linestyle="-.",label="Total Noise",linewidth=0.2)

plt.xlabel("Wavelength [{}]".format(all_wavelengths.unit))
plt.ylabel("Spectrum [{}]".format(full_obj_spectrum.unit))
# plt.ylim(1e0,5e1)
plt.grid()
plt.legend()<jupyter_output><empty_output><jupyter_text>We can also check the throughput of the system:<jupyter_code>th_sky = []
th_total = []
SR = []

fig = plt.figure(figsize=(30,10))

for hispec_filter in hispec.filters:    
    #Setup the instrument
    hispec.set_current_filter(hispec_filter)
    wavelengths = hispec.get_wavelength_range()
    hispec.set_observing_mode(3600,1,hispec_filter, wavelengths) 
    
    
    th_sky.append(keck.get_atmospheric_transmission(wavelengths))
    th_total.append(keck.get_telescope_throughput(wavelengths,hispec)*hispec.get_inst_throughput(wavelengths)*hispec.get_filter_transmission(wavelengths,hispec.current_filter))
    SR.append(hispec.compute_SR(wavelengths))            
    
    plt.plot(wavelengths,0.95*keck.get_atmospheric_transmission(wavelengths)*keck.get_telescope_throughput(wavelengths,hispec)*hispec.get_inst_throughput(wavelengths)*hispec.get_filter_transmission(wavelengths,hispec.current_filter),label=hispec_filter)
    
#     print(np.mean(hispec.current_dwvs))

th_sky = np.hstack(th_sky)
th_total = np.hstack(th_total)
SR = np.hstack(SR)

plt.xlim(0.9,2.4)
plt.ylim(0,0.15)
plt.grid()
plt.xlabel("Wavelength [micron]")
plt.ylabel("Throughput")

plt.legend()<jupyter_output><empty_output><jupyter_text>And the Strehl Ratio<jupyter_code>plt.figure(figsize=(30,10))

plt.plot(all_wavelengths,SR,'o')

# plt.ylim(1e-2,0.2)
plt.xlabel("Wavelength [{}]".format(all_wavelengths.unit))
plt.ylabel("Strehl Ratio")
# plt.ylim(1e-4,1e2)
plt.grid()
plt.ylim(0.3,0.9)
plt.xlim(0.9,2.4)<jupyter_output><empty_output><jupyter_text>Check the contrasts and mark the 0.4 arcsecond location for each wavelength<jupyter_code>separations = np.arange(0,1.,0.001)*u.arcsecond

plt.figure(figsize=(30,10))

for hispec_filter in hispec.filters:
    
    hispec.set_current_filter(hispec_filter)
    wavelengths = hispec.get_wavelength_range()
    hispec.set_observing_mode(3600,1,hispec_filter, wavelengths) 
    
#     print(wavelengths[int(wavelengths.shape[0]/2):int(wavelengths.shape[0]/2)+2])
    
    contrasts = hispec.get_speckle_noise(separations,host_properties["StarAOmag"],hispec.current_filter,
                                       wavelengths[int(wavelengths.shape[0]/2):int(wavelengths.shape[0]/2)+2],
                                         host_properties['StarSpT'],keck)[:,0]
    
    tsts = hispec.get_speckle_noise(np.array([0.4,0.5,0.6])*u.arcsecond,host_properties["StarAOmag"],
                                       hispec.current_filter,
                                       wavelengths[int(wavelengths.shape[0]/2):int(wavelengths.shape[0]/2)+2],
                                       host_properties['StarSpT'],keck)[0]
    
    plt.semilogy(separations/206265*keck.diameter/np.mean(wavelengths).to(u.m),contrasts,
                 label="{:.2f}".format(np.mean(wavelengths)))
    plt.axvline(0.4/206265*keck.diameter/np.mean(wavelengths).to(u.m))
    plt.plot([0.4/206265*keck.diameter/np.mean(wavelengths).to(u.m)],[tsts],'o')
    
# plt.xlim([0.01,100])
# plt.ylim([1e-7,1])
plt.xlabel('Angular separation ($\lambda$/D)')
plt.ylabel('Raw contrast')
plt.legend()

plt.grid(True)
plt.xlim(7.5,20)
plt.ylim(1e-5,1e-3)<jupyter_output><empty_output><jupyter_text>Cross Correlating Noisy Spectra with Model<jupyter_code>y = 0.94/1e5/3; j = 1.248 / 1e5 /3; h = 1.633 / 1e5 / 3; k = 2.200 / 1e5 / 3
exp_time = 3600 * u.s
delta_lb = np.array([y, j, h, k])*u.micron
samp = np.array([3,3,3,3])

## code to split the arrays for nicer plots. not necessary but splits plots into 4 filters
wavesplit = np.array([all_wavelengths[97:44205],all_wavelengths[44301:106705],all_wavelengths[106803:168167],all_wavelengths[169262::]]) 
flux_split = np.array([host_flux_at_obj[97:44205],host_flux_at_obj[44301:106705],host_flux_at_obj[106803:168167],host_flux_at_obj[169262::]]) 
thermal_flux_split = np.array([full_thermal_spec[97:44205],full_thermal_spec[44301:106705],full_thermal_spec[106803:168167],full_thermal_spec[169262::]]) 
obj_spec_split = np.array([full_obj_spectrum[97:44205].value,full_obj_spectrum[44301:106705].value,full_obj_spectrum[106803:168167].value,full_obj_spectrum[169262::].value]) 
ron_split = np.array([full_noise_components[1][97:44205],full_noise_components[1][44301:106705],full_noise_components[1][106803:168167],full_noise_components[1][169262::]]) 
sky_trans_split = np.array([sky_transmission_lsf[97:44205],sky_transmission_lsf[44301:106705],sky_transmission_lsf[106803:168167],sky_transmission_lsf[169262::]]) 
model_spec_split= np.array([full_obj_spectrum_nosky[97:44205],full_obj_spectrum_nosky[44301:106705],full_obj_spectrum_nosky[106803:168167],full_obj_spectrum_nosky[169262::]]) 
noisy_spec_filt, model_filt, ccf, ccf_model, velocity= signal.reduce_spec(wavesplit, delta_lb, samp, obj_spec_split, model_spec_split,flux_split, thermal_flux_split, ron_split, exp_time, sky_trans_split,cal, 301, 0.5)
for i in range(obj_spec_split.shape[0]):
    plt.figure(figsize=(30,10))
    plt.plot(wavesplit[i],model_filt[i], label='Model')
    plt.plot(wavesplit[i],noisy_spec_filt[i], label='Data',alpha=0.5)
    #plt.title(fac+'. Host: '+str(host_mag)+' mag (Vega, ' + host_filt +'), '+str(host_temp)+'K star')
    plt.xlabel('Wavelength (micron)')
    plt.ylabel('Noisy spectrum reduced, filtered')
    plt.legend()
    
    plt.figure(figsize=(30,10))
    plt.plot(velocity[i],ccf_model[i], label='Model ACF')
    plt.plot(velocity[i],ccf[i], label='CCF data with model')
    plt.xlim([-100,100])
    #plt.title(fac+'. Host: '+str(host_mag)+' mag (Vega, ' + host_filt +'), '+str(host_temp)+'K star')
    plt.xlabel('Velocity (km/s)')
    plt.ylabel('Correlation')
    plt.legend()<jupyter_output><empty_output><jupyter_text>Plot Spectra with Features<jupyter_code>for i in range(noisy_spec_filt.shape[0]):
    signal.plot_features(wavesplit[i],model_filt[i],noisy_spec_filt[i],features=['h2o','ch4','co','na','ca','vo','feh','k','h'])<jupyter_output><empty_output><jupyter_text>Compute PRV Accuracy<jupyter_code>#need to double check math
snr = full_obj_spectrum/full_total_noise
sigma_rv = signal.compute_prv_sigma(wavesplit, delta_lb, snr, model_spec_split, sky_trans_split, cal, 101)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11214223764093668, 'fraction_numerical': 0.027623590633130964, 'mean_word_length': 3.1558839430528023, 'pattern_counts': {'":': 12, '<': 101, '<?xml version=': 0, '>': 101, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 147, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '443391', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7969, 'n_tokens_neox': 7271, 'n_words': 1059, 'path': 'Tutorials/PSISIM HISPEC Tutorial.ipynb', 'repo_name': 'abgibbs/psisim', 'size': '23060'}