Source code for exoplasim.makestellarspec

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()