Line RT through a homogeneous, static, MgII sphere#
In this second tutorial, we again compute RT of radiation emitted by a central source through a homogeneous sphere filled with gas (and no dust initially). But this time, we look at what happens to an emission line when propagated through the sphere.
Setup and imports#
RascasDir = '../../../' # path to the rascas directory (could be an absolute path).
atomic_data_dir = '%s/ions_parameters/'%(RascasDir) # path to the rascas/ions_parameters directory.
ExperimentDir = './' # where the RASCAS run will happen (and where param file will be written).
# use a parallel environment or not
useMPI = True
# rascas
import sys
sys.path.append("%s/py/"%(RascasDir))
import write_param_files as wpf
import jphot as jp
import lya_utils as lya
# various things
import os
import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
from astropy.constants import c
1. Definition of parameters#
1.1. Model parameters#
1.1.1. Gas properties#
# Physical properties of the sphere
ColumnDensity_cgs = 1e14 # Column density of scatterers from center to border of sphere. [cm^-2]
Temperature = 1e4 # temperature of the gas. [K]
TurbulentVelocity_kms = 10.0 # turbulent velocity of the gas. [km/s]
# Geometrical parameters
boxsize_cm = 1e20 # [cm]. This could be anything, really...
Radius_boxUnits = 0.4 # Radius of sphere in units of box-size
# (0.5 means the sphere extends all the way to the edge).
# We should also specify the number and nature of the scatterer(s) (i.e. of the absorption line(s))
nscatterer = 1
scatterer_names = 'MgII-2796'
1.1.2. Properties of the radiation source#
NB: radiation is emitted from a point source at rest at the center of the computational box.
# Definition of the spectral shape
spec_type = 'Gauss' # We emit a Gaussian emission line
spec_gauss_l0_Ang = 2796.35 # centered at spec_gauss_l0_Ang [Angstrom]
spec_gauss_sigma_kms = 20. # with a width of spec_gauss_sigma_kms [km/s]
# Number of photon packets for the Monte Carlo experiment
nPhotonPackets = 10000
1.1.3. Numerical parameters#
These define how we will generate an adaptive mesh etc.
refine_lmax = 7 # no need for refinement with a homogeneous medium ...
refine_err_grad_d = 0.1 # will trigger refinement on density fluctuations (if above 10% level)
1.1.4. Parameters for mocks (with peeling off)#
# Here we need to specify the properties of the mocks we wish to construct, in terms of aim (coords
# of target and direction of observation) and in terms of observation (photometry, spectroscopy,
# imaging, 3D spectroscopy). This is done through parameters which need to be written to a file.
# We start by setting up one pointing and defining the observations for it.
target_position = np.array([0.5,0.5,0.5]) # aim for the center of the box (all in box units here... )
kobs = np.array([0,1,0]) # observe along the y direction (make sure k is normalised)
# now, let's add a spectroscopic observation
npix = 200 # the nb of wavelength bins
aperture = 1. # the radius of the circular aperture within which photon packets are collected (in box-size units)
lambda_min = spec_gauss_l0_Ang - 1. # min wavelength (Angstrom)
lambda_max = spec_gauss_l0_Ang + 1. # max wavelength (Angstrom)
spec = wpf.specParams(npix,aperture,lambda_min,lambda_max)
# let's make an image as well ...
npix = 200 # nb of pixels on a side
size = 1. # size of the image (side)
image = wpf.imageParams(npix,size)
# define the object containing all the above information
p = wpf.pointing(kobs,target_position,spec=spec,image=image)
# and write it into some file
mock_params_file = '%s/mock_param_file.txt'%ExperimentDir
f = open(mock_params_file,'w')
p.write(f)
f.close()
1.2. Setting up the parameter file#
Parameters = [] # initialise a list of parameter sections
### GEOMETRY
section = 'GenerateAMRmodel'
params = {
'DomDumpDir' : ExperimentDir,
'comput_dom_type' : 'sphere',
'comput_dom_pos' : '0.5,0.5,0.5',
'comput_dom_rsp' : '%f'%(Radius_boxUnits)
}
Parameters.append(wpf.param_section(section,params))
### MESH
section = 'mesh'
params = {
'verbose' : 'F',
'refine_lmax' : '%i'%(refine_lmax),
'refine_err_grad_d' : '%.2f'%(refine_err_grad_d),
'refine_err_grad_v' : '-1',
'refine_dv_over_vth' : 'F',
}
Parameters.append(wpf.param_section(section,params))
### GAS COMPOSITION
section = 'gas_composition'
params = {
'nscatterer' : '%i'%(nscatterer),
'scatterer_names' : scatterer_names,
'atomic_data_dir' : atomic_data_dir,
'ignoreDust': 'T'
}
Parameters.append(wpf.param_section(section,params))
### DUST (unused initially ... )
section = 'dust'
params = {
'albedo':0.32,
'g_dust':0.73,
'dust_model':'SMC'
}
Parameters.append(wpf.param_section(section,params))
### IDEALISED MODEL
section = 'IdealisedModel'
params ={
'ColumnDensity_cgs': '%e'%(ColumnDensity_cgs),
'boxsize_cm': '%e'%(boxsize_cm),
'Radius_boxUnits': '%e'%(Radius_boxUnits),
'Temperature':'%e'%(Temperature),
'TurbulentVelocity_kms':'%f'%(TurbulentVelocity_kms)
}
Parameters.append(wpf.param_section(section,params))
### SOURCE MODEL
section = 'PhotonsFromSourceModel'
params ={
'outputfile' : '%s/ppic.dat'%ExperimentDir,
'source_pos' : '0.5,0.5,0.5',
'source_vel' : '0,0,0', # cm/s
'source_type': 'pointlike',
'spec_type' : spec_type,
'spec_gauss_l0_Ang' : '%f'%(spec_gauss_l0_Ang),
'spec_gauss_sigma_kms' : '%f'%(spec_gauss_sigma_kms),
'nPhotonPackets' : '%i'%(nPhotonPackets)
}
Parameters.append(wpf.param_section(section,params))
# RASCAS
if useMPI:
section = 'rascas'
else:
section = 'rascas-serial'
nbundle = min(nPhotonPackets/100,1000)
params ={
'verbose':'T',
'DomDumpDir':ExperimentDir,
'PhotonICFile':'%s/ppic.dat'%ExperimentDir,
'fileout':'%s/result.dat'%ExperimentDir,
'nbundle':'%i'%nbundle
}
Parameters.append(wpf.param_section(section,params))
## WORKER (if MPI)
if useMPI:
section='worker'
params={'verbose':'T'}
Parameters.append(wpf.param_section(section,params))
### MOCK -- peeling off
section='mock'
params = {
'nDirections':'1',
'mock_parameter_file':'%s'%mock_params_file,
'mock_outputfilename':'%s/mock'%ExperimentDir
}
Parameters.append(wpf.param_section(section,params))
#### write parameter list to a file
wpf.write_parameter_file('%s/params.cfg'%(ExperimentDir),Parameters)
2. Run the codes#
2.1. Create mesh and visualize it#
Execute the following cell in order run the executable GenerateAMRmodel which generates the adaptive mesh containing the model. This will generate a few files in the ExperimentDir, including a log file GenerateAMRmodel.log, the .mesh file (which has the data) and some .dom files which describe different domains.
# run CreateDomDumpAMRModel
cmd = "./GenerateAMRmodel params.cfg > GenerateAMRmodel.log"
error = os.system(cmd)
if error != 0:
print('Something went wrong...')
print('-> check file GenerateAMRmodel.log or message below for more info')
2.2. Create the photon packets#
Execute the following cell in order run the executable PhotonsFromSourceModel which generates the initial conditions for photon packets. This will generate a couple files in the ExperimentDir: a log file PhotonsFromSourceModel.log and the ppic.dat file (which has the data).
# run PhotonsFromSourceModel
cmd = "./PhotonsFromSourceModel params.cfg > PhotonsFromSourceModel.log"
error = os.system(cmd)
if error != 0:
print('Something went wrong...')
print('-> check file PhotonsFromSourceModel.log or message below for more info')
2.3. Run the radiative transfer#
Execute the following cell in order run radiative transfer computation with rascas (or rascas-serial if you do not have MPI). This may take time depending on the experiment. If it is too slow, reduce the number of photon packets above.
# run RASCAS
if useMPI:
cmd = "time mpirun ./rascas params.cfg > rascas.log"
else:
cmd = "time ./rascas-serial params.cfg > rascas.log"
error = os.system(cmd)
if error != 0:
print('Something went wrong...')
print('-> check file rascas.log or message below for more info')
real 0m0.557s
user 0m2.991s
sys 0m0.238s
3. Plot the results#
3.1. Results from photon packets#
# load the results of the RASCAS experiment
p = jp.photonlist('ppic.dat','result.dat')
3.1.1. Spectrum (angle-averaged and spatially integrated)#
# get the observed spectrum
lobs,sobs = p.spectrum(frame='obs',nbins=200)
# get the intrinsic spectrum
lint,sint = p.spectrum(frame='ic',nbins=200)
# normalise spectra to the median value of the intrinsic spectrum (i.e. the continuum level)
continuum = np.median(sint)
sobs = sobs / continuum
sint = sint / continuum
# plot spectra
plt.plot(lobs,sobs,label='Escaping spectrum',lw=2,ls='-',c='k')
plt.plot(lint,sint,alpha=0.6,label='Intrinsic spectrum',lw=2,ls='-',c='firebrick')
plt.legend(fontsize=14)
plt.xlabel('wavelength (Angstrom)',fontsize=14)
plt.ylabel(r'$F_\lambda\ /\ F_{cont}$',fontsize=14)
plt.show()
3.1.2. Surface-brightness map and profile (angle-averaged)#
# compute projected distance of each photon packet to the source, in the plane perpendicular
# to the photon packet's propagation direction.
# (NB:Source is located at 0.5,0.5,0.5, i.e. center of the box.)
rp = p.projected_radius(xc=0.5,yc=0.5,zc=0.5)
/Users/leo/tmp/new_rascas_doc/test-new-doc/tutorials/IdealisedModels/tutorial_2/../../..//py/jphot.py:275: RuntimeWarning: invalid value encountered in divide
ux = np.nan_to_num((self.x-xc)/unorm)
/Users/leo/tmp/new_rascas_doc/test-new-doc/tutorials/IdealisedModels/tutorial_2/../../..//py/jphot.py:276: RuntimeWarning: invalid value encountered in divide
uy = np.nan_to_num((self.y-yc)/unorm)
/Users/leo/tmp/new_rascas_doc/test-new-doc/tutorials/IdealisedModels/tutorial_2/../../..//py/jphot.py:277: RuntimeWarning: invalid value encountered in divide
uz = np.nan_to_num((self.z-zc)/unorm)
# compute SB profile in arbitrary units (number of photon packets per annulus area, in box units)
h,edges = np.histogram(rp,bins=100)
bin_centers = 0.5*(edges[:-1]+edges[1:])
dr = bin_centers[1]-bin_centers[0]
h = h / (2.*np.pi*bin_centers*dr) # SB profile, arbitrary units
# plot SB profile
plt.plot(bin_centers,h, lw=2,ls='-',c='k')
plt.yscale('log')
plt.xlabel('radius (box units)',fontsize=14)
plt.ylabel('Surface Brightness (arbitrary units)',fontsize=14)
plt.show()
3.1.3. Spatially resolved spectrum (angle-averaged)#
# select and extract photon packets escaping from the inner region (point source -> take really small radius)
iiIN = np.where(rp < 1e-5)
pIN = p.extract_sample(iiIN)
lIN,sIN = pIN.spectrum(frame='obs',nbins=200,lmin=spec_gauss_l0_Ang-1.,lmax=spec_gauss_l0_Ang+1)
sIN /= continuum
plt.plot(lIN,sIN,alpha=0.6,label='Escaping spectrum from inner region',lw=3,ls='-',c='r')
# select and extract photon packets escaping from the outter region
iiIN = np.where(rp >= 1e-5)
pIN = p.extract_sample(iiIN)
lIN,sIN = pIN.spectrum(frame='obs',nbins=200,lmin=spec_gauss_l0_Ang-1.,lmax=spec_gauss_l0_Ang+1)
sIN /= continuum
plt.plot(lIN,sIN,alpha=0.6,label='Escaping spectrum from outter region',lw=3,ls='-',c='b')
# plot total as well
plt.plot(lobs,sobs,label='Total spectrum',lw=1,ls='--',c='k')
plt.xlabel('wavelength (Angstrom)',fontsize=14)
plt.ylabel(r'$F_\lambda\ /\ F_{cont}$',fontsize=14)
plt.legend()
plt.show()
3.2. Peeling off results#
from mocks import mockobs
m = mockobs('./','mock','./ppic.dat',load_spectrum=True,load_image=True)
plt.figure()
plt.plot(m.spec_lbda_Angstrom,m.spec/continuum,lw=4,ls='--',c='deepskyblue')
plt.plot(lobs,sobs,label='Escaping spectrum',lw=2,ls='-',c='k')
plt.plot(lint,sint,alpha=0.6,label='Intrinsic spectrum',lw=2,ls='-',c='firebrick')
plt.legend(fontsize=14)
plt.xlabel('wavelength (Angstrom)',fontsize=14)
plt.ylabel(r'$F_\lambda\ /\ F_{cont}$',fontsize=14)
plt.show()
m.show_image(smooth=False,vmin=1e-18,vmax=1e-14)
