Quatum states of Carbon Monoxide and Fortrat Diagram
====================================================

We here see the quantum states of Carbon Monoxide. Also, to see how the
band head appears, we would like to plot the Fortrat diagram with a
given quantum number and interval. To do so, we use
optional_quantum_states=True option in api.MdbExomol.

.. code:: ipython3

    from exojax.spec import api
    emf='CO/12C-16O/Li2015'   
    mdb = api.MdbExomol(emf, None, optional_quantum_states=True)


.. parsed-literal::

    /home/kawahara/exojax/src/exojax/utils/molname.py:133: FutureWarning: e2s will be replaced to exact_molname_exomol_to_simple_molname.
      warnings.warn(
    /home/kawahara/exojax/src/exojax/spec/api.py:153: UserWarning: nurange=None. Nonactive mode.
      warnings.warn("nurange=None. Nonactive mode.", UserWarning)


.. parsed-literal::

    HITRAN exact name= (12C)(16O)
    Background atmosphere:  H2
    Reading CO/12C-16O/Li2015/12C-16O__Li2015.trans.bz2
    DataFrame (self.df) available.


Check DataFrame. We see Li2015 contains the vibrational states for lower
and upper states, v_l, v_u.

.. code:: ipython3

    mdb.df[0:2]




.. raw:: html

    <table>
    <thead>
    <tr><th>#                            </th><th style="text-align: right;">  i_upper</th><th style="text-align: right;">  i_lower</th><th style="text-align: right;">        A</th><th style="text-align: right;">  nu_lines</th><th style="text-align: right;">  gup</th><th style="text-align: right;">  jlower</th><th style="text-align: right;">  jupper</th><th style="text-align: right;">  elower</th><th style="text-align: right;">  v_l</th><th style="text-align: right;">  v_u</th><th>kp_l  </th><th>kp_u  </th><th style="text-align: right;">        Sij0</th></tr>
    </thead>
    <tbody>
    <tr><td><i style='opacity: 0.6'>0</i></td><td style="text-align: right;">       84</td><td style="text-align: right;">       42</td><td style="text-align: right;">1.155e-06</td><td style="text-align: right;">   2.40559</td><td style="text-align: right;">    3</td><td style="text-align: right;">       0</td><td style="text-align: right;">       1</td><td style="text-align: right;"> 66960.7</td><td style="text-align: right;">   41</td><td style="text-align: right;">   41</td><td>e     </td><td>e     </td><td style="text-align: right;">3.81197e-164</td></tr>
    <tr><td><i style='opacity: 0.6'>1</i></td><td style="text-align: right;">       83</td><td style="text-align: right;">       41</td><td style="text-align: right;">1.161e-06</td><td style="text-align: right;">   2.44177</td><td style="text-align: right;">    3</td><td style="text-align: right;">       0</td><td style="text-align: right;">       1</td><td style="text-align: right;"> 65819.9</td><td style="text-align: right;">   40</td><td style="text-align: right;">   40</td><td>e     </td><td>e     </td><td style="text-align: right;">9.66303e-162</td></tr>
    </tbody>
    </table>



The Rovib transition changes both rotational and vibrational quantum
states. We here investigate the vibrational quantum state :math:`\nu`.
Let’s check how many :math:`\Delta \nu` Li2015 database contains:

.. code:: ipython3

    import numpy as np
    dv = mdb.df["v_u"]-mdb.df["v_l"]
    np.unique(dv.values)




.. parsed-literal::

    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])



So, we have 12 different :math:`\Delta \nu`. Let’s plot them.

.. code:: ipython3

    import matplotlib.pyplot as plt
    
    fig = plt.figure(figsize=(15, 4))
    ax = fig.add_subplot(111)
    for i, udv in enumerate(np.unique(dv.values)):
        mask = dv == udv
        mdf = mdb.df[mask]
        ax.plot(mdf["nu_lines"].values,
                 mdf["Sij0"].values,
                 ".",
                 alpha=0.3,
                 #alpha=0.01 + 0.005 * i,
                 color="gray")
        ax.text(
            np.sum(mdf["nu_lines"].values * mdf["Sij0"].values) /
            np.sum(mdf["Sij0"].values), 1.e2*np.max(mdf["Sij0"].values),"$\\Delta \\nu=$"+str(udv))
    
    for mic in [0.5,1,2,3,4,5,10,20]:
        x = 1.e4/mic
        plt.axvline(x,alpha=0.2,color="gray")
        #plt.text(x,1.e-210,str(mic)+" $\\mu$m",rotation="90")
        plt.text(x,1.e-39,str(mic)+" $\\mu$m",rotation="90")
    plt.yscale("log")
    plt.ylim(1.e-41,1.e-13)
    plt.tick_params(labelsize=14)
    plt.xlabel("wavenumber (cm-1)",fontsize=14)
    plt.ylabel("line strength",fontsize=14)
    plt.savefig("co_dnu.png", bbox_inches="tight", pad_inches=0.1)
    plt.show()



.. image:: Fortrat_files/Fortrat_8_0.png


Let’s go deeper! Expand this for :math:`\Delta \nu=2` (K-band feature).

.. code:: ipython3

    dv = mdb.df["v_u"] - mdb.df["v_l"]
    dJ = mdb.df["jupper"] - mdb.df["jlower"]
    
    fig = plt.figure(figsize=(10, 5))
    for i, vl in enumerate(np.unique(mdb.df["v_l"].values)):
        mask = (dv == 2) * (dJ == 1) * (mdb.df["v_l"] == vl)
        vdf = mdb.df[mask]
        plt.plot(vdf["nu_lines"].values, vdf["Sij0"].values, ".", color="black")
        if i < 10:
            plt.text(np.nanmean(vdf["nu_lines"].values),
                     8 * np.nanmax(vdf["Sij0"].values),
                     "$\\nu_{lower}=$" + str(vl),
                     fontsize=12)
        mask = (dv == 2) * (dJ == -1) * (mdb.df["v_l"] == vl)
        vdf = mdb.df[mask]
        plt.plot(vdf["nu_lines"].values, vdf["Sij0"].values, ".", color="gray")
    
    for mic in [2.3, 2.5, 2.7]:
        x = 1.e4 / mic
        plt.axvline(x, alpha=0.2, color="gray")
        #plt.text(x,1.e-210,str(mic)+" $\\mu$m",rotation="90")
        plt.text(x, 1.e-60, str(mic) + " $\\mu$m", rotation="90")
    
    plt.text(3800.0,
             1.e-25,
             "$\\Delta J$ = -1, P-branch",
             color="gray",
             fontsize=14)
    plt.text(4380.0,
             1.e-25,
             "$\\Delta J$ = 1, R-branch",
             color="black",
             fontsize=14)
    plt.yscale("log")
    plt.ylim(1.e-61, 1.e-13)
    plt.xlim(3500, 4620)
    plt.tick_params(labelsize=14)
    plt.xlabel("wavenumber (cm-1)", fontsize=14)
    plt.ylabel("line strength", fontsize=14)
    plt.savefig("co_dnu_expand.png", bbox_inches="tight", pad_inches=0.1)
    plt.show()



.. image:: Fortrat_files/Fortrat_10_0.png


Using DataFrame, we pick up the lines with :math:`\Delta \nu = 2`,
:math:`\Delta J = \pm 1` (R, P-branch), and :math:`\nu = 0` here.

.. code:: ipython3

    dv = mdb.df["v_u"]-mdb.df["v_l"]
    dJ = mdb.df["jupper"] - mdb.df["jlower"]
    vmask = mdb.df["v_l"] == 0
    mask_R = (dv == 2) * (dJ == 1) * vmask
    mask_P = (dv == 2) * (dJ == -1) * vmask
    
    df_R = mdb.df[mask_R]
    df_P = mdb.df[mask_P]


Let’s plot the Fortrat diagram. The y-axis of the Fortart diagram is
:math:`J_\mathrm{upper}` for R-branch and :math:`- J_\mathrm{lower}` for
P-branch.

.. code:: ipython3

    import matplotlib.pyplot as plt
    fig = plt.figure(figsize=(15,5))
    plt.plot(df_R["nu_lines"].values,df_R["jupper"].values,"+",alpha=0.8, color="black",label="R-branch")
    plt.plot(df_P["nu_lines"].values,- df_P["jupper"].values,"*",alpha=0.8, color="gray",label="P-branch")
    plt.tick_params(labelsize=16)
    plt.xlabel("wavenumber (cm-1)", fontsize=18)
    plt.ylabel("$\\mathcal{J}$", fontsize=18)
    plt.legend(fontsize=16)
    plt.title(emf+" $\\Delta \\nu = 2, \\nu_\\mathrm{lower} = 0$",fontsize=18)
    plt.xlim(4000.,4400)
    plt.ylim(-50,150)
    plt.savefig("fortrat.png",  bbox_inches="tight", pad_inches=0.1)
    plt.show()



.. image:: Fortrat_files/Fortrat_14_0.png