In [ ]:
# Test calc alpha ratio from spectra
from __future__ import print_function
import copy
import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt
from spectrum_overload import Spectrum
In [ ]:
wav_model = fits.getdata("/home/jneal/Phd/data/PHOENIX-ALL/PHOENIX/WAVE_PHOENIX-ACES-AGSS-COND-2011.fits")
wav_model /= 10 # nm
host = "/home/jneal/Phd/data/PHOENIX-ALL/PHOENIX/HD30501-lte05200-4.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits"
companion = "/home/jneal/Phd/data/PHOENIX-ALL/PHOENIX/HD30501b-lte02500-5.00-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits"
host = "/home/jneal/Phd/data/PHOENIX-ALL/PHOENIX/Z-0.0/lte05200-4.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits"
companion = "/home/jneal/Phd/data/PHOENIX-ALL/PHOENIX/Z-0.0/lte02600-5.00-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits"
host_f = fits.getdata(host)
comp_f = fits.getdata(companion)
plt.plot(wav_model, host_f, label="Host")
plt.plot(wav_model, comp_f, label="Companion")
plt.title("Phoenix spectra")
plt.xlabel("Wavelength (nm)")
plt.legend()
plt.show()
mask = (2000 < wav_model) & (wav_model < 2200)
wav_model = wav_model[mask]
host_f = host_f[mask]
comp_f = comp_f[mask]
plt.plot(wav_model, host_f, label="Host")
plt.plot(wav_model, comp_f, label="Companion")
plt.title("Phoenix spectra")
plt.legend()
plt.xlabel("Wavelength (nm)")
plt.show()
# Transform flux from /cm to /nm
host_f *= 1e-7
comp_f *= 1e-7
spec_host = Spectrum(xaxis=wav_model, flux=host_f, header=fits.getheader(host))
spec_comp = Spectrum(xaxis=wav_model, flux=comp_f, header=fits.getheader(companion))
print(spec_comp.flux)
normal_ratio = spec_host / spec_comp
print("normal_ratio", np.mean(normal_ratio.flux), np.std(normal_ratio.flux))
In [ ]:
def calc_alpha(model1, model2, chip=None, plot=False):
model1 = copy.copy(model1)
model2 = copy.copy(model2)
print("\nChip {}\n-------------------".format(chip))
"""Inherint flux ratio between the two.
Need to ransform from spectrum per area to surface area of each star.
For Phoenix Models - header parameters
PHXLUM - [W] Stellar luminosity
BUNIT - 'erg/s/cm^2/cm' Unit of flux
PHXREFF - [cm] Effective stellar radius
"""
def spec_area(spec):
# BUNIT 'erg/s/cm^2/cm' Unit of flux
# PHXREFF 67354000000.0 [cm] Effective stellar radius
radius = spec.header["PHXREFF"] * 1e-11 # in Gm gigametre
area = np.pi * radius ** 2
return area
area1 = spec_area(model1)
area2 = spec_area(model2)
# print(area1, area2)
area_ratio = area2 / area1
print("area_ratio =", area_ratio)
model1.flux = model1.flux* area1
model2.flux = model2.flux* area2
if chip in range(1,5):
ratio = np.nanmean(model2.flux / model1.flux)
print("spec ratio", ratio)
chip_limits = {1: [2111, 2124], 2: [2125, 2139], 3:[2140, 2152], 4: [2153,2169]}
model1.wav_select(*chip_limits[chip])
model2.wav_select(*chip_limits[chip])
chip_ratio = np.nanmean(model2.flux / model1.flux)
print("chip {} ratio = {}".format(chip, chip_ratio))
else:
model1.wav_select(2111, 2169)
model2.wav_select(2111, 2169)
full_ratio = np.nanmean(model2.flux / model1.flux)
print("full ratio", full_ratio)
lum_ratio = model2.header["PHXLUM"] / model1.header["PHXLUM"]
print("Lum ratio", lum_ratio)
if plot:
plt.figure
plt.plot(model1.xaxis, model1.flux, label="model1")
plt.plot(model2.xaxis, model2.flux, label="model2")
plt.legend()
plt.show()
model1.flux = model1.flux * area1
model2.flux = model2.flux * area2
plt.figure
plt.plot(model1.xaxis, model1.flux, label="model1")
plt.plot(model2.xaxis, model2.flux, label="model2")
plt.title("area scaled")
plt.legend()
plt.show()
plt.figure
plt.plot(full_ratio, label="model2/model1")
plt.plot(spec_ratio, label="model2/model1")
plt.plot(lum_ratio, label="model2/model1")
# plt.hlines(lum_ratio, 0, len(ratio), linestyle='--', label="luminosity ratio")
plt.legend()
plt.show()
return model1, model2, lum_ratio
In [ ]:
calc_alpha(spec_host, spec_comp, chip=None)
In [ ]:
calc_alpha(spec_host, spec_comp, chip=None)
calc_alpha(spec_host, spec_comp, chip=1)
calc_alpha(spec_host, spec_comp, chip=2)
calc_alpha(spec_host, spec_comp, chip=3)
calc_alpha(spec_host, spec_comp, chip=4)
In [ ]:
spec_host = Spectrum(xaxis=wav_model, flux=host_f/1e13, header=fits.getheader(host))
spec_comp = Spectrum(xaxis=wav_model, flux=comp_f/1e13, header=fits.getheader(companion))
spec_host2 = Spectrum(xaxis=wav_model, flux=host_f/1e13, header=fits.getheader(host))
spec_comp2 = Spectrum(xaxis=wav_model, flux=comp_f/1e13, header=fits.getheader(companion))
plt.plot(spec_host.xaxis, spec_host.flux, label="model1")
plt.plot(spec_host2.xaxis, spec_host2.flux, "--", label="model1b")
plt.plot(spec_comp.xaxis, spec_comp.flux, label="model2")
plt.plot(spec_comp2.xaxis, spec_comp2.flux, label="model2b")
plt.legend()
plt.show()
In [ ]:
# Get Radius and determine area
In [ ]:
r_host = spec_host.header["PHXREFF"] / 1e11
r_comp = spec_comp.header["PHXREFF"] / 1e11
print("r_host =", r_host, "\nr_comp =", r_comp)
print(type(r_host))
a_host = np.pi * r_host ** 2
a_comp = np.pi * r_comp ** 2
print("a_host =", a_host, "\na_comp =", a_comp)
print("a_ratio =", a_host/ a_comp, "or", a_comp / a_host)
In [ ]:
print(area1, area2, area2/area1, area1/area2)
spec_host.flux = spec_host.flux * a_host
spec_comp.flux = spec_comp.flux * a_comp
spec_ratio = spec_comp / spec_host
spec_ratio2 = spec_comp2 / spec_host2
plt.plot(spec_host.xaxis, spec_host.flux, label="model1")
#plt.plot(spec_host2.xaxis, spec_host2.flux / 1e4, "--", label="model1b")
plt.plot(spec_comp.xaxis, spec_comp.flux, label="model2")
plt.legend()
plt.show()
In [ ]:
plt.plot(spec_ratio.xaxis, spec_ratio.flux, label="with area")
plt.plot(spec_ratio2.xaxis, spec_ratio2.flux, label="without area")
plt.legend()
plt.show()
In [ ]:
np.mean(spec_ratio2.flux / spec_ratio.flux)
In [ ]: