Python model.Model类代码示例

def setup_model_shell(indir,outdir,outname,rin_shell=None,denser_wall=False,tsc=True,idl=False,plot=False,low_res=False,flat=True,scale=1.0):
    import numpy as np
    import astropy.constants as const
    import scipy as sci
    import matplotlib.pyplot as plt
    import matplotlib as mat
    import os
    from matplotlib.colors import LogNorm
    from scipy.optimize import fsolve
    from scipy.optimize import newton
    from scipy.integrate import nquad
    from envelope_func import func
    import hyperion as hp
    from hyperion.model import Model
    from plot_density import plot_density

    # Constants setup
    c         = const.c.cgs.value
    AU        = 1.49598e13     # Astronomical Unit       [cm]
    pc        = 3.08572e18     # Parsec                  [cm]
    MS        = 1.98892e33     # Solar mass              [g]
    LS        = 3.8525e33      # Solar luminosity        [erg/s]
    RS        = 6.96e10        # Solar radius            [cm]
    G         = 6.67259e-8     # Gravitational constant  [cm3/g/s^2]
    yr        = 60*60*24*365   # Years in seconds
    PI        = np.pi          # PI constant
    sigma     = const.sigma_sb.cgs.value  # Stefan-Boltzmann constant 


    m = Model()

    # Create dust properties

    # Hyperion needs nu, albedo, chi, g, p_lin_max
    from hyperion.dust import HenyeyGreensteinDust
    # Read in the dust opacity table used by RADMC-3D
    dust_radmc = dict()
    [dust_radmc['wl'], dust_radmc['abs'], dust_radmc['scat'], dust_radmc['g']] = np.genfromtxt('dustkappa_oh5_extended.inp',skip_header=2).T
    # opacity per mass of dust?
    dust_hy = dict()
    dust_hy['nu'] = c/dust_radmc['wl']*1e4
    ind = np.argsort(dust_hy['nu'])
    dust_hy['nu'] = dust_hy['nu'][ind]
    dust_hy['albedo'] = (dust_radmc['scat']/(dust_radmc['abs']+dust_radmc['scat']))[ind]
    dust_hy['chi'] = (dust_radmc['abs']+dust_radmc['scat'])[ind]
    dust_hy['g'] = dust_radmc['g'][ind]
    dust_hy['p_lin_max'] = 0*dust_radmc['wl'][ind]     # assume no polarization

    d = HenyeyGreensteinDust(dust_hy['nu'], dust_hy['albedo'], dust_hy['chi'], dust_hy['g'], dust_hy['p_lin_max'])
    # dust sublimation does not occur
    # d.set_sublimation_temperature(None)
    d.write(outdir+'oh5.hdf5')
    d.plot(outdir+'oh5.png')

    # Grids and Density
    # Calculation inherited from the script used for RADMC-3D

    # Grid Parameters
    nx        = 300L
    if low_res == True:
        nx    = 100L
    ny        = 400L
    nz        = 50L
    [nx, ny, nz] = [scale*nx, scale*ny, scale*nz]

    if tsc == False:
        # Parameters setup
        # Import the model parameters from another file 
        #
        params     = np.genfromtxt(indir+'/params.dat',dtype=None)
        tstar      = params[0][1]
        mstar      = params[1][1]*MS
        rstar      = params[2][1]*RS
        M_env_dot  = params[3][1]*MS/yr
        M_disk_dot = params[4][1]*MS/yr
        R_env_max  = params[5][1]*AU
        R_env_min  = params[6][1]*AU
        theta_cav  = params[7][1]
        R_disk_max = params[8][1]*AU
        R_disk_min = params[9][1]*AU
        R_cen      = R_disk_max
        M_disk     = params[10][1]*MS
        beta       = params[11][1]
        h100       = params[12][1]*AU
        rho_cav    = params[13][1]
        if denser_wall == True:
            wall       = params[14][1]*AU
            rho_wall   = params[15][1]
        rho_cav_center = params[16][1]
        rho_cav_edge   = params[17][1]*AU

        # Model Parameters
        #
        rin       = rstar
        rout      = R_env_max
        rcen      = R_cen

        # Star Parameters
        #
        mstar    = mstar
#.........这里部分代码省略.........

# A simple model to check what happens when a source is moving towards dust and
# we observe both the source and the dust. If we observe the source such that
# 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, 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.

def extract_hyperion(filename,indir=None,outdir=None,dstar=200.0,aperture=None,
                     save=True,filter_func=False,plot_all=False,clean=False,
                     exclude_wl=[],log=True,image=True,obj='BHR71',
                     print_data_w_aper=False,mag=1.5):
    """
    filename: The path to Hyperion output file
    indir: The path to the directory which contains observations data
    outdir: The path to the directory for storing extracted plots and ASCII files
    """
    def l_bol(wl,fv,dstar):
        import numpy as np
        import astropy.constants as const
        # wavelength unit: um
        # Flux density unit: Jy
        # constants setup
        #
        c = const.c.cgs.value
        pc = const.pc.cgs.value
        PI = np.pi
        SL = const.L_sun.cgs.value
        # Convert the unit from Jy to erg s-1 cm-2 Hz-1
        fv = np.array(fv)*1e-23
        freq = c/(1e-4*np.array(wl))

        diff_dum = freq[1:]-freq[0:-1]
        freq_interpol = np.hstack((freq[0:-1]+diff_dum/2.0,freq[0:-1]+diff_dum/2.0,freq[0],freq[-1]))
        freq_interpol = freq_interpol[np.argsort(freq_interpol)[::-1]]
        fv_interpol = np.empty(len(freq_interpol))
        # calculate the histogram style of spectrum
        #
        for i in range(0,len(fv)):
            if i == 0:
                fv_interpol[i] = fv[i]
            else:
                fv_interpol[2*i-1] = fv[i-1]
                fv_interpol[2*i] = fv[i]
        fv_interpol[-1] = fv[-1]

        dv = freq_interpol[0:-1]-freq_interpol[1:]
        dv = np.delete(dv,np.where(dv==0))

        fv = fv[np.argsort(freq)]
        freq = freq[np.argsort(freq)]

        return (np.trapz(fv,freq)*4.*PI*(dstar*pc)**2)/SL

    # function for properly calculating uncertainty of spectrophotometry value
    def unc_spectrophoto(wl, unc, trans):
        # adopting smiliar procedure as Trapezoidal rule
        # (b-a) * [ f(a) + f(b) ] / 2
        #
        return ( np.sum( trans[:-1]**2 * unc[:-1]**2 * (wl[1:]-wl[:-1])**2 ) / np.trapz(trans, x=wl)**2 )**0.5

    # to avoid X server error
    import matplotlib as mpl
    mpl.use('Agg')
    #
    import matplotlib.pyplot as plt
    import numpy as np
    import os
    from hyperion.model import ModelOutput, Model
    from scipy.interpolate import interp1d
    from hyperion.util.constants import pc, c, lsun, au
    from astropy.io import ascii
    import sys
    from phot_filter import phot_filter
    from get_obs import get_obs

    # Open the model
    m = ModelOutput(filename)

    # Read in the observation data and calculate the noise & variance
    if indir == None:
        indir = raw_input('Path to the observation data: ')
    if outdir == None:
        outdir = raw_input('Path for the output: ')

    # assign the file name from the input file
    print_name = os.path.splitext(os.path.basename(filename))[0]

    # use a canned function to extract observational data
    obs_data = get_obs(indir, obj=obj)        # unit in um, Jy
    wl_tot, flux_tot, unc_tot = obs_data['spec']
    flux_tot = flux_tot*1e-23    # convert unit from Jy to erg s-1 cm-2 Hz-1
    unc_tot = unc_tot*1e-23
    l_bol_obs = l_bol(wl_tot, flux_tot*1e23, dstar)

    wl_phot, flux_phot, flux_sig_phot = obs_data['phot']
    flux_phot = flux_phot*1e-23   # convert unit from Jy to erg s-1 cm-2 Hz-1
    flux_sig_phot = flux_sig_phot*1e-23

    if aperture == None:
        aperture = {'wave': [3.6, 4.5, 5.8, 8.0, 8.5, 9, 9.7, 10, 10.5, 11, 16, 20, 24, 35, 70, 100, 160, 250, 350, 500, 850],\
                    'aperture': [7.2, 7.2, 7.2, 7.2, 7.2, 7.2, 7.2, 7.2, 7.2, 7.2, 20.4, 20.4, 20.4, 20.4, 24.5, 24.5, 24.5, 24.5, 24.5, 24.5, 24.5]}
    # assign wl_aper and aper from dictionary of aperture
    wl_aper = aperture['wave']
    aper    = aperture['aperture']
    # create the non-repetitive aperture list and index array
    aper_reduced = list(set(aper))
    index_reduced = np.arange(1, len(aper_reduced)+1)  # '+1': the zeroth slice corresponds to infinite aperture
#.........这里部分代码省略.........

import numpy as np
from hyperion.model import Model
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.

def setup_model(cli):
	
    lsun_TRUST = 3.839e33
        
    #
    # Hyperion setup:
    #
    model = Model()


    if(cli.mode == "temperature"):
        #
        # Dust properties:
        #
        dust_properties = SphericalDust('dust_integrated_full_scattering.hdf5')
            
            
        #
        # Write dust properties:
        #
        dust_properties.write('dust_properties.hdf5')
        dust_properties.plot('dust_properties.png')
        
        
        #
        # Specify galaxy setup:
        #
        hR                     =  4000.0*pc             # [cm]
        Rmax                   =     5.0*hR             # [cm]
        hz_oldstars            =   350.0*pc             # [cm]
        hz_youngstars          =   200.0*pc             # [cm]
        hz_dust                =   200.0*pc             # [cm]
        zmax_oldstars          =     5.0*hz_oldstars    # [cm]
        zmax_youngstars        =     5.0*hz_youngstars  # [cm]
        zmax_dust              =     5.0*hz_dust        # [cm]
        zmax                   =  zmax_oldstars         # [cm]
        reff                   =  1600.0*pc             # [cm]
        n                      =     3.0
        q                      =     0.6
        bn                     = 2.0*n - 1.0/3.0 + 4.0/405.0/n + 46.0/25515.0/n/n + 131.0/1148175.0/n/n/n
        temperature_oldstars   =  3500.0                # [K]
        temperature_youngstars = 10000.0                # [K]
        temperature_bulge      =  3500.0                # [K]
        luminosity_oldstars    =     4.0e+10*lsun_TRUST # [ergs/s]
        luminosity_youngstars  =     1.0e+10*lsun_TRUST # [ergs/s]
        luminosity_bulge       =     3.0e+10*lsun_TRUST # [ergs/s]
        
        w_oldstars             =     0.25
        w_youngstars           =     0.75
        w_dust                 =     0.75
        phi0_oldstars          =     0.0
        phi0_youngstars        =    20.0 * pi/180.0
        phi0_dust              =    20.0 * pi/180.0
        modes                  =     2
        pitchangle             =    20.0 * pi/180.0
        
        
        
        #
        # Grid setup:
        #
        grid_wmin =  0.0
        grid_wmax =  Rmax
        grid_zmin = -zmax
        grid_zmax = +zmax
        grid_pmin =  0.0
        grid_pmax =  2.0*pi
        
        grid_dx = cli.resolution*pc
        grid_dw = grid_dx # uniform resolution
        grid_dz = grid_dx # uniform resolution
        grid_dp = grid_dx # resolution at characteristic radial disk spatial scale hR = 4000.0 pc
        
        grid_Nw   = int((grid_wmax - grid_wmin) / grid_dw) + 1
        grid_Nz   = int((grid_zmax - grid_zmin) / grid_dz) + 1
        if(cli.case == 1):
            grid_Np = 1
        if(cli.case == 2):
            grid_Np = int((grid_pmax - grid_pmin) * hR / grid_dp)
        
        if(cli.verbose):
            print("Grid setup:")
            print(" Grid resolution =",cli.resolution, "pc.")
            print(" grid_Nw =",grid_Nw)
            print(" grid_Nz =",grid_Nz)
            print(" grid_Np =",grid_Np)
        
        #grid_w      = np.logspace(np.log10(grid_wmin), np.log10(grid_wmax), grid_Nw)
        #grid_w      = np.hstack([0., grid_w]) # add innermost cell interface at w=0
        grid_w    = np.linspace(grid_wmin, grid_wmax, grid_Nw+1)
        grid_z    = np.linspace(grid_zmin, grid_zmax, grid_Nz+1)
        grid_p    = np.linspace(grid_pmin, grid_pmax, grid_Np+1)
        
        model.set_cylindrical_polar_grid(grid_w, grid_z, grid_p)
        
        #
        # Dust density and sources setup:
        #
        rho_oldstars   = np.zeros(model.grid.shape)
        rho_youngstars = np.zeros(model.grid.shape)
#.........这里部分代码省略.........

def extract_hyperion(filename,indir=None,outdir=None,dstar=178.0,wl_aper=None,save=True):
    def l_bol(wl,fv,dist=178.0):
        import numpy as np
        import astropy.constants as const
        # wavelength unit: um
        # Flux density unit: Jy
        #
        # constants setup
        #
        c = const.c.cgs.value
        pc = const.pc.cgs.value
        PI = np.pi
        SL = const.L_sun.cgs.value
        # Convert the unit from Jy to erg s-1 cm-2 Hz-1
        fv = np.array(fv)*1e-23
        freq = c/(1e-4*np.array(wl))
        
        diff_dum = freq[1:]-freq[0:-1]
        freq_interpol = np.hstack((freq[0:-1]+diff_dum/2.0,freq[0:-1]+diff_dum/2.0,freq[0],freq[-1]))
        freq_interpol = freq_interpol[np.argsort(freq_interpol)[::-1]]
        fv_interpol = np.empty(len(freq_interpol))
        # calculate the histogram style of spectrum
        #
        for i in range(0,len(fv)):
            if i == 0:
                fv_interpol[i] = fv[i]
            else:
                fv_interpol[2*i-1] = fv[i-1]
                fv_interpol[2*i] = fv[i]
        fv_interpol[-1] = fv[-1]
        
        dv = freq_interpol[0:-1]-freq_interpol[1:]
        dv = np.delete(dv,np.where(dv==0))

        fv = fv[np.argsort(freq)]
        freq = freq[np.argsort(freq)]

        return (np.trapz(fv,freq)*4.*PI*(dist*pc)**2)/SL


    import matplotlib.pyplot as plt
    import numpy as np
    import os
    from hyperion.model import ModelOutput
    from hyperion.model import Model
    from scipy.interpolate import interp1d
    from hyperion.util.constants import pc, c, lsun

    # Read in the observation data and calculate the noise & variance
    if indir == None:
        indir = '/Users/yaolun/bhr71/'
    if outdir == None:
        outdir = '/Users/yaolun/bhr71/hyperion/'

    # assign the file name from the input file
    print_name = os.path.splitext(os.path.basename(filename))[0]
    #
    [wl_pacs,flux_pacs,unc_pacs] = np.genfromtxt(indir+'BHR71_centralSpaxel_PointSourceCorrected_CorrectedYES_trim_continuum.txt',\
                                        dtype='float',skip_header=1).T
    # Convert the unit from Jy to erg cm-2 Hz-1
    flux_pacs = flux_pacs*1e-23
    [wl_spire,flux_spire] = np.genfromtxt(indir+'BHR71_spire_corrected_continuum.txt',dtype='float',skip_header=1).T
    flux_spire = flux_spire*1e-23 
    wl_obs = np.hstack((wl_pacs,wl_spire))
    flux_obs = np.hstack((flux_pacs,flux_spire))

    [wl_pacs_data,flux_pacs_data,unc_pacs_data] = np.genfromtxt(indir+'BHR71_centralSpaxel_PointSourceCorrected_CorrectedYES_trim.txt',\
                                                  dtype='float').T
    [wl_spire_data,flux_spire_data] = np.genfromtxt(indir+'BHR71_spire_corrected.txt',\
                                                    dtype='float').T

    [wl_pacs_flat,flux_pacs_flat,unc_pacs_flat] = np.genfromtxt(indir+'BHR71_centralSpaxel_PointSourceCorrected_CorrectedYES_trim_flat_spectrum.txt',\
                                        dtype='float',skip_header=1).T
    [wl_spire_flat,flux_spire_flat] = np.genfromtxt(indir+'BHR71_spire_corrected_flat_spectrum.txt',dtype='float',skip_header=1).T

    # Convert the unit from Jy to erg cm-2 Hz-1
    flux_pacs_flat = flux_pacs_flat*1e-23 
    flux_spire_flat = flux_spire_flat*1e-23
    flux_pacs_data = flux_pacs_data*1e-23
    flux_spire_data = flux_spire_data*1e-23


    wl_pacs_noise = wl_pacs_data
    flux_pacs_noise = flux_pacs_data-flux_pacs-flux_pacs_flat
    wl_spire_noise = wl_spire_data
    flux_spire_noise = flux_spire_data-flux_spire-flux_spire_flat

    # Read in the Spitzer IRS spectrum
    [wl_irs, flux_irs]= (np.genfromtxt(indir+'bhr71_spitzer_irs.txt',skip_header=2,dtype='float').T)[0:2]
    # Convert the unit from Jy to erg cm-2 Hz-1
    flux_irs = flux_irs*1e-23
    # Remove points with zero or negative flux 
    ind = flux_irs > 0
    wl_irs = wl_irs[ind]
    flux_irs = flux_irs[ind]
    # Calculate the local variance (for spire), use the instrument uncertainty for pacs
    #
    wl_noise_5 = wl_spire_noise[(wl_spire_noise > 194)*(wl_spire_noise <= 304)]
    flux_noise_5 = flux_spire_noise[(wl_spire_noise > 194)*(wl_spire_noise <= 304)]
    wl_noise_6 = wl_spire_noise[wl_spire_noise > 304]
#.........这里部分代码省略.........

import numpy as np

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)

import os

import numpy as np

from hyperion.model import Model
from hyperion.dust import SphericalDust
from hyperion.util.constants import pc, au, sigma, pi, rsun

NPHOTONS = 1e7

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

# TODO: remove dust around source

m = Model()

x = np.linspace(0., 60. * au, 256)
y = np.linspace(0., 60. * au, 256)
z = np.linspace(0., 60. * au, 256)

x = np.hstack([-10 * au, x])
y = np.hstack([-10 * au, y])
z = np.hstack([-10 * au, z])

m.set_cartesian_grid(x, y, z)

# Grain Properties:

d = SphericalDust('integrated_hg_scattering.hdf5')
chi_v = d.optical_properties.interp_chi_wav(0.55)

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'

import numpy as np
from hyperion.model import Model
from hyperion.dust import SphericalDust
from hyperion.util.constants import au

m = Model.read('bm2_eff_vor_temperature.rtout', only_initial=False)

m.set_n_initial_iterations(0)

del m.n_photons['initial']
del m.n_photons['last']

i = m.add_peeled_images()
i.set_viewing_angles([0., 90., 90., 90., 90., 180.], [0., 0., 90., 180., 270., 0.])
i.set_image_limits(-60 * au, 60 * au, -60 * au, 60 * au)
i.set_image_size(300, 300)

# Set up monochromatic mode
m.set_monochromatic(True, wavelengths=[0.10019, 0.55165, 2.00293, 10.03850, 101.15800])

# Use raytracing
m.set_raytracing(True)

# Set up number of photons
m.set_n_photons(imaging_sources=1e7, imaging_dust=1e7,
                raytracing_sources=1, raytracing_dust=1e7)

# Write out and run
m.write('bm2_eff_images.rtin', overwrite=True)
m.run('bm2_eff_images.rtout', mpi=True)

import numpy as np
from hyperion.model import Model
from hyperion.dust import SphericalDust
from hyperion.util.constants import pc

for tau_v in [0.1, 1.0, 20.0]:

    m = Model()

    # Global geometry:
    #
    # * slab
    # * system size = 10x10x10 pc
    # * system coordinates (x,y,z min/max) = -5 to +5 pc
    # * slab z extent = -2 to -5 pc
    # * slab xy extend = -5 pc to 5 pc
    # * z optical depth @0.55um in slab = 0.1, 1, 20
    # * optical depth outside slab = 0

    x = np.linspace(-5 * pc, 5 * pc, 100)
    y = np.linspace(-5 * pc, 5 * pc, 100)
    z = np.hstack([np.linspace(-5 * pc, -2 * pc, 100), 5 * pc])

    m.set_cartesian_grid(x, y, z)

    # Grain Properties:

    d = SphericalDust('integrated_hg_scattering.hdf5')
    chi_v = d.optical_properties.interp_chi_wav(0.55)

    # Determine density in slab

def setup_model(outdir,record_dir,outname,params,dust_file,tsc=True,idl=False,plot=False,\
                low_res=True,flat=True,scale=1,radmc=False,mono=False,record=True,dstar=178.,\
                aperture=None,dyn_cav=False,fix_params=None,alma=False,power=2,better_im=False,ellipsoid=False,\
                TSC_dir='~/programs/misc/TSC/', IDL_path='/Applications/exelis/idl83/bin/idl',auto_disk=0.25):
    """
    params = dictionary of the model parameters
    alma keyword is obsoleted 
    outdir: The directory for storing Hyperion input files
    record_dir: The directory contains "model_list.txt" for recording parameters
    TSC_dir: Path the TSC-related IDL routines
    IDL_path: The IDL executable 
    """
    import numpy as np
    import astropy.constants as const
    import scipy as sci
    # to avoid X server error
    import matplotlib as mpl
    mpl.use('Agg')
    #
    import matplotlib.pyplot as plt
    import os
    from matplotlib.colors import LogNorm
    from scipy.integrate import nquad
    from hyperion.model import Model
    from record_hyperion import record_hyperion
    from outflow_inner_edge import outflow_inner_edge
    from pprint import pprint
    # import pdb
    # pdb.set_trace()

    # Constants setup
    c         = const.c.cgs.value
    AU        = 1.49598e13     # Astronomical Unit       [cm]
    pc        = 3.08572e18     # Parsec                  [cm]
    MS        = 1.98892e33     # Solar mass              [g]
    LS        = 3.8525e33      # Solar luminosity        [erg/s]
    RS        = 6.96e10        # Solar radius            [cm]
    G         = 6.67259e-8     # Gravitational constant  [cm3/g/s^2]
    yr        = 60*60*24*365   # Years in seconds
    PI        = np.pi          # PI constant
    sigma     = const.sigma_sb.cgs.value  # Stefan-Boltzmann constant 
    mh        = const.m_p.cgs.value + const.m_e.cgs.value
    g2d       = 100.
    mmw       = 2.37   # Kauffmann 2008


    m = Model()

    # Create dust properties

    # Hyperion needs nu, albedo, chi, g, p_lin_max
    from hyperion.dust import HenyeyGreensteinDust
    # Read in the dust opacity table used by RADMC-3D
    dust = dict()
    # [dust_radmc['wl'], dust_radmc['abs'], dust_radmc['scat'], dust_radmc['g']] = np.genfromtxt(dust_file,skip_header=2).T
    [dust['nu'], dust['albedo'], dust['chi'], dust['g']] = np.genfromtxt(dust_file).T
    # opacity per mass of dust?
    # dust_hy = dict()
    # dust_hy['nu'] = c/dust_radmc['wl']*1e4
    # ind = np.argsort(dust_hy['nu'])
    # dust_hy['nu'] = dust_hy['nu'][ind]
    # dust_hy['albedo'] = (dust_radmc['scat']/(dust_radmc['abs']+dust_radmc['scat']))[ind]
    # dust_hy['chi'] = (dust_radmc['abs']+dust_radmc['scat'])[ind]
    # dust_hy['g'] = dust_radmc['g'][ind]
    # dust_hy['p_lin_max'] = 0*dust_radmc['wl'][ind]     # assume no polarization

    # d = HenyeyGreensteinDust(dust_hy['nu'], dust_hy['albedo'], dust_hy['chi'], dust_hy['g'], dust_hy['p_lin_max'])
    d = HenyeyGreensteinDust(dust['nu'], dust['albedo'], dust['chi'], dust['g'], dust['g']*0)
    # dust sublimation option
    d.set_sublimation_temperature('slow', temperature=1600.0)
    d.set_lte_emissivities(n_temp=3000,
                       temp_min=0.1,
                       temp_max=2000.)
    # try to solve the freq. problem
    d.optical_properties.extrapolate_nu(3.28e15, 4e15)
    #
    d.write(outdir+os.path.basename(dust_file).split('.')[0]+'.hdf5')
    d.plot(outdir+os.path.basename(dust_file).split('.')[0]+'.png')
    plt.clf()

    # Grids and Density
    # Calculation inherited from the script used for RADMC-3D

    # Grid Parameters
    nx        = 300L
    if low_res == True:
        nx    = 100L
    ny        = 400L
    nz        = 50L
    [nx, ny, nz] = [int(scale*nx), int(scale*ny), int(scale*nz)]

    # TSC model input setting
    # params    = np.genfromtxt(indir+'/tsc_params.dat', dtype=None)
    dict_params = params # input_reader(params_file)
    # TSC model parameter
    cs        = dict_params['Cs']*1e5
    t         = dict_params['age']  # year
    omega     = dict_params['Omega0']
    # calculate related parameters
    M_env_dot = 0.975*cs**3/G
#.........这里部分代码省略.........

import glob

import numpy as np

from hyperion.model import Model
from hyperion.dust import SphericalDust
from hyperion.util.constants import pc

import yaml
settings = yaml.load(open('settings.yml'))
 
WAV = np.logspace(-1, 3, 45)

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

    m = Model.read(model_path, only_initial=False)

    m.set_n_initial_iterations(0)

    del m.n_photons['initial']
    del m.n_photons['last']

    i = m.add_peeled_images(sed=True, image=False)
    i.set_viewing_angles([0., 30., 60., 90., 120., 150., 180.],
                         [0., 0., 0., 0., 0., 0., 0.])
    i.set_track_origin('basic')

    i = m.add_peeled_images(sed=True, image=False)
    i.set_viewing_angles([0., 30., 60., 90., 120., 150., 180.],
                         [0., 0., 0., 0., 0., 0., 0.])
    i.set_track_origin('basic')

# A simple model to check what happens when a source is moving towards dust and
# we observe both the source and the dust. If we observe the source such that
# 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

def setup_model(outdir, record_dir, outname, params, dust_file, wav_range, aperture,
                tsc=True, idl=False, plot=False, low_res=True, max_rCell=100,
                scale=1, radmc=False, mono_wave=None, norecord=False,
                dstar=200., dyn_cav=False, fix_params=None,
                power=2, mc_photons=1e6, im_photons=1e6, ellipsoid=False,
                TSC_dir='~/programs/misc/TSC/',
                IDL_path='/Applications/exelis/idl83/bin/idl', auto_disk=0.25,
                fast_plot=False, image_only=False, ulrich=False):
    """
    params = dictionary of the model parameters
    'alma' keyword is obsoleted
    outdir: The directory for storing Hyperion input files
    record_dir: The directory contains "model_list.txt" for recording parameters
    TSC_dir: Path the TSC-related IDL routines
    IDL_path: The IDL executable
    fast_plot: Do not plot the polar plot of the density because the rendering
               takes quite a lot of time.
    mono: monochromatic radiative transfer mode (need to specify the wavelength
          or a list of wavelength with 'mono_wave')
    image_only: only run for images
    """
    import numpy as np
    import astropy.constants as const
    import scipy as sci
    # to avoid X server error
    import matplotlib as mpl
    mpl.use('Agg')
    #
    import matplotlib.pyplot as plt
    import os
    from matplotlib.colors import LogNorm
    from scipy.integrate import nquad
    from hyperion.model import Model
    from record_hyperion import record_hyperion
    from pprint import pprint

    # Constants setup
    c         = const.c.cgs.value
    AU        = const.au.cgs.value     # Astronomical Unit       [cm]
    pc        = const.pc.cgs.value     # Parsec                  [cm]
    MS        = const.M_sun.cgs.value  # Solar mass              [g]
    LS        = const.L_sun.cgs.value  # Solar luminosity        [erg/s]
    RS        = const.R_sun.cgs.value  # Solar radius            [cm]
    G         = const.G.cgs.value      # Gravitational constant  [cm3/g/s^2]
    yr        = 60*60*24*365           # Years in seconds
    PI        = np.pi                  # PI constant
    sigma     = const.sigma_sb.cgs.value  # Stefan-Boltzmann constant
    mh        = const.m_p.cgs.value + const.m_e.cgs.value
    g2d       = 100.
    mmw       = 2.37                   # Kauffmann 2008

    m = Model()

    # min and max wavelength to compute (need to define them first for checking dust properties)
    wav_min, wav_max, wav_num = wav_range

    # Create dust properties
    # Hyperion needs nu, albedo, chi, g, p_lin_max
    from hyperion.dust import HenyeyGreensteinDust
    dust = dict()
    [dust['nu'], dust['albedo'], dust['chi'], dust['g']] = np.genfromtxt(dust_file).T
    d = HenyeyGreensteinDust(dust['nu'], dust['albedo'], dust['chi'], dust['g'], dust['g']*0)
    # dust sublimation option
    # dust sublimation temperture specified here
    T_sub = 1600.0
    d.set_sublimation_temperature('slow', temperature=T_sub)
    d.set_lte_emissivities(n_temp=3000,
                           temp_min=0.1,
                           temp_max=2000.)
    # if the min and/or max wavelength fall out of range
    if c/wav_min/1e-4 > dust['nu'].max():
        d.optical_properties.extrapolate_nu(dust['nu'].min(), c/wav_min/1e-4)
        print('minimum wavelength is out of dust model.  The dust model is extrapolated.')
    if c/wav_max/1e-4 < dust['nu'].min():
        d.optical_properties.extrapolate_nu(c/wav_max/1e-4, dust['nu'].max())
        print('maximum wavelength is out of dust model.  The dust model is extrapolated.')

    # try to solve the freq. problem
    d.optical_properties.extrapolate_nu(3.28e15, 5e15)
    #
    d.write(outdir+os.path.basename(dust_file).split('.')[0]+'.hdf5')
    d.plot(outdir+os.path.basename(dust_file).split('.')[0]+'.png')
    plt.clf()

    # Grids and Density

    # Grid Parameters
    nx        = 300L
    if low_res == True:
        nx    = 100L
    ny        = 400L
    nz        = 50L
    [nx, ny, nz] = [int(scale*nx), int(scale*ny), int(scale*nz)]

    # TSC model input setting
    dict_params = params
    # TSC model parameter
    cs        = dict_params['Cs']*1e5
    t         = dict_params['age']  # year
    omega     = dict_params['Omega0']
#.........这里部分代码省略.........