import numpy as np
import sys
from pathlib import Path
import argparse as ag
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()
def _processArgs():
parser = ag.ArgumentParser(description="""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.""")
parser.add_argument("spectrumfile",help="Path to the spectrum file",)
parser.add_argument("name",help="Name of the star")
parser.add_argument("-p","--plot",action="store_true",help="Display plots of the converted spectrum")
parser.add_argument('-w',"--numwavelengths",default=2048,type=int,help="Number of wavelengths in the hi-res spectrum")
parser.add_argument('-n',"--normalize",action="store_true",help="Normalize the output spectra")
args = parser.parse_args()
return args
[docs]
def convert(spectrumfile, name, plot=False, numwavelengths=2048, normalize=False):
"""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.
Parameters
----------------
spectrumfile : str
Path to the file to convert
name : str
Identifying name for the star to use in filenames
plot : bool, optional
If set, will plot the spectra. Requires matplotlib
numwavelengths : int, optional
Number of wavelengths to use for the hi-res spectrum. Default 2048
normalize : bool, optional
If set, output spectra will be normalized.
Yields
------
mystarsname.dat
965 wavelengths
mystarsname_hr.dat
2048 wavelengths
"""
sfile = spectrumfile
#wave0 = float(sys.argv[3]) #microns
#wave1 = float(sys.argv[4]) #microns
numw = numwavelengths
norm = normalize
w,f,u = readspec(sfile)
print(u,w.min(),w.max())
if plot:
import matplotlib.pyplot as plt
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))])
if plot:
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:]
if plot:
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)
if plot:
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()
[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 [spectrumfile] [name] --plot --numwavelengths --normalize
Parameters
----------------
-p, --plot : bool, optional
If set, will plot the spectra. Requires matplotlib
-w, --numwavelengths : int, optional
Number of wavelengths to use for the hi-res spectrum. Default 2048
-n, --normalize : bool, optional
If set, output spectra will be normalized.
Yields
------
mystarsname.dat
965 wavelengths
mystarsname_hr.dat
2048 wavelengths
"""
args = _processArgs()
if args.plot:
import matplotlib.pyplot as plt
sfile = args.spectrumfile
name = args.name
#wave0 = float(sys.argv[3]) #microns
#wave1 = float(sys.argv[4]) #microns
numw = args.numwavelengths
norm = args.normalize
w,f,u = readspec(sfile)
print(u,w.min(),w.max())
if args.plot:
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))])
if args.plot:
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:]
if args.plot:
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)
if args.plot:
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()