.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/volcanology/plot_10_wetherill_concordia.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_volcanology_plot_10_wetherill_concordia.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_volcanology_plot_10_wetherill_concordia.py:


Wetherill's concordia (1956)
============================

George Wetherill plotted the two uranium-lead clocks against each other.
Minerals that stayed closed lie on a curve, the concordia, where both give the
same age. Zircons that lost part of their lead in a later event fall on a
straight chord, the discordia, whose two intersections with the concordia date
the crystallization and the lead loss.

.. GENERATED FROM PYTHON SOURCE LINES 11-47



.. image-sg:: /api/gallery/volcanology/images/sphx_glr_plot_10_wetherill_concordia_001.png
   :alt: Discordia intercepts: 2.7 Ga and 0.5 Ga
   :srcset: /api/gallery/volcanology/images/sphx_glr_plot_10_wetherill_concordia_001.png
   :class: sphx-glr-single-img





.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.optimize import brentq

    from earthsciencekit.units import convert
    from earthsciencekit.volcanology import concordia_ratios

    formed, reset = convert(2.7, "Ga", "s"), convert(0.5, "Ga", "s")
    upper, lower = concordia_ratios(formed), concordia_ratios(reset)
    loss = np.array([0.1, 0.3, 0.5, 0.7])
    zircons = (1 - loss)[:, None] * upper + loss[:, None] * lower
    slope, intercept = np.polyfit(zircons[:, 0], zircons[:, 1], 1)


    def gap(age):
        x, y = concordia_ratios(age)
        return y - (intercept + slope * x)


    ages = [brentq(gap, convert(a, "Ga", "s"), convert(b, "Ga", "s")) for a, b in [(0.1, 1.5), (1.5, 4.0)]]
    assert np.allclose(convert(np.array(ages), "s", "Ga"), [0.5, 2.7])

    curve = concordia_ratios(convert(np.linspace(0, 3.2, 200), "Ga", "s"))
    ticks = concordia_ratios(convert(np.arange(0.5, 3.5, 0.5), "Ga", "s"))
    fig, ax = plt.subplots()
    ax.plot(*curve.T, label="concordia")
    ax.plot(*ticks.T, "k.")
    for age, point in zip(np.arange(0.5, 3.5, 0.5), ticks):
        ax.annotate(f"{age:g} Ga", point, textcoords="offset points", xytext=(5, -8), fontsize=8)
    ax.plot([lower[0], upper[0]], [lower[1], upper[1]], "--", color="0.5", label="discordia")
    ax.plot(*zircons.T, "o", color="C3", label="zircons with lead loss")
    ax.set(xlabel="207Pb* / 235U", ylabel="206Pb* / 238U", title="Discordia intercepts: 2.7 Ga and 0.5 Ga")
    ax.legend()

    plt.show()


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

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


.. _sphx_glr_download_api_gallery_volcanology_plot_10_wetherill_concordia.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/volcanology/plot_10_wetherill_concordia.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

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

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