Python Model.add_density_grid方法代码示例

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]
import numpy as np

from hyperion.model import Model
from hyperion.util.constants import pc, lsun

# Initialize model
m = Model()

# Set one-cell cartesian grid
w = np.linspace(-pc, pc, 32)
m.set_cartesian_grid(w, w, w)

# Add density grid with constant density
m.add_density_grid(np.ones(m.grid.shape) * 4.e-20, 'kmh_lite.hdf5')

# Add a point source in the center
s = m.add_point_source()
s.luminosity = 1000 * lsun
s.temperature = 6000.

# Add 10 SEDs for different viewing angles
image = m.add_peeled_images(sed=True, image=False)
image.set_wavelength_range(250, 0.01, 5000.)
image.set_viewing_angles(np.linspace(0., 90., 10), np.repeat(20., 10))
image.set_track_origin('basic')

# Add multi-wavelength image for a single viewing angle
image = m.add_peeled_images(sed=False, image=True)
image.set_wavelength_range(30, 1., 1000.)
image.set_viewing_angles([30.], [20.])
image.set_image_size(200, 200)

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]
from hyperion.util.constants import au, lsun, rsun
from hyperion.dust import SphericalDust


# Model
m = Model()
dist = 20000 * au
x = np.linspace(-dist, dist, 101)
y = np.linspace(-dist, dist, 101)
z = np.linspace(-dist, dist, 101)
m.set_cartesian_grid(x,y,z)

# Dust
d = SphericalDust('kmh.hdf5')
d.set_sublimation_temperature('fast', temperature=1600.)
m.add_density_grid(np.ones((100,100,100)) * 1.e-18,'kmh.hdf5')

# Alpha centauri
sourceA = m.add_spherical_source()
sourceA.luminosity = 1.519 * lsun
sourceA.radius = 1.227 * rsun
sourceA.temperature = 5790.
sourceA.position = (0., 0., 0.)

# Beta centauri
sourceB = m.add_spherical_source()
sourceB.luminosity = 0.5 * lsun
sourceB.radius = 0.865 * rsun
sourceB.temperature = 5260.
sourceB.position = (-11.2 * au, 0., 0.)

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]

#.........这里部分代码省略.........
            return y0
        rho_env_copy = np.array(rho_env_tsc)
        for ithetac in range(0, len(thetac)):
            rho_dum = np.log10(rho_env_copy[(rc > 1.1*R_inf) & (np.isnan(rho_env_copy[:,ithetac]) == False),ithetac])
            rc_dum = np.log10(rc[(rc > 1.1*R_inf) & (np.isnan(rho_env_copy[:,ithetac]) == False)])
            rc_dum_nan = np.log10(rc[(rc > 1.1*R_inf) & (np.isnan(rho_env_copy[:,ithetac]) == True)])
            for i in range(0, len(rc_dum_nan)):
                rho_extrapol = poly(rc_dum, rho_dum, rc_dum_nan[i])
                rho_env_copy[(np.log10(rc) == rc_dum_nan[i]),ithetac] = 10**rho_extrapol
        rho_env2d = rho_env_copy
        rho_env = np.empty((nx,ny,nz))
        for i in range(0, nz):
            rho_env[:,:,i] = rho_env2d
        # create the array of density of disk and the whole structure
        #
        rho      = np.zeros([len(rc),len(thetac),len(phic)])
        # The function for calculating the normalization of disk using the total disk mass
        #
        for ir in range(0,len(rc)):
            for itheta in range(0,len(thetac)):
                for iphi in range(0,len(phic)):
                    if rc[ir] > rin_shell:
                        # Envelope profile
                        rho[ir,itheta,iphi] = rho_env[ir,itheta,iphi]
                    else:
                        rho[ir,itheta,iphi] = 1e-25
        rho_env  = rho_env  + 1e-40
        rho      = rho      + 1e-40

    # Call function to plot the density
    plot_density(rho, rc, thetac,'/Users/yaolun/bhr71/hyperion/', plotname='shell')
    # Insert the calculated grid and dust density profile into hyperion
    m.set_spherical_polar_grid(ri, thetai, phii)
    m.add_density_grid(rho.T, outdir+'oh5.hdf5')    # numpy read the array in reverse order

    # Define the luminsoity source
    source = m.add_spherical_source()
    source.luminosity = (4*PI*rstar**2)*sigma*(tstar**4)  # [ergs/s]
    source.radius = rstar  # [cm]
    source.temperature = tstar  # [K]
    source.position = (0., 0., 0.)
    print 'L_center =  % 5.2f L_sun' % ((4*PI*rstar**2)*sigma*(tstar**4)/LS)

    # Setting up the wavelength for monochromatic radiative transfer
    lambda0 = 0.1
    lambda1 = 2.0
    lambda2 = 50.0
    lambda3 = 95.0
    lambda4 = 200.0
    lambda5 = 314.0
    lambda6 = 670.0
    n01     = 10.0
    n12     = 20.0
    n23     = (lambda3-lambda2)/0.02
    n34     = (lambda4-lambda3)/0.03
    n45     = (lambda5-lambda4)/0.1
    n56     = (lambda6-lambda5)/0.1

    lam01   = lambda0 * (lambda1/lambda0)**(np.arange(n01)/n01)
    lam12   = lambda1 * (lambda2/lambda1)**(np.arange(n12)/n12)
    lam23   = lambda2 * (lambda3/lambda2)**(np.arange(n23)/n23)
    lam34   = lambda3 * (lambda4/lambda3)**(np.arange(n34)/n34)
    lam45   = lambda4 * (lambda5/lambda4)**(np.arange(n45)/n45)
    lam56   = lambda5 * (lambda6/lambda5)**(np.arange(n56+1)/n56)

    lam     = np.concatenate([lam01,lam12,lam23,lam34,lam45,lam56])

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]
in_box = (x > 0.) & (x < 60 * au) & (y > 0.) & (y < 60 * au) & (z > 0.) & (z < 60 * au)
density[in_box] = rho_0

# Set up sphere 1
in_sphere_1 = (x - 10 * au) ** 2 + (y - 15 * au) ** 2 + (z - 20 * au) ** 2 < r_1 ** 2
density[in_sphere_1] = rho_1

# Set up sphere 2
in_sphere_2 = (x - 26.666667 * au) ** 2 + (y - 31.666667 * au) ** 2 + (z - 28.333333 * au) ** 2 < r_2 ** 2
density[in_sphere_2] = rho_2

# Remove dust close to source
in_rsub = np.sqrt(x * x + y * y + z * z) < RSUB
density[in_rsub] = 0.

m.add_density_grid(density, d)

# m.set_propagation_check_frequency(1.0)

# Set up illuminating source:
s = m.add_spherical_source()
s.radius = 6.6 * rsun
s.temperature = 33000.
s.luminosity = 4 * pi * s.radius ** 2 * sigma * s.temperature ** 4

# Set up number of photons
m.set_n_photons(initial=NPHOTONS, imaging=0)

# Write out and run
m.write(os.path.join('models', 'bm2_eff_vor_temperature.rtin'), overwrite=True)
m.run(os.path.join('models', 'bm2_eff_vor_temperature.rtout'), mpi=True, overwrite=True)

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]
from hyperion.util.constants import au, lsun, rsun
from hyperion.dust import SphericalDust


# Model
m = Model()
dist = 20000 * au
x = np.linspace(-dist, dist, 101)
y = np.linspace(-dist, dist, 101)
z = np.linspace(-dist, dist, 101)
m.set_cartesian_grid(x, y, z)

# Dust
d = SphericalDust("kmh.hdf5")
d.set_sublimation_temperature("fast", temperature=1600.0)
m.add_density_grid(np.ones((100, 100, 100)) * 1.0e-18, "kmh.hdf5")

# Alpha centauri
sourceA = m.add_spherical_source()
sourceA.luminosity = 1.519 * lsun
sourceA.radius = 1.227 * rsun
sourceA.temperature = 5790.0
sourceA.position = (0.0, 0.0, 0.0)

# Beta centauri
sourceB = m.add_spherical_source()
sourceB.luminosity = 0.5 * lsun
sourceB.radius = 0.865 * rsun
sourceB.temperature = 5260.0
sourceB.position = (-11.2 * au, 0.0, 0.0)

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]

#.........这里部分代码省略.........
			for j in range(0, grid_Nz):
				for i in range(0, grid_Nw):
				
					t = model.grid.gz[k,j,i] / (10*pc)
				
					if(cli.filament == "linear"):
						filament_center_x  = 0
						filament_center_y  = 0
					elif(cli.filament == "spiraling"):
						filament_center_x  = - math.sin(2*pi*t)*pc
						filament_center_y  = + math.cos(2*pi*t)*pc
				
					spherical_grid_r   = model.grid.gw[k,j,i]
					spherical_grid_phi = model.grid.gp[k,j,i]
				
					cartesian_grid_x   = spherical_grid_r * math.cos(spherical_grid_phi)
					cartesian_grid_y   = spherical_grid_r * math.sin(spherical_grid_phi)
				
					rsquared = (
								(cartesian_grid_x - filament_center_x)**2
								+
								(cartesian_grid_y - filament_center_y)**2
								)
				
					rho[k,j,i] = nC / (1.0 + (rsquared / (RC*RC)))
				
					if rsquared**0.5 > 3*pc:
						rho[k,j,i] = 0

		rho[model.grid.gw > grid_wmax] = 0
		rho[model.grid.gz < grid_zmin] = 0
		rho[model.grid.gz > grid_zmax] = 0

		model.add_density_grid(rho, 'dust_properties.hdf5')


		#
		# Check optical depth through the filament:
		#
		#  (y,z = 0, 2.5 pc goes through the filament center in all setups)
		
		#
		# Determine index of closest grid cell to z = 2.5 pc:
		#
		dz_last = 2*abs(grid_zmax-grid_zmin)
		for j in range(0, grid_Nz):
			dz = abs(model.grid.gz[0,j,0] - 2.5*pc)
			if(dz > dz_last):
				j=j-1
				break
			else:
				dz_last = dz

		#
		# Opacity at 1.0 micron (per gram dust):
		#
		chi = dust_properties.optical_properties.interp_chi_wav(1.0)

		tau_max = 0
		for k in range(0, grid_Np):
			tau = 0
			for i in range(0, grid_Nw):
				dr = model.grid.widths[0,k,j,i]
				dtau = dr * rho[k,j,i] * chi
				tau += dtau
			tau_max = max(tau_max, tau)

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [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")

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]

#.........这里部分代码省略.........
                            else:
                                # condition for the outer ellipsoid
                                cav_con = (2*(w/b_out)**2 + ((abs(z)-z_out)/a_out)**2) < 1

                            if cav_con:
                                # open cavity
                                if ellipsoid == False:
                                    if (rc[ir] <= rho_cav_edge) & (rc[ir] >= R_env_min):
                                        rho_dum = g2d * rho_cav_center
                                    else:
                                        rho_dum = g2d * rho_cav_center*discont*(rho_cav_edge/rc[ir])**power
                                else:
                                    # condition for the inner ellipsoid
                                    if (2*(w/b_in)**2 + ((abs(z)-z_in)/a_in)**2) > 1:
                                        rho_dum = rho_cav_out
                                    else:
                                        rho_dum = rho_cav_in

                        rho[ir, itheta, iphi] = rho_env[ir, itheta, iphi] + rho_dum

                    else:
                        rho[ir,itheta,iphi] = 1e-40

                    # add the dust mass into the total count
                    cell_mass = rho[ir, itheta, iphi] * (1/3.)*(ri[ir+1]**3 - ri[ir]**3) * (phii[iphi+1]-phii[iphi]) * -(np.cos(thetai[itheta+1])-np.cos(thetai[itheta]))
                    total_mass = total_mass + cell_mass
    # apply gas-to-dust ratio of 100
    rho_dust = rho/g2d
    total_mass_dust = total_mass/MS/g2d
    print('Total dust mass = %f Solar mass' % total_mass_dust)

    # Insert the calculated grid and dust density profile into hyperion
    m.set_spherical_polar_grid(ri, thetai, phii)
    m.add_density_grid(rho_dust.T, d)

    # Define the luminsoity source
    source = m.add_spherical_source()
    source.luminosity = (4*PI*rstar**2)*sigma*(tstar**4)  # [ergs/s]
    source.radius = rstar  # [cm]
    source.temperature = tstar  # [K]
    source.position = (0., 0., 0.)
    print('L_center =  % 5.2f L_sun' % ((4*PI*rstar**2)*sigma*(tstar**4)/LS))

    # radiative transfer settigs
    m.set_raytracing(True)

    # determine the number of photons for imaging
    # the case of monochromatic
    if mono_wave != None:
        if (type(mono_wave) == int) or (type(mono_wave) == float) or (type(mono_wave) == str):
            mono_wave = float(mono_wave)
            mono_wave = [mono_wave]

        # Monochromatic radiative transfer setting
        m.set_monochromatic(True, wavelengths=mono_wave)
        m.set_n_photons(initial=mc_photons, imaging_sources=im_photon,
                        imaging_dust=im_photon, raytracing_sources=im_photon,
                        raytracing_dust=im_photon)
    # regular SED
    else:
        m.set_n_photons(initial=mc_photons, imaging=im_photon * wav_num,
                        raytracing_sources=im_photon,
                        raytracing_dust=im_photon)
    # number of iteration to compute dust specific energy (temperature)
    m.set_n_initial_iterations(20)
    m.set_convergence(True, percentile=dict_params['percentile'],

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]
# the dust is directly behind, and the source is moving towards the dust, we
# should see red-shifted emission from the source and blue-shifted scattered
# light emission.

import numpy as np
from hyperion.model import Model
from hyperion.util.constants import c

m = Model()

m.set_cartesian_grid([-1.0, 0, 1], [-1.0, 1.0], [-1.0, 1])

density = np.zeros(m.grid.shape)
density[:, :, 0] = 1.0

m.add_density_grid(density, "kmh_lite.hdf5")

# narrow emission line spectrum at 1 micron
wav = np.array([0.9999, 1.0001])
fnu = np.array([1.0, 1.0])
nu = c / (wav * 1.0e-4)

s = m.add_spherical_source()
s.position = 0.5, 0.0, 0.0
s.velocity = -1e8, 0.0, 0.0
s.spectrum = nu[::-1], fnu[::-1]
s.luminosity = 1
s.radius = 0.1

# Set up images

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]

#.........这里部分代码省略.........
                        #                 else:
                        #                     mu_o_dum = roots[imu]
                        #         if mu_o_dum == -0.5:
                        #             print 'Problem with cubic solving, roots are: ', roots
                        #     mu_o = mu_o_dum.real
                        #     rho_env[ir,itheta,iphi] = M_env_dot/(4*PI*(G*mstar*rcen**3)**0.5)*(rc[ir]/rcen)**(-3./2)*(1+mu/mu_o)**(-0.5)*(mu/mu_o+2*mu_o**2*rcen/rc[ir])**(-1)
                        # # Disk profile
                        # if ((w >= R_disk_min) and (w <= R_disk_max)) == True:
                        #     h = ((w/(100*AU))**beta)*h100
                        #     rho_disk[ir,itheta,iphi] = rho_0*(1-np.sqrt(rstar/w))*(rstar/w)**(beta+1)*np.exp(-0.5*(z/h)**2)
                        # # Combine envelope and disk
                        # rho[ir,itheta,iphi] = rho_disk[ir,itheta,iphi] + rho_env[ir,itheta,iphi]
                    else:
                        rho[ir,itheta,iphi] = 1e-30
        rho_env  = rho_env  + 1e-40
        rho_disk = rho_disk + 1e-40
        rho      = rho      + 1e-40
    else:
        for ir in range(0,len(rc)):
            for itheta in range(0,len(thetac)):
                for iphi in range(0,len(phic)):
                    # Envelope profile
                    w = abs(rc[ir]*np.cos(thetac[itheta]))
                    z = rc[ir]*np.sin(thetac[itheta])
                    z_cav = c*abs(w)**1.5
                    z_cav_wall = c*abs(w-wall)**1.5
                    if z_cav == 0:
                        z_cav = R_env_max
                    if abs(z) > abs(z_cav):
                        # rho_env[ir,itheta,iphi] = rho_cav
                        # Modification for using density gradient in the cavity
                        if rc[ir] <= 20*AU:
                            rho_env[ir,itheta,iphi] = rho_cav_center*((rc[ir]/AU)**2)
                        else:
                            rho_env[ir,itheta,iphi] = rho_cav_center*discont*(20*AU/rc[ir])**2
                        i += 1
                    elif (abs(z) > abs(z_cav_wall)) and (abs(z) < abs(z_cav)):
                        rho_env[ir,itheta,iphi] = rho_wall
                    else:
                        j += 1
                        mu = abs(np.cos(thetac[itheta]))
                        mu_o = np.abs(fsolve(func,[0.5,0.5,0.5],args=(rc[ir],rcen,mu))[0])
                        rho_env[ir,itheta,iphi] = M_env_dot/(4*PI*(G*mstar*rcen**3)**0.5)*(rc[ir]/rcen)**(-3./2)*(1+mu/mu_o)**(-0.5)*(mu/mu_o+2*mu_o**2*rcen/rc[ir])**(-1)
                    # Disk profile
                    if ((w >= R_disk_min) and (w <= R_disk_max)) == True:
                        h = ((w/(100*AU))**beta)*h100
                        rho_disk[ir,itheta,iphi] = rho_0*(1-np.sqrt(rstar/w))*(rstar/w)**(beta+1)*np.exp(-0.5*(z/h)**2)
                    # Combine envelope and disk
                    rho[ir,itheta,iphi] = rho_disk[ir,itheta,iphi] + rho_env[ir,itheta,iphi]
        rho_env  = rho_env  + 1e-40
        rho_disk = rho_disk + 1e-40
        rho      = rho      + 1e-40

    # Insert the calculated grid and dust density profile into hyperion
    m.set_spherical_polar_grid(ri, thetai, phii)
    m.add_density_grid(rho.T, outdir+'oh5.hdf5')    # numpy read the array in reverse order

    # Define the luminsoity source
    source = m.add_spherical_source()
    source.luminosity = (4*PI*rstar**2)*sigma*(tstar**4)  # [ergs/s]
    source.radius = rstar  # [cm]
    source.temperature = tstar  # [K]
    source.position = (0., 0., 0.)
    print 'L_center =  % 5.2f L_sun' % ((4*PI*rstar**2)*sigma*(tstar**4)/LS)

    # Setting up images and SEDs
    image = m.add_peeled_images()
    image.set_wavelength_range(300, 2.0, 670.0)
    # pixel number
    image.set_image_size(300, 300)
    image.set_image_limits(-R_env_max, R_env_max, -R_env_max, R_env_max)
    image.set_viewing_angles([82.0], [0.0])
    image.set_uncertainties(True)
    # output as 64-bit
    image.set_output_bytes(8)

    # Radiative transfer setting

    # number of photons for temp and image
    m.set_raytracing(True)
    m.set_n_photons(initial=1000000, imaging=1000000, raytracing_sources=1000000, raytracing_dust=1000000)
    # number of iteration to compute dust specific energy (temperature)
    m.set_n_initial_iterations(5)
    m.set_convergence(True, percentile=99., absolute=1.5, relative=1.02)
    m.set_mrw(True)   # Gamma = 1 by default

    # Output setting
    # Density
    m.conf.output.output_density = 'last'

    # Density difference (shows where dust was destroyed)
    m.conf.output.output_density_diff = 'none'

    # Energy absorbed (using pathlengths)
    m.conf.output.output_specific_energy = 'last'

    # Number of unique photons that passed through the cell
    m.conf.output.output_n_photons = 'last'

    m.write(outdir+'old_setup2.rtin')

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]
from hyperion.model import Model
from hyperion.util.constants import pc, lsun

# Initialize model
m = Model()

# Set up 64x64x64 cartesian grid
w = np.linspace(-pc, pc, 64)
m.set_cartesian_grid(w, w, w)

# Add density grid with constant density and add a higher density cube inside to
# cause a shadow.
density = np.ones(m.grid.shape) * 1e-21
density[26:38, 26:38, 26:38] = 1.e-18
m.add_density_grid(density, 'kmh_lite.hdf5')

# Add a point source in the center
s = m.add_point_source()
s.position = (0.4 * pc, 0., 0.)
s.luminosity = 1000 * lsun
s.temperature = 6000.

# Add multi-wavelength image for a single viewing angle
image = m.add_peeled_images(sed=False, image=True)
image.set_wavelength_range(1, 190., 210.)
image.set_viewing_angles(np.repeat(45., 36), np.linspace(5., 355., 36))
image.set_image_size(400, 400)
image.set_image_limits(-1.5 * pc, 1.5 * pc, -1.5 * pc, 1.5 * pc)

# Set runtime parameters. We turn off scattering for the imaging since it is not

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]
import random
random.seed('hyperion')  # ensure that random numbers are the same every time

import numpy as np
from hyperion.model import Model
from hyperion.util.constants import pc, lsun

# Define cell walls
x = np.linspace(-10., 10., 101) * pc
y = np.linspace(-10., 10., 101) * pc
z = np.linspace(-10., 10., 101) * pc

# Initialize model and set up density grid
m = Model()
m.set_cartesian_grid(x, y, z)
m.add_density_grid(np.ones((100, 100, 100)) * 1.e-20, 'kmh_lite.hdf5')

# Generate random sources
for i in range(100):
    s = m.add_point_source()
    xs = random.uniform(-10., 10.) * pc
    ys = random.uniform(-10., 10.) * pc
    zs = random.uniform(-10., 10.) * pc
    s.position = (xs, ys, zs)
    s.luminosity = 10. ** random.uniform(0., 3.) * lsun
    s.temperature = random.uniform(3000., 8000.)

# Specify that the specific energy and density are needed
m.conf.output.output_specific_energy = 'last'
m.conf.output.output_density = 'last'

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]
import numpy as np

from hyperion.model import Model
from hyperion.util.constants import pc, lsun

# Initialize model
m = Model()

# Set one-cell cartesian grid
w = np.linspace(-pc, pc, 32)
m.set_cartesian_grid(w, w, w)

# Add density grid with constant density
m.add_density_grid(np.ones(m.grid.shape) * 4.0e-20, "kmh_lite.hdf5")

# Add a point source in the center
s = m.add_point_source()
s.luminosity = 1000 * lsun
s.temperature = 6000.0

# Add 10 SEDs for different viewing angles
image = m.add_peeled_images(sed=True, image=False)
image.set_wavelength_range(250, 0.01, 5000.0)
image.set_viewing_angles(np.linspace(0.0, 90.0, 10), np.repeat(20.0, 10))
image.set_track_origin("basic")

# Add multi-wavelength image for a single viewing angle
image = m.add_peeled_images(sed=False, image=True)
image.set_wavelength_range(30, 1.0, 1000.0)
image.set_viewing_angles([30.0], [20.0])
image.set_image_size(200, 200)

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]

#.........这里部分代码省略.........
        [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)

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

        ax.set_ylim([0,15])
        fig.gca().set_xlim(left=np.log10(0.05))
        # ax.set_xlim([np.log10(0.8),np.log10(10000)])

        # subplot shows the radial density profile along the midplane
        ax_mid = plt.axes([0.2,0.2,0.2,0.2], frameon=True)
        ax_mid.plot(np.log10(rc/AU), np.log10(rho2d[:,199]/g2d/mmw/mh),'o',color='b',linewidth=1, markersize=2)
        ax_mid.plot(np.log10(rc/AU), np.log10(rho_env_tsc2d[:,199]/mmw/mh),'-',color='r',linewidth=1, markersize=2)
        # ax_mid.set_ylim([0,10])
        # ax_mid.set_xlim([np.log10(0.8),np.log10(10000)])
        ax_mid.set_ylim([0,15])
        fig.savefig(outdir+outname+'_gas_radial.pdf',format='pdf',dpi=300,bbox_inches='tight')
        fig.clf()

    # Insert the calculated grid and dust density profile into hyperion
    m.set_spherical_polar_grid(ri, thetai, phii)
    # temperary for comparing full TSC and infall-only TSC model
    # import sys
    # sys.path.append(os.path.expanduser('~')+'/programs/misc/')
    # from tsc_comparison import tsc_com
    # rho_tsc, rho_ulrich = tsc_com()
    m.add_density_grid(rho_dust.T, d)
    # m.add_density_grid(rho.T, outdir+'oh5.hdf5')    # numpy read the array in reverse order

    # Define the luminsoity source
    source = m.add_spherical_source()
    source.luminosity = (4*PI*rstar**2)*sigma*(tstar**4)  # [ergs/s]
    source.radius = rstar  # [cm]
    source.temperature = tstar  # [K]
    source.position = (0., 0., 0.)
    print 'L_center =  % 5.2f L_sun' % ((4*PI*rstar**2)*sigma*(tstar**4)/LS)

    # # add an infrared source at the center
    # L_IR = 0.04
    # ir_source = m.add_spherical_source()
    # ir_source.luminosity = L_IR*LS
    # ir_source.radius = rstar      # [cm]
    # ir_source.temperature = 500 # [K]  peak at 10 um
    # ir_source.position = (0., 0., 0.)
    # print 'Additional IR source, L_IR = %5.2f L_sun' % L_IR

    # Setting up the wavelength for monochromatic radiative transfer
    lambda0 = 0.1
    lambda1 = 2.0
    lambda2 = 50.0
    lambda3 = 95.0
    lambda4 = 200.0
    lambda5 = 314.0
    lambda6 = 1000.0
    n01     = 10.0
    n12     = 20.0
    n23     = 50.0

    lam01   = lambda0 * (lambda1/lambda0)**(np.arange(n01)/n01)

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import add_density_grid [as 别名]
import numpy as np
from hyperion.model import Model
from hyperion.util.constants import c

m = Model()

m.set_cartesian_grid([-1.,0, 1], [-1., 1.], [-1., 1])

density = np.zeros(m.grid.shape)
density[:,:,0] = 1.

vx = np.ones(m.grid.shape) * -1e8
vy = np.zeros(m.grid.shape)
vz = np.zeros(m.grid.shape)

m.add_density_grid(density, 'kmh_lite.hdf5', velocity=(vx, vy, vz))

# narrow emission line spectrum at 1 micron
wav = np.array([0.9999, 1.0001])
fnu = np.array([1., 1.])
nu = c / (wav * 1.e-4)

s = m.add_spherical_source()
s.position = 0.5, 0., 0.
s.velocity = -1e8, 0., 0.
s.spectrum = nu[::-1], fnu[::-1]
s.luminosity = 1
s.radius = 0.1

# Set up images