<jupyter_start><jupyter_text>20190529 Looking at LRS2 data. Rough here, will probably turn into a script to make coadded spectra.<jupyter_code>runs = glob.glob('programs/ut*')
runs
run3 = runs[0]
spectra = glob.glob(run3+'/spectrum*')
hdul = fits.open(spectra[0])
hdul.info()
hdul[0].header
stars  = []
for spectrum in spectra:
    star = fits.getheader(spectrum,0)['QOBJECT']
    stars.append(star)
stars = np.unique(stars)
stars
#star = stars[1]
uvspec = glob.glob(run3+'/*uv.fits')
ospec = glob.glob(run3+'/*orange.fits')
def spectra_adder(fluxes, errors):
    """
    combines the flux at each overlap
    """
    weight_f = np.average(fluxes, axis =0, weights=(1/errors))
    weight_e = np.average((weight_f - fluxes)**2, axis=0, weights = (1/errors))**0.5
    return weight_f, weight_e
stars[1]
fs = []
es = []
for sp in ospec:
      if fits.getheader(sp,0)['QOBJECT'] == stars[1]:
        #print(sp)
        data = fits.getdata(sp,0)
        fs.append(data[1])
        es.append(data[3]) #2 is sky
        w = data[0] #fixed w scale, no need to interpolate
f, e = spectra_adder(np.array(fs),np.array(es))

plt.plot(w,f)        
ow, of, oe, = w, f, e
fs = []
es = []
for sp in uvspec:
      if fits.getheader(sp,0)['QOBJECT'] == stars[1]:
        #print(sp)
        data = fits.getdata(sp,0)
        fs.append(data[1])
        es.append(data[3]) #2 is sky
        w = data[0] #fixed w scale, no need to interpolate
f, e = spectra_adder(np.array(fs),np.array(es))

plt.plot(w,f)        
uw, uf, ue = w, f, e
plt.plot(uw, uf)
plt.plot(ow, of)
#plt.xlim(4600, 4800)<jupyter_output><empty_output><jupyter_text>(Whispers) how is the absolute flux calibration?<jupyter_code>phottab = 'gj649_phot.vot'
from astropy.io.votable import parse
votable = parse(phottab)
table = votable.get_first_table()
data = table.array
c = const.c.value
masked_data = data.data
w = []
f = []
e = []
filters = []
for i in range(len(masked_data['sed_freq'])):
    w.append(10.*c/masked_data['sed_freq'][i])
    f.append(2.99792458E-05*masked_data['sed_flux'][i]/(w[i]**2))
    e.append((masked_data['sed_eflux'][i]*f[i])/masked_data['sed_flux'][i])
    filters.append(masked_data['sed_filter'][i])
#print(w,f,e)
w, f ,e = np.array(w), np.array(f), np.array(e)
arr1inds = w.argsort()
w = w[arr1inds]
f = f[arr1inds]
e = e[arr1inds]


plt.figure('sed', figsize=(10,10))
for i in range(len(w)):
    if e[i] !=0.0:
        plt.errorbar(w[i], f[i], yerr=e[i], marker='o', ls='none')
    else:
        plt.errorbar(w[i], f[i], yerr=f[i]*0.1, marker='o', ls='none', uplims=True)
plt.yscale('log')
plt.xscale('log')
mask = (w > uw[0]) & (w < ow[-1]) 
w, f, e = w[mask], f[mask], e[mask]
plt.errorbar(w,f, yerr=e)
plt.plot(uw, uf)
plt.plot(ow, of)
plt.scatter(w,f, c='C2', zorder=10)<jupyter_output><empty_output><jupyter_text>It's...it's alright. Photometry is all over the place, that's more of a problem. What about the STIS spectrum? (Is there one yet)? Yes.<jupyter_code>stisccd = '../MegaMUSCLES/GJ_649/HST/STIS/odlm12010_sx1.fits'
data = fits.getdata(stisccd,1)
ws, fs = data['WAVELENGTH'], data['FLUX']
plt.step(ws, fs)
mask = ws > 3500
ws, fs = ws[mask], fs[mask]
plt.plot(uw, uf)
plt.plot(ow, of)
plt.plot(ws, fs)<jupyter_output><empty_output><jupyter_text>Bit concerned about the relative flux at the blue end...<jupyter_code>plt.plot(uw, uf)
plt.plot(ow, of)
plt.plot(ws, fs)
plt.yscale('log')
plt.plot(uw[:-1], np.diff(uw))
plt.plot(ow[:-1], np.diff(ow))
for wi, fi in zip([uw, ow],[uf, of]):
    mask = (ws > wi[0]) & (ws < wi[-1])
    newf = interpolate.interp1d(ws[mask], fs[mask], fill_value ='extrapolate')(wi[wi < ws[-1]])
    scale = fi[wi < ws[-1]]/newf
    plt.plot(wi[wi < ws[-1]], scale/np.mean(scale))
plt.xlim(3600,5700)
plt.axhline(1, ls ='--', c='r')<jupyter_output><empty_output><jupyter_text>Great, looks line we can get away with a fixed normfac<jupyter_code>for wi, fi in zip([uw, ow],[uf, of]):
    smask = (ws > wi[0]) & (ws < wi[-1])
    lmask = (wi < ws[-1])
    lf = np.trapz(fi[lmask], wi[lmask])
    sf = np.trapz(fs[smask], ws[smask])
    normfac = sf/lf
    print(normfac)
    plt.plot(wi, fi*normfac)

plt.step(ws[:-1], fs[:-1])<jupyter_output><empty_output><jupyter_text>Perfect.What about the overlap?<jupyter_code>for wi, fi in zip([uw, ow],[uf, of]):
    smask = (ws > wi[0]) & (ws < wi[-1])
    lmask = (wi < ws[-1])
    lf = np.trapz(fi[lmask], wi[lmask])
    sf = np.trapz(fs[smask], ws[smask])
    normfac = sf/lf
    print(normfac)
    plt.plot(wi, fi*normfac)

plt.step(ws[:-1], fs[:-1])
plt.xlim(4500, 4800)<jupyter_output>1.299439270267638
1.2362533847921378<jupyter_text>It's very small. Could just pick one or the other, easer than averaging. UV has higher res? Going to turn this into a script, LRS2_add.py  Trying the other star, causing problems in the script...<jupyter_code>fs = []
es = []
for sp in uvspec:
      if fits.getheader(sp,0)['QOBJECT'] == stars[0]:
        #print(sp)
        data = fits.getdata(sp,0)
        fs.append(data[1])
        es.append(data[3]) #2 is sky
        w = data[0] #fixed w scale, no need to interpolate
f, e = spectra_adder(np.array(fs),np.array(es))

plt.plot(w,f)<jupyter_output><empty_output><jupyter_text>Ok it exists, must be a script issue. Fixed.New issue, older spectra don't have QOBJECT.<jupyter_code>runs
hdul.close()
hdul = fits.open(runs[1]+'/'+os.listdir(runs[1])[0])
hdul[0].header<jupyter_output><empty_output><jupyter_text>Ok, sadly it looks like it will have to be 'OBJECT' To do - move all data to one directory, no point splitting it up by program. Put stuff on github. Sort spectra by date. 20190531 Got some new data. GJ649 has a "W" and "E" observation, wonder how they're different. LHS 2686 is flaring strongly, check sub-exposures.<jupyter_code>datafiles = glob.glob('spectra/*')
gj649_spec = glob.glob('spectra/GJ649*')
gj649_spec
for s in gj649_spec:
    data = Table.read(s)
    plt.plot(data['WAVELENGTH'], data['FLUX'])<jupyter_output><empty_output><jupyter_text>Slight offset in flux calibration, coadd should be fine.<jupyter_code>from lrs2_add import stars_finder as stars_finder
stars = stars_finder('data/')
stars
star = stars[-1]
for sp in glob.glob('data/spectrum*'):
    if fits.getheader(sp,0)['OBJECT'] == star:
        data = fits.getdata(sp,0)
        plt.plot(data[0], data[1])
plt.xlim(3600, 7100)<jupyter_output><empty_output><jupyter_text>No real difference :(<jupyter_code>print(normfac)<jupyter_output>1.2362533847921378
{'fraction_non_alphanumeric': 0.17982592477463474, 'fraction_numerical': 0.036369288156667706, 'mean_word_length': 4.305028854080791, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '227395', 'license': 'MIT', 'n_tokens_mistral': 2844, 'n_tokens_neox': 2636, 'n_words': 580, 'path': 'lrs2_rough.ipynb', 'repo_name': 'davidjwilson/lrs2_fun', 'size': '6434'}