.. _mesa-web-profile: ********************** MESA-Web Profile Files ********************** Preliminaries ============= In the previous labs, you've explored the :ref:`basics of Python and Jupyter notebooks`, undertaken a :ref:`MESA-Web calculation `, and learned how to :ref:`plot data from a MESA-Web history file `. In this lab, you'll use Jupyter to plot data from the profile files created during your calculation. Exploring a Profile File ======================== Using the Profiles Index ------------------------ Profile files provide a snapshot of the full internal structure of the model star --- i.e., variables such as radial coordinate :math:`r`, interior mass :math:`m`, temperature :math:`T` --- at a single timestep. Because they can be quite large, MESA-Web only writes them for selected timesteps. An index of the profile files, and the model numbers they correspond to, is provided in the MESA-Web output folder as the :file:`profiles.index` file. Go ahead and take a look at this index file for your :ref:`calculation `; you can do this using Notepad (on Windows), TextEdit (on MacOS), or from the file browser that Jupyter provides. You should see something like this: .. literalinclude:: profiles.index :language: console The first column of this file lists model numbers for which corresponding profile files exist; and the third column gives the profile number (don't worry about the second column for now). Let's figure out if there is a profile file for the present-day sun. In :ref:`Lab III ` we found that the present-day Sun (defined by the point where :math:`L_{\rm H} = L_{\odot}`) occurs in our calculation at model number 298. The index file shows that there isn't a specific profile file for this model number; however, there *is* a file for model number 297 (which should have a very similar structure), and it has a profile number of 8. Creating a Notebook ------------------- Fire up Jupyter in your web browser and create a new notebook in the ``astro_310`` folder, naming it ``profile-notebook``. Then, import modules: .. code:: # Import the numpy module to provide numerical functionality import numpy as np # Import the matplotlib.pyplot module to provide plotting functionality import matplotlib.pyplot as plt # Tell matplotlib.pyplot to do inline plots %matplotlib inline # Import the mesa-web module to simplify reading MESA-Web files import mesa_web as mw (as usual, cut and paste this code into a notebook cell, and then execute the cell). Then, customize the following code to read in the file ``profileP.data`` from the MESA-Web folder (where ``P`` is the profile number determined above). .. code:: # Read history data. Be sure to replace the NNNNNNNNNN with the # specific digits of your folder, and P with the profile number prof_data = mw.read_profile('MESA-Web_Job_NNNNNNNNNN/profileP.data') # Inspect the prof_data variable print(type(prof_data)) print(prof_data.keys()) As was the case with ``hist_data``, we see that ``prof_data`` is a dict containing many key/value pairs, and we can use ``help()`` command to access the documentation describing these data: .. code:: # Print out documentation for the read_profile function help(mw.read_profile) Let's verify that we read in the correct profile, by checking the model number: .. code:: # Check the model number print(prof_data['model_number']) Plotting Profile Data --------------------- We're now in a position to make some plots. First, a simple plot showing how the pressure, density and temperature vary as a function of fractional radius :math:`r/R` throughout the star: .. code:: # Define constants RSUN = 6.957E10 # Solar radius, in cm MSUN = 1.989E33 # Solar mass, in g # Plot pressure, density and temperature versus fractional radius P = prof_data['pressure'] rho = 10**prof_data['logRho'] T = 10**prof_data['logT'] r = prof_data['radius']*RSUN # Convert radii to cm R = prof_data['photosphere_r']*RSUN # Same plt.figure() plt.plot(r/R, P, color='r', label='P (dyn/cm^2)') plt.plot(r/R, rho, color='g', label='rho (g/cm^3)') plt.plot(r/R, T, color='b', label='T (K)') plt.xlabel('r/R') plt.yscale('log') plt.legend() Note how we make the vertical axis be logarithmic by using the ``plt.yscale('log')`` command. We're not limited to having radius on the horizontal axis; we can use any variable we like. Here a plot of the hydrogen, helium, carbon and nitrogen mass fractions versus fractional interior mass :math:`m/M`: .. code:: # Plot abundances versus fractional interior mass X = prof_data['h1'] Y = prof_data['he4'] Z_C12 = prof_data['c12'] Z_N14 = prof_data['n14'] m = prof_data['mass']*MSUN # Convert masses to g M = prof_data['star_mass']*MSUN # Same plt.figure() plt.plot(m/M, X, color='r', label='X') plt.plot(m/M, Y, color='g', label='Y') plt.plot(m/M, Z_C12, color='b', label='Z_C12') plt.plot(m/M, Z_N14, color='m', label='Z_N14') plt.xlabel('m/M') plt.yscale('log') plt.legend() Analyzing Profile Data ---------------------- To finish the lab, let's do a little analysis on the profile data. Specifically, let's see if we can confirm that the model for the present-day Sun satisfies the hydrostatic equilibrium equation .. math:: \frac{{\rm d}P}{{\rm d}r} = -\frac{G m}{r^{2}} \rho. First, let's calculate the left-hand and right-hand sides of this equation. To evaluate the pressure gradient, we can use the ``np.gradient`` function, which uses `finite differences `__ to evaluate a numerical gradient estimate: .. code:: # Define constants G = 6.674E-8 # Gravitational constant, in cm^3/g/s^2 # Evaluate lhs and rhs of the hydrostatic equilibrium equation lhs = np.gradient(P, r) rhs = -G*m/r**2 * rho Now, create a plot comparing the two: .. code:: # Plot the lhs and rhs of the hydrostatic equilibrium equation plt.figure() plt.plot(r/R, lhs, color='b', label='lhs') plt.plot(r/R, rhs, color='r', label='rhs', dashes=(2,1)) plt.xlabel('r/R') plt.legend() The correspondence between the two is reasonable, but the pressure gradient data are very noisy. This is a consequence of using finite differences to evaluate the gradient. Internally, MESA uses a more sophisticated approach, assuring that hydrostatic balance is achieved to a high degree of accuracy.