Simple Disk Envelope (SiDE) fit framework
4.6 Example of basic usage 87
# Set image parameters
impar = [{'npix':512,'wav':[1100.,3000.],'sizeau':11000,'incl':67.}]
• Model image options In visdata and impar have to be define wav which is the wavelength in microns, npix is the number of pixels on the rectangular images, sizeau is the diameter of the image in au and incl is the inclination angle of the source. In this framework the input original image size can be re-scaled to speed-up the model image time computing. This was implemented in the previous analysis of Chapter 3, where we required a large image size to model the observation, but only a small percent from extended emission contributed to the modelled uv-range.
To avoid do the raytracing on pixels that doesn’t contribute much to the overall visibilitysizeau_extis used. The new image conserves the original image resolution (in terms of au/px), the old image size is specified insizeau_ext andsizeau covers the new image size. The original image is stored in the center of the new image, while the extended image size is filled out with zero values. It is recommended to check the data-sets withGALARIO, using get_image_size, before to implementing this change. In case of a typical model using for example the data-sets of Chapter 3, the computation time is dominated by the raytracing step. To fully resolve the models on the same scale as the observations, large synthetic images are needed. Typical runtime of the thermal radiation transport is 30 – 40 s, forGALARIO runtime is 6 – 10 s.
• Model options In the paramfile is possible to setup the stellar properties, dust opacity and density profiles for different components. The disk model is given by a flaring density profile in which is possible to change parameters as the height of the disk, flaring index, density power law, inner disk radius among others. The envelope profile has four options: powerlaw, Tafalla2004, Ulrich1976andUlrich1976_2.
In Table 4.1 are the parameters that have to be setup for each envelope density pro-file. In the case of powerlaw and Tafalla2004density models the parameterr0Env refers to the radius at which rho0Envis defined, while in the case of theUlrich1976 density profile r0Env is the centrifugal radius.
The block of dust opacity contains the option to change the number of grain pop-ulations, ngpop. This version of the framework accepts up to three different dust populations, one for the disk and two for the envelope. If Ulrich1976_2 is setup, Envcuthas to be defined in theparamfile. Envcutis a limit in units of au where one dust population is in the inner part of Envcut and another dust population outside Envcut.
Additionally, if the model contains an outflow cavity then it is possible to setup dif-ferent cavity shapes with modeCav. The available cavity shapes are: Sheehan2017, edgecenter and gridcenter, see Fig. 4.2. Sheehan2017 is a outflow cavity model from Sheehan & Eisner (2017) characterized by:
z ą1 au`rζ, (4.2)
Table 4.1: Model setup parameters Parameter
Envelope model rho0Env r0Enva rTrunEnv prhoEnv ngpop
powerlaw 3 3 3 3 1
Tafalla2004 3 3 3 3 1
Ulrich1976 3 3 3 7 1
Ulrich1976_2b 3 3 3 7 2
(a)r0Env is the centrifugal radius in the Ulrich profile. (b) To setup 2 different dust populations in the Ulrich profile Envcut has to be define in units of au.
where 1 au is the height above the midplane where the cavity starts and ζ charac-terizes the shape and opening angle of the cavity. More details in Eisner (2012).
In the gridcenter option the cavity start at the center of the model grid (0,0,0), while in theedgecenter the cavity start at the outer radius of the disk following the equation:
arctan
ˆrsinθ´rdisk rcosθ
˙
ăθcav, (4.3)
where rdisk is the outer radius of the disk in au, and θcav is the open angle onf the cavity in degrees.
• Fitting arguments The next step is to define a dictionary that contains keyword arguments for the logarithm of the posterior probability function lnpostfn(). An example of kwargs is given by:
kwargs = {'dpc': 125., 'incl': 67., 'impar': impar, 'verbose': True, 'idisk':True, 'ienv':True, 'icav':False,
'cleanModel': True, 'binary': True, 'chi2_only':True, 'galario_check':False, 'time':True }
where dpc is the distace to the source in unit of parsec, impar collects the image parameters, verbose is an option to print a summary of the model parameters and the inclusion of the model profiles withidisk (disk profile),ienv(envelope profile), icav (outflow cavity) and islab (slab model). cleanModel is an option to delete the RADMC-3D model folder from disk after the posterior probability estimation, therefore the model files are not stored. If binary is True then RADMC3D will use binary format otherwise the file format is ascii. The use of binary files improves the computation speed and reduces disk space usage when models are kept.
The keyword chi2_onlycomputes and store the χ2 between the data and model. If it is True the synthetic visibility is computed but not stored (a zero value array is stored instead). Theχ2 is still computed and stored. Set this to True when running MCMC fitting in order to improve speed.
Figure 4.2: Three examples of 2D density maps for different shapes of outflow cavity. In gridcenter the outflow cavity starts at the center of the grid, edgecentercavity starts atrTrunEnv, the inner radius of the envelope density profile and the Sheehan2017is a cavity described in Sheehan & Eisner (2017).
Finally, galario_check and time, return the Nyquist criterion for computing the synthetic visibilities in the (u,v) locations provided by GALARIO and the function run-time information, respectively.
The following step is to choose the fitting parameters, the range of the space ex-ploration and the initial guess for the parameters. The fitted parameter names are listed in parname. These can be parameters from the model setup file and some special parameters as: source distance, inclination, position angle and offset in right ascension (RA) and declination (DEC). p_rangesare the limits for each parameter.
# Choose fitting parameters and initial values
parname = ['mdisk','rho0Env','gsmax_disk','gsmax_env']
p_ranges = [[-10., -2.], # log disk dust mass [solar mass]
[-23., -19.], # log envelope dust density [g/cm**3]
[-6., 0.], # log disk grain size [cm]
[-6., 0.]] # log envelope grain size [cm]
Additionally, the initial guess for the parameters is assigned to p0. The mode of exploration of the initial guess parameter can be defined through the lists: p_form, p_formpriorand p_sigma.
# initial guess for the parameters p0 = [-5, -20, -4., -4.]
p_form = ['log', 'lin', 'log', 'log']
p_formprior = ['uniform', 'normal', 'uniform', 'uniform']
p_sigma = [0., 0.5, 0., 0., 0., 0.]
Thep_formsets whether the prior probabilityp[i]is logarithmic (value = 10**p[i]) or linear (value = p[i]). p_formpriorsets the functional form of the prior proba-bility distribution. It should be set to normal oruniform for a Gaussian or rectan-gular distribution, respectively. p_sigma is the width of the Gaussian function.
• SiDE on computer clusters. SiDE has the option use_mpi where it is possible to use MPI pools instead of Python threads. This alternative is useful for running MCMC procedures on computer clusters using multiple nodes. To implement this, first it is necessary to add the runtime parameters as number of walkers (nwalkers) and the number of walker steps in the main run (nsteps) in the fitting Python script.
# Runtime parameters nwalkers = 240 nburnin = 0