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

.. only:: html

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

        :ref:`Go to the end <sphx_glr_download_api_gallery_climate_plot_05_milankovic_orbital_insolation.py>`
        to download the full example code or to run this example in your browser via JupyterLite.

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

.. _sphx_glr_api_gallery_climate_plot_05_milankovic_orbital_insolation.py:


Milanković's orbital insolation (1920)
======================================

Milutin Milanković computed the sunlight each latitude receives through the
year and how it changes as the tilt and orbit of the Earth vary. He proposed
that ice ages begin when high northern summers receive little sunlight. The
daily mean at the top of the atmosphere depends on latitude, the Sun's
declination and the Earth-Sun distance. Whatever the season, the global mean
is a quarter of the solar constant.

.. GENERATED FROM PYTHON SOURCE LINES 12-40



.. image-sg:: /api/gallery/climate/images/sphx_glr_plot_05_milankovic_orbital_insolation_001.png
   :alt: Sunlight through the year, 65°N at the summer solstice
   :srcset: /api/gallery/climate/images/sphx_glr_plot_05_milankovic_orbital_insolation_001.png
   :class: sphx-glr-single-img





.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np

    from earthsciencekit.climate import daily_insolation

    obliquity = np.deg2rad(23.44)
    day = np.arange(365)
    declination = np.arcsin(np.sin(obliquity) * np.sin(2 * np.pi * (day - 80) / 365.25))
    latitude = np.deg2rad(np.linspace(-90, 90, 181))
    flux = daily_insolation(1361, latitude[:, None], declination[None, :], 1)
    global_mean = np.trapezoid(flux * np.cos(latitude)[:, None], latitude, axis=0) / 2
    np.testing.assert_allclose(global_mean, 1361 / 4, rtol=1e-3)
    assert np.isclose(daily_insolation(1361, np.pi / 2, obliquity, 1), 1361 * np.sin(obliquity))

    tilt = np.deg2rad(np.linspace(22.1, 24.5, 50))
    summer_65n = daily_insolation(1361, np.deg2rad(65), tilt, 1)
    assert np.all(np.diff(summer_65n) > 0)

    fig, (left, right) = plt.subplots(1, 2, figsize=(10, 4))
    image = left.contourf(day, np.rad2deg(latitude), flux, levels=20, cmap="magma")
    fig.colorbar(image, ax=left, label="daily-mean insolation [W/m²]")
    left.set(xlabel="day of year", ylabel="latitude [degree]", title="Sunlight through the year")
    right.plot(np.rad2deg(tilt), summer_65n)
    right.set(xlabel="obliquity [degree]", ylabel="insolation [W/m²]", title="65°N at the summer solstice")
    fig.tight_layout()

    plt.show()


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

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


.. _sphx_glr_download_api_gallery_climate_plot_05_milankovic_orbital_insolation.py:

.. only:: html

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

    .. container:: lite-badge

      .. image:: images/jupyterlite_badge_logo.svg
        :target: ../../../lite/lab/index.html?path=api/gallery/climate/plot_05_milankovic_orbital_insolation.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_