.. _tutorial-non-lte:

myXCLASSFit with non-LTE
========================

The purpose of this tutorial is to describe how to model the contribution
of a molecule in Non-LTE, see Sect. ":ref:`myxclass-NonLTE-molecules`"
using the *myXCLASSFit* function, see Sect. ":ref:`myxclassfit-short`".

All files used within this tutorial can be downloaded
`here <https://xclass.astro.uni-koeln.de/sites/xclass/files/software/Tutorial/XCLASS-Tutorial-NLTE.zip>`_.



.. _tutorial-non-lte-start:

Introduction
------------

Here, we want to model the contribution of CH\ :math:`_3` \ OH\ :math:`_{v=0}` A
and take the contribution of CH\ :math:`_3` \ OH\ :math:`_{v=0}` E into
account as well. We use observational data of the LMC.


.. _tutorial-non-lte-obs-xml:

Observational xml file
----------------------

For the *myXCLASSFit* function we have to define the frequency range,
(here from 241.45 to 241.75 GHz), the size of the interferometric
beam (here 0.14 arcsec), the background and dust parameters, respectively.

**Observational xml file** ``tutorial-nlte__obs__LMC.xml`` used for the
*myXCLASSFit* function:

::

    <?xml version="1.0" encoding="UTF-8"?>
    <ExpFiles>


       <!-- define number of observation files -->
       <NumberExpFiles>1</NumberExpFiles>


       <!-- ********************************************************************** -->
       <!-- define parameters for 1st observation data file -->
       <file>


          <!-- define path and name of first data file -->
          <FileNamesExpFiles>files/data/CH3OH_WS_px_1050_1039.dat</FileNamesExpFiles>


          <!-- define import filter -->
          <ImportFilter>xclassASCII</ImportFilter>


          <!-- define number of data ranges -->
          <NumberExpRanges>1</NumberExpRanges>


          <!-- define parameters for each data ranges -->
          <FrequencyRange>
                <MinExpRange>2.414511573992605845e+05</MinExpRange>
                <MaxExpRange>2.417543671618033550e+05</MaxExpRange>
                <StepFrequency>1.2128390501893591</StepFrequency>


                <!-- define background temperature and temperature slope -->
                <t_back_flag>True</t_back_flag>
                <BackgroundTemperature>0.0</BackgroundTemperature>
                <TemperatureSlope>0.0</TemperatureSlope>


                <!-- define hydrogen column density, beta for dust, and kappa -->
                <HydrogenColumnDensity>1.4e+19</HydrogenColumnDensity>
                <DustBeta>1.4</DustBeta>
                <Kappa>0.4</Kappa>
          </FrequencyRange>


          <!-- define size of telescope -->
          <TelescopeSize>0.14</TelescopeSize>


          <!-- define if interferometric observation is modeled -->
          <Inter_Flag>True</Inter_Flag>


          <!-- define local standard of rest (vLSR) -->
          <GlobalvLSR>233.0</GlobalvLSR>


          <!-- define parameters for ASCII file import -->
          <ErrorY>no</ErrorY>
          <NumberHeaderLines>0</NumberHeaderLines>
          <SeparatorColumns> </SeparatorColumns>
       </file>


       <!-- *********************************** -->
       <!-- taken local-overlap into account? -->
       <LocalOverlap_Flag>False</LocalOverlap_Flag>


       <!-- *************************************************************** -->
       <!-- parameters for isotopologues -->
       <iso_flag>True</iso_flag>
       <IsoTableFileName>files/tutorial-nlte__iso-ratios.dat</IsoTableFileName>


       <!-- *************************************************************************************** -->
       <!-- define file describing collision partners -->
       <CollisionPartnerFileName>files/tutorial-nlte__collision-partners.dat</CollisionPartnerFileName>
    </ExpFiles>


In contrast to the LTE description, the Non-LTE description requires
the collision rates for the corresponding molecules and the density of
the collision partners, which are described in the collision partner
file, see Sect. ":ref:`myxclass-collpart`". In the aforementioned obs.
xml file, the path and name of this collision partner file is described
by tag ``<CollisionPartnerFileName>``.



.. _tutorial-non-lte-molfit:

Molfit file
-----------

Methanol is the simplest asymmetric-top molecule with one hindered torsion
and eleven normal vibration modes. The rotational and vibration-rotational
spectra of methanol are complicated by various interactions among torsion,
rotation, and vibrations, which have been extensively studied for many
years in microwave, far-infrared, and infrared frequency regions.

It has a symmetric group and a locally rotated axis (i.e., CH\ :math:`_3` \
axis) of symmetry, it rotates around this methyl group axis with internal
rotation or torsion of the CH\ :math:`_3` \ group relative to the OH group. The
identical nuclei in the CH\ :math:`_3` \ group give rise to the existence
of a specific type of nuclear spin isomers in methanol. It is from the
three spin-1/2 hydrogen nuclei of the CH\ :math:`_3` \ group that the
ortho-CH\ :math:`_3`\ OH (I = 3/2) and para-CH\ :math:`_3`\ OH
(I = 1/2) are formed in methanol, which are distinguished respectively
by the symmetry quantum numbers \ :math:`\sigma`\  = 0 and ± 1.
The ortho-CH\ :math:`_3`\ OH isomer corresponds to the symmetry quantum
number \ :math:`\sigma`\ = 0. Accordingly for the para-CH\ :math:`_3`\ OH
isomer to \ :math:`\sigma`\ = ±1. Here each \ :math:`\sigma`\ value
corresponds to a torsional-symmetry species of A-species
(\ :math:`\sigma`\ = 0, CH\ :math:`_3`\ OH\ :math:`_{v=0}` A) and
E-species (\ :math:`\sigma`\ = ±1, CH\ :math:`_3`\ OH\ :math:`_{v=0}` E).

The relative abundance of the ortho and para nuclear spin isomers in
astrophysical environments can provide insights into the temperature
and physical conditions of the observed sources. Ortho-methanol is the
more abundant nuclear spin isomer at low temperatures due to its lower
energy state. Para-methanol is the more abundant nuclear spin isomer
at high temperatures due to its higher energy state.

In order to describe a component of a certain molecule in Non-LTE, the
keyword ``Non-LTE`` has to be added on the right side of the column
describing the ordering (distance) parameter of the corresponding
component.

The **molfit file** ``tutorial-non-lte__Band_6.molfit`` used for the
*myXCLASSFit* function:

::

   %=======================================================================================================================================================================================================================================================================
   %
   % define database parameters:
   %
   %=======================================================================================================================================================================================================================================================================
   %%MinNumTransitionsSQL = 1                % (min. number of transitions)
   %%MaxNumTransitionsSQL = 90000            % (max. number of transitions)
   %%MaxElowSQL = 8.000e+03                  % (max. lower energy, i.e. upper limit of the lower energy of a transition)
   %%MingASQL = 1.000e-09                    % (minimal intensity, i.e. lower limit of gA of a transition)
   %%TransOrderSQL = 1                       % (order of transitions: (=1): by lower energy, (=2) by gA, (=3) gA/E_low^2, else trans. freq.)



   %=======================================================================================================================================================================================================================================================================
   %
   % define molecules and their components:
   %
   %=======================================================================================================================================================================================================================================================================
   CH3OH;v=0;A   1
               n        1.00     1000.00         500.00000    y       10.00      200.00          30.00000    y    1.00e+12    1.00e+21       8.00000e+15    y        0.10        9.00           2.00000    y     -250.00      100.00           0.00000   c      Non-LTE



.. _tutorial-non-lte-iso-ratio:

Iso ratio file
--------------

Here, the iso ratio file is used to describe the ratio between
CH\ :math:`_3`\ OH\ :math:`_{v=0}` A and CH\ :math:`_3`\ OH\ :math:`_{v=0}` E.

The **iso-ratio file** ``tutorial-nlte__iso-ratios.dat`` used for the
*myXCLASSFit* function:

::

   % Isotopologue:      Iso-master:        Iso-ratio:      Lower-limit:      Upper-limit:       collision rate file:
   CH3OH;v=0;E          CH3OH;v=0;A               1.0               1.0               1.0       files/e-ch3oh.dat

Here, the last column describes the path and name of the file describing
the collision rates of CH\ :math:`_3`\ OH\ :math:`_{v=0}` E.



.. _tutorial-non-lte-coll-part:

Collision partner file
----------------------

The collision partner file is used, among other things, to describe the name
and density of the collision partners, see Sect. ":ref:`myxclass-collpart`".


The **collision partner file** ``tutorial-nlte__collision-partners.dat`` used
for the *myXCLASSFit* function:

::

   %---------------------------------------------------------------------------------------------------------------------------------------------------------
   %
   % input file for collision partners
   %
   % densities can be fitted by declaring lower and upper limits, if both limits
   % are identical, corresponding density is not fitted.
   %
   % comments are marked with "%" sign
   %
   %
   %---------------------------------------------------------------------------------------------------------------------------------------------------------
   1                                                                                       % Geometry (=1: sphere, =2: LVG, =3: slab)
   1                                                                                       % number of molecules
   CH3OH;v=0;A                                                                             % molecule
   files/tutorial-nlte__coll-rates__a-ch3oh.dat                                            % path and name of file containing collision rates
   1                                                                                       % number of collision partners
   H2                                                                                      % first collision partner
   5.2E+03            1.e2            5.e5                                                 % density of first collision partner

Here, the initial density of the collision partner H\ :math:`_2` \ is set to
:math:`5.2 \cdot 10^{3}` cm\ :math:`^{-3}`.



.. _tutorial-non-lte-alg-xml:

Algorithm xml file
------------------

In order to fit the parameters described in the molfit, iso ratio, and collision
partner files, we have to define the optimization algorithm(s) which is (are)
used by the *myXCLASSFit* function. Here, we use the **trf** algorithm, where we
want to limit the number of iterations to 10 by using the tag
``"<number_iterations>"``. Furthermore, we define another stopping criterion
by setting tag the upper limit of the :math:`\chi^2` function to
:math:`1 \cdot 10^{-15}`, using the tag ``"<limit_of_chi2>"``. Additionally, we
want to use 14 processors cores, which is described by tag ``"<NumberProcessors>"``.

**Algorithm xml file** ``tutorial-nlte__algorithm.xml`` used for the
*myXCLASSFit* function:

::

   <?xml version="1.0" encoding="UTF-8"?>
   <FitControl>
      <!-- settings for fit process -->


      <!-- set number of used algorithms -->
      <NumberOfFitAlgorithms>1</NumberOfFitAlgorithms>


      <!-- define algorithm -->
      <algorithm>
         <FitAlgorithm>trf</FitAlgorithm>


         <!-- define value of the variation -->
         <VariationValue>1.e-4</VariationValue>


         <!-- set max. number of iterations -->
         <number_iterations>10</number_iterations>


         <!-- set max. number of processors -->
         <NumberProcessors>14</NumberProcessors>


         <!-- set path and name of host file -->
         <MPIHostFileName> </MPIHostFileName>


         <!-- settings for chi^2 -->
         <limit_of_chi2>1e-15</limit_of_chi2>
         <RenormalizedChi2>yes</RenormalizedChi2>
         <DeterminationChi2>default</DeterminationChi2>
         <SaveChi2>yes</SaveChi2>


         <!-- set plot options -->
         <PlotAxisX>Rest Frequency [MHz]</PlotAxisX>
         <PlotAxisY>T_mb [K]</PlotAxisY>
         <PlotIteration>no</PlotIteration>
      </algorithm>
   </FitControl>



.. _tutorial-non-lte-call:

Call of myXCLASSFit function
----------------------------

Now everything is prepared to start the *myXCLASSFit* function.

Start **myXCLASSFit** function

::

    >>> from xclass import task_myXCLASSFit
    >>> import os

    ## get path of current directory
    >>> LocalPath = os.getcwd() + "/"

    # import XCLASS packages
    >>> from xclass import task_myXCLASSFit

    # define path and name of molfit file
    >>> MolfitsFileName = LocalPath + "files/tutorial-nlte__LMC.molfit"

    # define path and name of obs. xml file
    >>> ObsXMLFileName = LocalPath + "files/tutorial-nlte__obs__LMC.xml"

    # define path and name of algorithm xml file
    >>> AlgorithmXMLFileName = LocalPath + "files/tutorial-nlte__algorithm.xml"

    # set optimization package
    >>> Optimizer = "scipy"

    # define if results of myXCLASSFit function are return as dictionary
    >>> DictOutFlag = True

    # call myXCLASSFit function
    >>> OutputDictionary = task_myXCLASSFit.myXCLASSFitCore( \
                                 MolfitsFileName = MolfitsFileName, \
                                 AlgorithmXMLFileName = AlgorithmXMLFileName, \
                                 ObsXMLFileName = ObsXMLFileName, \
                                 Optimizer = Optimizer, \
                                 DictOutFlag = DictOutFlag)


.. _tutorial-non-lte-results:

Results
-------

.. figure:: ../figures/tutorial-nlte__final-plot.png
   :align: center
   :width: 100%
   :name: fig-tutorial-non-lte-results_plot

   Fitted spectrum (left panel) and corresponding :math:`\chi^2` function
   (right panel) created by the myXCLASSFit function.

After finishing we get the following molfit file

::

   %=========================================================================================================================================================
   %
   % define database parameters:
   %
   %=========================================================================================================================================================
   %%MinNumTransitionsSQL = 1                % (min. number of transitions)
   %%MaxNumTransitionsSQL = 90000            % (max. number of transitions)
   %%MaxElowSQL = 8.000e+03                  % (max. lower energy, i.e. upper limit of the lower energy of a transition)
   %%MingASQL = 1.000e-09                    % (minimal intensity, i.e. lower limit of gA of a transition)
   %%TransOrderSQL = 1                       % (order of transitions: (=1): by Elow, (=2) by gA, (=3) gA/E_low^2, else trans. freq.)


   %=========================================================================================================================================================
   %
   % define molecules and their components:
   %
   %=========================================================================================================================================================
   CH3OH;v=0;A  1
      n        1.00     1000.00         500.00000    y       10.00      200.00          11.49368    y    1.00e+12    1.00e+21       9.90429e+15    y        0.10        9.00           1.96329    y     -250.00      100.00           0.62725      c      Non-LTE


Additionally, we get an density of :math:`8.53 \cdot 10^{4}` cm\ :math:`^{-3}`
for the collision partner H\ :math:`_2`.



.. ----------------------------------------------------------------------------------------



.. Footnotes
.. ---------

.. Footnotes