Using Wakeflow with MCFOST
This tutorial demonstrates the use of wakeflow models with the radiative transfer code MCFOST to create synthetic observations. It is assumed that you have read the Quickstart Tutorial.
Prerequisites
Firstly, we will need working installations of both `MCFOST <https://github.com/cpinte/mcfost>`__ and `pymcfost <https://github.com/cpinte/pymcfost>`__. Please see their linked Github pages and documentation for instructions on installing. Also ensure that you read their usage and citation guidelines.
Let’s check that our MCFOST installation is up to date and working by running the following shell command:
[1]:
!mcfost -v
You are running MCFOST 3.0.44
Git SHA = 8686df4719b1cadb7b1d82257fc576ec8072e0c0
Binary compiled the Jul 18 2022 at 11:56:15
with INTEL compiler version 2021
Checking last version ...
MCFOST is up-to-date
And checking the pymcfost installation by importing it:
[2]:
import pymcfost
WARNING: mpl_scatter_density is not present
WARNING: progressbar is not present
Provided both of these installations are working, and wakeflow is installed, you should be ready to proceed with the rest of the tutorial.
Required Parameters
In order for a wakeflow model to be compatible with MCFOST, it must be run on the correct grid geometry. wakeflow achieves this by calling MCFOST to generate the geometry. We therefore have to make sure that wakeflow is configured correctly so that the resultant model is readable by MCFOST.
The crucial parameters that you must specify are:
cw_rotation=False(this is default)grid_type='mcfost'dimensionless=False(this is default)write_FITS=True
The mcfost grid type is cylindrical, with logarithmically placed radii and a flared \(z\) coordinate. See the MCFOST docs for more information. We also have to specify for the disk to rotate anticlockwise, if you wish for your disk model to rotate clockwise you may achieve this through transformations with MCFOST later.
Running Wakeflow
Let us generate a wakeflow model identical to that in the Quickstart Tutorial, except with a planet mass of \(4.0 \, \mathrm{M_J}\) and using the required parameters for MCFOST.
[3]:
from wakeflow import WakeflowModel
hd163_model = WakeflowModel()
hd163_model.configure(
name = "mcfost_tutorial",
system = "HD_163296",
m_star = 1.9,
m_planet = 4.0,
r_outer = 500,
r_planet = 256,
r_ref = 256,
q = 0.35,
p = 2.15,
hr = 0.09,
cw_rotation = False,
grid_type = "mcfost",
n_r = 200,
n_phi = 280,
n_z = 40,
show_midplane_plots = False,
write_FITS = True
)
hd163_model.run()
Model initialised.
Model configured.
* Performing checks on model parameters:
M_thermal = 0.967 M_Jup
M_planet = 4.135 M_th
WARNING: M_planet > M_thermal. This may break the solution.
Parameters Ok - continuing
Generating MCFOST grid data...
You are running MCFOST 3.0.44
Git SHA = 8686df4719b1cadb7b1d82257fc576ec8072e0c0
WARNING: polarization is turned off in ray-traced SEDs
it can be turned back on with -rt2
Input file read successfully
Computation of disk structure
Creating directory ././data_disk
Parallelized code on 10 processors
Tue 23 Aug 2022 22:27:27 AEST
Forcing 3D mode
Using ray-tracing method 1
Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc
s/mcfost_utils/Dust/Draine_Si_sUV.dat
Number of regions detected: 1
zone 1 --> region= 1 : R=100.00 to 500.00 AU
Using 4.480000 million cells
Total gas mass in model: 0.1000000 Msun
Total dust mass in model: 1.0000000E-03 Msun
Using scattering method 2
Writing disk structure files in data_disk ...
Exiting
Done
* Creating 4.0 Mj model:
Generating unperturbed background disk
Extracting linear perturbations nearby planet
Propagating outer wake...
Completed in 0.17 s
Propagating inner wake...
Completed in 0.17 s
* Mapping to physical coords: 100%|██████████| 200/200 [00:05<00:00, 37.37it/s]
* Saving results:
Perturbations saved to HD_163296/mcfost_tutorial/4.0Mj
Total saved to HD_163296/mcfost_tutorial/4.0Mj
Saved FITS file formatted for MCFOST.
* Done!
Lets see what wakeflow has made for us:
[4]:
!ls HD_163296/mcfost_tutorial
4.0Mj mcfost_tutorial_config.yaml
mcfost
So far, this is the same as usual, except there is a new directory called mcfost. This directory contains the grid geometry data. Looking in our results directory
[5]:
!ls HD_163296/mcfost_tutorial/4.0Mj
delta_PHI.npy delta_v_r.npy total_Z.npy vr_z0.pdf
delta_R.npy mcfost.para total_rho.npy wakeflow_model.fits
delta_Z.npy rho_z0.pdf total_v_phi.npy
delta_rho.npy total_PHI.npy total_v_r.npy
delta_v_phi.npy total_R.npy vphi_z0.pdf
We see two files of interest for us. mcfost.para and wakeflow_model.fits. These are the two files we need to run MCFOST. The .fits file contains the results from the wakeflow model, while the .para file is the MCFOST parameter file. Most likely you will need to change some of the values in the .para file for your purposes (see MCFOST docs), but you should NOT change anything in the Grid Geometry or Density Structure settings, as wakeflow has already
chosen the appropriate values.
Running MCFOST
Once you have configured the .para file to your liking, we are ready to run MCFOST. First we move to the results directory:
[6]:
%cd HD_163296/mcfost_tutorial/4.0Mj
/Users/tomhilder/Documents/Research/Protoplanetary_Discs/wakeflow/docs/tutorials/HD_163296/mcfost_tutorial/4.0Mj
We can now run MCFOST on our model easily using the following:
[7]:
!mcfost mcfost.para -df wakeflow_model.fits -mol
You are running MCFOST 3.0.44
Git SHA = 8686df4719b1cadb7b1d82257fc576ec8072e0c0
WARNING: polarization is turned off in ray-traced SEDs
it can be turned back on with -rt2
Input file read successfully
Thermal equilibrium calculation
Temperature calculation under LTE approximation
Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc
s/mcfost_utils/Molecules/co@xpol.dat
Molecular line transfer under LTE approximation
Parallelized code on 10 processors
Tue 23 Aug 2022 22:28:11 AEST
Creating directory ././data_th
Creating directory ././data_CO
Forcing 3D mode
Using ray-tracing method 1
Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc
s/mcfost_utils/Dust/Draine_Si_sUV.dat
Number of regions detected: 1
zone 1 --> region= 1 : R=100.00 to 500.00 AU
Using 4.480000 million cells
Reading density file : wakeflow_model.fits
read_gas_density = 1
read_gas_velocity = 2
Reading dust density ...
No grain size found
The density file only has positive z, making it symmetric
Dust density range: 8.2326818E-10 11.77502
Using grain size distribution from parameter file
Reading gas density ...
Updating gas-to-dust ratio to 100.0000
Gas density range: 8.2326818E-10 11.77502
Reading gas velocity ...
Velocity field is in cylindrical coordinates.
Velocity range: -842.0961 4076.110
Constant spatial distribution
Total gas mass in model: 0.1000000 Msun
Total dust mass in model: 1.0000000E-03 Msun
Using scattering method 2
Trying to find appropriate stellar spectra ...
Star # 1 --> lte9200-4.0.NextGen.fits.gz
Done
Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc
s/mcfost_utils/Stellar_Spectra/lte9200-4.0.NextGen.fits.gz
Computing dust properties ... Done
lambda = 0.871355626679430
Integ tau dans plan eq. = 81.53822
Column density (g/cm�) = 0.1843978
Integ tau (i =75.0 deg) = 0.10330E+01
Column density (g/cm�) = 2.3361719E-03
No dark zone
Initializing thermal properties ... Done
Initialization complete in 45.00s
Computing temperature structure ...
100% |==================================================|
Max. temperature = 96.65446
Temperature calculation complete in 0h 1m 12s
Source fct time 0.0000000E+00 s
RT time 0.0000000E+00 s
-------------------------------------------------------
Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc
s/mcfost_utils/Molecules/co@xpol.dat
-------------------------------------------------------
CO molecular file read successfully
WARNING : memory size if lots of pixels
Setting constant abundance
Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc
s/mcfost_utils/Dust/Draine_Si_sUV.dat
Computing dust properties for 40 wavelength
Trying to find appropriate stellar spectra ...
Star # 1 --> lte9200-4.0.NextGen.fits.gz
Done
Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc
s/mcfost_utils/Stellar_Spectra/lte9200-4.0.NextGen.fits.gz
-------------------------------
Transition J= 2 - 1
tau_mol = 1157.638
tau_dust= 0.4059193
Vertical Tau_mol=1 (for r=100 au) at z= 13.05515 au
Processing complete in 0h 2m 51s
CPU time used 0h 10m 41s
[8]:
!ls
_dust_prop_th.tmp delta_Z.npy rho_z0.pdf total_v_phi.npy
data_CO delta_rho.npy total_PHI.npy total_v_r.npy
data_th delta_v_phi.npy total_R.npy vphi_z0.pdf
delta_PHI.npy delta_v_r.npy total_Z.npy vr_z0.pdf
delta_R.npy mcfost.para total_rho.npy wakeflow_model.fits
We see that MCFOST has generated two files from our configuration, data_CO and data_th. data_CO contains the synthetic line emission cube, that we can go on to read with astropy or plot with pymcfost, and is beyond the scope of this tutorial.
[10]:
%cd ../../..
/Users/tomhilder/Documents/Research/Protoplanetary_Discs/wakeflow/docs/tutorials
Extra Integrations
You can also ask wakeflow to run MCFOST automatically, instead of doing it manually as above. There are extra parameters that you can give wakeflow that will be written to the .para file, they are
inclinationPAPAprun_mcfosttemp_stardistancev_maxn_v
For example, we could generate models of HD 163296 with a range of planet masses, and run MCFOST on each one with only a few lines of code:
[15]:
hd163_model_2 = WakeflowModel()
hd163_model_2.configure(
name = "mcfost_tutorial_2",
system = "HD_163296",
m_star = 1.9,
m_planet = [1.0, 2.0, 3.0],
r_outer = 500,
r_planet = 256,
r_ref = 256,
q = 0.35,
p = 2.15,
hr = 0.09,
cw_rotation = False,
grid_type = "mcfost",
n_r = 200,
n_phi = 280,
n_z = 40,
show_midplane_plots = False,
write_FITS = True,
run_mcfost = True,
inclination = 225,
PA = 45,
PAp = 125,
temp_star = 9250,
distance = 101.5,
v_max = 4.0,
n_v = 36
)
hd163_model_2.run()
Model initialised.
Model configured.
* Performing checks on model parameters:
M_thermal = 0.967 M_Jup
M_planets = [1.034 2.068 3.101] M_th
At least one planet mass exceeds thermal mass. This may break the solution.
Parameters Ok - continuing
Generating MCFOST grid data...
You are running MCFOST 3.0.44
Git SHA = 8686df4719b1cadb7b1d82257fc576ec8072e0c0
WARNING: polarization is turned off in ray-traced SEDs
it can be turned back on with -rt2
Input file read successfully
Computation of disk structure
Creating directory ././data_disk
Parallelized code on 10 processors
Tue 23 Aug 2022 22:33:19 AEST
Forcing 3D mode
Using ray-tracing method 1
Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc
s/mcfost_utils/Dust/Draine_Si_sUV.dat
Number of regions detected: 1
zone 1 --> region= 1 : R=100.00 to 500.00 AU
Using 4.480000 million cells
Total gas mass in model: 0.1000000 Msun
Total dust mass in model: 1.0000000E-03 Msun
Using scattering method 2
Writing disk structure files in data_disk ...
Exiting
Done
* Creating 1.0 Mj model:
Generating unperturbed background disk
Extracting linear perturbations nearby planet
Propagating outer wake...
Completed in 0.68 s
Propagating inner wake...
Completed in 0.71 s
* Mapping to physical coords: 100%|██████████| 200/200 [00:05<00:00, 33.71it/s]
* Saving results:
Perturbations saved to HD_163296/mcfost_tutorial_2/1.0Mj
Total saved to HD_163296/mcfost_tutorial_2/1.0Mj
Saved FITS file formatted for MCFOST.
* Creating 2.0 Mj model:
Generating unperturbed background disk
Extracting linear perturbations nearby planet
Propagating outer wake...
Completed in 0.33 s
Propagating inner wake...
Completed in 0.35 s
* Mapping to physical coords: 100%|██████████| 200/200 [00:05<00:00, 35.19it/s]
* Saving results:
Perturbations saved to HD_163296/mcfost_tutorial_2/2.0Mj
Total saved to HD_163296/mcfost_tutorial_2/2.0Mj
Saved FITS file formatted for MCFOST.
* Creating 3.0 Mj model:
Generating unperturbed background disk
Extracting linear perturbations nearby planet
Propagating outer wake...
Completed in 0.22 s
Propagating inner wake...
Completed in 0.23 s
* Mapping to physical coords: 100%|██████████| 200/200 [00:05<00:00, 36.14it/s]
* Saving results:
Perturbations saved to HD_163296/mcfost_tutorial_2/3.0Mj
Total saved to HD_163296/mcfost_tutorial_2/3.0Mj
Saved FITS file formatted for MCFOST.
* Running MCFOST on each model:
1.0 Mj...
2.0 Mj...
3.0 Mj...
* Done!