<jupyter_start><jupyter_text>Lyman alpha fittingThe equation for the cross section is (Diplas \& Savage 1994)$ \sigma(\lambda) = 4.26e^{-20} \text{cm}^2 / [6.04e^{-10}+ (\lambda - \lambda_0)^2] $where $\lambda_0 = 1215.67 \AA$. We might need a velocity doppler shift parameter for $\lambda_0$. The model for the spectrum, with $N = N(\text{H I})$ as a parameter, is written as$ f_m(\lambda; N) = f_c(\lambda) \exp(-N \sigma(\lambda)) $.The function $f_c(\lambda)$ describes the continuum, the flux without any lines. It is usually a constant or 1st order approximation, determined by fitting or sometimes by eye. $N(\text{HI})$ can be found by minimizing a $\chi^2$ function$ \chi^2(N) = \frac{1}{n - 1} \sum_i \sigma_c^{-1} (f_i - f_m(\lambda_i; N))^2$. Here, a suitable noise model $\sigma_c$ needs to be chosen. In DS94, $\sigma_c = \text{constant} \times \exp(\sigma(\lambda_i) N)$. This is not an error on the measurement, but rather a factor to treat the other features of the spectrum as noise. For the constant, they use the RMS deviation with respect to $f_c$, over a wavelength range outside the line, deemed to have similar noise as the line wings. The workflow goes as follows:1. Find a suitable approximation for the continuum.2. Find a suitable noise level with respect to the continuum.2. Mask wavelengths that have troublesome features.3. Optimize the $\chi^2$ function for the remaining data points $(\lambda_i, f_i)$, minimizing the differences between $f_i$ and $f_m(\lambda_i)$.<jupyter_code>def sigma(l):
    l0 = 1215.67
    return 4.26e-20 / (6.04e-10 + np.square(l - l0))
    
def fmodel(fcontinuum, NHI):
    return lambda l : fcontinuum(l) * np.exp(-NHI * sigma(l))

def chi2(NHI, fc, sigma_c, wavs, fs):
    extinctions = np.exp(NHI * sigma(l))
    deltas = fc(wavs) - fs * extinctions
    sigmas = sigma_c * extinctions
    chi2 = np.square((deltas / sigmas)).sum() / (len(deltas) - 1)
    return chi2

from matplotlib import pyplot as plt
import numpy as np
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 226

%matplotlib inline

def prepare_axes():
    plt.xlabel('wavelength')
    plt.ylabel('flux (erg cm$^{-2}$s$^{-1}\\AA^{-1}$)')
    
def plot_iue_order(table, index):
    data_column_name = 'ABS_CAL'
    flux = table[index][data_column_name]
    wavs = [table[index]['WAVELENGTH'] + table[index]['DELTAW'] * i for i in range(len(flux))]
#    dq = table[index]['QUALITY']
 #   good = dq == -2  
  #  plt.plot(wavs[good], flux[good])
    plt.plot(wavs, flux)
    
def plot_stis_order(table, index):
    wavs = table[index]['WAVELENGTH']
    data_column_name = 'FLUX'
    flux = table[index][data_column_name]
    plt.plot(wavs, flux)
    
def plot_function_at_lya(f, ylabel='f(x)', *args):
    x = np.linspace(1100, 1300, 1000)
    y = f(x, *args)
    plt.xlabel('wavelength')<jupyter_output><empty_output><jupyter_text>Data exploration<jupyter_code># open some data
from astropy.io import fits
from astropy.table import Table
from pathlib import Path
def sort_by_wav(wavs, flux):
    idxs = np.argsort(wavs)
    return wavs[idxs], flux[idxs]<jupyter_output><empty_output><jupyter_text>STIS data looks like this<jupyter_code>spectrum_file = '/Users/dvandeputte/Projects/FUSE H2/plotting_git/data/HD094493/mastDownload/HST/o54306010/o54306010_x1d.fits'
fits.info(spectrum_file)
t = Table.read(spectrum_file)
t
t['ERROR'][20]
import numpy as np
def stis_merge_orders(t, extra_columns=None):
    """
    Put the data from all spectral orders of a set of STIS data into one big array
    
    t: Table containing STIS data
    
    extra_columns: list of strings, indicating columns other than 'WAVELENGTH' and 'FLUX' to be merged and sorted
    
    Returns
    -------
        allwavs, allflux: ndarray, wavelengths and fluxes, sorted by wavelength
        
        list of arrays: if extra_columns is not None, the requested columns merged and sorted in the same way
    """
    allwavs = np.concatenate(t['WAVELENGTH'])
    allflux = np.concatenate(t['FLUX'])
    idxs = np.argsort(allwavs)
    
    if extra_columns is None:
        return allwavs[idxs], allflux[idxs]
    else:    
        extras = []
        for name in extra_columns:
            if not isinstance(t[name], np.ndarray):
                print("This STIS data column does not have per-wavelength data: ", name)
                raise
            extras.append(np.concatenate(t[name]))
        return allwavs[idxs], allflux[idxs], [array[idxs] for array in extras]

def stis_binned(t, numbins):
    wavmin = 1150
    wavmax = 1280
    wavbins = np.linspace(wavmin, wavmax, numbins, endpoint=True)

    wavs, flux, (errs,) = stis_merge_orders(t, extra_columns=['ERROR'])
    
    # np.digitize returns list of indices. 1 means that the data point is between wav[0] (first) and wav[1]
    # n-1 means between wav[n-2] and wav[n-1] (last). 0 and n mean out of range.
    bs = np.digitize(wavs, wavbins)
    newwavs = np.zeros(len(wavbins) - 1)
    newflux = np.zeros(len(wavbins) - 1)
    for i in range(0, len(wavbins) - 1): 
        b = i + 1 # from 1 to n-1
        in_bin = bs == i+1
        
        weights = 1 / np.square(errs[in_bin])
        newwavs[i] = np.average(wavs[in_bin], weights=weights)
        newflux[i] = np.average(flux[in_bin], weights=weights)
        
    return newwavs, newflux
    
# try plotting this
x, y = stis_merge_orders(t)
xbin, ybin = stis_binned(t, 1000)
plt.plot(x, y)
plt.plot(xbin, ybin)<jupyter_output>/Users/dvandeputte/opt/anaconda3/lib/python3.8/site-packages/numpy/lib/function_base.py:412: RuntimeWarning: invalid value encountered in true_divide
  avg = np.multiply(a, wgt, dtype=result_dtype).sum(axis)/scl
<ipython-input-32-5a40e05ef52c>:48: UserWarning: Warning: converting a masked element to nan.
  newwavs[i] = np.average(wavs[in_bin], weights=weights)
<ipython-input-32-5a40e05ef52c>:49: UserWarning: Warning: converting a masked element to nan.
  newflux[i] = np.average(flux[in_bin], weights=weights)<jupyter_text>IUE high res like this<jupyter_code>spectrum_file = '/Users/dvandeputte/Projects/FUSE H2/plotting_git/data/HD197770/swp49267.mxhi.gz'
fits.info(spectrum_file)
t = Table.read(spectrum_file)
t
# check the contents of one row
row = 1
print(t[row]['WAVELENGTH'])
t[row]['NET'].shape
dq = t[row]['QUALITY']
dq == -2
def iue_h_merge_orders(t):
    
    def iue_wavs(i):
        return t[i]['WAVELENGTH'] + t[i]['DELTAW'] * np.arange(t[i]['NPOINTS'])
    
    def iue_flux(i):
        pixrange = slice(t[i]['STARTPIX'], t[i]['STARTPIX'] + t[i]['NPOINTS'])
        return t[i]['ABS_CAL'][pixrange]
    
    allwavs = np.concatenate([iue_wavs(i) for i in range(len(t))])
    allflux = np.concatenate([iue_flux(i) for i in range(len(t))])
    
    return sort_by_wav(allwavs, allflux)
allwavs, allflux = iue_h_merge_orders(t)
plt.plot(allwavs, allflux)
# separate orders
#for i in range(1, 30):
#    plot_iue_order(t, i)
    
plot_iue_order(t, 14)<jupyter_output><empty_output><jupyter_text>IUE low res<jupyter_code>spectrum_file = "/Users/dvandeputte/Projects/FUSE H2/plotting_git/data/BD+52d3210/swp34153mxlo_vo.fits"
t = Table.read(spectrum_file)
t<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13264306702112774, 'fraction_numerical': 0.02700433748425913, 'mean_word_length': 3.8005372733378104, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '354943', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2826, 'n_tokens_neox': 2632, 'n_words': 682, 'path': 'lyafitting.ipynb', 'repo_name': 'drvdputt/dust_fuse_h2', 'size': '7147'}