import numpy as np
import matplotlib.pyplot as plt
import sys
from pathlib import Path
cc = 299792458.0 #Speed of light
[docs]def readspec(sfile,cgs=False):
"""Read a Phoenix stellar spectrum and return wavelengths, fluxes, and units
Takes as input a spectrum file generated by the Phoenix stellar spectrum web simulator,
https://phoenix.ens-lyon.fr/simulator-jsf22-26.
Parameters
----------
sfile : str
Path to the spectrum file
cgs : bool, optional
Whether or not we should try to use CGS units (irrelevant for exoplasim)
Returns
-------
numpy.ndarray, numpy.ndarray, str
Returns wavelengths, fluxes, and the units
"""
f=open(sfile,"r")
ftxt = f.read()
if ftxt[0]=='#':
newstyle=True
else:
newstyle=False
f.close()
if newstyle:
lines = ftxt.split('\n')[:-2] #chop off timing line
nhl = 0
for entry in lines:
if entry[0]=='#':
nhl+=1
else:
break
units = lines[nhl-2].split()[-1]
if lines[nhl-2].split()[1]=="(ANGSTROM)":
cgs=False
nhi = 0
for entry in lines[nhl:]:
if float(entry.split()[0])<0.1: #angstroms
nhi+=1
else:
break
print("Cutting at %d: "%(nhl+nhi),lines[nhl+nhi])
lines = lines[nhl+nhi:]
wvls = np.zeros(len(lines))
fluxes = np.zeros(len(lines))
n=0
for l in lines:
wvls[n] = float(l.split()[0])
fluxes[n] = float(l.split()[1])
n+=1
#print l.split()
else:
lines = ftxt.split('\n')[:-1]
wvls = np.zeros(len(lines))
fluxes = np.zeros(len(lines))
n=0
for l in lines:
wvls[n] = float(l.split()[0])
fluxes[n] = 10**(float(l.split()[1].replace('D','E'))-8.0)
n+=1
units='ergs/sec/cm^2/angstrom'
if not cgs:
units='W/m^2/micron'
wvls[:] = wvls[:]*1.0e-4 #microns
fluxes[:] = fluxes[:]*10.0 #W/m^2/micron
return wvls,fluxes,units
def _coarsen(wvls,fluxes,w1,w2,num=5000):
inds = np.arange(len(wvls))
nlow = (inds[wvls>w1])[0]
nhigh= (inds[wvls>w2])[0]
inc = int(nhigh-nlow)//num
waves = np.zeros(num)
energ = wvls*fluxes
flxs2 = np.zeros(num)
for n in range(0,num):
v0 = nlow+n*inc-1
v1 = nlow+n*inc
intg = 0.0
for v in range(0,inc):
intg+=0.5*(wvls[v1+v]-wvls[v0+v])*(fluxes[v1+v]+fluxes[v0+v])
dnu = wvls[v1+inc]-wvls[v0]
waves[n] = 0.5*(wvls[v1+inc]+wvls[v0])
flxs2[n] = intg/dnu
return waves,flxs2
def _integrate_trap(x,y):
net = 0
for n in range(1,len(y)):
net += 0.5*(x[n]-x[n-1])*(y[n]+y[n-1])
return net
def _normalize(wvls,fluxes):
netf = np.trapz(fluxes,x=wvls)
factor = 1366.941858/netf
return fluxes*factor
[docs]def writedat(wvls,fluxes,name):
"""Write wavelengths and fluxes to a .dat file ExoPlaSim can read.
Parameters
----------
wvls : array-like
An array-like object containing a monotonic list of wavelengths going from short to long
fluxes : array-like
An array-like object containing a flux for each wavelength in ``wvls``.
name : str
A name to use when generating output files.
Yields
------
name.dat
The provided spectrum in a format ExoPlaSim can read.
"""
f=open(name+".dat","w")
sdat = ' Wavelength Flux \n'
for n in range(0,len(wvls)):
sdat+=str(wvls[n])+' '+str(fluxes[n])+'\n'
f.write(sdat)
f.close()
[docs]def main():
"""Convert spectral files to exoplasim-compliant formats, including resampling to the necessary resolutions.
Must give as input a spectrum file generated by the Phoenix stellar spectrum web simulator,
https://phoenix.ens-lyon.fr/simulator-jsf22-26. **Do not use as an imported function; only call directly as a command-line program.**
**Usage**
>>> python makestellarspec.py spectrum.txt mystarsname
Yields
------
mystarsname.dat
965 wavelengths
mystarsname_hr.dat
2048 wavelengths
"""
sfile = sys.argv[1]
name = sys.argv[2]
#wave0 = float(sys.argv[3]) #microns
#wave1 = float(sys.argv[4]) #microns
try:
numw = int(sys.argv[3])
except:
numw = int(2048)
try:
norm = sys.argv[4]
except:
norm = False
w,f,u = readspec(sfile)
print(u,w.min(),w.max())
plt.plot(w,f)
plt.xscale('log')
plt.yscale('log')
plt.show()
#wc,fc = _coarsen(w,f,wave0,wave1,num=numw*4)
#plt.plot(wc,fc)
#plt.xscale('log')
#plt.yscale('log')
#plt.xlabel("$\lambda$ [$\mu$m]")
#plt.ylabel("$F_\lambda$ [W/m$^2$/$\mu$m]")
#plt.show()
w2 = np.concatenate([np.geomspace(0.2,0.75,num=int(numw/2)+1)[:-1],np.geomspace(0.75,100.0,num=int(numw/2))])
plt.plot(w2)
plt.yscale('log')
plt.show()
f2 = np.interp(w2,w,f)
E0 = np.trapz(f,x=w)
E1 = np.trapz(f2,x=w2)
print(abs(E0-E1)/E0)
if norm:
factor = 1366.941858/E1
f2*=factor
print(np.trapz(f2,x=w2))
print(np.trapz(_normalize(w,f),x=w))
#print f2[np.argwhere(w2>40.0)],w[-10:]
plt.plot(w2,f2)
#plt.plot(wc,fc)
plt.xscale('log')
plt.yscale('log')
plt.xlabel("$\lambda$ [$\mu$m]")
plt.ylabel("$F_\lambda$ [W/m$^2$/$\mu$m]")
plt.show()
wvref = np.loadtxt(Path(__file__).parent.resolve()+"/wvref.txt")
f3 = np.interp(wvref,w2,f2)
writedat(w2,f2/cc,name+"_hr")
writedat(wvref,f3/cc,name)
plt.plot(wvref,f3)
plt.xscale('log')
plt.yscale('log')
plt.xlabel("$\lambda$ [$\mu$m]")
plt.ylabel("$\lambda F_\lambda$ [W/m$^2$]")
plt.show()
if __name__=="__main__" and (Path(sys.argv[0]).name!="sphinx-build" and
Path(sys.argv[0]).name!="build.py"):
main()