本文整理汇总了Python中hyperion.model.ModelOutput.get_quantities方法的典型用法代码示例。如果您正苦于以下问题：Python ModelOutput.get_quantities方法的具体用法？Python ModelOutput.get_quantities怎么用？Python ModelOutput.get_quantities使用的例子？那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类hyperion.model.ModelOutput的用法示例。
在下文中一共展示了ModelOutput.get_quantities方法的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: getRadialDensity
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
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,:]示例2: plot_results
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
def plot_results(cli):
file = filename(cli, "plot")
file += ".rtout"
#
# Read in the model:
#
model = ModelOutput(file)
if(cli.mode == "images"):
#
# Extract the quantities
#
g = model.get_quantities()
#
# Get the wall positions:
#
ww = g.w_wall / pc
zw = g.z_wall / pc
pw = g.p_wall
grid_Nw = len(ww) - 1
grid_Nz = len(zw) - 1
grid_Np = len(pw) - 1
#
# Graphics:
#
fig = plt.figure()
los = [0 for i in range(3)]
los[0] = 'x'
los[1] = 'y'
los[2] = 'z'
#Imaxp = [0 for i in range(4)]
##Imaxp[0] = 1e-4
#Imaxp[1] = 1e-5
#Imaxp[2] = 1e-7
#Imaxp[3] = 1e-8
for k in range(0, 3):
if(cli.verbose):
print("Group: ", k)
image = model.get_image(distance=1*pc, units='MJy/sr', inclination=0, component='total', group=k)
source_emit = model.get_image(distance=1*pc, units='MJy/sr', inclination=0, component='source_emit', group=k)
dust_emit = model.get_image(distance=1*pc, units='MJy/sr', inclination=0, component='dust_emit' , group=k)
source_scat = model.get_image(distance=1*pc, units='MJy/sr', inclination=0, component='source_scat', group=k)
dust_scat = model.get_image(distance=1*pc, units='MJy/sr', inclination=0, component='dust_scat' , group=k)
if(cli.verbose):
print(" Data cube: ", image.val.shape)
print(" Wavelengths =", image.wav)
print(" Uncertainties =", image.unc)
image_Nx=image.val.shape[0]
image_Ny=image.val.shape[1]
Nwavelength=image.val.shape[2]
if(cli.verbose):
print(" Image Nx =", image_Nx)
print(" Image Ny =", image_Ny)
print(" Nwavelength =", Nwavelength)
for i in range(0, Nwavelength):
if(cli.verbose):
print(" Image #", i,":")
print(" Wavelength =", image.wav[i])
Imin = np.min(image.val[:, :, i])
Imax = np.max(image.val[:, :, i])
# TODO: compute the mean value as well and use this for specifying the maximum value/color?!
if(cli.verbose):
print(" Intensity min =", Imin)
print(" Intensity max =", Imax)
#Imax=Imaxp[i]
#ax = fig.add_subplot(2, 1, 2)
ax = fig.add_subplot(1, 1, 1)
if(image.wav[i] < 10.0):
ax.imshow(source_scat.val[:, :, i] + dust_scat.val[:, :, i], vmin=Imin, vmax=Imax/10, cmap=plt.cm.gist_heat, origin='lower')
else:
ax.imshow(image.val[:, :, i], vmin=Imin, vmax=Imax/10, cmap=plt.cm.gist_heat, origin='lower')
ax.set_xticks([0,100,200,300], minor=False)
ax.set_yticks([0,100,200,300], minor=False)
ax.set_xlabel('x (pixel)')
ax.set_ylabel('y (pixel)')
ax.set_title(str(image.wav[i]) + ' microns' + '\n' + los[k] + '-direction', y=0.88, x=0.5, color='white')
#ax = fig.add_subplot(2, 1, 1)
#ax.imshow([np.logspace(np.log10(Imin+1e-10),np.log10(Imax/10),100),np.logspace(np.log10(Imin+1e-10),np.log10(Imax/10),100)], vmin=Imin, vmax=Imax/10, cmap=plt.cm.gist_heat)
#ax.set_xticks(np.logspace(np.log10(Imin+1e-10),np.log10(Imax/10),1), minor=False)
##ax.set_xticks(np.linspace(np.log10(Imin+1e-10),np.log10(Imax/10),10), minor=False)
#.........这里部分代码省略.........
示例3: temp_hyperion
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
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
#.........这里部分代码省略.........
示例4: range
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
ax_top = grid[i].twiny()
ax_top.set_xticks(r_ticks)
ax_top.set_xticklabels([])
grid[i].tick_params('both',labelsize=14)
for i in range(4,8):
# get the H-band simulated image
m = ModelOutput(filename[i-4])
image = m.get_image(group=0, inclination=0, distance=178 * pc, units='MJy/sr')
# Find the closest wavelength
iwav = np.argmin(np.abs(wave - image.wav))
# Calculate the image width in arcseconds given the distance used above
# get the max radius
rmax = max(m.get_quantities().r_wall)
w = np.degrees(rmax / image.distance) * 3600.
# Image in the unit of MJy/sr
# Change it into erg/s/cm2/Hz/sr
factor = 1e-23*1e6
# avoid zero in log
# flip the image, because the setup of inclination is upside down
val = image.val[::-1, :, iwav] * factor + 1e-30
if size != 'full':
pix_e2c = (w-size/2.)/w * len(val[:,0])/2
val = val[pix_e2c:-pix_e2c, pix_e2c:-pix_e2c]
w = size/2.
if convolve:示例5: extract_hyperion
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
#.........这里部分代码省略.........
dum.use_sources(filename)
L_cen = dum.sources[0].luminosity/lsun
# legend
lg_data = ax_sed.legend([irs, photometry, aper, aper_obs],
[r'$\rm{observation}$',
r'$\rm{photometry}$',r'$\rm{F_{aper,sim}}$',r'$\rm{F_{aper,obs}}$'],
loc='upper left',fontsize=14*mag,numpoints=1,framealpha=0.3)
if clean == False:
lg_sim = ax_sed.legend([sim],[r'$\rm{L_{bol,sim}=%5.2f\,L_{\odot},\,L_{center}=%5.2f\,L_{\odot}}$' % (l_bol_sim, L_cen)], \
loc='lower right',fontsize=mag*16)
plt.gca().add_artist(lg_data)
# plot setting
ax_sed.set_xlabel(r'$\rm{log\,\lambda\,[{\mu}m]}$',fontsize=mag*20)
ax_sed.set_ylabel(r'$\rm{log\,\nu S_{\nu}\,[erg\,s^{-1}\,cm^{-2}]}$',fontsize=mag*20)
[ax_sed.spines[axis].set_linewidth(1.5*mag) for axis in ['top','bottom','left','right']]
ax_sed.minorticks_on()
ax_sed.tick_params('both',labelsize=mag*18,width=1.5*mag,which='major',pad=15,length=5*mag)
ax_sed.tick_params('both',labelsize=mag*18,width=1.5*mag,which='minor',pad=15,length=2.5*mag)
# fix the tick label font
ticks_font = mpl.font_manager.FontProperties(family='STIXGeneral',size=mag*18)
for label in ax_sed.get_xticklabels():
label.set_fontproperties(ticks_font)
for label in ax_sed.get_yticklabels():
label.set_fontproperties(ticks_font)
# Write out the plot
fig.savefig(outdir+print_name+'_sed.pdf',format='pdf',dpi=300,bbox_inches='tight')
fig.clf()
# option for suppress image plotting (for speed)
if image:
# Package for matching the colorbar
from mpl_toolkits.axes_grid1 import make_axes_locatable, ImageGrid
# Users may change the unit: mJy, Jy, MJy/sr, ergs/cm^2/s, ergs/cm^2/s/Hz
# !!!
image = m.get_image(group=len(aper_reduced)+1, inclination=0,
distance=dstar*pc, units='MJy/sr')
# Open figure and create axes
fig = plt.figure(figsize=(12,12))
grid = ImageGrid(fig, 111,nrows_ncols=(3,3),direction='row',
add_all=True,label_mode='1',share_all=True,
cbar_location='right',cbar_mode='single',
cbar_size='3%',cbar_pad=0)
for i, wav in enumerate([3.6, 8.0, 9.7, 24, 40, 100, 250, 500, 1000]):
ax = grid[i]
# Find the closest wavelength
iwav = np.argmin(np.abs(wav - image.wav))
# Calculate the image width in arcseconds given the distance used above
# get the max radius
rmax = max(m.get_quantities().r_wall)
w = np.degrees(rmax / image.distance) * 3600.
# Image in the unit of MJy/sr
# Change it into erg/s/cm2/Hz/sr
factor = 1e-23*1e6
# avoid zero in log
# flip the image, because the setup of inclination is upside down
val = image.val[::-1, :, iwav] * factor + 1e-30
# This is the command to show the image. The parameters vmin and vmax are
# the min and max levels for the colorscale (remove for default values).
cmap = plt.cm.CMRmap
im = ax.imshow(np.log10(val), vmin= -22, vmax= -12,
cmap=cmap, origin='lower', extent=[-w, w, -w, w], aspect=1)
ax.set_xlabel(r'$\rm{RA\,Offset\,[arcsec]}$', fontsize=14)
ax.set_ylabel(r'$\rm{Dec\,Offset\,[arcsec]}$', fontsize=14)
# fix the tick label font
ticks_font = mpl.font_manager.FontProperties(family='STIXGeneral',size=14)
for label in ax.get_xticklabels():
label.set_fontproperties(ticks_font)
for label in ax.get_yticklabels():
label.set_fontproperties(ticks_font)
# Colorbar setting
cb = ax.cax.colorbar(im)
cb.solids.set_edgecolor('face')
cb.ax.minorticks_on()
cb.ax.set_ylabel(r'$\rm{log(I_{\nu})\,[erg\,s^{-1}\,cm^{-2}\,Hz^{-1}\,sr^{-1}]}$',fontsize=18)
cb_obj = plt.getp(cb.ax.axes, 'yticklabels')
plt.setp(cb_obj,fontsize=18)
ticks_font = mpl.font_manager.FontProperties(family='STIXGeneral',size=18)
for label in cb.ax.get_yticklabels():
label.set_fontproperties(ticks_font)
ax.tick_params(axis='both', which='major', labelsize=16)
ax.text(0.7,0.88,str(wav) + r'$\rm{\,\mu m}$',fontsize=16,color='white', transform=ax.transAxes)
fig.savefig(outdir+print_name+'_image_gridplot.pdf', format='pdf', dpi=300, bbox_inches='tight')
fig.clf()示例6: extract_hyperion
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
#.........这里部分代码省略.........
ax_sed.minorticks_on()
ax_sed.tick_params('both',labelsize=mag*18,width=1.5*mag,which='major',pad=15,length=5*mag)
ax_sed.tick_params('both',labelsize=mag*18,width=1.5*mag,which='minor',pad=15,length=2.5*mag)
ax_sed.set_ylim([-13,-7.5])
ax_sed.set_xlim([0,3])
# lg_data = ax_sed.legend([sim, aper], [r'$\rm{w/o~aperture}$', r'$\rm{w/~aperture}$'], \
# loc='upper left', fontsize=14*mag, framealpha=0.3, numpoints=1)
lg_data = ax_sed.legend([irs, photometry, aper, aper_obs],\
[r'$\rm{observation}$',\
r'$\rm{photometry}$',r'$\rm{F_{aper,sim}}$',r'$\rm{F_{aper,obs}}$'],\
loc='upper left',fontsize=14*mag,numpoints=1,framealpha=0.3)
# plt.gca().add_artist(lg_sim)
# Write out the plot
fig.savefig(outdir+print_name+'_sed.pdf',format='pdf',dpi=300,bbox_inches='tight')
fig.clf()
# Package for matching the colorbar
from mpl_toolkits.axes_grid1 import make_axes_locatable
# 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.
image = m.get_image(group=len(wl_aper)+1, inclination=0, distance=dstar * pc, units='MJy/sr')
# image = m.get_image(group=14, inclination=0, distance=dstar * pc, units='MJy/sr')
# Open figure and create axes
# fig = plt.figure(figsize=(8, 8))
fig, axarr = plt.subplots(3, 3, sharex='col', sharey='row',figsize=(13.5,12))
# Pre-set maximum for colorscales
VMAX = {}
# VMAX[3.6] = 10.
# VMAX[24] = 100.
# VMAX[160] = 2000.
# VMAX[500] = 2000.
VMAX[100] = 10.
VMAX[250] = 100.
VMAX[500] = 2000.
VMAX[1000] = 2000.
# We will now show four sub-plots, each one for a different wavelength
# for i, wav in enumerate([3.6, 24, 160, 500]):
# for i, wav in enumerate([100, 250, 500, 1000]):
# for i, wav in enumerate([4.5, 9.7, 24, 40, 70, 100, 250, 500, 1000]):
for i, wav in enumerate([3.6, 8.0, 9.7, 24, 40, 100, 250, 500, 1000]):
# ax = fig.add_subplot(3, 3, i + 1)
ax = axarr[i/3, i%3]
# Find the closest wavelength
iwav = np.argmin(np.abs(wav - image.wav))
# Calculate the image width in arcseconds given the distance used above
rmax = max(m.get_quantities().r_wall)
w = np.degrees(rmax / image.distance) * 3600.
# w = np.degrees((1.5 * pc) / image.distance) * 60.
# Image in the unit of MJy/sr
# Change it into erg/s/cm2/Hz/sr
factor = 1e-23*1e6
# avoid zero in log
val = image.val[:, :, iwav] * factor + 1e-30
# This is the command to show the image. The parameters vmin and vmax are
# the min and max levels for the colorscale (remove for default values).
im = ax.imshow(np.log10(val), vmin= -22, vmax= -12,
cmap=plt.cm.jet, origin='lower', extent=[-w, w, -w, w], aspect=1)
# Colorbar setting
# create an axes on the right side of ax. The width of cax will be 5%
# of ax and the padding between cax and ax will be fixed at 0.05 inch.
if (i+1) % 3 == 0:
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cb = fig.colorbar(im, cax=cax)
cb.solids.set_edgecolor("face")
cb.ax.minorticks_on()
cb.ax.set_ylabel(r'$\rm{log(I_{\nu})\,[erg\,s^{-2}\,cm^{-2}\,Hz^{-1}\,sr^{-1}]}$',fontsize=12)
cb_obj = plt.getp(cb.ax.axes, 'yticklabels')
plt.setp(cb_obj,fontsize=12)
if (i+1) == 7:
# Finalize the plot
ax.set_xlabel('RA Offset (arcsec)', fontsize=14)
ax.set_ylabel('Dec Offset (arcsec)', fontsize=14)
ax.tick_params(axis='both', which='major', labelsize=16)
ax.set_adjustable('box-forced')
ax.text(0.7,0.88,str(wav) + r'$\rm{\,\mu m}$',fontsize=18,color='white',weight='bold',transform=ax.transAxes)
fig.subplots_adjust(hspace=0,wspace=-0.2)
# Adjust the spaces between the subplots
# plt.tight_layout()
fig.savefig(outdir+print_name+'_cube_plot.png', format='png', dpi=300, bbox_inches='tight')
fig.clf()示例7: ModelOutput
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
import numpy as np
from hyperion.model import ModelOutput
import matplotlib.pyplot as plt
from yt.mods import write_bitmap, ColorTransferFunction
plt.rcParams['font.family'] = 'Arial'
# Read in model from Hyperion
m = ModelOutput('pla704850_lev7_129.rtout')
grid = m.get_quantities()
# Convert quantities to yt
pf = grid.to_yt()
# Instantiate the ColorTransferfunction.
tmin, tmax = 1.3, 2.3
tf_temp = ColorTransferFunction((tmin, tmax))
dmin, dmax = -20, -16
tf_dens = ColorTransferFunction((dmin, dmax))
# Set up the camera parameters: center, looking direction, width, resolution
c = (pf.domain_right_edge + pf.domain_left_edge) / 2.0
L = np.array([1.0, 1.0, 1.0])
W = 0.7 / pf["unitary"]
N = 512
# Create camera objects
cam_temp = pf.h.camera(c, L, W, N, tf_temp,示例8: YSOModelSim
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
#.........这里部分代码省略.........
ax.set_ylabel(r'$\lambda F_\lambda$ [ergs/cm$^2/s$]')
self.plotData(ax,sourcename)
ax.set_xscale('log')
ax.set_yscale('log')
#ax.set_ylabel(r'$F_{Jy}$ [Jy]')
#plt.legend(loc=4)
ax=fig.add_subplot(2,3,2)
sed = self.mo.get_sed(aperture=-1, inclination=self.inc, distance=self.dist)
ext=np.exp(-tau)
ax.loglog(sed.wav, sed.val.transpose()*ext.T, lw=3,color='black',label='source_total')
ax.set_xlim(sed.wav.min(), 1300)
ax.set_ylim(1e-13, 1e-7) ### for lamFlam
sed = self.mo.get_sed(aperture=-1, inclination=self.inc, distance=self.dist,component='source_emit')
ax.loglog(sed.wav, sed.val.transpose()*ext.T, color='blue',label='source_emit')
sed = self.mo.get_sed(aperture=-1, inclination=self.inc, distance=self.dist,component='source_scat')
ax.loglog(sed.wav, sed.val.transpose()*ext.T, color='teal',label='source_scat')
sed = self.mo.get_sed(aperture=-1, inclination=self.inc, distance=self.dist,component='dust_emit')
ax.loglog(sed.wav, sed.val.transpose()*ext.T, color='red',label='dust_emit')
sed = self.mo.get_sed(aperture=-1, inclination=self.inc, distance=self.dist,component='dust_scat')
ax.loglog(sed.wav, sed.val.transpose()*ext.T, color='orange',label='dust_scat')
self.plotData(ax,sourcename)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_title('seds_inc=inc')
ax.set_xlabel(r'$\lambda$ [$\mu$m]')
ax.set_ylabel(r'$\lambda F_\lambda$ [ergs/cm$^2/s$]')
#ax.set_ylabel(r'$F_{Jy}$ [Jy]')
leg = ax.legend(loc=4,fontsize='small')
#leg = plt.gca().get_legend()
#plt.setp(leg.get_text(),fontsize='small')
# Extract the quantities
g = self.mo.get_quantities()
# Get the wall positions for r and theta
rw, tw = g.r_wall / au, g.t_wall
# Make a 2-d grid of the wall positions (used by pcolormesh)
R, T = np.meshgrid(rw, tw)
# Calculate the position of the cell walls in cartesian coordinates
X, Z = R * np.sin(T), R * np.cos(T)
# Make a plot in (x, z) space for different zooms
from matplotlib.colors import LogNorm,PowerNorm
# Make a plot in (r, theta) space
ax = fig.add_subplot(2, 3, 3)
if g.shape[-1]==2:
c = ax.pcolormesh(X, Z, g['temperature'][0].array[0, :, :]+g['temperature'][1].array[0, :, :],norm=PowerNorm(gamma=0.5,vmin=1,vmax=500))
else :
c = ax.pcolormesh(X, Z, g['temperature'][0].array[0, :, :],norm=PowerNorm(gamma=0.5,vmin=1,vmax=500))
#ax.set_xscale('log')
#ax.set_yscale('log')
ax.set_xlim(X.min(), X.max()/5.)
ax.set_ylim(Z.min()/10., Z.max()/10.)
ax.set_xlabel('x (au)')
ax.set_ylabel('z (au)')
#ax.set_yticks([np.pi, np.pi * 0.75, np.pi * 0.5, np.pi * 0.25, 0.])
#ax.set_yticklabels([r'$\pi$', r'$3\pi/4$', r'$\pi/2$', r'$\pi/4$', r'$0$'])
cb = fig.colorbar(c)
ax.set_title('Temperature structure')
cb.set_label('Temperature (K)')
#fig.savefig(modelname+'_temperature_spherical_rt.png', bbox_inches='tight')
示例9: zeros
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
m.set_spherical_polar_grid(r, t, p)
dens = zeros((nr-1,nt-1,np-1)) + 1.0e-17
m.add_density_grid(dens, d)
source = m.add_spherical_source()
source.luminosity = lsun
source.radius = rsun
source.temperature = 4000.
m.set_n_photons(initial=1000000, imaging=0)
m.set_convergence(True, percentile=99., absolute=2., relative=1.02)
m.write("test_spherical.rtin")
m.run("test_spherical.rtout", mpi=False)
n = ModelOutput('test_spherical.rtout')
grid = n.get_quantities()
temp = grid.quantities['temperature'][0]
for i in range(9):
plt.imshow(temp[i,:,:],origin="lower",interpolation="nearest", \
vmin=temp.min(),vmax=temp.max())
plt.colorbar()
plt.show()
示例10: run_thermal_hyperion
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
def run_thermal_hyperion(self, nphot=1e6, mrw=False, pda=False, \
niterations=20, percentile=99., absolute=2.0, relative=1.02, \
max_interactions=1e8, mpi=False, nprocesses=None):
d = []
for i in range(len(self.grid.dust)):
d.append(IsotropicDust( \
self.grid.dust[i].nu[::-1].astype(numpy.float64), \
self.grid.dust[i].albedo[::-1].astype(numpy.float64), \
self.grid.dust[i].kext[::-1].astype(numpy.float64)))
m = HypModel()
if (self.grid.coordsystem == "cartesian"):
m.set_cartesian_grid(self.grid.w1*AU, self.grid.w2*AU, \
self.grid.w3*AU)
elif (self.grid.coordsystem == "cylindrical"):
m.set_cylindrical_polar_grid(self.grid.w1*AU, self.grid.w3*AU, \
self.grid.w2)
elif (self.grid.coordsystem == "spherical"):
m.set_spherical_polar_grid(self.grid.w1*AU, self.grid.w2, \
self.grid.w3)
for i in range(len(self.grid.density)):
if (self.grid.coordsystem == "cartesian"):
m.add_density_grid(numpy.transpose(self.grid.density[i], \
axes=(2,1,0)), d[i])
if (self.grid.coordsystem == "cylindrical"):
m.add_density_grid(numpy.transpose(self.grid.density[i], \
axes=(1,2,0)), d[i])
if (self.grid.coordsystem == "spherical"):
m.add_density_grid(numpy.transpose(self.grid.density[i], \
axes=(2,1,0)), d[i])
sources = []
for i in range(len(self.grid.stars)):
sources.append(m.add_spherical_source())
sources[i].luminosity = self.grid.stars[i].luminosity * L_sun
sources[i].radius = self.grid.stars[i].radius * R_sun
sources[i].temperature = self.grid.stars[i].temperature
m.set_mrw(mrw)
m.set_pda(pda)
m.set_max_interactions(max_interactions)
m.set_n_initial_iterations(niterations)
m.set_n_photons(initial=nphot, imaging=0)
m.set_convergence(True, percentile=percentile, absolute=absolute, \
relative=relative)
m.write("temp.rtin")
m.run("temp.rtout", mpi=mpi, n_processes=nprocesses)
n = ModelOutput("temp.rtout")
grid = n.get_quantities()
self.grid.temperature = []
temperature = grid.quantities['temperature']
for i in range(len(temperature)):
if (self.grid.coordsystem == "cartesian"):
self.grid.temperature.append(numpy.transpose(temperature[i], \
axes=(2,1,0)))
if (self.grid.coordsystem == "cylindrical"):
self.grid.temperature.append(numpy.transpose(temperature[i], \
axes=(2,0,1)))
if (self.grid.coordsystem == "spherical"):
self.grid.temperature.append(numpy.transpose(temperature[i], \
axes=(2,1,0)))
os.system("rm temp.rtin temp.rtout")示例11: azimuthal_simulation
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
def azimuthal_simulation(rtout, beam_size, wave, dist=200., group=22):
"""
rtout: the filepath to the output file of Hyperion
beam_size: the beam size used for the width of annulus
dist: the physical distance to the source
group: the group which contains image
"""
import numpy as np
import matplotlib.pyplot as plt
import astropy.constants as const
from hyperion.model import ModelOutput
# constant setup
pc = const.pc.cgs.value
au = const.au.cgs.value
output = {'wave': wave, 'annuli': [], 'flux_annuli': []}
# Read in the Hyperion output file
m = ModelOutput(rtout)
# get image
image = m.get_image(group=5, inclination=0, distance=dist*pc, units='Jy')
# Calculate the image width in arcsec given the distance to the source
rmax = max(m.get_quantities().r_wall)
w = np.degrees(rmax / image.distance) * 3600
# grid of radii of annulus
annuli = np.linspace(beam_size/2., np.floor((w-beam_size/2.)/beam_size)*beam_size, np.floor((w-beam_size/2.)/beam_size)) # plot
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
# iternate through wavelength
if type(wave) == int or type(wave) == float:
wave = [wave]
color_list = plt.cm.viridis(np.linspace(0, 1, len(wave)+1))
for i in range(len(wave)):
wav = wave[i]
# Find the closest wavelength
iwav = np.argmin(np.abs(wav - image.wav))
# avoid zero when log, and flip the image
val = image.val[::-1, :, iwav]
# determine the center of the image
npix = len(val[:,0])
center = np.array([npix/2. + 0.5, npix/2. + 0.5])
scale = 2*rmax/npix
# create index array of the image
x = np.empty_like(val)
for j in range(len(val[0,:])):
x[:,j] = j
flux_annuli = np.empty_like(annuli)
for k in range(len(annuli)):
flux_annuli[k] = np.sum(val[(((x-center[0])**2+(x.T-center[1])**2)**0.5*2*w/npix >= annuli[k]-beam_size/2.) & \
(((x-center[0])**2+(x.T-center[1])**2)**0.5*2*w/npix < annuli[k]+beam_size/2.)])
output['annuli'].append(np.array(annuli))
output['flux_annuli'].append(flux_annuli)
flux_annuli = flux_annuli/np.nanmax(flux_annuli)
ax.plot(np.log10(annuli*dist), np.log10(flux_annuli), 'o-', color=color_list[i], \
markersize=3, mec='None', label=r'$\rm{'+str(wav)+'\,\mu m}$')
ax.axvline(np.log10((w-beam_size/2.)*dist), linestyle='--', color='k')
ax.axvline(np.log10(w*dist), linestyle='-', color='k')
ax.legend(loc='best', fontsize=12, numpoints=1, ncol=2)
ax.set_xlabel(r'$\rm{log(Radius)\,[AU]}$', fontsize=18)
ax.set_ylabel(r'${\rm log(}F/F_{\rm max})$', fontsize=18)
fig.gca().set_ylim(top=0.1)
[ax.spines[axis].set_linewidth(1.5) for axis in ['top','bottom','left','right']]
ax.minorticks_on()
ax.tick_params('both',labelsize=18,width=1.5,which='major',pad=15,length=5)
ax.tick_params('both',labelsize=18,width=1.5,which='minor',pad=15,length=2.5)
fig.savefig('/Users/yaolun/test/annuli_profile.pdf', format='pdf', dpi=300, bbox_inches='tight')
fig.clf()
return output示例12: hyperion_image
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
def hyperion_image(rtout, wave, plotdir, printname, dstar=178., group=0, marker=0,
size='full', convolve=False, unit=None):
# to avoid X server error
import matplotlib as mpl
mpl.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import astropy.constants as const
from hyperion.model import ModelOutput
# Package for matching the colorbar
from mpl_toolkits.axes_grid1 import make_axes_locatable
pc = const.pc.cgs.value
if unit == None:
unit = r'$\rm{log(I_{\nu})\,[erg\,s^{-1}\,cm^{-2}\,Hz^{-1}\,sr^{-1}]}$'
m = ModelOutput(rtout)
# Extract the image.
image = m.get_image(group=group, inclination=0, distance=dstar * pc, units='mJy')
print np.shape(image.val)
# Open figure and create axes
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)
# Find the closest wavelength
iwav = np.argmin(np.abs(wave - image.wav))
# Calculate the image width in arcseconds given the distance used above
# get the max radius
rmax = max(m.get_quantities().r_wall)
w = np.degrees(rmax / image.distance) * 3600.
# Image in the unit of MJy/sr
# Change it into erg/s/cm2/Hz/sr
# factor = 1e-23*1e6
factor = 1
# avoid zero in log
# flip the image, because the setup of inclination is upside down
val = image.val[::-1, :, iwav] * factor + 1e-30
if convolve:
from astropy.convolution import convolve, Gaussian2DKernel
img_res = 2*w/len(val[:,0])
kernel = Gaussian2DKernel(0.27/2.354/img_res)
val = convolve(val, kernel)
if size != 'full':
pix_e2c = (w-size/2.)/w * len(val[:,0])/2
val = val[pix_e2c:-pix_e2c, pix_e2c:-pix_e2c]
w = size/2.
# This is the command to show the image. The parameters vmin and vmax are
# the min and max levels for the colorscale (remove for default values).
# cmap = sns.cubehelix_palette(start=0.1, rot=-0.7, gamma=0.2, as_cmap=True)
cmap = plt.cm.CMRmap
# im = ax.imshow(np.log10(val), vmin= -20, vmax= -15,
# cmap=cmap, origin='lower', extent=[-w, w, -w, w], aspect=1)
im = ax.imshow(val,
cmap=cmap, origin='lower', extent=[-w, w, -w, w], aspect=1)
print val.max()
# plot the marker for center position by default or user input offset
ax.plot([0],[-marker], '+', color='ForestGreen', markersize=10, mew=2)
ax.set_xlim([-w,w])
ax.set_ylim([-w,w])
# ax.plot([0],[-10], '+', color='m', markersize=10, mew=2)
# fix the tick label font
ticks_font = mpl.font_manager.FontProperties(family='STIXGeneral',size=14)
for label in ax.get_xticklabels():
label.set_fontproperties(ticks_font)
for label in ax.get_yticklabels():
label.set_fontproperties(ticks_font)
# Colorbar setting
# create an axes on the right side of ax. The width of cax will be 5%
# of ax and the padding between cax and ax will be fixed at 0.05 inch.
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cb = fig.colorbar(im, cax=cax)
cb.solids.set_edgecolor("face")
cb.ax.minorticks_on()
cb.ax.set_ylabel(unit,fontsize=18)
cb_obj = plt.getp(cb.ax.axes, 'yticklabels')
plt.setp(cb_obj,fontsize=14)
# fix the tick label font
ticks_font = mpl.font_manager.FontProperties(family='STIXGeneral',size=14)
for label in cb.ax.get_yticklabels():
label.set_fontproperties(ticks_font)
ax.set_xlabel(r'$\rm{RA\,Offset\,(arcsec)}$', fontsize=18)
ax.set_ylabel(r'$\rm{Dec\,Offset\,(arcsec)}$', fontsize=18)
ax.tick_params(axis='both', which='major', labelsize=18)
ax.text(0.7,0.88,str(wave) + r'$\rm{\,\mu m}$',fontsize=20,color='white', transform=ax.transAxes)
#.........这里部分代码省略.........
示例13: azimuthal_avg_radial_intensity
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
#.........这里部分代码省略.........
linestyle_dum = '-'
rrange_dum = rrange
annulus_width_dum = annulus_width
else:
label_dum = o['label']
color_dum = o['plot_color']
linestyle_dum = o['plot_linestyle']
rrange_dum = o['rrange']
annulus_width_dum = o['annulus_width']
r_dum, I_dum, I_err_dum = ExtractIntensityObs(rrange_dum, annulus_width_dum, o)
# determine the label
I_obs.append({'imgpath':o['imgpath'], 'r':r_dum, 'I':I_dum, 'I_err':I_err_dum, 'label': label_dum,
'plot_color':color_dum, 'plot_linestyle':linestyle_dum})
# The first image should be the one to be compared primarily, and written out
I = I_obs[0]['I']
I_err = I_obs[0]['I_err']
imgpath = I_obs[0]['imgpath']
#
# read in from RTout
rtout = ModelOutput(rtout)
im = rtout.get_image(group=group, inclination=0, distance=dstar*pc, units='Jy', uncertainties=True)
factor = 1
# Find the closest wavelength
iwav = np.argmin(np.abs(wave - im.wav))
# avoid zero when log, and flip the image
val = im.val[::-1, :, iwav]
unc = im.unc[::-1, :, iwav]
w = np.degrees(max(rtout.get_quantities().r_wall) / im.distance) * 3600
npix = len(val[:,0])
pix2arcsec = 2*w/npix
I_sim = np.empty_like(r[:-1])
I_sim_hi = np.empty_like(r[:-1])
I_sim_low = np.empty_like(r[:-1])
I_sim_err = np.empty_like(r[:-1])
# for calculating the uncertainty from the variation within each annulus
# construct the x- and y-matrix
grid_x, grid_y = np.meshgrid(np.linspace(0,npix-1,npix),
np.linspace(0,npix-1,npix))
dist_x = abs(grid_x - ((npix-1)/2.))
dist_y = abs(grid_y - ((npix-1)/2.))
grid_dist = (dist_x**2+dist_y**2)**0.5
# iteration
for ir in range(len(r)-1):
aperture = CircularAnnulus((npix/2.+0.5, npix/2.+0.5),
r_in=r[ir]/pix2arcsec, r_out=r[ir+1]/pix2arcsec)
phot = ap(val, aperture, error=unc)
I_sim[ir] = phot['aperture_sum'].data / aperture.area()
# uncertainty
im_dum = np.where((grid_dist < r[ir+1]/pix2arcsec) & (grid_dist >= r[ir]/pix2arcsec), val, np.nan)
# I_sim_err[ir] = phot['aperture_sum_err'].data / aperture.area()
# I_sim_err[ir] = (np.nanstd(im_dum)**2+phot['aperture_sum_err'].data**2)**0.5 * factor / aperture.area()
offset = -1
aperture = CircularAnnulus((npix/2.+0.5, npix/2.+0.5),示例14: ModelOutput
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
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()示例15: hyperion_image
# 需要导入模块: from hyperion.model import ModelOutput [as 别名]
# 或者: from hyperion.model.ModelOutput import get_quantities [as 别名]
def hyperion_image(rtout, wave, plotdir, printname, dstar=200., group=0, marker=0,
size='full', convolve=False, unit=None, scalebar=None):
# to avoid X server error
import matplotlib as mpl
mpl.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import astropy.constants as const
from hyperion.model import ModelOutput
# Package for matching the colorbar
from mpl_toolkits.axes_grid1 import make_axes_locatable
pc = const.pc.cgs.value
if unit == None:
unit = 'erg\,s^{-1}\,cm^{-2}\,Hz^{-1}\,sr^{-1}'
m = ModelOutput(rtout)
# Extract the image.
image = m.get_image(group=group, inclination=0, distance=dstar * pc, units='MJy/sr')
# print np.shape(image.val)
# Open figure and create axes
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)
# Find the closest wavelength
iwav = np.argmin(np.abs(wave - image.wav))
# Calculate the image width in arcseconds given the distance used above
# get the max radius
rmax = max(m.get_quantities().r_wall)
w = np.degrees(rmax / image.distance) * 3600.
# Image in the unit of MJy/sr
# Change it into erg/s/cm2/Hz/sr
# factor = 1e-23*1e6
factor = 1
# avoid zero in log
# flip the image, because the setup of inclination is upside down
val = image.val[::-1, :, iwav] * factor + 1e-30
if convolve:
from astropy.convolution import convolve, Gaussian2DKernel
img_res = 2*w/len(val[:,0])
kernel = Gaussian2DKernel(0.27/2.354/img_res)
val = convolve(val, kernel)
if size != 'full':
pix_e2c = (w-size/2.)/w * len(val[:,0])/2
val = val[pix_e2c:-pix_e2c, pix_e2c:-pix_e2c]
w = size/2.
# This is the command to show the image. The parameters vmin and vmax are
# the min and max levels for the colorscale (remove for default values).
# cmap = sns.cubehelix_palette(start=0.1, rot=-0.7, gamma=0.2, as_cmap=True)
cmap = plt.cm.CMRmap
im = ax.imshow(val,
# norm=mpl.colors.LogNorm(vmin=1.515e-01, vmax=4.118e+01),
norm=mpl.colors.LogNorm(vmin=1e-04, vmax=1e+01),
cmap=cmap, origin='lower', extent=[-w, w, -w, w], aspect=1)
# draw the flux extraction regions
# x = 100
# y = 100
# area = x*y / 4.25e10
# offset = 50
#
# pos_n = (len(val[0,:])/2.-1,len(val[0,:])/2.-1 + offset*len(val[0,:])/2/w)
# pos_s = (len(val[0,:])/2.-1,len(val[0,:])/2.-1 - offset*len(val[0,:])/2/w)
#
# import matplotlib.patches as patches
# ax.add_patch(patches.Rectangle((-x/2, -y), x, y, fill=False, edgecolor='lime'))
# ax.add_patch(patches.Rectangle((-x/2, 0), x, y, fill=False, edgecolor='lime'))
# plot the marker for center position by default or user input offset
ax.plot([0],[-marker], '+', color='lime', markersize=10, mew=2)
ax.set_xlim([-w,w])
ax.set_ylim([-w,w])
# ax.plot([0],[-10], '+', color='m', markersize=10, mew=2)
print(w)
# plot scalebar
if scalebar != None:
ax.plot([0.85*w-scalebar, 0.85*w], [-0.8*w, -0.8*w], color='w', linewidth=3)
# add text
ax.text(0.85*w-scalebar/2, -0.9*w, r'$\rm{'+str(scalebar)+"\,arcsec}$",
color='w', fontsize=18, fontweight='bold', ha='center')
# fix the tick label font
ticks_font = mpl.font_manager.FontProperties(family='STIXGeneral',size=16)
for label in ax.get_xticklabels():
label.set_fontproperties(ticks_font)
for label in ax.get_yticklabels():
label.set_fontproperties(ticks_font)
# Colorbar setting
#.........这里部分代码省略.........