Get started with *mascaret*
===========================

The **mascaret** module is a linear oscillation code fully written in
Python and dedicated to computing eigenfunctions of stellar modes of
oscillations considering an equatorial-channel approximation of
rotation. The following tutorial is aimed at introducing the main
functionalities of the module.

.. code:: ipython3

    import mascaret as em
    import matplotlib.pyplot as plt
    em.__version__




.. parsed-literal::

    '0.2'



After importing the module, we can start by setting some parameters for
the mode exploration we will carry out.

.. code:: ipython3

    # Solver to use, that will be passed to scipy.integrate.solve_ivp
    em.solver_kwargs["method"] = "RK23" #default solver is RK45
    em.solver_kwargs["max_step"] = 1e-3
    # Number of nodes to use for the parallelised step
    nodes = 8
    # Grid search method, frequency interval to explore, and number of
    # grid points to consider
    method = "inverse"
    omega_0, omega_1 = 20, 30 #(muHz)
    n_grid = 300
    # Angular degree and azimuthal number of the modes to look for
    l_m = {1:[1]}

The module includes an example stellar profile of a 1
:math:`\rm M_\odot` star that was computed with
`MESA <https://docs.mesastar.org/en/latest/installation.html>`__.

.. code:: ipython3

    model = em.load_example_profile ()

In the first steps of this tutorial, we will compute the properties of
dipolar (:math:`\ell=1`) gravity modes (g modes) in the example stellar
profile. The Brunt-Väisälä frequency, :math:`N`, is the key physical
quantity to consider when dealing with these modes, which will be
propagative only when :math:`N>0` (radiative stratified regions) and
evanescent otherwise (convective regions).

.. code:: ipython3

    fig, ax = plt.subplots (1, 1, figsize=(5,4))
    
    ax.plot (model.x, model.N, 
             lw=2, color="darkorange")
    ax.set_xlabel (r"$r/R_\star$")
    ax.set_ylabel (r"$N$ ($\mu$Hz)")
    ax.set_xlim (0, 1)
    ax.set_ylim (0, 500)




.. parsed-literal::

    (0.0, 500.0)




.. image:: quickstart_files/quickstart_8_1.png


A case without rotation
-----------------------

We start by considering a case without rotation.

.. code:: ipython3

    rot_solar_rate = 0

Depending on the chosen search method (simple or double shooting), it is
possible to define a :math:`\Delta (\omega)` function (with
:math:`\omega` the wave frequency) for which the eigenmodes of the
problem will correspond to :math:`\Delta (\omega) = 0`. The frequency
grid is thus scanned to search for occurrences where
:math:`\Delta (\omega)` changes sign, then the eigenfrequencies and
eigenfunctions are refined by performing a dichotomy search around the
found intervals.

.. code:: ipython3

    results = em.mode_search (model, 
                            l_m, 
                            omega_0=omega_0, 
                            omega_1=omega_1,
                            n_grid=n_grid, 
                            rot_solar_rate=rot_solar_rate, 
                            method=method, 
                            parallelise=True, 
                            nodes=nodes,
                            double_shooting=True, 
                            output_dir="quickstart",
                            start_full_span="stitch")


.. parsed-literal::

    Looking for l=1, m=1 modes.
    Iterating on grid chunks.


.. parsed-literal::

    9it [00:04,  1.96it/s]                                                                                                                                          

.. parsed-literal::

    11 zeros found


.. parsed-literal::

    100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 11/11 [00:09<00:00,  1.11it/s]


.. code:: ipython3

    results




.. raw:: html

    <div><i>Table length=11</i>
    <table id="table5796150352" class="table-striped table-bordered table-condensed">
    <thead><tr><th>n</th><th>n_g</th><th>n_p</th><th>l</th><th>m</th><th>nu</th><th>delta</th><th>transmission</th><th>rot_solar_rate</th></tr></thead>
    <thead><tr><th>int32</th><th>int32</th><th>int32</th><th>int32</th><th>int32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th></tr></thead>
    <tr><td>-33</td><td>-33</td><td>0</td><td>1</td><td>1</td><td>20.371603</td><td>0.0</td><td>-0.12487564</td><td>0.0</td></tr>
    <tr><td>-32</td><td>-32</td><td>0</td><td>1</td><td>1</td><td>20.995756</td><td>0.0</td><td>-0.13381477</td><td>0.0</td></tr>
    <tr><td>-31</td><td>-31</td><td>0</td><td>1</td><td>1</td><td>21.66339</td><td>0.0</td><td>-0.14416665</td><td>0.0</td></tr>
    <tr><td>-30</td><td>-30</td><td>0</td><td>1</td><td>1</td><td>22.372078</td><td>0.0</td><td>-0.15531035</td><td>0.0</td></tr>
    <tr><td>-29</td><td>-29</td><td>0</td><td>1</td><td>1</td><td>23.126274</td><td>0.0</td><td>-0.16792531</td><td>0.0</td></tr>
    <tr><td>-28</td><td>-28</td><td>0</td><td>1</td><td>1</td><td>23.937809</td><td>0.0</td><td>-0.181721</td><td>0.0</td></tr>
    <tr><td>-27</td><td>-27</td><td>0</td><td>1</td><td>1</td><td>24.80276</td><td>0.0</td><td>-0.19798696</td><td>0.0</td></tr>
    <tr><td>-26</td><td>-26</td><td>0</td><td>1</td><td>1</td><td>25.736795</td><td>0.0</td><td>-0.21623734</td><td>0.0</td></tr>
    <tr><td>-25</td><td>-25</td><td>0</td><td>1</td><td>1</td><td>26.738888</td><td>0.0</td><td>-0.23738325</td><td>0.0</td></tr>
    <tr><td>-24</td><td>-24</td><td>0</td><td>1</td><td>1</td><td>27.823713</td><td>0.0</td><td>-0.26373035</td><td>0.0</td></tr>
    <tr><td>-23</td><td>-23</td><td>0</td><td>1</td><td>1</td><td>28.997177</td><td>0.0</td><td>-0.2918132</td><td>0.0</td></tr>
    </table></div>



An important property of the gravity modes we are considering here is
that their period spacing is asymptotically constant (:math:`|n|\gg1`).

.. code:: ipython3

    delta_p = - np.diff (1e6/(results["nu"]*60))
    
    fig, ax = plt.subplots (1, 1)
    
    ax.scatter (results["nu"][1:], delta_p, marker="o", 
                color="darkorange", edgecolor="black")
    
    ax.set_xlabel (r"$\nu$ ($\mu$Hz)")
    ax.set_ylabel (r"$\Delta P$ (min)")
    
    ax.set_ylim (0, 30)




.. parsed-literal::

    (0.0, 30.0)




.. image:: quickstart_files/quickstart_16_1.png


Including rotation
------------------

Let us now consider the rotating case. Note that, in **mascaret**, the
stellar rotation rate has to be provided in unity of solar rotation
frequency (which is set to ``4.14e-7`` Hz, and can be accessed through
the ``mascaret.solar_rotation_frequency`` global variable). It is
possible to set a uniform rotation rate by passing ``rot_solar_rate`` as
a float; radial rotation profile are also accepted, in this case, the
variable has to be provided as an array.

.. code:: ipython3

    rot_solar_rate = 50

.. code:: ipython3

    results_rot = em.mode_search (model, 
                                l_m, 
                                omega_0=omega_0, 
                                omega_1=omega_1,
                                n_grid=n_grid, 
                                rot_solar_rate=rot_solar_rate, 
                                method=method, 
                                parallelise=True, 
                                nodes=nodes,
                                double_shooting=True, 
                                output_dir="quickstart",
                                start_full_span="stitch")


.. parsed-literal::

    Looking for l=1, m=1 modes.
    Iterating on grid chunks.


.. parsed-literal::

    9it [00:07,  1.27it/s]                                                                                                                                          

.. parsed-literal::

    12 zeros found


.. parsed-literal::

    
    100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 12/12 [00:11<00:00,  1.07it/s]


.. code:: ipython3

    results_rot




.. raw:: html

    <div><i>Table length=12</i>
    <table id="table5798433488" class="table-striped table-bordered table-condensed">
    <thead><tr><th>n</th><th>n_g</th><th>n_p</th><th>l</th><th>m</th><th>nu</th><th>delta</th><th>transmission</th><th>rot_solar_rate</th></tr></thead>
    <thead><tr><th>int32</th><th>int32</th><th>int32</th><th>int32</th><th>int32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th></tr></thead>
    <tr><td>-34</td><td>-34</td><td>0</td><td>1</td><td>1</td><td>20.181133</td><td>0.0</td><td>-0.17386039</td><td>50.0</td></tr>
    <tr><td>-33</td><td>-33</td><td>0</td><td>1</td><td>1</td><td>20.754906</td><td>0.0</td><td>-0.18290704</td><td>50.0</td></tr>
    <tr><td>-32</td><td>-32</td><td>0</td><td>1</td><td>1</td><td>21.36104</td><td>0.0</td><td>-0.19318701</td><td>50.0</td></tr>
    <tr><td>-31</td><td>-31</td><td>0</td><td>1</td><td>1</td><td>22.011572</td><td>0.0</td><td>-0.2050534</td><td>50.0</td></tr>
    <tr><td>-30</td><td>-30</td><td>0</td><td>1</td><td>1</td><td>22.706648</td><td>0.0</td><td>-0.21858346</td><td>50.0</td></tr>
    <tr><td>-29</td><td>-29</td><td>0</td><td>1</td><td>1</td><td>23.4481</td><td>0.0</td><td>-0.23327672</td><td>50.0</td></tr>
    <tr><td>-28</td><td>-28</td><td>0</td><td>1</td><td>1</td><td>24.246815</td><td>0.0</td><td>-0.24978416</td><td>50.0</td></tr>
    <tr><td>-27</td><td>-27</td><td>0</td><td>1</td><td>1</td><td>25.100496</td><td>0.0</td><td>-0.27011144</td><td>50.0</td></tr>
    <tr><td>-26</td><td>-26</td><td>0</td><td>1</td><td>1</td><td>26.025555</td><td>0.0</td><td>-0.2911211</td><td>50.0</td></tr>
    <tr><td>-25</td><td>-25</td><td>0</td><td>1</td><td>1</td><td>27.017998</td><td>0.0</td><td>-0.31664476</td><td>50.0</td></tr>
    <tr><td>-24</td><td>-24</td><td>0</td><td>1</td><td>1</td><td>28.09541</td><td>0.0</td><td>-0.34629923</td><td>50.0</td></tr>
    <tr><td>-23</td><td>-23</td><td>0</td><td>1</td><td>1</td><td>29.261377</td><td>0.0</td><td>-0.38139284</td><td>50.0</td></tr>
    </table></div>



.. code:: ipython3

    delta_p_rot = - np.diff (1e6/(results_rot["nu"]*60))
    
    fig, ax = plt.subplots (1, 1)
    
    ax.scatter (results["nu"][1:], delta_p, marker="o", 
                color="darkorange", edgecolor="black")
    ax.scatter (results_rot["nu"][1:], delta_p_rot, marker="o", 
                color="blue", edgecolor="black")
    
    ax.set_xlabel (r"$\nu$ ($\mu$Hz)")
    ax.set_ylabel (r"$\Delta P$ (min)")
    
    ax.set_ylim (0, 30)




.. parsed-literal::

    (0.0, 30.0)




.. image:: quickstart_files/quickstart_22_1.png


Comparing the eigenfunctions
----------------------------

.. code:: ipython3

    n = results["n"][0]
    mode = em.load_eigenfunction ("quickstart/nabs{}_n{}_l1_m1_0_solar_rot.h5".format (str (np.abs (n)).zfill (4), 
                                                                                         str (n).zfill (4)))
    mode_rot = em.load_eigenfunction ("quickstart/nabs{}_n{}_l1_m1_{}_solar_rot.h5".format (str (np.abs (n)).zfill (4), 
                                                                                         str (n).zfill (4), 
                                                                                             rot_solar_rate))

.. code:: ipython3

    fig = em.plot_eigenfunction (mode["x"], 
                                 mode["xi_r"]/model.R, 
                                 mode["W"], 
                                 axs=None, model=model, omega=None,
                                 rot_solar_rate=None,
                                 c1="darkgrey", c2="darkgrey") 
    axs = fig.get_axes ()
    _ = em.plot_eigenfunction (mode_rot["x"], 
                               mode_rot["xi_r"]/model.R, 
                               mode_rot["W"], 
                               axs=axs, model=None, omega=None,
                               rot_solar_rate=None,
                               c1="darkorange", c2="blue") 
    
    axs[0].set_ylim (-50, 50)




.. parsed-literal::

    (-50.0, 50.0)




.. image:: quickstart_files/quickstart_25_1.png


Eigenfunction of a truncated shell
----------------------------------

It is also possible to compute the eigenfunction considering an inner
boundary located upper than the centre. We can for example explore the
behaviour of the isolated convective envelope in order to take a look at
the properties of pure inertial modes. For simplicity, we choose to
force the buoyancy of the medium to remain neutral (:math:`N^2 = 0`).

.. code:: ipython3

    model = em.load_example_profile (convective_neutral=False,
                                     remove_atmosphere=True)

.. code:: ipython3

    rot_solar_rate = 1
    # Considering only the convective envelope
    x_span = (model.x[model.index2+1], 1)
    omega_0, omega_1 = 0.02, 0.1
    xf = 0.85
    method = "linear"

.. code:: ipython3

    results_truncated = em.mode_search (model, 
                                        l_m, 
                                        omega_0=omega_0, 
                                        omega_1=omega_1,
                                        n_grid=n_grid, 
                                        rot_solar_rate=rot_solar_rate, 
                                        method=method, 
                                        parallelise=True, 
                                        nodes=nodes,
                                        double_shooting=True, 
                                        output_dir="quickstart",
                                        start_full_span="stitch",
                                        x_span=x_span, 
                                        xf=xf)


.. parsed-literal::

    Looking for l=1, m=1 modes.
    Iterating on grid chunks.


.. parsed-literal::

    9it [00:04,  1.95it/s]                                                                                                                                          

.. parsed-literal::

    2 zeros found


.. parsed-literal::

    
    100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:01<00:00,  1.47it/s]


.. code:: ipython3

    results_truncated




.. raw:: html

    <div><i>Table length=2</i>
    <table id="table5799082832" class="table-striped table-bordered table-condensed">
    <thead><tr><th>n</th><th>n_g</th><th>n_p</th><th>l</th><th>m</th><th>nu</th><th>delta</th><th>transmission</th><th>rot_solar_rate</th></tr></thead>
    <thead><tr><th>int32</th><th>int32</th><th>int32</th><th>int32</th><th>int32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th></tr></thead>
    <tr><td>-3</td><td>-3</td><td>0</td><td>1</td><td>1</td><td>0.027354332</td><td>0.0</td><td>inf</td><td>1.0</td></tr>
    <tr><td>-2</td><td>-2</td><td>0</td><td>1</td><td>1</td><td>0.0662436</td><td>0.0</td><td>-inf</td><td>1.0</td></tr>
    </table></div>



.. code:: ipython3

    filename_1 = "quickstart/nabs0003_n-003_l1_m1_{}_solar_rot.h5".format (rot_solar_rate)
    filename_2 = "quickstart/nabs0002_n-002_l1_m1_{}_solar_rot.h5".format (rot_solar_rate)
    mode_truncated_1 = em.load_eigenfunction (filename_1)
    mode_truncated_2 = em.load_eigenfunction (filename_2)
    
    fig = em.plot_eigenfunction (mode_truncated_1["x"], 
                               mode_truncated_1["xi_r"]/np.abs (mode_truncated_1["xi_r"]).max(), 
                               mode_truncated_1["W"]/np.abs (mode_truncated_1["W"]).max(), 
                               model=None, omega=None,
                               rot_solar_rate=None,
                               c1="darkorange", c2="blue") 
    axs = fig.get_axes ()
    _ = em.plot_eigenfunction (mode_truncated_2["x"], 
                               mode_truncated_2["xi_r"]/np.abs (mode_truncated_2["xi_r"]).max(), 
                               mode_truncated_2["W"]/np.abs (mode_truncated_2["W"]).max(), 
                               axs=axs,
                               model=None, omega=None,
                               rot_solar_rate=None,
                               c1="darkorange", c2="blue", ls="--") 



.. image:: quickstart_files/quickstart_31_0.png