.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/plot_03_custom_fitting.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_plot_03_custom_fitting.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_plot_03_custom_fitting.py:


Custom fitting example
======================

This is a short example for fitting data, which is not covered by the fitting widget.
It relies on calling lmfit manually. Therefore one needs to provide a fitting function, which calculates the numerical residual between measurement and model.
This notebook is about fitting multiple datasets from different angles of incidents simultaneously, but can be adapted to a wide range of different use cases.

.. GENERATED FROM PYTHON SOURCE LINES 12-20

.. code-block:: Python

    import numpy as np
    import elli
    from elli.fitting import ParamsHist
    from lmfit import minimize, fit_report
    import matplotlib.pyplot as plt










.. GENERATED FROM PYTHON SOURCE LINES 22-25

Data import
-----------
Read Ψ (psi) and Δ (delta) ellipsometric spectra for different angles from a NeXus file, limited to a wavelength range of 210–800 nm, to be compatible with a tabulated Silicon dispersion.

.. GENERATED FROM PYTHON SOURCE LINES 25-31

.. code-block:: Python

    data = elli.read_nexus_psi_delta("SiO2onSi.ellips.nxs").loc[
        (slice(None), slice(210, 800)), :
    ]
    lbda = data.loc[50].index.get_level_values("Wavelength").to_numpy()
    data






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th></th>
          <th>Ψ</th>
          <th>Δ</th>
        </tr>
        <tr>
          <th>Angle of Incidence</th>
          <th>Wavelength</th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th rowspan="5" valign="top">50</th>
          <th>210.0</th>
          <td>40.028156</td>
          <td>146.370560</td>
        </tr>
        <tr>
          <th>211.0</th>
          <td>40.012478</td>
          <td>146.607407</td>
        </tr>
        <tr>
          <th>212.0</th>
          <td>40.007420</td>
          <td>146.799438</td>
        </tr>
        <tr>
          <th>213.0</th>
          <td>40.005928</td>
          <td>147.055832</td>
        </tr>
        <tr>
          <th>214.0</th>
          <td>40.001190</td>
          <td>147.294281</td>
        </tr>
        <tr>
          <th>...</th>
          <th>...</th>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th rowspan="5" valign="top">70</th>
          <th>796.0</th>
          <td>9.051100</td>
          <td>174.423874</td>
        </tr>
        <tr>
          <th>797.0</th>
          <td>9.054414</td>
          <td>174.494553</td>
        </tr>
        <tr>
          <th>798.0</th>
          <td>9.044289</td>
          <td>174.405136</td>
        </tr>
        <tr>
          <th>799.0</th>
          <td>9.048944</td>
          <td>174.376541</td>
        </tr>
        <tr>
          <th>800.0</th>
          <td>9.021362</td>
          <td>174.418228</td>
        </tr>
      </tbody>
    </table>
    <p>1773 rows × 2 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 32-35

Setting up invariant materials and fitting parameters
-----------------------------------------------------
Set up the optical constants for the substrate from RII and define initial fit parameters for the Cauchy SiO₂ layer.

.. GENERATED FROM PYTHON SOURCE LINES 35-48

.. code-block:: Python

    rii_db = elli.db.RII()
    Si = rii_db.get_mat("Si", "Aspnes")

    params = ParamsHist()
    params.add("SiO2_n0", value=1.452, min=-100, max=100, vary=True)
    params.add("SiO2_n1", value=36.0, min=-40000, max=40000, vary=True)
    params.add("SiO2_n2", value=0, min=-40000, max=40000, vary=False)
    params.add("SiO2_k0", value=0, min=-100, max=100, vary=False)
    params.add("SiO2_k1", value=0, min=-40000, max=40000, vary=False)
    params.add("SiO2_k2", value=0, min=-40000, max=40000, vary=False)
    params.add("SiO2_d", value=20, min=0, max=40000, vary=True)









.. GENERATED FROM PYTHON SOURCE LINES 49-52

Model helper function
---------------------
This model function is not strictly needed, but simplifies the fit function, as the model only needs to be defined once.

.. GENERATED FROM PYTHON SOURCE LINES 52-71

.. code-block:: Python

    def model(lbda, angle, params):
        SiO2 = elli.Cauchy(
            params["SiO2_n0"],
            params["SiO2_n1"],
            params["SiO2_n2"],
            params["SiO2_k0"],
            params["SiO2_k1"],
            params["SiO2_k2"],
        ).get_mat()

        structure = elli.Structure(
            elli.AIR,
            [elli.Layer(SiO2, params["SiO2_d"])],
            Si,
        )

        return structure.evaluate(lbda, angle, solver=elli.Solver2x2)









.. GENERATED FROM PYTHON SOURCE LINES 72-76

Defining the fit function
-------------------------
The fit function follows the protocol defined by the lmfit package and needs the parameters dictionary as first argument. It has to return a residual value, which will be minimized.
Here psi and delta are used across all angles to calculate the residual, but could be changed to any other measured quantity like transmission or reflection data.

.. GENERATED FROM PYTHON SOURCE LINES 76-93

.. code-block:: Python



    def fit_function(params, lbda, data):
        residual = []

        for phi_i in [50, 60, 70]:
            model_result = model(lbda, phi_i, params)

            resid_psi = data.loc[(phi_i, "Ψ")].to_numpy() - model_result.psi
            resid_delta = data.loc[(phi_i, "Δ")].to_numpy() - model_result.delta

            residual.append(resid_psi)
            residual.append(resid_delta)

        return np.concatenate(residual)









.. GENERATED FROM PYTHON SOURCE LINES 94-99

Running the fit
---------------
The fitting is performed by calling the minimize function with the fit_function and the needed arguments.
It is possible to change the underlying algorithm by providing the method kwarg.
The fit_report function provides fitted values, uncertainties, and goodness-of-fit statistics.

.. GENERATED FROM PYTHON SOURCE LINES 99-103

.. code-block:: Python


    out = minimize(fit_function, params, args=(lbda, data), method="leastsq")
    print(fit_report(out))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    [[Fit Statistics]]
        # fitting method   = leastsq
        # function evals   = 119
        # data points      = 3546
        # variables        = 3
        chi-square         = 157.748841
        reduced chi-square = 0.04452409
        Akaike info crit   = -11031.1779
        Bayesian info crit = -11012.6572
    [[Variables]]
        SiO2_n0:  1.63549687 +/- 0.03012997 (1.84%) (init = 1.452)
        SiO2_n1:  27.6408771 +/- 8.75248298 (31.66%) (init = 36)
        SiO2_n2:  0 (fixed)
        SiO2_k0:  0 (fixed)
        SiO2_k1:  0 (fixed)
        SiO2_k2:  0 (fixed)
        SiO2_d:   1.75831329 +/- 0.02404272 (1.37%) (init = 20)
    [[Correlations]] (unreported correlations are < 0.100)
        C(SiO2_n0, SiO2_d)  = -0.9916
        C(SiO2_n0, SiO2_n1) = -0.9596
        C(SiO2_n1, SiO2_d)  = +0.9259




.. GENERATED FROM PYTHON SOURCE LINES 104-107

Plotting the results
----------------------------------------------
The measurent results are displayed as scatter plot and the PyElli results are overlayed as solid lines.

.. GENERATED FROM PYTHON SOURCE LINES 107-137

.. code-block:: Python


    fit_50 = model(lbda, 50, out.params)
    fit_60 = model(lbda, 60, out.params)
    fit_70 = model(lbda, 70, out.params)

    fig = plt.figure(dpi=100)
    ax = fig.add_subplot(1, 1, 1)
    ax.scatter(lbda, data.loc[(50, "Ψ")], s=20, alpha=0.1, label="50° Measurement")
    ax.scatter(lbda, data.loc[(60, "Ψ")], s=20, alpha=0.1, label="Psi 60° Measurement")
    ax.scatter(lbda, data.loc[(70, "Ψ")], s=20, alpha=0.1, label="Psi 70° Measurement")
    (psi50,) = ax.plot(lbda, fit_50.psi, c="tab:blue", label="Psi 50°")
    (psi60,) = ax.plot(lbda, fit_60.psi, c="tab:orange", label="Psi 60°")
    (psi70,) = ax.plot(lbda, fit_70.psi, c="tab:green", label="Psi 70°")
    ax.set_xlabel("wavelenth / nm")
    ax.set_ylabel("psi / degree")
    ax.legend(handles=[psi50, psi60, psi70], loc="lower left")
    fig.canvas.draw()

    fig = plt.figure(dpi=100)
    ax = fig.add_subplot(1, 1, 1)
    ax.scatter(lbda, data.loc[(50, "Δ")], s=20, alpha=0.1, label="Delta 50° Measurement")
    ax.scatter(lbda, data.loc[(60, "Δ")], s=20, alpha=0.1, label="Delta 60° Measurement")
    ax.scatter(lbda, data.loc[(70, "Δ")], s=20, alpha=0.1, label="Delta 70° Measurement")
    (delta50,) = ax.plot(lbda, fit_50.delta, c="tab:blue", label="Delta 50°")
    (delta60,) = ax.plot(lbda, fit_60.delta, c="tab:orange", label="Delta 60°")
    (delta70,) = ax.plot(lbda, fit_70.delta, c="tab:green", label="Delta 70°")
    ax.set_xlabel("wavelenth / nm")
    ax.set_ylabel("delta / degree")
    ax.legend(handles=[delta50, delta60, delta70], loc="lower right")
    fig.canvas.draw()



.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /auto_examples/images/sphx_glr_plot_03_custom_fitting_001.png
         :alt: plot 03 custom fitting
         :srcset: /auto_examples/images/sphx_glr_plot_03_custom_fitting_001.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_examples/images/sphx_glr_plot_03_custom_fitting_002.png
         :alt: plot 03 custom fitting
         :srcset: /auto_examples/images/sphx_glr_plot_03_custom_fitting_002.png
         :class: sphx-glr-multi-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.646 seconds)


.. _sphx_glr_download_auto_examples_plot_03_custom_fitting.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_03_custom_fitting.ipynb <plot_03_custom_fitting.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_03_custom_fitting.py <plot_03_custom_fitting.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_03_custom_fitting.zip <plot_03_custom_fitting.zip>`