<jupyter_start><jupyter_text>An introduction to HRS data reduction<NAME>  In this notebook, we describe how to use `hrs` package to reduce HRS data.To install `hrs` package, type `pip install hrs` in your terminal.<jupyter_code>import os
import numpy as np
import matplotlib.pyplot as plt

import ccdproc
import hrs

# from astropy.table import Table, Column, Row
os.chdir('/town/HRS/20161110test')
# %matplotlib 
ls<jupyter_output>2016111000100.fits  20161110002.fits   201611100054.fits  201611100079.fits
2016111000101.fits  201611100030.fits  201611100055.fits  20161110007.fits
2016111000103.fits  201611100031.fits  201611100056.fits  201611100081.fits
2016111000104.fits  201611100032.fits  201611100057.fits  201611100082.fits
2016111000105.fits  201611100034.fits  201611100058.fits  201611100083.fits
2016111000106.fits  201611100035.fits  201611100059.fits  201611100084.fits
2016111000107.fits  201611100036.fits  20161110005.fits   201611100085.fits
201611100010.fits   201611100037.fits  201611100060.fits  201611100086.fits
201611100012.fits   201611100038.fits  201611100061.fits  201611100087.fits
201611100013.fits   201611100039.fits  201611100062.fits  201611100088.fits
201611100014.fits   20161110003.fits   201611100063.fits  201611100089.fits
201611100015.fits   201611100040.fits  201611100064.fits  20161110008.fits
201611100016.fits   201611100041.fits  201611100065.fits  201611100090.fits
2[...]<jupyter_text>1. Read log table<jupyter_code>logt = hrs.LogTable.read_log('./20161110.txt')
logt.pprint(1000, 1000)<jupyter_output>file      obj    exp_start exp_time      RA         Dec     epoch note      filename     
----------- -------- --------- -------- ----------- ----------- ----- ---- ------------------
20161110001     flat  19:39:54      9.0                                      20161110001.fits
20161110002     flat  19:39:54      9.0                                      20161110002.fits
20161110003     flat  19:39:54      9.0                                      20161110003.fits
20161110004     flat  19:39:54      9.0                                      20161110004.fits
20161110005     flat  19:39:54      9.0                                      20161110005.fits
20161110006     flat  19:39:54      9.0                                      20161110006.fits
20161110007     flat  19:39:54      9.0                                      20161110007.fits
20161110008     flat  19:39:54      9.0                                      20161110008.fits
20161110009     flat  19:39:54      9.0                     [...]<jupyter_text>2. Process BIAS<jupyter_code>ind_bias = logt['obj']=='bias'
bias = ccdproc.combine(','.join(logt['filename'][ind_bias]), unit='adu', method='median')<jupyter_output>INFO: splitting each image into 1 chunks to limit memory usage to 16000000000.0 bytes. [ccdproc.combiner]<jupyter_text>3. Process FLAT  bias substraction & trim<jupyter_code>ind_flat1_1 = (logt['obj']=='flat')*(logt['exp_time']==8)
ind_flat1_2 = (logt['obj']=='flat')*(logt['exp_time']==9)
ind_flat2_1 = (logt['obj']=='flat')*(logt['exp_time']==32)
ind_flat2_2 = (logt['obj']=='flat')*(logt['exp_time']==36)
ind_flat3 = (logt['obj']=='flat')*(logt['exp_time']==144)

flat1_1 = ccdproc.combine(','.join(logt['filename'][ind_flat1_1]), unit='adu', method='median')
flat1_2 = ccdproc.combine(','.join(logt['filename'][ind_flat1_2]), unit='adu', method='median')
flat2_1 = ccdproc.combine(','.join(logt['filename'][ind_flat2_1]), unit='adu', method='median')
flat2_2 = ccdproc.combine(','.join(logt['filename'][ind_flat2_2]), unit='adu', method='median')
flat3 = ccdproc.combine(','.join(logt['filename'][ind_flat3]), unit='adu', method='median')
flat1 = ccdproc.combine([flat1_1.multiply(9.).divide(8.), flat1_2], method='average')
flat2 = ccdproc.combine([flat2_1.multiply(36.).divide(32.), flat2_2], method='average')

flat1_bias = ccdproc.subtract_bias(flat1, bias)
flat2_bias = ccdproc.subtract_bias(flat2, bias)
flat3_bias = ccdproc.subtract_bias(flat3, bias)

flat1_bias_trim = ccdproc.trim_image(flat1_bias[:, :4096])
flat2_bias_trim = ccdproc.trim_image(flat2_bias[:, :4096])
flat3_bias_trim = ccdproc.trim_image(flat3_bias[:, :4096])<jupyter_output>INFO: splitting each image into 1 chunks to limit memory usage to 16000000000.0 bytes. [ccdproc.combiner]
INFO: splitting each image into 1 chunks to limit memory usage to 16000000000.0 bytes. [ccdproc.combiner]
INFO: splitting each image into 1 chunks to limit memory usage to 16000000000.0 bytes. [ccdproc.combiner]
INFO: splitting each image into 1 chunks to limit memory usage to 16000000000.0 bytes. [ccdproc.combiner]
INFO: splitting each image into 1 chunks to limit memory usage to 16000000000.0 bytes. [ccdproc.combiner]
INFO: splitting each image into 1 chunks to limit memory usage to 16000000000.0 bytes. [ccdproc.combiner]
INFO: splitting each image into 1 chunks to limit memory usage to 16000000000.0 bytes. [ccdproc.combiner]<jupyter_text>combine FLAT<jupyter_code>flat_list = [flat1_bias_trim, flat2_bias_trim, flat3_bias_trim]

# find & combine & group apertures
ap_comb = hrs.combine_apertures(flat_list, n_jobs=10)
cheb_coefs, ap_uorder_interp = hrs.group_apertures(ap_comb, start_col=2100, order_dist=10)

# combine flat
flat_comb, flat_origin = hrs.combine_flat(flat_list, ap_uorder_interp, sat_count=45000, p=95)
flat_comb = ccdproc.CCDData(flat_comb, unit='adu')

# scattered light substraction
flat_comb_sl = hrs.substract_scattered_light(flat_comb, ap_uorder_interp, ap_width=10, shrink=.85)
flat1d = hrs.extract_1dspec(flat_comb_sl, ap_uorder_interp, ap_width=7)[0]<jupyter_output>@Cham: 143 possible apertures found (exciting)!<jupyter_text>4. Process ThAr  bias substraction<jupyter_code>ind_thar1 = (logt['obj']=='thar')*(logt['exp_time']==30)
ind_thar2 = (logt['obj']=='thar')*(logt['exp_time']==60)
ind_thar3 = (logt['obj']=='thar')*(logt['exp_time']==120)

thar1 = ccdproc.combine(','.join(logt['filename'][ind_thar1]), unit='adu', method='average')
thar2 = ccdproc.combine(','.join(logt['filename'][ind_thar2]), unit='adu', method='average')
thar3 = ccdproc.CCDData.read(','.join(logt['filename'][ind_thar3]), unit='adu', method='average')

thar1_bias = ccdproc.subtract_bias(thar1, bias)
thar2_bias = ccdproc.subtract_bias(thar2, bias)
thar3_bias = ccdproc.subtract_bias(thar3, bias)

thar1_bias_trim = ccdproc.trim_image(thar1_bias[:, :4096])
thar2_bias_trim = ccdproc.trim_image(thar2_bias[:, :4096])
thar3_bias_trim = ccdproc.trim_image(thar3_bias[:, :4096])<jupyter_output>INFO: splitting each image into 1 chunks to limit memory usage to 16000000000.0 bytes. [ccdproc.combiner]
INFO: splitting each image into 1 chunks to limit memory usage to 16000000000.0 bytes. [ccdproc.combiner]<jupyter_text>combine & trim<jupyter_code>""" combine thar """
thar_list = [thar1_bias_trim, thar2_bias_trim, thar3_bias_trim]
thar_comb, thar_origin = hrs.combine_flat(thar_list, ap_uorder_interp, sat_count=45000, p=100)
thar_comb = ccdproc.CCDData(thar_comb, unit='adu')

""" extract 1d thar """
#thar_comb_sl = substract_scattered_light(thar_comb, ap_uorder_interp, ap_width=10, shrink=.5)
thar1d = hrs.extract_1dspec(thar_comb, ap_uorder_interp, ap_width=7)[0]
thar1d_scaled = thar1d/flat1d<jupyter_output>@Cham: marking pixels of order [0]
@Cham: marking pixels of order [1]
@Cham: marking pixels of order [2]
@Cham: marking pixels of order [3]
@Cham: marking pixels of order [4]
@Cham: marking pixels of order [5]
@Cham: marking pixels of order [6]
@Cham: marking pixels of order [7]
@Cham: marking pixels of order [8]
@Cham: marking pixels of order [9]
@Cham: marking pixels of order [10]
@Cham: marking pixels of order [11]
@Cham: marking pixels of order [12]
@Cham: marking pixels of order [13]
@Cham: marking pixels of order [14]
@Cham: marking pixels of order [15]
@Cham: marking pixels of order [16]
@Cham: marking pixels of order [17]
@Cham: marking pixels of order [18]
@Cham: marking pixels of order [19]
@Cham: marking pixels of order [20]
@Cham: marking pixels of order [21]
@Cham: marking pixels of order [22]
@Cham: marking pixels of order [23]
@Cham: marking pixels of order [24]
@Cham: marking pixels of order [25]
@Cham: marking pixels of order [26]
@Cham: marking pixels of order [27]
@C[...]<jupyter_text>5. Wavelength calibration  initial estimation<jupyter_code># laod template thar
thar_temp_path = '/home/cham/PycharmProjects/hrs/hrs/calibration/thar_template/thar_temp_w20160120022t.fits'
wave_temp, thar_temp, order_temp = hrs.load_thar_temp(thar_temp_path)

# fix thar
thar1d_fixed = hrs.fix_thar_sat_neg(thar1d, arm=30, sat_count=500000)
thar_temp_fixed = hrs.fix_thar_sat_neg(thar_temp, arm=30, sat_count=500000)

""" initial estimation of wavelength """
shift, corr2d = hrs.thar_corr2d(thar1d_fixed, thar_temp, ytrim=(30, 70), y_shiftmax = 3, x_shiftmax=20)
print ("@Cham: the shift is ", shift)

wave_init = hrs.interpolate_wavelength_shift(wave_temp, shift, thar_temp, thar1d_fixed)
order_init = hrs.interpolate_order_shift(order_temp, shift, thar1d_fixed)<jupyter_output>@Cham: computing 2D cross-correlation...
('@Cham: the shift is ', (array([-2]), array([0])))<jupyter_text>refine wavelength<jupyter_code>thar_list = np.loadtxt('/home/cham/PycharmProjects/hrs/hrs/calibration/iraf/thar.dat')
poly_order = (3, 5)

# refine thar positions
lc_coord, lc_order, lc_thar, popt, pcov = hrs.calibration.refine_thar_positions(
        wave_init, order_init, thar1d_fixed, thar_list, n_jobs=20, verbose=10)

# fit grating function    
x_mini_lsq, ind_good_thar, scaler_coord, scaler_order, scaler_ml = hrs.calibration.fit_grating_equation(
    lc_coord, lc_order, lc_thar, popt, pcov, poly_order=(3, 5), max_dev_threshold=1., iter=False)
    
# recover the wavelength (grid)
grid_coord = np.repeat(np.arange(wave_init.shape[1]).reshape(1,-1), 104, axis=0)
grid_order = order_init
wave_fitted =  hrs.calibration.grating_equation_predict(
    grid_coord, grid_order, x_mini_lsq, poly_order, scaler_coord, scaler_order, scaler_ml)

# recover the wavelength (thar)
lc_fitted = hrs.grating_equation_predict(
    lc_coord, lc_order, x_mini_lsq, poly_order, scaler_coord, scaler_order, scaler_ml)
lc_fitted_diff = lc_fitted - lc_thar    
""" residual check zoom in """
%matplotlib inline
fig = plt.figure(figsize=(16,8));
ax = fig.add_subplot(221)
ax.plot(lc_order,lc_fitted-lc_thar, 'b+', alpha=.8)
ax.plot(lc_order[ind_good_thar],lc_fitted[ind_good_thar]-lc_thar[ind_good_thar], 'r+', alpha=.8)
plt.xlabel("Order")
plt.axis([61, 161, -0.1, 0.1])

ax = fig.add_subplot(222)
ax.plot(lc_coord,lc_fitted-lc_thar, 'b+', alpha=.8)
ax.plot(lc_coord[ind_good_thar],lc_fitted[ind_good_thar]-lc_thar[ind_good_thar], 'r+', alpha=.8)
plt.xlabel("CCD X Coordinate")
plt.xticks(np.arange(5)*1024)
plt.axis([0, 4096, -0.1, 0.1])

ax = fig.add_subplot(212)
ax.plot(lc_thar,lc_fitted-lc_thar, 'b+', alpha=.8)
ax.plot(lc_thar[ind_good_thar],lc_fitted[ind_good_thar]-lc_thar[ind_good_thar], 'r+', alpha=.8)
plt.xlabel("True ThAr wavelength")
plt.axis([3700, 10000, -0.1, 0.1])

fig.tight_layout()
fig.show()<jupyter_output><empty_output><jupyter_text>6. Process Science frames<jupyter_code>"""
Master variables:
    bias
    ap_uorder_interp
    flat1d
    wave_fitted
"""

# 1. read
sci = ccdproc.CCDData.read(logt['filename'][41], unit='adu')

# 2. BIAS correction
sci_bias = ccdproc.subtract_bias(sci, bias)

# 3. trim 
sci_bias_trim = ccdproc.trim_image(sci_bias[:, :4096])

# 4. scattered light correction
sci_bias_trim_sl = hrs.substract_scattered_light(sci_bias_trim, ap_uorder_interp, ap_width=10, shrink=.85)

# 5. extract 1D spectra
spec1d = hrs.extract_1dspec(sci_bias_trim_sl, ap_uorder_interp, ap_width=7)[0]

# 6. de-blaze
spec1d_db = spec1d/flat1d

# 7. roughly scale to 1.
spec1d_nm = spec1d_db/np.percentile(spec1d_db, 95, axis=1)[:, None]
plt.figure(figsize=(16, 5))
plt.plot(wave_fitted.T, spec1d_nm.T);
plt.xlim(5900, 6000)
plt.ylim(-0.1, 1.2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13459119496855346, 'fraction_numerical': 0.11932914046121593, 'mean_word_length': 4.9959778783308195, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 10, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '391191', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5578, 'n_tokens_neox': 4587, 'n_words': 909, 'path': 'example/example_basic.ipynb', 'repo_name': 'hypergravity/hrs', 'size': '11927'}