.. container::
   :name: toc

Forward modeling of the emission spectrum using VALD3
=====================================================

| Tako Ishikawa, Hajime Kawahara
| created: : 2021/07/20,　　last update: 2022/10/22

.. raw:: html

   <!-- 
   written with reference to :  
   "exojax/examples/tutorial/Forward\ modeling.ipynb"  
   "ghR/exojax_0/examples/testlines/line_strength_CO.py"  

   cd ~/work

   -->

This example provides how to use VALD3 for forward modeling of the
emission spectrum. Currenty, we use exojax.spec.moldb as API for VALD3
because we have not implemented the VALD3 API in radis.api yet. Someday,
we (or someone who is reading this document.) would include the VALD3
API in radis.api!

.. code:: ipython3

    from exojax.utils.grids import wavenumber_grid
    from exojax.spec.rtransfer import pressure_layer 
    from exojax.spec import moldb, molinfo, contdb
    from exojax.spec import atomll
    from exojax.spec.exomol import gamma_exomol
    from exojax.spec import SijT, doppler_sigma
    from exojax.spec import planck
    import matplotlib.pyplot as plt
    import jax.numpy as jnp
    from jax import vmap, jit
    import numpy as np

T-P profile
-----------

.. code:: ipython3

    #Assume ATMOSPHERE                                                                     
    NP=100
    T0=3000. #10000. #3000. #1295.0 #K
    Parr, dParr, k=pressure_layer(NP=NP)
    H_He_HH_VMR = [0.0, 0.16, 0.84] #typical quasi-"solar-fraction"
    Tarr = T0*(Parr)**0.1
    
    PH = Parr* H_He_HH_VMR[0]
    PHe = Parr* H_He_HH_VMR[1]
    PHH = Parr* H_He_HH_VMR[2]
    
    fig=plt.figure(figsize=(6,4))
    plt.plot(Tarr,Parr)
    plt.plot(Tarr, PH, '--'); plt.plot(Tarr, PHH, '--'); plt.plot(Tarr, PHe, '--')
    plt.plot(Tarr[80],Parr[80], marker='*', markersize=15)
    plt.yscale("log")
    plt.xlabel("temperature (K)")
    plt.ylabel("pressure (bar)")
    plt.gca().invert_yaxis()
    plt.show()



.. image:: Forward_modeling_for_metal_line_files/Forward_modeling_for_metal_line_6_0.png


Wavenumber
----------

.. code:: ipython3

    #We set a wavenumber grid using wavenumber_grid.
    nus,wav,res = wavenumber_grid(10380, 10430, 4500, unit="AA") 


.. parsed-literal::

    xsmode assumes ESLOG in wavenumber space: mode=lpf


Load a database of atomic lines from VALD3
------------------------------------------

.. code:: ipython3

    #Loading a database of a few atomic lines from VALD3  #BU: CO and CIA (H2-H2)... 
    """
        valdlines:  fullpath to the input line list obtained from VALD3 (http://vald.astro.uu.se/):
                VALD data access is free but requires registration through the Contact form (http://vald.astro.uu.se/~vald/php/vald.php?docpage=contact.html). 
                After the registration, you can login and select one of the following modes depending on your purpose: "Extract All", "Extract Stellar", or "Extract Element".
            For a example in this notebook, the request form of "Extract All" mode was filled as:
              Extract All
                Starting wavelength :    10380
                Ending wavelength :    10430
                Extraction format :    Long format
                Retrieve data via :    FTP
                (Hyperfine structure:    N/A)
                (Require lines to have a known value of :    N/A)
                Linelist configuration :    Default
                Unit selection:    Energy unit: eV - Medium: vacuum - Wavelength unit: angstrom - VdW syntax: default
            Please assign the fullpath of the output file sent by VALD ([user_name_at_VALD].[request_number_at_VALD].gz;  "vald2600.gz" in the code below) to the variable "valdlines".
            Note that the number of spectral lines that can be extracted in a single request is limited to 1000 in VALD (https://www.astro.uu.se/valdwiki/Restrictions%20on%20extraction%20size).
    """
    
    valdlines = '.database/HiroyukiIshikawa.4214450.gz'
    adbFe = moldb.AdbVald(valdlines, nus)



.. parsed-literal::

    Reading VALD file


Relative partition function
---------------------------

.. code:: ipython3

    #Computing the relative partition function,
    
    qt_284 = vmap(adbFe.QT_interp_284)(Tarr)
    qt = np.zeros([len(adbFe.QTmask), len(Tarr)])
    for i, mask in enumerate(adbFe.QTmask):
        qt[i] = qt_284[:, mask]  #e.g., qt_284[:,76] #Fe I
    qt = jnp.array(qt)


Pressure and Natural broadenings (Lorentzian width)
---------------------------------------------------

.. code:: ipython3

    gammaLMP = jit(vmap(atomll.gamma_vald3,(0,0,0,0,None,None,None,None,None,None,None,None,None,None,None)))\
            (Tarr, PH, PHH, PHe, adbFe.ielem, adbFe.iion, \
                    adbFe.dev_nu_lines, adbFe.elower, adbFe.eupper, adbFe.atomicmass, adbFe.ionE, \
                    adbFe.gamRad, adbFe.gamSta, adbFe.vdWdamp, 1.0)  

Doppler broadening
------------------

.. code:: ipython3

    sigmaDM=jit(vmap(doppler_sigma,(None,0,None)))\
        (adbFe.nu_lines, Tarr, adbFe.atomicmass)

Line strength
-------------

.. code:: ipython3

    SijM=jit(vmap(SijT,(0,None,None,None,0)))\
        (Tarr, adbFe.logsij0, adbFe.nu_lines, adbFe.elower, qt.T)

nu matrix
---------

.. code:: ipython3

    from exojax.spec.initspec import init_lpf
    numatrix=init_lpf(adbFe.nu_lines,nus)

Compute dtau for each atomic species (or ion) in a SEPARATE array
-----------------------------------------------------------------

Separate species

.. code:: ipython3

    def get_unique_list(seq):
        seen = []
        return [x for x in seq if x not in seen and not seen.append(x)]
    
    uspecies = get_unique_list(jnp.vstack([adbFe.ielem, adbFe.iion]).T.tolist())

Set the stellar/planetary parameters

.. code:: ipython3

    #Parameters of Objects
    Rp = 0.36*10 #R_sun*10    #Rp=0.88 #[R_jup]
    Mp = 0.37*1e3 #M_sun*1e3    #Mp=33.2 #[M_jup]
    g = 2478.57730044555*Mp/Rp**2
    print('logg: '+str(np.log10(g))) #check


.. parsed-literal::

    logg: 4.849799190511717


Calculate delta tau

.. code:: ipython3

    #For now, ASSUME all atoms exist as neutral atoms.
    #In fact, we can't ignore the effect of molecular formation e.g. TiO (」゜□゜)」
    
    from exojax.spec.lpf import xsmatrix
    from exojax.spec.rtransfer import dtauM
    from exojax.spec.atomllapi import load_atomicdata
    
    ipccd = load_atomicdata()
    ieleml = jnp.array(ipccd['ielem'])
    Narr = jnp.array(10**(12 + ipccd['solarA']))  #number density
    massarr = jnp.array(ipccd['mass'])  #mass of each neutral atom
    Nmassarr = Narr * massarr  #mass of each neutral species
    
    dtaual = np.zeros([len(uspecies), len(Tarr), len(nus)])
    maskl = np.zeros(len(uspecies)).tolist()
    
    for i, sp in enumerate(uspecies):
        maskl[i] = (adbFe.ielem==sp[0])\
                        *(adbFe.iion==sp[1])
    
        #Currently not dealing with ionized species yet... (#tako %\\\\20210814)
        if sp[1] > 1:
            continue
    
        #Providing numatrix, thermal broadening, gamma, and line strength, we can compute cross section.
        xsm = xsmatrix(numatrix[maskl[i]], sigmaDM.T[maskl[i]].T,
                       gammaLMP.T[maskl[i]].T, SijM.T[maskl[i]].T)
        #Computing delta tau for atomic absorption
        MMR_X_I = Nmassarr[jnp.where(ieleml == sp[0])[0][0]] / jnp.sum(Nmassarr)
        mass_X_I = massarr[jnp.where(ieleml == sp[0])[0][
            0]]  #MMR and mass of neutral atom X (if all elemental species are neutral)
        dtaual[i] = dtauM(dParr, xsm, MMR_X_I * np.ones_like(Tarr), mass_X_I, g)


compute delta tau for CIA

.. code:: ipython3

    cdbH2H2=contdb.CdbCIA('.database/H2-H2_2011.cia', nus)
    
    from exojax.spec.rtransfer import dtauCIA
    mmw=2.33 #mean molecular weight
    mmrH2=0.74
    molmassH2=molinfo.molmass("H2")
    vmrH2=(mmrH2*mmw/molmassH2) #VMR
    dtaucH2H2=dtauCIA(nus,Tarr,Parr,dParr,vmrH2,vmrH2,\
                mmw,g,cdbH2H2.nucia,cdbH2H2.tcia,cdbH2H2.logac)


.. parsed-literal::

    H2-H2


Total delta tau
---------------

.. code:: ipython3

    dtau = np.sum(dtaual, axis=0) + dtaucH2H2

Plot contribution function
--------------------------

.. code:: ipython3

    from exojax.plot.atmplot import plotcf
    plotcf(nus,dtau,Tarr,Parr,dParr)
    plt.show()



.. image:: Forward_modeling_for_metal_line_files/Forward_modeling_for_metal_line_33_0.png


Radiative transfer
------------------

.. code:: ipython3

    from exojax.spec import planck
    from exojax.spec.rtransfer import rtrun
    sourcef = planck.piBarr(Tarr, nus)
    F0=rtrun(dtau, sourcef)

.. code:: ipython3

    fig=plt.figure(figsize=(5, 3))
    plt.plot(wav[::-1],F0)
    plt.show()



.. image:: Forward_modeling_for_metal_line_files/Forward_modeling_for_metal_line_36_0.png


.. code:: ipython3

    #Check line species
    print(np.unique(adbFe.ielem))


.. parsed-literal::

    [12 13 14 17 18 20 21 22 24 25 26 27 28 29 32 38 59 64 65 66 70 90]


Rotational & instrumental broadening
------------------------------------

.. code:: ipython3

    from exojax.spec import response
    from exojax.utils.constants import c #[km/s]
    import jax.numpy as jnp
    
    wavd=jnp.linspace(10380, 10450,500) #observational wavelength grid
    nusd = 1.e8/wavd[::-1]
    
    RV=10.0 #RV km/s
    vsini=20.0 #Vsini km/s
    u1=0.0 #limb darkening u1
    u2=0.0 #limb darkening u2
    
    R=100000.
    beta=c/(2.0*np.sqrt(2.0*np.log(2.0))*R) #IP sigma need check 
    
    Frot=response.rigidrot(nus,F0,vsini,u1,u2)
    F=response.ipgauss_sampling(nusd,nus,Frot,beta,RV)

.. code:: ipython3

    fig=plt.figure(figsize=(5, 3))
    plt.plot(wav[::-1],F0, label='F0')
    plt.plot(wavd[::-1],F, label='F')
    plt.legend()
    plt.show()



.. image:: Forward_modeling_for_metal_line_files/Forward_modeling_for_metal_line_40_0.png