Merge pull request VirtualPlanetaryLaboratory#48 from VirtualPlanetar…

…yLaboratory/dev

Dev

INSTALL

@@ -9,20 +9,34 @@ simulate the examples in this repository.
         N.B. We recommend installing python via anaconda
     C. matplotlib v1.4.0 or higher
 
-2. Clone the repository from GitHub.
-    Compile the code: > make opt
-
-3. Pip install or clone vplot.
-    > pip install vplot
-    OR
-    > git clone https://github.com/VirtualPlanetaryLaboratory/vplot.git
-    > cd vplot
-    > python setup.py develop
-    We recommend installing as administrator, but it does work for regular users
-      by changing the last command to
-    > python setup.py install --user
-
-4. If you encountered no errors, then everything should be setup. To quickly
+2. Clone the repository from GitHub:
+      > git clone https://github.com/VirtualPlanetaryLaboratory/vplanet.git
+   Compile the code:
+      > make opt
+
+3. At this point you should be able to run the code on one of the examples. To
+   verify it's working, execute the following commands:
+   > cd examples/EarthInterior
+   > vplanet vpl.in
+
+   This should create files called earth.earth.forward earth.log that contain
+   results from the simulation that you can check. This example recreates the
+   thermal and magnetic evolution of Earth and should require about 1 second.
+
+4. VPLanet has its own plotting package which can streamline your research. This
+   package is not required, but you can only reproduce the plots in this
+   repository exactly if you install it. There are two ways to install it, via
+   pip, or cloning from GitHub:
+   > pip install vplot
+   OR
+   > git clone https://github.com/VirtualPlanetaryLaboratory/vplot.git
+   > cd vplot
+   > python setup.py install --user
+   We recommend installing as administrator, which can be completed by changing
+      the last command to
+   > sudo python setup.py develop
+
+5. If you encountered no errors, then everything should be setup. To quickly
     determine if it worked, execute the following from the vplanet directory:
 
     > cd examples/EarthInterior
@@ -32,7 +46,7 @@ simulate the examples in this repository.
     These commands require less than 1 minute and should create 2 png figures
       identical to those in the examples/EarthInterior directory.
 
-5. If you did encounter errors, it is probably because either a) vplanet is not
+6. If you did encounter errors, it is probably because either a) vplanet is not
     in your PATH, b) the proper version of python is not in your PATH, c) you do
     not have the proper version of python installed, or d) you need to update
     matplotlib. If you have followed these instructions exactly and still cannot

examples/MagneticBraking/README.rst

@@ -7,7 +7,7 @@ Overview
 Rotational evolution of stars due to magnetic effects.
 
 ===================   ============
-**Date**              10/21/18
+**Date**              06/19/19
 **Author**            David Fleming
 **Modules**           STELLAR
 **Approx. runtime**   2 minutes
@@ -18,6 +18,9 @@ evolution and magnetic braking. This example shows how the three available
 magnetic braking laws (`Reiners & Mohanty [2012] <https://ui.adsabs.harvard.edu/abs/2012ApJ...746...43R/abstract>`_, `Repetto & Nelemans [2014] <https://ui.adsabs.harvard.edu/abs/2014MNRAS.444..542R/abstract>`_,
 and `Matt et al. [2015] <https://ui.adsabs.harvard.edu/abs/2015ApJ...799L..23M/abstract>`_) impact the rotation period.
 
+This example also demonstrates how the VPLanet implementation of the `Matt et al. [2015] <https://ui.adsabs.harvard.edu/abs/2015ApJ...799L..23M/abstract>`_
+magnetic braking model can be used to reproduce the upper enveloped of the Kepler
+field stellar rotation period distribution, a main result from `Matt et al. [2015] <https://ui.adsabs.harvard.edu/abs/2015ApJ...799L..23M/abstract>`_.
 
 To run this example
 -------------------
@@ -37,3 +40,19 @@ Expected output
 
 Rotation period evolution of 0.1 and 1 solar-mass stars due to stellar evolution
 (`Baraffe et al. 2015 <https://ui.adsabs.harvard.edu/abs/2015A%26A...577A..42B/abstract>`_) and the three magnetic braking laws.
+
+.. figure:: kepler.png
+   :width: 600px
+   :align: center
+
+Rotation period distribution of a ~4 Gyr-old synthetic cluster of stars
+simulated using STELLAR with the `Matt et al. [2015] <https://ui.adsabs.harvard.edu/abs/2015ApJ...799L..23M/abstract>`_ magnetic braking model
+(black, adapting VPLanet simulations from `Fleming et al. [2019] <https://ui.adsabs.harvard.edu/abs/2019arXiv190305686F/abstract>`_). Following Fig. (3) in
+Matt et al. [2015], we compare the `Fleming et al. [2019] <https://ui.adsabs.harvard.edu/abs/2019arXiv190305686F/abstract>`_ simulated distribution to
+the rotation distribution of Kepler field stars (red) measured by
+`McQuillan et al. [2014] <https://ui.adsabs.harvard.edu/abs/2014ApJS..211...24M/abstract>`_. For reference, we plot the modern solar rotation period
+as a blue star. Using STELLAR, `Fleming et al. [2019] <https://ui.adsabs.harvard.edu/abs/2019arXiv190305686F/abstract>`_
+recover the `Matt et al. [2015] <https://ui.adsabs.harvard.edu/abs/2015ApJ...799L..23M/abstract>`_
+result that the upper envelope of the Kepler stellar rotation period
+distribution is well-matched by a 4 Gyr-old synthetic cluster, validating the
+STELLAR implementation of the `Matt et al. [2015] <https://ui.adsabs.harvard.edu/abs/2015ApJ...799L..23M/abstract>`_ magnetic braking model.

examples/MagneticBraking/makeplot.py

@@ -1,7 +1,6 @@
 """
 This script produces a figure comparing VPLanet's various magnetic braking
 implementations.
-
 David P. Fleming, University of Washington, 2018
 """
 
@@ -12,6 +11,12 @@
 import numpy as np
 import sys
 import vplot as vpl
+import pandas as pd
+from matplotlib.patches import Patch
+from matplotlib.lines import Line2D
+
+mpl.rcParams['figure.figsize'] = (9,8)
+mpl.rcParams['font.size'] = 18.0
 
 # Check correct number of arguments
 if (len(sys.argv) != 2):
@@ -23,9 +28,7 @@
     print('Options are: pdf, png')
     exit(1)
 
-# Make the plot!
-mpl.rcParams['figure.figsize'] = (9,8)
-mpl.rcParams['font.size'] = 18.0
+### Magnetic braking validation figure ###
 
 ms = ["a_matt", "a_reiners", "a_sk"]
 gs = ["b_matt", "b_reiners", "b_sk"]
@@ -66,3 +69,49 @@
     plt.savefig('MagneticBraking.pdf', bbox_inches="tight", dpi=600)
 if (sys.argv[1] == 'png'):
     plt.savefig('MagneticBraking.png', bbox_inches="tight", dpi=600)
+
+### Now make Kepler comparison figure ###
+
+#Typical plot parameters that make for pretty plots
+mpl.rcParams['font.size'] = 17
+
+## for Palatino and other serif fonts use:
+mpl.rc('font',**{'family':'serif'})
+mpl.rc('text', usetex=True)
+
+# Load data
+kep = pd.read_csv("mcSingleMarch27.csv")
+mcq = pd.read_csv("mcquillan2014.tsv", delimiter="\t", comment="#", header=0)
+
+# Plot distribution within 100 Myr of age ~ 4 Gyr
+mask = np.fabs(kep["Age"].values - 4.0e9) < 1.0e8
+massesUp = kep["Pri_dMass"].values[mask]
+protsUp = kep["Pri_ProtAge"].values[mask]
+
+fig, ax = plt.subplots(figsize=(6,5))
+
+ax.scatter(massesUp, protsUp, color="k", s=10, edgecolor=None, zorder=10)
+ax.scatter(mcq["Mass"], mcq["Prot"], color="r", s=2, edgecolor=None, zorder=0, alpha=0.2)
+ax.scatter([1.0], [26.3], marker="*", color="C0", s=100, edgecolor=None, zorder=20)
+
+# Format
+ax.set_xlabel(r"Stellar Mass [$M_{\odot}$]")
+ax.set_ylabel("Rotation Period [d]")
+ax.set_xlim(0.1, 1.025)
+ax.set_ylim(0.1, 60)
+ax.set_yscale("log")
+legend_elements = [Line2D([0], [0], marker='o', color='w', label='VPLanet Age = 4 Gyr',
+                          markerfacecolor='k', markersize=10),
+                   Line2D([0], [0], marker='o', color='w', label='McQuillan et al. (2014)',
+                          markerfacecolor='red', markersize=10),
+                   Line2D([0], [0], marker='*', color='w', label='Sun',
+                          markerfacecolor='C0', markersize=15)]
+ax.legend(handles=legend_elements, loc='best', fontsize=12)
+
+# Save!
+if (sys.argv[1] == 'pdf'):
+    fig.savefig("kepler.pdf", bbox_inches="tight", dpi=200)
+if (sys.argv[1] == 'png'):
+    fig.savefig("kepler.png", bbox_inches="tight", dpi=200)
+
+# Done!