Python model.ModelOutput类代码示例

def export_to_fits(cli):
	
	#
	# Read in the model:
	#
	file = filename(cli, "plot")
	file += ".rtout"
	model = ModelOutput(file)
	
	
	#
	# Write fits file:
	#
	if(cli.mode == "images"):
		
		los = [0 for i in range(3)]
		los[0] = 'x'
		los[1] = 'y'
		los[2] = 'z'

		for k in range(0, 3):
			image = model.get_image(distance=1*pc, units='MJy/sr', inclination=0, component='total', group=k)
			Nwavelength=image.val.shape[2]
			for i in range(0, Nwavelength):
				file = filename(cli, "fits")
				file += "_wavelength=" + str(image.wav[i]) + "micron_los=" + los[k] + ".fits"
				fits.writeto(file, image.val[:, :, i], clobber=True)
				if(cli.verbose):
					print("  The fits file was written to", file)

	else:
		print("ERROR: The specified mode", mode, "is not available. Use 'images' only.")

	def calcChi2(self,dist_pc=140,extinction=0, sourcename='Oph.1'):
		self.dist=dist_pc*pc
		self.extinction=extinction
		chi = np.loadtxt('kmh94_3.1_full.chi')
		wav = np.loadtxt('kmh94_3.1_full.wav')
		Chi = interp1d(wav,chi,kind='linear')
		modelname = self.folder+self.name
		self.mo = ModelOutput(modelname+'.rtout')
		
		# get the sed of all inclination
		sed = self.mo.get_sed(aperture=-1, inclination='all', distance=self.dist,units='Jy')
				
		# calculate the optical depth at all wavelengths
		tau = self.extinction*Chi(sed.wav)/Chi(0.550)/1.086
		
		# calculate extinction values
		ext = np.array([np.exp(-tau) for i in range(sed.val.shape[0])])
		
		# apply extinction to model
		extinct_values = np.log10(sed.val.transpose()*ext.T)
		
		# data points and errors
		folder_export="/n/a2/mrizzo/Dropbox/SOFIA/Processed_Data/"
		sourcetable = pickle.load(open(folder_export+"totsourcetable_fits.data","r"))
		TwoMASS = ['j','h','ks']
		uTwoMASS = ["e_"+col for col in TwoMASS]
		wlTwoMASS = [1.3,1.6,2.2]
		labelTwoMASS = '2MASS'
		Spitzer = ['i1','i2','i3','i4']
		uSpitzer = ["e_"+col for col in Spitzer]
		wlSpitzer = [3.6,4.5,5.8,8.]
		labelSpitzer = 'Spitzer'
		SOFIA = ['F11','F19','F31','F37']
		uSOFIA = ["e_"+col for col in SOFIA]
		wlSOFIA = [11.1,19.7,31.5,37.1]
		labelSOFIA = 'SOFIA'
		sources = sourcetable.group_by('SOFIA_name')
		for key,source in zip(sources.groups.keys,sources.groups):
			if sourcename == source['SOFIA_name'][0]:	
				datapoints = source[TwoMASS+Spitzer+SOFIA]
				dataerrors = source[uTwoMASS+uSpitzer+uSOFIA]
				print p.nptable(datapoints),p.nptable(dataerrors)
				
				# calculate log10 of quantities required for chi squared
				logFnu = np.log10(p.nptable(datapoints))-0.5*(1./np.log(10.))*p.nptable(dataerrors)**2/p.nptable(datapoints)**2
				varlogFnu = (1./np.log(10)/p.nptable(datapoints))**2*p.nptable(dataerrors)**2
				print extinct_values,extinct_values.shape
				
				# for each inclination, calculate chi squared; need to interpolate to get model at required wavelengths
				Ninc = extinct_values.shape[1]
				chi2 = np.zeros(Ninc)
				wl=wlTwoMASS+wlSpitzer+wlSOFIA
				N = len(wl)
				for j in range(Ninc):
					interp_func = interp1d(sed.wav,extinct_values[:,j],kind='linear')
					interp_vals = interp_func(wl)
					chi2[j] = 1./N * np.sum((logFnu - interp_vals)**2/varlogFnu)
					
				print chi2

def extract(model):

    # Extract model name
    model_name = os.path.basename(model).replace('.rtout', '').replace('external_', '')

    m = ModelOutput(model)

    wav, flux = m.get_image(group=0, units='MJy/sr', distance=1000. * kpc)  # distance should not matter as long as it is large
    flux = flux[0, :, :, :]

    # Convolve with filters
    flux_conv = np.zeros((len(filters), flux.shape[0], flux.shape[1]))
    for i, filtname in enumerate(filters):
        transmission = rebin_filter(filtname, c / (wav * 1.e-4))
        flux_conv[i, :, :] = np.sum(transmission[np.newaxis, np.newaxis:] * flux, axis=2)

    pyfits.writeto('models/external/external_%s.fits' % model_name, flux, clobber=True)
    pyfits.writeto('models/external/external_%s_conv.fits' % model_name, flux_conv, clobber=True)

def getRadialDensity(rtout, angle, plotdir):
    """
    """
    import numpy as np
    from hyperion.model import ModelOutput


    m = ModelOutput(rtout)
    q = m.get_quantities()
    r_wall = q.r_wall; theta_wall = q.t_wall; phi_wall = q.p_wall
    # get the cell coordinates
    rc = r_wall[0:len(r_wall)-1] + 0.5*(r_wall[1:len(r_wall)]-r_wall[0:len(r_wall)-1])
    thetac = theta_wall[0:len(theta_wall)-1] + \
             0.5*(theta_wall[1:len(theta_wall)]-theta_wall[0:len(theta_wall)-1])
    phic = phi_wall[0:len(phi_wall)-1] + \
           0.5*(phi_wall[1:len(phi_wall)]-phi_wall[0:len(phi_wall)-1])
    #
    rho = q['density'].array[0]

    # find the closest angle in the thetac grid
    ind = np.argsort(abs(thetac-angle*np.pi/180.))[0]

    return rc, rho[0,ind,:]

def temp_hyperion(rtout,outdir, bb_dust=False):
    import numpy as np
    import matplotlib as mpl
    mpl.use('Agg')
    import matplotlib.pyplot as plt
    import os
    from hyperion.model import ModelOutput
    import astropy.constants as const
    from matplotlib.colors import LogNorm

    # seaborn colormap
    import seaborn.apionly as sns

    # constants setup
    AU = const.au.cgs.value

    # misc variable setup
    print_name = os.path.splitext(os.path.basename(rtout))[0]

    m = ModelOutput(rtout)
    q = m.get_quantities()

    # get the grid info
    ri, thetai = q.r_wall, q.t_wall
    rc     = 0.5*( ri[0:len(ri)-1]     + ri[1:len(ri)] )
    thetac = 0.5*( thetai[0:len(thetai)-1] + thetai[1:len(thetai)] )

    # get the temperature profile
    # and average across azimuthal angle
    # temperature array in [phi, theta, r]
    temp = q['temperature'][0].array.T
    temp2d = np.sum(temp**2, axis=2)/np.sum(temp, axis=2)
    temp2d_exp = np.hstack((temp2d,temp2d,temp2d[:,0:1]))
    thetac_exp = np.hstack((thetac-np.pi/2, thetac+np.pi/2, thetac[0]-np.pi/2))

    mag = 1
    fig = plt.figure(figsize=(mag*8,mag*6))
    ax = fig.add_subplot(111, projection='polar')

    # cmap = sns.cubehelix_palette(light=1, as_cmap=True)
    cmap = plt.cm.CMRmap
    im = ax.pcolormesh(thetac_exp, rc/AU, temp2d_exp, cmap=cmap, norm=LogNorm(vmin=5, vmax=100))
    #
    # cmap = plt.cm.RdBu_r
    # im = ax.pcolormesh(thetac_exp, np.log10(rc/AU), temp2d_exp/10, cmap=cmap, norm=LogNorm(vmin=0.1, vmax=10))
    #
    print temp2d_exp.min(), temp2d_exp.max()
    im.set_edgecolor('face')

    ax.set_xlabel(r'$\rm{Polar\,angle\,(Degree)}$',fontsize=20)
    # ax.set_ylabel(r'$\rm{Radius\,(AU)}$',fontsize=20, labelpad=-140, color='grey')
    # ax.set_ylabel('',fontsize=20, labelpad=-140, color='grey')
    ax.tick_params(labelsize=16)
    ax.tick_params(axis='y', colors='grey')
    ax.set_yticks(np.hstack((np.arange(0,(int(max(rc)/AU/10000.)+1)*10000, 10000),max(rc)/AU)))
    #
    # ax.set_yticks(np.log10(np.array([1, 10, 100, 1000, 10000, max(rc)/AU])))
    #
    ax.set_yticklabels([])
    ax.grid(True, color='LightGray', linewidth=1.5)
    # ax.grid(True, color='k', linewidth=1)

    ax.set_xticklabels([r'$\rm{90^{\circ}}$',r'$\rm{45^{\circ}}$',r'$\rm{0^{\circ}}$',r'$\rm{-45^{\circ}}$',\
                            r'$\rm{-90^{\circ}}$',r'$\rm{-135^{\circ}}$',r'$\rm{180^{\circ}}$',r'$\rm{135^{\circ}}$'])
    cb = fig.colorbar(im, pad=0.1)
    cb.ax.set_ylabel(r'$\rm{Averaged\,Temperature\,(K)}$',fontsize=20)
    cb.set_ticks([5,10,20,30,40,50,60,70,80,90,100])
    cb.set_ticklabels([r'$\rm{5}$',r'$\rm{10}$',r'$\rm{20}$',r'$\rm{30}$',r'$\rm{40}$',r'$\rm{50}$',r'$\rm{60}$',r'$\rm{70}$',r'$\rm{80}$',r'$\rm{90}$',r'$\rm{>100}$'])
    #
    # cb.ax.set_ylabel(r'$\rm{log(T/10)}$',fontsize=20)
    # cb.set_ticks([0.1, 10**-0.5, 1, 10**0.5, 10])
    # cb.set_ticklabels([r'$\rm{-1}$',r'$\rm{-0.5}$',r'$\rm{0}$',r'$\rm{0.5}$',r'$\rm{\geq 1}$'])
    #
    cb_obj = plt.getp(cb.ax.axes, 'yticklabels')
    plt.setp(cb_obj,fontsize=20)

    # fix the tick label font
    ticks_font = mpl.font_manager.FontProperties(family='STIXGeneral',size=20)
    for label in ax.get_yticklabels():
        label.set_fontproperties(ticks_font)

    fig.savefig(outdir+print_name+'_temperature.png', format='png', dpi=300, bbox_inches='tight')
    fig.clf()

    # Plot the radial temperature profile
    fig = plt.figure(figsize=(12,9))
    ax = fig.add_subplot(111)

    plot_grid = [0,99,199]
    label_grid = [r'$\rm{outflow}$', r'$\rm{45^{\circ}}$', r'$\rm{midplane}$']
    alpha = np.linspace(0.3,1.0,len(plot_grid))
    color_list = [[0.8507598215729224, 0.6322174528970308, 0.6702243543099417],\
                  [0.5687505862870377, 0.3322661256969763, 0.516976691731939],\
                  [0.1750865648952205, 0.11840023306916837, 0.24215989137836502]]

    for i in plot_grid:
        temp_rad, = ax.plot(np.log10(rc/AU), np.log10(temp2d[:,i]),'-',color=color_list[plot_grid.index(i)],\
                            linewidth=2, markersize=3,label=label_grid[plot_grid.index(i)])

    # plot the theoretical prediction for black body dust without considering the extinction
#.........这里部分代码省略.........

import matplotlib.pyplot as plt

from hyperion.model import ModelOutput
from hyperion.util.constants import pc

mo = ModelOutput('class1_example.rtout')
sed = mo.get_sed(aperture=-1, distance=140. * pc)

fig = plt.figure(figsize=(5, 4))
ax = fig.add_subplot(1, 1, 1)
ax.loglog(sed.wav, sed.val.transpose(), color='black')
ax.set_xlim(0.03, 2000.)
ax.set_ylim(2.e-15, 1e-8)
ax.set_xlabel(r'$\lambda$ [$\mu$m]')
ax.set_ylabel(r'$\lambda F_\lambda$ [ergs/cm$^2/s$]')
fig.savefig('class1_example_sed.png', bbox_inches='tight')

name = ['model']

angles=np.arccos(np.linspace(0,1.,20))*180./np.pi
inclinations=angles[::-1]

d = SphericalDust()
d.read('d03_5.5_3.0_A.hdf5')
chi = d.optical_properties.chi
chi = chi[::-1]
wav = d.optical_properties.wav
wav = wav[::-1]
Chi = interp1d(wav,chi,kind='linear')

sorted_grid = pickle.load(open(folder[0]+name[0]+"_"+target+".grid.dat",'r'))
best_model_fname = folder[0]+sorted_grid['name'][0]+'.rtout'
best_model = ModelOutput(fname)
inc = int(np.argwhere(inclinations==sorted_grid['inc'][0]))
sed = best_model.get_sed(aperture=-1, inclination=inc, distance=dist,units='Jy')
N = len(sed.wav)
vec = np.zeros(N,len(target_list)+1)
vec[:,0] = sed.wav

for i in range(len(target_list)):
	target = target_list[i]
	sorted_grid = pickle.load(open(folder[0]+name[0]+"_"+target+".grid.dat",'r'))
	best_model_fname = folder[0]+sorted_grid['name'][0]+'.rtout'
	best_model = ModelOutput(fname)
	extinction = sorted_grid['ext'][0]
				
	# get inclination
	inc = int(np.argwhere(inclinations==sorted_grid['inc'][0]))

import pyfits

from hyperion.model import ModelOutput
from hyperion.util.constants import pc

# Open the model - we specify the name without the .rtout extension
m = ModelOutput('tutorial_model.rtout')

# Extract the image for the first inclination, and scale to 300pc. We
# have to specify group=1 as there is no image in group 0
wav, nufnu = m.get_image(group=1, inclination=0, distance=300 * pc)

# The image extracted above is a 3D array. We can write it out to FITS.
# We need to swap some of the directions around so as to be able to use
# the ds9 slider to change the wavelength of the image.
pyfits.writeto('image_cube.fits', nufnu.swapaxes(0, 2).swapaxes(1, 2), \
               clobber=True)

# We can also just output one of the wavelengths
pyfits.writeto('image_slice.fits', nufnu[:, :, 0], clobber=True)

import matplotlib.pyplot as plt
from hyperion.model import ModelOutput
from hyperion.util.constants import pc

mo = ModelOutput('pure_scattering.rtout')

wav, fnu = mo.get_image(inclination=0, units='MJy/sr', distance=300. * pc)
wav, pol = mo.get_image(inclination=0, stokes='linpol')

fig = plt.figure(figsize=(8, 8))

# Make total intensity sub-plot

ax = fig.add_axes([0.1, 0.3, 0.4, 0.4])
ax.imshow(fnu[:, :, 0], extent=[-13, 13, -13, 13],
          interpolation='none', cmap=plt.cm.gist_heat,
          origin='lower', vmin=0., vmax=4e9)
ax.set_xlim(-13., 13.)
ax.set_ylim(-13., 13.)
ax.set_xlabel("x (solar radii)")
ax.set_ylabel("y (solar radii)")
ax.set_title("Surface brightness")

# Make linear polarization sub-plot

ax = fig.add_axes([0.51, 0.3, 0.4, 0.4])
im = ax.imshow(pol[:, :, 0] * 100., extent=[-13, 13, -13, 13],
               interpolation='none', cmap=plt.cm.gist_heat,
               origin='lower', vmin=0., vmax=100.)
ax.set_xlim(-13., 13.)
ax.set_ylim(-13., 13.)

wav = np.loadtxt('kmh94_3.1_full.wav')
Chi = interp1d(wav,chi,kind='linear')

# now load up the grid
name = ['IRAS20050']
folder = ['/cardini3/mrizzo/2012SOFIA/SED_Models/hyperion/IRAS20050_new/']
filename = folder[0]+name[0]+"_"+target+".grid.dat"
if os.path.exists(filename):
	grid = pickle.load(open(filename,'r'))
	# for each fit that is desired, search for the right line in the grid
	fitnames = grid.group_by('name')
	for fitkey,fitname in zip(fitnames.groups.keys,fitnames.groups):
		name = fitkey['name']
		#print name
		fname = folder[0]+name+'.rtout'
		mo = ModelOutput(fname)
		sed = mo.get_sed(aperture=-1, inclination='all', distance=dist,units='Jy')
		
		if name in plotlist:
			# sort table according to chi2
			fitname.sort('chi2')
			
			# the first line is then the best fit. let's select the extinction
			extinction = fitname['ext'][0]
			print "extinction = ",extinction
			
			# inclination
			angles=np.arccos(np.linspace(0,1.,20))*180./np.pi
			incs=angles[::-1]
			#incs = [0.,10.,20.,30.,40.,50.,60.,70.,80.,90.]

import numpy as np
from hyperion.model import ModelOutput
from hyperion.util.constants import au, lsun

RES = 256

mo = ModelOutput('bm2_eff_vor_temperature.rtout')

g = mo.get_quantities()

from scipy.spatial import cKDTree

sites = np.array([g.x, g.y, g.z]).transpose()

tree = cKDTree(sites)

ymin, ymax = 0 * au, 60 * au
zmin, zmax = 0 * au, 60 * au

y = np.linspace(ymin, ymax, RES)
z = np.linspace(zmin, zmax, RES)

Y, Z = np.meshgrid(y, z)
YR = Y.ravel()
ZR = Z.ravel()

for x_cut in [10 * au, 26.666667 * au]:

    XR = np.ones(YR.shape) * x_cut

    map_sites = np.array([XR, YR, ZR]).transpose()

def hyperion_sedcom(modellist, outdir, plotname, obs_data=None, labellist=None, lbol=False, legend=True, mag=1.5,\
					obs_preset='sh', dstar=1, aper=[3.6, 4.5, 5.8, 8.0, 10, 20, 24, 70, 160, 250, 350, 500, 850]):
	"""
	obs_data: dictionary which obs_data['spec'] is spectrum and obs_data['phot'] is photometry
			  obs_data['label'] = (wave, Fv, err) in um and Jy by default
	"""

	import numpy as np
	import os
	import matplotlib.pyplot as plt
	import astropy.constants as const
	from hyperion.model import ModelOutput
	from scipy.interpolate import interp1d
	from l_bol import l_bol
	import seaborn as sb
	# from seaborn import color_palette
	# from seaborn_color import seaborn_color

	# constant setup
	c = const.c.cgs.value
	pc = const.pc.cgs.value

	if labellist == None:
		if legend == True:
			print 'Model labels are not provided.  Use their filename instead.'
		labellist = []
		for i in range(0, len(modellist)):
			labellist.append(r'$\mathrm{'+os.path.splitext(os.path.basename(modellist[i]))[0]+'}$')

	# cm = seaborn_color('colorblind',len(modellist))
	sb.set(style="white")
	cm = sb.color_palette('husl', len(modellist))

	# create figure object
	fig = plt.figure(figsize=(8*mag,6*mag))
	ax = fig.add_subplot(111)
	# sb.set_style('ticks')

	print 'plotting with aperture at ', aper, 'um'

	# if the obs_data is provided than plot the observation first.  In this way, models won't be blocked by data
	if obs_data != None:
		if 'spec' in obs_data.keys():
			(wave, fv, err) = obs_data['spec']
			vfv = c/(wave*1e-4)*fv*1e-23
			l_bol_obs = l_bol(wave, fv, dstar)
			if legend == True:
				ax.text(0.75,0.9,r'$\mathrm{L_{bol}= %5.2f L_{\odot}}$' % l_bol_obs,fontsize=mag*16,transform=ax.transAxes)

			# general plotting scheme
			if obs_preset == None:
				spec, = ax.plot(np.log10(wave),np.log10(vfv),'-',color='k',linewidth=1.5*mag, label=r'$\mathrm{observations}$')
			# plot spitzer, Herschel pacs and spire in different colors
			elif obs_preset == 'sh':
				# spitzer
				spitz, = ax.plot(np.log10(wave[wave < 50]),np.log10(vfv[wave < 50]),'-',color='b',linewidth=1*mag,\
									label=r'$\mathrm{\it Spitzer}$')
				# herschel
				pacs, = ax.plot(np.log10(wave[(wave < 190.31) & (wave > 50)]),np.log10(vfv[(wave < 190.31) & (wave > 50)]),'-',\
									color='Green',linewidth=1*mag, label=r'$\mathrm{{\it Herschel}-PACS}$')
				spire, = ax.plot(np.log10(wave[wave >= 190.31]),np.log10(vfv[wave >= 190.31]),'-',color='k',linewidth=1*mag,\
									label=r'$\mathrm{{\it Herschel}-SPIRE}$')
				spec = [spitz, pacs, spire]

		if 'phot' in obs_data.keys():
			(wave_p, fv_p, err_p) = obs_data['phot']
			vfv_p = c/(wave_p*1e-4)*fv_p*1e-23
			vfv_p_err = c/(wave_p*1e-4)*err_p*1e-23
			phot, = ax.plot(np.log10(wave_p),np.log10(vfv_p),'s',mfc='DimGray',mec='k',markersize=8)
			ax.errorbar(np.log10(wave_p),np.log10(vfv_p),yerr=[np.log10(vfv_p)-np.log10(vfv_p-vfv_p_err), np.log10(vfv_p+vfv_p_err)-np.log10(vfv_p)],\
						fmt='s',mfc='DimGray',mec='k',markersize=8)

	modplot = dict()
	for imod in range(0, len(modellist)):
		m = ModelOutput(modellist[imod])
		# if not specified, distance of the star will be taken as 1 pc. 
		if aper == None:
			sed_dum = m.get_sed(group=0, inclination=0, aperture=-1, distance=dstar * pc)
			modplot['mod'+str(imod+1)], = ax_sed.plot(np.log10(sed_dum.wav), np.log10(sed_dum.val), '-', color='GoldenRod', linewidth=1.5*mag)
		else:
			vfv_aper = np.empty_like(aper)
			for i in range(0, len(aper)):
				sed_dum = m.get_sed(group=i+1, inclination=0, aperture=-1, distance=dstar * pc)
				f = interp1d(sed_dum.wav, sed_dum.val)
				vfv_aper[i] = f(aper[i])
			modplot['mod'+str(imod+1)], = ax.plot(np.log10(aper),np.log10(vfv_aper),'o',mfc='None',mec=cm[imod],markersize=12,\
													markeredgewidth=3, label=labellist[imod], linestyle='-',color=cm[imod],linewidth=1.5*mag)

	# plot fine tune
	ax.set_xlabel(r'$\mathrm{log~\lambda~({\mu}m)}$',fontsize=mag*20)
	ax.set_ylabel(r'$\mathrm{log~\nu S_{\nu}~(erg/cm^{2}/s)}$',fontsize=mag*20)
	[ax.spines[axis].set_linewidth(1.5*mag) for axis in ['top','bottom','left','right']]
	ax.minorticks_on()
	ax.tick_params('both',labelsize=mag*18,width=1.5*mag,which='major',pad=15,length=5*mag)
	ax.tick_params('both',labelsize=mag*18,width=1.5*mag,which='minor',pad=15,length=2.5*mag)

	if obs_preset == 'sh':
		ax.set_ylim([-14,-7])
		ax.set_xlim([0,3])

#.........这里部分代码省略.........

import os

import numpy as np
import matplotlib.pyplot as plt

from hyperion.model import ModelOutput
from hyperion.util.constants import pc

# Create output directory if it does not already exist
if not os.path.exists('frames'):
    os.mkdir('frames')

# Open model
m = ModelOutput('flyaround_cube.rtout')

# Read image from model
image = m.get_image(distance=300 * pc, units='MJy/sr')

# image.val is now an array with four dimensions (n_view, n_y, n_x, n_wav)

for iview in range(image.val.shape[0]):

    # Open figure and create axes
    fig = plt.figure(figsize=(3, 3))
    ax = fig.add_subplot(1, 1, 1)

    # This is the command to show the image. The parameters vmin and vmax are
    # the min and max levels for the grayscale (remove for default values).
    # The colormap is set here to be a heat map. Other possible heat maps
    # include plt.cm.gray (grayscale), plt.cm.gist_yarg (inverted grayscale),
    # plt.cm.jet (default, colorful). The np.sqrt() is used to plot the

import matplotlib.pyplot as plt

from hyperion.model import ModelOutput
from hyperion.util.constants import pc

m = ModelOutput('class2_sed.rtout')

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

# Extract all SEDs
sed = m.get_sed(inclination='all', aperture=-1, distance=300 * pc)

# Plot SED for each inclination
for i in range(sed.val.shape[0]):
    ax.loglog(sed.wav, sed.val[i, :], color='black')

ax.set_xlabel(r'$\lambda$ [$\mu$m]')
ax.set_ylabel(r'$\lambda F_\lambda$ [ergs/s/cm$^2$]')
ax.set_xlim(0.1, 2000.)
ax.set_ylim(2.e-16, 2.e-9)
fig.savefig('class2_sed_plot_incl.png')

matplotlib.use('Agg')

import numpy as np
from astropy.io import fits

from hyperion.model import ModelOutput
from hyperion.util.constants import kpc

import matplotlib.pyplot as plt

if not os.path.exists('seds'):
    os.mkdir('seds')

for model_path in glob.glob(os.path.join('models', '*_seds.rtout')):

    m = ModelOutput(model_path)

    model_name = os.path.basename(model_path).replace('_seds.rtout', '')

    for iincl, theta in enumerate([0, 30, 60, 90, 120, 150, 180]):

        sed_total = m.get_sed(inclination=iincl, units='Jy', distance=10. * kpc, aperture=-1, group=0, component='total')
        sed_semit = m.get_sed(inclination=iincl, units='Jy', distance=10. * kpc, aperture=-1, group=0, component='source_emit')
        sed_sscat = m.get_sed(inclination=iincl, units='Jy', distance=10. * kpc, aperture=-1, group=0, component='source_scat')
        sed_demit = m.get_sed(inclination=iincl, units='Jy', distance=10. * kpc, aperture=-1, group=0, component='dust_emit')
        sed_dscat = m.get_sed(inclination=iincl, units='Jy', distance=10. * kpc, aperture=-1, group=0, component='dust_scat')

        sed_trans = m.get_sed(inclination=iincl, units='Jy', distance=10. * kpc, aperture=-1, group=1, component='source_emit')

        output_file = 'seds/{name}_i{theta:03d}a000.sed'.format(name=model_name, theta=theta)