{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Using Wakeflow with MCFOST" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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_." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prerequisites" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "Let's check that our MCFOST installation is up to date and working by running the following shell command:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " You are running MCFOST 3.0.44\n", " Git SHA = 8686df4719b1cadb7b1d82257fc576ec8072e0c0\n", " Binary compiled the Jul 18 2022 at 11:56:15\n", " with INTEL compiler version 2021\n", " \n", " Checking last version ...\n", " MCFOST is up-to-date\n" ] } ], "source": [ "!mcfost -v" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And checking the `pymcfost` installation by importing it:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "WARNING: mpl_scatter_density is not present\n", "WARNING: progressbar is not present\n" ] } ], "source": [ "import pymcfost" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Provided both of these installations are working, and `wakeflow` is installed, you should be ready to proceed with the rest of the tutorial." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Required Parameters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The _crucial_ parameters that you _must_ specify are:\n", "\n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running Wakeflow" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Model initialised.\n", "Model configured.\n", "\n", "* Performing checks on model parameters:\n", "M_thermal = 0.967 M_Jup\n", "M_planet = 4.135 M_th\n", "WARNING: M_planet > M_thermal. This may break the solution.\n", "Parameters Ok - continuing\n", "Generating MCFOST grid data...\n", " You are running MCFOST 3.0.44\n", " Git SHA = 8686df4719b1cadb7b1d82257fc576ec8072e0c0\n", " WARNING: polarization is turned off in ray-traced SEDs\n", " it can be turned back on with -rt2\n", " Input file read successfully\n", " Computation of disk structure\n", " Creating directory ././data_disk\n", " Parallelized code on 10 processors\n", " \n", "Tue 23 Aug 2022 22:27:27 AEST\n", " Forcing 3D mode\n", " Using ray-tracing method 1\n", " Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc\n", " s/mcfost_utils/Dust/Draine_Si_sUV.dat\n", " Number of regions detected: 1\n", " zone 1 --> region= 1 : R=100.00 to 500.00 AU\n", " Using 4.480000 million cells\n", " Total gas mass in model: 0.1000000 Msun\n", " Total dust mass in model: 1.0000000E-03 Msun\n", " Using scattering method 2\n", " Writing disk structure files in data_disk ...\n", " Exiting\n", "Done\n", "\n", "* Creating 4.0 Mj model:\n", "Generating unperturbed background disk\n", "Extracting linear perturbations nearby planet\n", "Propagating outer wake... \n", "Completed in 0.17 s\n", "Propagating inner wake... \n", "Completed in 0.17 s\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "* Mapping to physical coords: 100%|██████████| 200/200 [00:05<00:00, 37.37it/s] \n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "* Saving results:\n", "Perturbations saved to HD_163296/mcfost_tutorial/4.0Mj\n", "Total saved to HD_163296/mcfost_tutorial/4.0Mj\n", "Saved FITS file formatted for MCFOST.\n", "\n", "* Done!\n" ] } ], "source": [ "from wakeflow import WakeflowModel\n", "\n", "hd163_model = WakeflowModel()\n", "\n", "hd163_model.configure(\n", " name = \"mcfost_tutorial\",\n", " system = \"HD_163296\",\n", " m_star = 1.9,\n", " m_planet = 4.0,\n", " r_outer = 500,\n", " r_planet = 256,\n", " r_ref = 256,\n", " q = 0.35,\n", " p = 2.15,\n", " hr = 0.09,\n", " cw_rotation = False,\n", " grid_type = \"mcfost\",\n", " n_r = 200,\n", " n_phi = 280,\n", " n_z = 40,\n", " show_midplane_plots = False,\n", " write_FITS = True\n", ")\n", "\n", "hd163_model.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets see what `wakeflow` has made for us:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[34m4.0Mj\u001b[m\u001b[m mcfost_tutorial_config.yaml\n", "\u001b[1m\u001b[34mmcfost\u001b[m\u001b[m\n" ] } ], "source": [ "!ls HD_163296/mcfost_tutorial" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "delta_PHI.npy delta_v_r.npy total_Z.npy vr_z0.pdf\n", "delta_R.npy mcfost.para total_rho.npy wakeflow_model.fits\n", "delta_Z.npy rho_z0.pdf total_v_phi.npy\n", "delta_rho.npy total_PHI.npy total_v_r.npy\n", "delta_v_phi.npy total_R.npy vphi_z0.pdf\n" ] } ], "source": [ "!ls HD_163296/mcfost_tutorial/4.0Mj" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running MCFOST" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once you have configured the `.para` file to your liking, we are ready to run `MCFOST`. First we move to the results directory:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/tomhilder/Documents/Research/Protoplanetary_Discs/wakeflow/docs/tutorials/HD_163296/mcfost_tutorial/4.0Mj\n" ] } ], "source": [ "%cd HD_163296/mcfost_tutorial/4.0Mj" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now run `MCFOST` on our model easily using the following:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " You are running MCFOST 3.0.44\n", " Git SHA = 8686df4719b1cadb7b1d82257fc576ec8072e0c0\n", " WARNING: polarization is turned off in ray-traced SEDs\n", " it can be turned back on with -rt2\n", " Input file read successfully\n", " Thermal equilibrium calculation\n", " Temperature calculation under LTE approximation\n", " Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc\n", " s/mcfost_utils/Molecules/co@xpol.dat\n", " Molecular line transfer under LTE approximation\n", " Parallelized code on 10 processors\n", " \n", "Tue 23 Aug 2022 22:28:11 AEST\n", " Creating directory ././data_th\n", " Creating directory ././data_CO\n", " Forcing 3D mode\n", " Using ray-tracing method 1\n", " Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc\n", " s/mcfost_utils/Dust/Draine_Si_sUV.dat\n", " Number of regions detected: 1\n", " zone 1 --> region= 1 : R=100.00 to 500.00 AU\n", " Using 4.480000 million cells\n", " Reading density file : wakeflow_model.fits\n", " read_gas_density = 1\n", " read_gas_velocity = 2\n", " Reading dust density ...\n", " No grain size found\n", " The density file only has positive z, making it symmetric\n", " Dust density range: 8.2326818E-10 11.77502 \n", " Using grain size distribution from parameter file\n", " Reading gas density ...\n", " Updating gas-to-dust ratio to 100.0000 \n", " Gas density range: 8.2326818E-10 11.77502 \n", " Reading gas velocity ...\n", " Velocity field is in cylindrical coordinates.\n", " Velocity range: -842.0961 4076.110 \n", " Constant spatial distribution\n", " Total gas mass in model: 0.1000000 Msun\n", " Total dust mass in model: 1.0000000E-03 Msun\n", " Using scattering method 2\n", " Trying to find appropriate stellar spectra ...\n", " Star # 1 --> lte9200-4.0.NextGen.fits.gz\n", " Done\n", " Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc\n", " s/mcfost_utils/Stellar_Spectra/lte9200-4.0.NextGen.fits.gz\n", " Computing dust properties ... Done\n", " lambda = 0.871355626679430 \n", " Integ tau dans plan eq. = 81.53822 \n", " Column density (g/cm�) = 0.1843978 \n", " Integ tau (i =75.0 deg) = 0.10330E+01\n", " Column density (g/cm�) = 2.3361719E-03\n", " No dark zone\n", " Initializing thermal properties ... Done\n", " Initialization complete in 45.00s\n", " Computing temperature structure ...\n", " 100% |==================================================|\n", " Max. temperature = 96.65446 \n", " Temperature calculation complete in 0h 1m 12s\n", " Source fct time 0.0000000E+00 s\n", " RT time 0.0000000E+00 s\n", " \n", " -------------------------------------------------------\n", " Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc\n", " s/mcfost_utils/Molecules/co@xpol.dat\n", " \n", " -------------------------------------------------------\n", " CO molecular file read successfully\n", " WARNING : memory size if lots of pixels\n", " Setting constant abundance\n", " Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc\n", " s/mcfost_utils/Dust/Draine_Si_sUV.dat\n", " Computing dust properties for 40 wavelength\n", " Trying to find appropriate stellar spectra ...\n", " Star # 1 --> lte9200-4.0.NextGen.fits.gz\n", " Done\n", " Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc\n", " s/mcfost_utils/Stellar_Spectra/lte9200-4.0.NextGen.fits.gz\n", " -------------------------------\n", " Transition J= 2 - 1\n", " tau_mol = 1157.638 \n", " tau_dust= 0.4059193 \n", " Vertical Tau_mol=1 (for r=100 au) at z= 13.05515 au\n", " Processing complete in 0h 2m 51s\n", " CPU time used 0h 10m 41s\n" ] } ], "source": [ "!mcfost mcfost.para -df wakeflow_model.fits -mol" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "_dust_prop_th.tmp delta_Z.npy rho_z0.pdf total_v_phi.npy\n", "\u001b[1m\u001b[34mdata_CO\u001b[m\u001b[m delta_rho.npy total_PHI.npy total_v_r.npy\n", "\u001b[1m\u001b[34mdata_th\u001b[m\u001b[m delta_v_phi.npy total_R.npy vphi_z0.pdf\n", "delta_PHI.npy delta_v_r.npy total_Z.npy vr_z0.pdf\n", "delta_R.npy mcfost.para total_rho.npy wakeflow_model.fits\n" ] } ], "source": [ "!ls" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/tomhilder/Documents/Research/Protoplanetary_Discs/wakeflow/docs/tutorials\n" ] } ], "source": [ "%cd ../../.." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Extra Integrations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Model initialised.\n", "Model configured.\n", "\n", "* Performing checks on model parameters:\n", "M_thermal = 0.967 M_Jup\n", "M_planets = [1.034 2.068 3.101] M_th\n", "At least one planet mass exceeds thermal mass. This may break the solution. \n", "Parameters Ok - continuing\n", "Generating MCFOST grid data...\n", " You are running MCFOST 3.0.44\n", " Git SHA = 8686df4719b1cadb7b1d82257fc576ec8072e0c0\n", " WARNING: polarization is turned off in ray-traced SEDs\n", " it can be turned back on with -rt2\n", " Input file read successfully\n", " Computation of disk structure\n", " Creating directory ././data_disk\n", " Parallelized code on 10 processors\n", " \n", "Tue 23 Aug 2022 22:33:19 AEST\n", " Forcing 3D mode\n", " Using ray-tracing method 1\n", " Reading /Users/tomhilder/Users/tomhilder/Documents/Research/Protoplanetary_Disc\n", " s/mcfost_utils/Dust/Draine_Si_sUV.dat\n", " Number of regions detected: 1\n", " zone 1 --> region= 1 : R=100.00 to 500.00 AU\n", " Using 4.480000 million cells\n", " Total gas mass in model: 0.1000000 Msun\n", " Total dust mass in model: 1.0000000E-03 Msun\n", " Using scattering method 2\n", " Writing disk structure files in data_disk ...\n", " Exiting\n", "Done\n", "\n", "* Creating 1.0 Mj model:\n", "Generating unperturbed background disk\n", "Extracting linear perturbations nearby planet\n", "Propagating outer wake... \n", "Completed in 0.68 s\n", "Propagating inner wake... \n", "Completed in 0.71 s\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "* Mapping to physical coords: 100%|██████████| 200/200 [00:05<00:00, 33.71it/s] \n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "* Saving results:\n", "Perturbations saved to HD_163296/mcfost_tutorial_2/1.0Mj\n", "Total saved to HD_163296/mcfost_tutorial_2/1.0Mj\n", "Saved FITS file formatted for MCFOST.\n", "\n", "* Creating 2.0 Mj model:\n", "Generating unperturbed background disk\n", "Extracting linear perturbations nearby planet\n", "Propagating outer wake... \n", "Completed in 0.33 s\n", "Propagating inner wake... \n", "Completed in 0.35 s\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "* Mapping to physical coords: 100%|██████████| 200/200 [00:05<00:00, 35.19it/s] \n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "* Saving results:\n", "Perturbations saved to HD_163296/mcfost_tutorial_2/2.0Mj\n", "Total saved to HD_163296/mcfost_tutorial_2/2.0Mj\n", "Saved FITS file formatted for MCFOST.\n", "\n", "* Creating 3.0 Mj model:\n", "Generating unperturbed background disk\n", "Extracting linear perturbations nearby planet\n", "Propagating outer wake... \n", "Completed in 0.22 s\n", "Propagating inner wake... \n", "Completed in 0.23 s\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "* Mapping to physical coords: 100%|██████████| 200/200 [00:05<00:00, 36.14it/s] \n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "* Saving results:\n", "Perturbations saved to HD_163296/mcfost_tutorial_2/3.0Mj\n", "Total saved to HD_163296/mcfost_tutorial_2/3.0Mj\n", "Saved FITS file formatted for MCFOST.\n", "\n", "* Running MCFOST on each model:\n", "1.0 Mj...\n", "2.0 Mj...\n", "3.0 Mj...\n", "\n", "* Done!\n" ] } ], "source": [ "hd163_model_2 = WakeflowModel()\n", "\n", "hd163_model_2.configure(\n", " name = \"mcfost_tutorial_2\",\n", " system = \"HD_163296\",\n", " m_star = 1.9,\n", " m_planet = [1.0, 2.0, 3.0],\n", " r_outer = 500,\n", " r_planet = 256,\n", " r_ref = 256,\n", " q = 0.35,\n", " p = 2.15,\n", " hr = 0.09,\n", " cw_rotation = False,\n", " grid_type = \"mcfost\",\n", " n_r = 200,\n", " n_phi = 280,\n", " n_z = 40,\n", " show_midplane_plots = False,\n", " write_FITS = True,\n", " run_mcfost = True,\n", " inclination = 225,\n", " PA = 45,\n", " PAp = 125,\n", " temp_star = 9250,\n", " distance = 101.5,\n", " v_max = 4.0,\n", " n_v = 36\n", ")\n", "\n", "hd163_model_2.run()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.9.13 64-bit ('base')", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" }, "orig_nbformat": 4, "vscode": { "interpreter": { "hash": "2ac3d5d3a4336f34c662130573d1ea696fb4dfcbf6baf7c9925b84262e1f6ae2" } } }, "nbformat": 4, "nbformat_minor": 2 }