ProDiMo-MCFOST Interfaces¶
MCFOST to ProDiMo¶
WARNING: This document is likely outdated. Especially on the MCFOST side a lot has changed. It is currently unclear if the interface is still working.
ProDiMo can run in tandem with the code MCFOST (multidimensional Monte-Carlo radiative transfer code). You need to know how to run MCFOST, the ProDiMo developers do not maintain nor support the code MCFOST, please enquire Christophe Pinte. See also the MCFOST documentation.
The idea is to first model the dust by running MCFOST with the -prodimo option, which creates a interface data file for ProDiMo, and then to model the gas by running ProDiMo on top of that.
19/08/2013: The new MCFOST-ProDiMo interface transfers both the actual GAS density structure and the actual DUST density structure. Therefore the user has to re-run MCFOST each time that he/she changes the dust-to-gas mass ratio. The dust-to-gas is now a parameter in MCFOST that can be ProDiMo-zone dependent.
Step 1¶
When you run MCFOST in this modus, make sure to use the ProDiMo.lambda wavelength file in /mcfost/utils/Lambda/, and set the flag to use this (i.e. not the default MCFOST wavelength range). The corresponding section in the MCFOST input file should look like this:
#Grid geometry and size
1 1 = cylindrical, 2 = spherical
100 80 1 20 n_rad, nz (or n_theta), n_az, n_rad_in
#Wavelength
70 0.1 3000.0 n_lambda, min wl, max wl
T T F compute temperature?, compute sed?, use default wavelength grid ?
ProDiMo.lambda. wavelength file (if previous parameter is F)
F F separation of different contributions?, stokes parameters?You can also use ProDiMo.lambda_UV1_3, ProDiMo.lambda_UV3_9, or ProDiMo.lambda_UV6_15, attached to this document, but then you need to link it to ProDiMo.lambda, because the -prodimo switch overrules your settings. It must be one of these grid files, otherwise ProDiMo will reject and stop, because we need to make sure that the lambda boundary values match with the definitions of the UV bands in ProDiMo.
The ProDiMo.lambda_UV1_3 can be safely used for a powerlaw-UV, however, if you have a steep or spiky, line-dominated UV, I'd recommend to use ProDiMo.lambda_UV6_15. Note that the file mcfost/utils/Lambda/ProDiMo.lambda has recently been changed. It should look like this:
0.0942359
0.1006141
0.1074241
0.1168222
0.1293989
0.1433296
0.1587600
0.1758515
0.1947831
0.2267423
0.2773893
0.3393494
0.4151493
0.5078805
0.6213250
0.7601094
0.9298938
1.1376028
1.3917073
1.7025708
2.0828713
2.5481190
3.1172884
3.8135922
4.4361089
4.9066020
5.4269956
6.0025820
6.6392150
7.3433693
8.1222060
8.9836461
9.9364504
10.9903090
12.1559397
13.4451970
14.8711928
16.4484296
18.1929478
20.1224894
22.2566779
24.6172182
27.2281171
30.1159275
33.3100188
36.8428749
40.7504252
45.0724099
52.4298639
64.1410401
78.4681233
95.9954246
117.4377716
143.6696620
175.7609284
215.0203705
263.0491324
321.8060034
393.6873042
481.6246182
589.2043516
720.8140008
881.8210904
1078.7920804
1319.7601706
1614.5529240
1975.1930712
2416.3888407
2956.1338153
3616.4407760The number of photon packages should be quite large for the first MCFOST-phase, because we need smooth Tdust temperatures. I recommend
#Number of photon packages
128 nbr_parallel_loop
50000 nbr_photons_eq_th : T computationNext, when you run MCFOST, make sure not to forget the -prodimo switch
./mcfost my_disc.para -prodimoThis will generate a file called forProDiMo.fits.gz in the data_th directory.
'''NOTE: You cannot combine -prodimo with -rt ! The ray-tracing SED results from MCFOST are incorrect if you use -prodimo. If you need -rt, you must run MCFOST twice, once with -rt to get the proper ray-traced SED, and once with -prodimo to get the proper forProDiMo.fits.gz.
Step 2¶
Go to the directory where you want to run ProDiMo, and copy forProDiMo.fits.gz here. Also copy sed2.fits.gz to that directory if you want to compare MCFOST's SED with the one calculated by ProDiMo based on the passed temperatures, mean intensities, and opacities (good safty-check).
Next, ensure that all the required *.in files needed by ProDiMo are present, typically Elements.in, LineTransferList.in, Parameter.in, Reactions.in and Species.in. At the moment, I would recommend to use the 3 input files provided in ~/ProDiMo/input_examples/Xrays+PAHs/developer.
In the Parameter.in file, specify the following input as follows:
0.25 ! Rin : inner disk radius [AU]
300.0 ! Rout : outer disk radius [AU]
0.01 ! dust_to_gas : dust/gas mass ratio
0.1 ! fPAH : PAH abundance with respect to ISM
70 ! NXX : number of radial grid points
70 ! NZZ : number of vertical grid points
.false. ! solve_diskstruc : solve the vertical hydrostatic eq.?
.true. ! solve_temp : solve the heating/cooling balance?
.false. ! conserve_pressure : conserve pgas instead of n<H>?
.false. ! restart : restart from checkpoint file?
.false. ! radtrans : calculate dust continuum rad. transfer?
.true. ! chi_from_RT : calculate chi from UV rad. transfer?
.true. ! Td_from_RT : calculate dust temp. from rad. transfer?
.true. ! Jback_from_RT : calculate background Jnu from rad. transfer?
.true. ! UVpumping : use large model atoms?
.true. ! PAH_from_RT : PAH heating from cross-sections?
.true. ! Rphoto_from_RT : calculate photorates from rad. transfer?
.false. ! Rphoto_bandint : use band-integrated photo-rates?
------ mcfost ------
.false. ! use_MCFOST_rgrid
.true. ! readMCFOST
forProDiMo.fits.gz
------ image and SED parameters ------
.true. ! calcSED : calculate SED and images?
.true. ! monoSED : use monochromatic mode?
300 ! NlamSED : number of lambda gridpoints
0.1 ! lminSED : minimum lambda[mic]
3000.0 ! lmaxSED : minimum lambda[mic]
140.0 ! dist [pc] : distance
30.0 ! incl [deg] : inclination (0=face-on)
------ line transfer? ------
.true. ! line_transfer : calculate line transfer?Most of the other physical parameters, e.g. the stellar parameters, the dust size and opacity parameters, the dust settling and the disk shape parameters are meaningless, because they will be retrieved from forProDiMo.fits.gz and overwritten.
I recommend 30x30 grid for tests, 70x70 for serious models, and 100x100 for publications. The switch use_MCFOST_rgrid should be .false.. If this parameter is switched on, the r-gridpoints of MCFOST are taken over and the settings for NXX are ignored, however, Nzz still specifies the number of vertical grid points, creating a 100x70 grid in this case. However, the MCFOST grid is not very suitable for ProDiMo, in particular there are serious problems with mutli-zone models and use_MCFOST_rgrid=.true.
Rin and Rout should be the same than in the MCFOST parameter file.
radtrans=.false.: ProDiMo will not use it's own radiative transfer, instread the disk internal radiation field J_nu(r,z) and the dust temperature structure Tdust(r,z) are taken over from MCFOST. The radiation field can be noisy, in that case increase the number of photons in the MCFOST parameter file. The noise is stronger for optically thick disks. For test purposes, it is interesting to run with radtrans=.true., in which case ProDiMo will recompute the disk internal radiation field using MCFOST's dust opacities, and use Tdust and Jnu as initial guesses, however, remember that ProDiMo can only do isotropic scattering.
Rphoto_bandint=.false. debatable. If you have line-dominated UV input, one should use .true., however, if you model e.g. an A_star without excess UV, the interstellar hard UV will start to dominate at some radius, and then .true. will produce wrong results. If you have a strong powerlaw-UV, it shouldn't matter. Problem is that the spectral shape of the original stellar spectrum (+UV) should be the same on both sides MCFOST and ProDiMo, and the way to average bands should be the same, too, but this is not guaranteed at the moment (would need more testing). Therefore, I'd recommend Rphoto_bandint=.false. at the moment.
calcSED=.true. will re-calculate the SED using ProDiMo's raytracing code. This takes about 40sec, but provides an inportant sanity check. The SEDs of MCFOST and ProDiMo should look very similar in that case, in particular if you specify "-isotropic" for your MCFOST run. If you see large differences, it means that something went wrong - probably in the way the gridpoints are set.
The figure shows the kind of SED accuracy we are looking for. If you can achieve that quality, you can be quite sure that the dust opacities, temperatures, mean intensities are properly passed on to ProDiMo, and that ProDiMo produced a suitable grid.
Multi-zone models¶
When you run a multi-zone model, make sure to specify dust/gas and fPAH in each zone, noting that the outer zone in ProDiMo is zone 1, the innermost zone in ProDiMo is zone 2 etc.
zone2 ... [zone 3 ... [zone 4 ...]] zone 1Where "[]" means you can omit. In a 2-zone-model, you should specify in Parameter.in
0.005 ! dust_to_gas in outer zone
0.01 ! d2ust_to_gas in inner zone
0.01 ! fPAH in outer zone
0.001 ! f2PAH in inner zoneThe re-adjustment of a suitable radial ProDiMo grid from MCFOST's exchange file forProDiMo.fits.gz should work automatically, including the determination of the number of zones NZONES, and gaps, but it is wise to study the stdout from ProDiMo carefully:
READING FROM MCFOST: forProDiMo.fits.gz ...
MCFOST version 2.14.1
MCFOST to ProDiMo version 2
n_rad, nz, n_lambda = 120 70 70
Mstar = 1.900000
Teff = 6200.000
Rstar = 3.600000
Lstar = 17.1870152705656
fUV = 0.0000000E+00
pUV = 2.200000
Mdust = 2.0001001E-04
Rin = 0.2500000
Rout = 30.00000
H0 = 1.100000
r0 = 10.00000
epsi = 0.5000000
beta = 1.000000
amin = 0.1000000
amax = 3.000000
apow = 3.500000
rho_gr = 3.500000
settle = 0.0000000E+00
asettle = 1.000000
... detected 1 holes and 2 zones:
zone 2 from 0.25000002335467 AU to 28.26718910431741 AU
zone 1 from 140.00000197610575 AU to 300.00000000000000 AU
disk from 0.25000002335467 AU to 300.00000000000000 AU
zone 2: dust/gas = 1.000E-02 fPAH = 1.000E-02
zone 1: dust/gas = 1.000E-02 fPAH = 1.000E-02
beta_max= 30.9426872605686
detected Nzone = 2
detected Ngaps = 1
...
taking over lambda-grid & incident Int.
NUV1,NUV= 3 9
...
calculating column densities
1 0.250000023354674 5.908288995499093E+021
2 0.251530556370283 5.890286025001251E+021
3 0.256122272190482 5.837247221078691E+021
4 0.286339433536542 5.509062128063851E+021
5 0.301110609204870 5.383549741286824E+021
6 0.334781252912953 5.096120290735470E+021
7 0.366987008073658 4.858816950382459E+021
8 0.384902504346483 4.761638167286310E+021
9 0.435174136721880 4.478166331312444E+021
10 0.492011683824119 4.211570268093337E+021
11 0.556272711524112 3.960845298886777E+021
12 0.628926791293445 3.725046596281422E+021
13 0.711070129115120 3.503285578347285E+021
14 0.803942105057949 3.294726523238409E+021
15 0.908943973063969 3.098583443040433E+021
16 1.14004427285060 2.761037101408729E+021
17 1.31363325208360 2.577476374780669E+021
18 1.48520524026011 2.424032960794197E+021
19 1.67918603019324 2.279724429288384E+021
20 2.14646342083968 2.016368999370176E+021
21 2.42681031071742 1.896329665410184E+021
22 3.08705753519752 1.680318860612417E+021
23 3.50729947623887 1.577412685935551E+021
24 4.45804311171118 1.398143427450189E+021
25 5.06885501586194 1.312129867894473E+021
26 6.47939637506116 1.160551678880762E+021
27 7.32566219277675 1.091461226944690E+021
28 9.36421912014953 9.653748279905089E+020
29 11.9700577808020 8.538540334777842E+020
30 13.5334519844599 8.030220369311671E+020
31 17.2994886331865 7.102563455621818E+020
32 22.1135233873364 6.282070182259526E+020
33 28.2671891043174 5.556360907750159E+020
34 28.5498609953606 25604635691649.7
35 36.3094386203707 32563729406387.8
36 46.1779948119784 41414237859992.9
37 58.7287296603579 52670229386798.2
38 74.6906335271345 66985491150081.9
39 94.9908293427696 85191503375954.2
40 120.808423133123 108345734618776.
41 133.920397965329 120105068187205.
42 138.600001956345 124301920682935.
43 140.000001976106 2.887202287376513E+023
44 140.129504749818 2.884533409562567E+023
45 140.518022951484 2.876555751826006E+023
46 142.072095758149 2.845083468930210E+023
47 144.171453014339 2.803648161775390E+023
48 147.016568287865 2.749384571347583E+023
49 152.737489711649 2.646340929800989E+023
50 155.893245257319 2.592825318925674E+023
51 158.969681291275 2.542658969101490E+023
52 165.305884907291 2.445202110423491E+023
53 168.568072310692 2.397884030098453E+023
54 171.894636530810 2.351487013651214E+023
55 178.587746278871 2.263335794889140E+023
56 182.273420221898 2.217612375856819E+023
57 185.870449182822 2.174693098458693E+023
58 192.594573399793 2.098599447401853E+023
59 197.093075005053 2.050848155181887E+023
60 200.982558715447 2.011153616892906E+023
61 204.948798463733 1.972229746473382E+023
62 212.450945976816 1.902462264592444E+023
63 217.323351212632 1.859919634771975E+023
64 221.612064219967 1.823925626531229E+023
65 225.985411755329 1.788625183899020E+023
66 233.715409165855 1.729330516426738E+023
67 239.259237361201 1.689363168157565E+023
68 244.359049470396 1.654163299603391E+023
69 249.181291664232 1.622151247351913E+023
70 254.098697183619 1.590756097798077E+023
71 262.376974363645 1.540383579644759E+023
72 268.397378807194 1.505881714815123E+023
73 274.516345789186 1.472398078111046E+023
74 280.180216424611 1.442660825499823E+023
75 285.709362426977 1.414742452906459E+023
76 290.126807337815 1.393105748760227E+023
77 293.887285002163 1.375240336365060E+023
78 296.902248084777 1.361395808031002E+023
79 299.035739889000 1.351981547435153E+023
80 300.000000000000 1.347838923338465E+023
checking dust and gas masses ...
dust_to_gas_c = 1.000000000000000E-002
ProDiMo: Mgas/Msun = 1.999E-02 Mdust/Msun = 2.224E-04
MCFOST: Mgas/Msun = 2.140E-02 Mdust/Msun = 2.371E-04
end reading data from MCFOST-ProDiMo fileThe lower part of the output shows the final choice of ProDiMo's radial gridpoints and associated column densities. There should be a point with large column density right at the extracted value of the outer radius of the inner zone R2out, next grid point should have very low column density situated shortly thereafter in the gap, then a point with large column density at the inner radius of the outer zone Rin and the previous point, still in the gap, should have very low column density and small distance to Rin, see multiple_zones for further details.
Peter Woitke, 24.10.2012
ProDiMo to MCFOST¶
This interface allows to run the MCFOST line radiative transfer on top of a ProDiMo model. Where ProDiMo provides the molecular concentrations, the gas temperature and the population levels.
TODO: How To do this? MCFOST still seems to have the option.