Source code for exoplasim.makestellarspec

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, 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 = 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, **Do not use as an imported function; only call directly as a command-line program.** **Usage** >>> python 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') #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]") 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') 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]") 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$]")
if __name__=="__main__" and (Path(sys.argv[0]).name!="sphinx-build" and Path(sys.argv[0]).name!=""): main()