Python Model.set_raytracing方法代码示例

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

#.........这里部分代码省略.........
                        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])
    nlam    = len(lam)

    # Create camera wavelength points
    n12     = 70.0
    n23     = 70.0
    n34     = 70.0
    n45     = 50.0
    n56     = 50.0
    
    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_cam = np.concatenate([lam12,lam23,lam34,lam45,lam56])
    n_lam_cam = len(lam_cam)

    # Radiative transfer setting

    # number of photons for temp and image
    m.set_raytracing(True)
    m.set_monochromatic(True, wavelengths=[3.6, 4.5, 5.8, 8.0, 24, 70, 100, 160, 250, 350, 500])
    m.set_n_photons(initial=1000000, imaging_sources=1000000, imaging_dust=1000000,raytracing_sources=1000000, raytracing_dust=1000000)
    # imaging=100000, raytracing_sources=100000, raytracing_dust=100000
    # 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
    # m.set_forced_first_scattering(forced_first_scattering=True)
    # Setting up images and SEDs
    image = m.add_peeled_images()
    # image.set_wavelength_range(300, 2.0, 670.0)
    # use the index of wavelength array used by the monochromatic radiative transfer
    image.set_wavelength_index_range(2,12)
    # 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)

    # 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+outname+'.rtin')

# 需要导入模块: from hyperion.model import Model [as 别名]
# 或者: from hyperion.model.Model import set_raytracing [as 别名]
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)
image.set_image_limits(-1.5 * pc, 1.5 * pc, -1.5 * pc, 1.5 * pc)

# Add a fly-around at 500 microns
image = m.add_peeled_images(sed=False, image=True)
image.set_wavelength_range(1, 499., 501.)
image.set_viewing_angles(np.repeat(45., 36), np.linspace(5., 355., 36))
image.set_image_size(200, 200)
image.set_image_limits(-1.5 * pc, 1.5 * pc, -1.5 * pc, 1.5 * pc)

# Set runtime parameters
m.set_n_initial_iterations(5)
m.set_raytracing(True)
m.set_n_photons(initial=1e5, imaging=1e6,
                raytracing_sources=1e5, raytracing_dust=1e5)

# Write out input file
m.write('tutorial_model.rtin')

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

#.........这里部分代码省略.........
        except IOError:
            print("ERROR: File '", file, "' cannot be found. \nERROR: This file, containing the specific energy density, has to be computed first via calling hyperion.")
            exit(2)
        
		#
		# To compute total photon numbers:
		#
        grid_Nw = len(model.grid.gw[0,0,:])
        grid_Nz = len(model.grid.gw[0,:,0])
        grid_Np = len(model.grid.gw[:,0,0])
        grid_N = grid_Nw * grid_Nz * grid_Np
        if(cli.verbose):
            print("Grid setup:")
            print(" grid_Nw =",grid_Nw)
            print(" grid_Nz =",grid_Nz)
            print(" grid_Np =",grid_Np)
            print("Radiation setup:")
            print(" photons_temperature / cell =", cli.photons_temperature)
            print(" photons_temperature total  =", grid_N * cli.photons_temperature)
            print(" photons_raytracing / cell  =", cli.photons_raytracing)
            print(" photons_raytracing total   =", grid_N * cli.photons_raytracing)
            print(" photons_imaging / cell     =", cli.photons_imaging)
            print(" photons_imaging total      =", grid_N * cli.photons_imaging)
        
        file = filename(cli, "")
        file += ".rtin"


    ##
    ## Temperature, Images, and SEDs:
    ##
    if(cli.mode == "temperature"):
    
        model.set_raytracing(True)
        model.set_n_photons(
            initial            = grid_N * cli.photons_temperature,
            raytracing_sources = grid_N * cli.photons_raytracing,
            raytracing_dust    = grid_N * cli.photons_raytracing,
            imaging            = grid_N * cli.photons_imaging
        )
        
    elif(cli.mode == "images"):
        
        model.set_n_initial_iterations(0)
        model.set_raytracing(True)
        # old setup: model.set_monochromatic(True, wavelengths=[0.4, 1.0, 10.0, 100.0, 500.0])
        model.set_monochromatic(True, wavelengths=[0.45483, 1.2520, 26.114, 242.29])
        model.set_n_photons(
            raytracing_sources = grid_N * cli.photons_raytracing,
            raytracing_dust    = grid_N * cli.photons_raytracing,
            imaging_sources    = grid_N * cli.photons_imaging,
            imaging_dust       = grid_N * cli.photons_imaging
        )
    
        # group = 0
        image1 = model.add_peeled_images(sed=False, image=True)
        image1.set_image_size(501, 501)
        image1.set_image_limits(-12500.0*pc, +12500.0*pc, -12500.0*pc, +12500.0*pc)
        image1.set_viewing_angles([30],[0])
        image1.set_uncertainties(True)
        image1.set_output_bytes(8)
        image1.set_track_origin('basic')
    
        # group = 1
        image2 = model.add_peeled_images(sed=False, image=True)
        image2.set_image_size(501, 501)

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

#.........这里部分代码省略.........
                                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'],
                            absolute=dict_params['absolute'],
                            relative=dict_params['relative'])
    m.set_mrw(True)   # Gamma = 1 by default

    # Setting up images and SEDs
    if not image_only:
        # SED setting
        # Infinite aperture
        syn_inf = m.add_peeled_images(image=False)
        # use the index of wavelength array used by the monochromatic radiative transfer
        if mono_wave == None:

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

#.........这里部分代码省略.........
    lambda6 = 1000.0
    n01     = 10.0
    n12     = 20.0
    n23     = 50.0

    lam01   = lambda0 * (lambda1/lambda0)**(np.arange(n01)/n01)
    lam12   = lambda1 * (lambda2/lambda1)**(np.arange(n12)/n12)
    lam23   = lambda2 * (lambda6/lambda2)**(np.arange(n23+1)/n23)

    lam      = np.concatenate([lam01,lam12,lam23])
    nlam    = len(lam)

    # Create camera wavelength points
    n12     = 70.0
    n23     = 70.0
    n34     = 70.0
    n45     = 50.0
    n56     = 50.0
    
    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_cam = np.concatenate([lam12,lam23,lam34,lam45,lam56])
    n_lam_cam = len(lam_cam)

    # Radiative transfer setting

    # number of photons for temp and image
    lam_list = lam.tolist()
    # print lam_list
    m.set_raytracing(True)
    # option of using more photons for imaging
    if better_im == False:
        im_photon = 1e6
    else:
        im_photon = 5e7

    if mono == True:
        # Monechromatic radiative transfer setting
        m.set_monochromatic(True, wavelengths=lam_list)
        m.set_n_photons(initial=1000000, imaging_sources=im_photon, imaging_dust=im_photon,raytracing_sources=1000000, raytracing_dust=1000000)
    else:
        # regular wavelength grid setting
        m.set_n_photons(initial=1000000, imaging=im_photon,raytracing_sources=1000000, raytracing_dust=1000000)    
    # number of iteration to compute dust specific energy (temperature)
    m.set_n_initial_iterations(20)
    # m.set_convergence(True, percentile=95., absolute=1.5, relative=1.02)
    m.set_convergence(True, percentile=dict_params['percentile'], absolute=dict_params['absolute'], relative=dict_params['relative'])
    m.set_mrw(True)   # Gamma = 1 by default
    # m.set_forced_first_scattering(forced_first_scattering=True)

    # Setting up images and SEDs
    # SED setting

    # Infinite aperture
    syn_inf = m.add_peeled_images(image=False)
    # use the index of wavelength array used by the monochromatic radiative transfer
    if mono == False:
        syn_inf.set_wavelength_range(1400, 2.0, 1400.0)
    syn_inf.set_viewing_angles([dict_params['view_angle']], [0.0])
    syn_inf.set_uncertainties(True)
    syn_inf.set_output_bytes(8)