.. _spectra-example:

Spectra: comparison between plan-parallel and spherical
=======================================================

.. highlight:: python



In this example, we present how you can calculate a reflectance spectrum with
**htrdr-planets** from a set of 1D data. Here, we compare between plan-parallel
and spherical geometries.

Imports
-------

We start by importing the necessary libraries.::

    import numpy as np
    import matplotlib.pyplot as plt
    import htrdrPy as htrdr
    import scipy.constants as cst
    import json

Data
----

The first step is to generate an :py:class:`.Data` object that will contain the
optical properties. This object aims at handling the data in order to generate
the input files. In this example, the optical properties are loaded from an
**.npz** file that you can find `here`_. We first create an instance of the
:py:class:`.Data` object with the correct planet radius::

    radius = 2.575e6 # m

    dataSp= htrdr.Data(radius = radius, nTheta=80, nPhi=80, name="Sphere")
    dataPP = htrdr.Data(radius = 1e9, name="PP")

.. note::
   In plan parallel, the radius is used to define the x and y extension of the
   atmosphere. Make sure this value is large enough, but also note that a too
   large value can result in issues within **htrdr**. For this example, we found
   that any value larger the 1e9 would result in inconsistent results.

.. note::
   At the creation of each Data instance, an **inputs_{name}/** repository is
   created and will contain all the input files necessary for htrdr. A
   **outputs_{name}/** repository is also created and will contain the files
   created by **htrdr**.

.. _here: https://www.dropbox.com/scl/fo/qk20qvsvocy9w3zuik254/ADNYNbijPepUqOOgQGoQH78?rlkey=r95p5zzsmdkylihtyacfslgoa&st=pbtkmwxa&dl=0

Then, we read the data and convert it into a python disctionnary.  The
:py:class:`.Data` object will work with the dictionnary, which therefore needs
to have the correct keys. In this example, the keys used to generate the
**.npz** file are already the right ones. Each corresponds to a table with a
given shape.  All these informations (keys name and array shapes) are provided
further down in this documentation (c.f. :ref:`data`). ::

    data = np.load("Example_spectra_data/1D_spetral_data.npz")
    data = dict(data)
    for key, array in data.items():
        print(key, ", shape : ", array.shape)

Additional information, such as the dimensions (number of angles in the phase
function, number of wavelength, etc.) are also required to be present in the
dictionnary. We therefore add the necessary information::

    nWavelength, nLevel, nCoeff = data["absorption (m-1)"].shape
    nAngle = len(data["angles (°)"])
    wavelengths = data["wavelength"]

    data.update({
            'nLevel' : nLevel,
            'nWavelength' : nWavelength,
            'nCoeff' : nCoeff,
            'nAngle' : nAngle,
            })

    print(data.keys())

The next step is to provide those data to the :py:class:`.Data` object.
Depending on the context, different methods exist to pass the information. In
this example, we use the :py:meth:`.Data.makeMixture` method, which handles the
case where the data provided are those of the atmosphere mixture comprising all
the gases and aerosols. The **dim** parameter indicates the dimension of the
provided data: 1 means that we provide a column, 2 for a slice along altitude
latitude and 3 for a fully 3D heteorgeneous atmosphere. Obviously, the shapes of
the tables provided through the dictionnary change accordingly. **dim=0** sets a
plan-parallel atmosphere. ::

    dataSp.makeMixture(data, dim=1)
    dataPP.makeMixture(data, dim=0)

Then, we generate the surface by providing the temperature map (a single float
in 1D) and the **brdf** dictionnary to the :py:class:`.Data` instance through
the :py:meth:`.Data.makeGround` methods. The method to use depends on the
dimension of the input data. ::

    brdf = {
            "kind": 'lambertian',
            "albedo": data['surface albedo'],
            "wavelengths": wavelengths
            }
    dataPP.makeGroundFrom1D_PP(91, brdf)
    dataSp.makeGroundFrom1D(91, brdf)

Generating input files
----------------------

We can now generate the input files with the command
:py:meth:`.Data.writeInputs`::

    dataPP.writeInputs()
    dataSp.writeInputs()

We can also generate the VTK files with the command
:py:meth:`.Data.writeVTKfiles`. VTK files are not necessary and only serves to
visualize the data that will be send to **htrdr-planets**. ::

    dataPP.writeVTKfiles()
    dataSp.writeVTKfiles()

Observation geometry
--------------------

Once the input files are generated (this can actually be done before), we have
to define the observation geometry. These informations concerns the camera (such
as its position or the target point observed), the image (the (x,y) definition
and the number of sample per pixel) and the source (size, radius, distance,
longitude, latitude, etc.). We start by creating an instance of
:py:class:`.Geometry` (two actually, one for the plan-parallel
calculation and one for the spherical calulation)::

    geomPP = htrdr.Geometry(case='PP')
    geomSp = htrdr.Geometry(case='Sphere')

.. note::
   A geometry can be used multiple times, with different
   :py:class:`.Data` and :py:class:`.Script` object. Here, we have
   to create one insatnce for each because the target won't be at the same
   position in both case: in the spherical scenario, the planet radius adds a
   shift to the scene.

In the current scenario, we aim to reproduce some observation with constraints
being the incidence, emergence, phase and azimut angles. To calculate the camera
and source position, the :py:class:`.Geometry` objects has a method
:py:meth:`.Geometry.makeGeomFromAPIE` (a routine that makes you happy!) that
will automatically generate the geometry from the constraints on the observation
angles. Those information are read from a file (contained in the data repository
downloaded earlier). Additional information are requirted by the routine, such
as the distance between the camera and target point and the source properties
(distance, size and temperature). ::

    with open("Example_spectra_data/observation.json", 'r') as f: 
        obs = json.loads(f.read())
    cameraDist = 1000 * cst.kilo    # m
    solDist = 10 * cst.au   # m
    solRad = 7e8    # m
    solTemp = 5800  # K

    geomPP.setImage([64,64], 10)
    geomPP.makeGeomFromAPIE(obs, 0, cameraDist, solDist, solRad, srcTemp=solTemp)
    geomPP.exportGeometry()

    geomSp.setImage([64,64], 10)
    geomSp.makeGeomFromAPIE(obs, radius, cameraDist, solDist, solRad, srcTemp=solTemp)
    geomSp.exportGeometry()

.. note::
   The source temperature will be used by htrdr to determine the spectral
   distribution of the incoming stellar radiation, but note that is it also
   pausible to directly use a spectrum (c.f.
   :py:meth:`.Geometry.setSource` documentation).

.. note::
   The **image** information must be provided separately through
   :py:meth:`.Geometry.setImage` or when creating the
   :py:class:`.Geometry` instance.

.. note::
   :py:meth:`.Geometry.exportGeometry` will produce a file within the
   **geometries** repository, containing the source, camera and image
   parameters.

Scripts
-------

Now, we create an instance of :py:class:`.Script` that handles the last details
of the calulation and the call to **htrdr**. The different methods of the
:py:class:`.Script` object simplifies the call to **htrdr** via a bunch of
predefined scirpts (c.f. :ref:`script` for a full review of the
possibilities). Here, we start a reflectance spectrum (I/F) and we therefore use
the :py:meth:`.Script.reflectanceSpectrum` method, which requires the instance
of :py:class:`.Geometry` previously created, the type of calculation ("sw" or
"lw") and the list of wavlengths::

    scriptPP = htrdr.Script(case='Spectrum',
            MPIcmd="mpirun -np 4 --map-by socket:PE=9")
    scriptPP.reflectanceSpectrum(geomPP, "sw", wavelengths)

    scriptSp = htrdr.Script(case='Spectrum',
            MPIcmd="mpirun -np 4 --map-by socket:PE=9")
    scriptSp.reflectanceSpectrum(geomSp, "sw", wavelengths)

Then, we can start the calculation by calling the instance of
:py:class:`.Script` with the previously created instance of :py:class:`.Data`::

    scriptPP(dataPP)
    scriptSp(dataSp)

Post-process
------------

Finally, a call to :py:class:`.Postprocess` will treat the raw output of
**htrdr** to generate the required results, here it will be spectra, stored in
the **results_Sphere/** and the **results_PP/** repositories::

    htrdr.Postprocess(scriptPP)
    htrdr.Postprocess(scriptSp)

We can now recover and plot the calculated spectra along with the observed
specrtum as well as additonal data, all provided in the **spectrum.txt** file
previously downloaded::

    with open("results_PP/reflectance_spectrum_Spectrum.json", 'r') as f:
        resPP = json.loads(f.read())

    with open("results_Sphere/reflectance_spectrum_Spectrum.json", 'r') as f:
        resSp= json.loads(f.read())

    with open("Example_spectra_data/spectrum.txt", 'r') as f:
        f.readline()
        wvl = []
        IF_obs   = []
        IF_SHDOM = []
        IF_MCC   = []
        for line in f:
            l = line.split()
            wvl     .append(l[0])
            IF_obs  .append(l[2])
            IF_SHDOM.append(l[3])
            IF_MCC  .append(l[4])

    wvl      = np.array(wvl     , dtype=float)
    IF_obs   = np.array(IF_obs  , dtype=float)
    IF_SHDOM = np.array(IF_SHDOM, dtype=float)
    IF_MCC   = np.array(IF_MCC  , dtype=float)

    print(wvl)
    print(IF_obs)

    fig, ax = plt.subplots()

    ax.plot(wvl, IF_obs  , label="Observed", ls="", marker="o")
    ax.plot(wvl, IF_SHDOM, label="SHDOMPP", ls="", marker="o")
    ax.plot(wvl, IF_MCC  , label="MCCSPHE", ls="", marker="o")

    ax.errorbar(np.array(resPP['wavelength'])/cst.micro,
                np.array(resPP['reflectance spectrum']),
                yerr = 3 * np.array(resPP['reflectance std deviation']),
                label = "plan-parallel")

    ax.errorbar(np.array(resSp['wavelength'])/cst.micro,
                np.array(resSp['reflectance spectrum']),
                yerr = 3 * np.array(resSp['reflectance std deviation']),
                label = "spherical")

    ax.set_ylabel(r'$\frac{I}{F}$', rotation='horizontal', fontsize=15)
    ax.set_xlabel('Wavlength [µm]', fontsize=15)

    ax.legend()

    fig.savefig('Spectrum.png')

.. image:: /plot_examples/spectrum.png