Skip to content

Instantly share code, notes, and snippets.

@YannBouyeron

YannBouyeron/chrono.md

Last active February 28, 2023 22:21
Show Gist options
  • Save YannBouyeron/e4d0315fbf3b2b67859aba84c8d9c352 to your computer and use it in GitHub Desktop.
Save YannBouyeron/e4d0315fbf3b2b67859aba84c8d9c352 to your computer and use it in GitHub Desktop.
Download ZIP
Radio chronologie 14C 40K/40Ar U/Pb Rb/Sr
Raw
chrono.md

Datation absolue avec chrono.py

Principe de la datation absolue ou radio-chronologie

Télécharger chrono.py

Télécharger au format zip depuis Github

Télécharger avec git:

git clone https://gist.github.com/e4d0315fbf3b2b67859aba84c8d9c352.git

Télécharger avec wget:

wget https://gist.github.com/YannBouyeron/e4d0315fbf3b2b67859aba84c8d9c352/archive/15f7064b79d8f640a0c5e223133bdfa4d406a139.zip
unzip 15f7064b79d8f640a0c5e223133bdfa4d406a139.zip

Le fichier chrono.py doit ensuite être placé dans votre repertoire de travail pour pouvoir être importé avec python.

Importer matplotlib.pyplot

Si vous êtes dans jupyter:

import matplotlib.pyplot as plt
%matplotlib inline

Si vous travaillez en console linux sans interface graphique:

#import matplotlib
#import os

#if "DISPLAY" not in os.environ:
    
    #matplotlib.use("Agg")

#import matplotlib.pyplot as plt

Modification de l´apparence des graphiques pour permettre l´export depuis Jupyter en html / markdown dark theme, tout en conservant les axes visibles.

plt.style.use('seaborn-white') # blanc

#plt.style.use('dark_background') # noir

Importer chrono.py

Attention le fichier chrono.py doit se trouver dans le même repertoire que celui à partir duquel vous avez lancé python.

from chrono import Chrono

c = Chrono()

Datation au 14C

On mesure une concentration en 14C de 3 dpm dans un échantillon.

On utlilise alors la méthode c14(dpm) pour retourner l'âge de l'échantillon et la courbe de désintégration du 14C avec la position de la mesure.

age = c.c14(3)

png

age
12467.041271439173

On peut aussi simplement afficher la courbe de désintégration du 14C:

c.c14()

png

L'argument facultatif xmax permet de restreindre l'axe du temps en fixant un t max pour la représentation graphique:

c.c14(3, xmax=20000)
12467.041271439173

png

Datation par K/Ar

On utilise la méthode kar(k, Ar) avec k la concentration en 40K et Ar la concentration en 40Ar mesurées sur l'objet à dater.

age = c.kar(0.008, 0.0002)

png

Attention, sur le graphique le temps est bien en années , ici en x 109

age
385910162.54356444

L'argument facultatif ymax (rapport k/Ar max) permet d'ajuster la représentation graphique:

c.kar(0.008, 0.0002, ymax=0.05)
385910162.54356444

png

La datation par U/Pb

Fiche synthèse datation U/Pb

Premier exemple:

On commence par tracer la concordia sur l'intervalle de temps présumé avec la méthode concordia:

help(c.concordia)
Help on method concordia in module chrono:

concordia(ti=0, tf=1000000000, incre=10000000, dotlabel=True, labelincre=3) method of chrono.Chrono instance
    Arguments:
    
            ti: age présumé inférieur à réouverture
    
            tf: age présumé superieur à age roche
    
            incre: incrément de temps pour la construction de la concordia en années
    
            dotlabel: si True (defaut): indication age sur concordia
    
            labelincre: icrémentation des labels des points de la concordia pour ne pas surcharger trop
    
    Return: array des points de la concordia.

c.concordia(ti=0, tf=800000000.0, incre=100000000)
(array([0.        , 0.01562075, 0.0314855 , 0.04759808, 0.06396234,
        0.08058223, 0.09746174, 0.11460491]),
 array([0.        , 0.10349785, 0.21770751, 0.34373762, 0.48281158,
        0.63627939, 0.80563079, 0.9925097 ]))

png

On ajoute ensuite nos mesures effectuées sur les zircons de l'échantillon de la roche:

pb206U = [0.06231448, 0.06247915, 0.06264385, 0.06280857, 0.06297332, 0.06313809,
 0.06330289, 0.06346772, 0.06363257, 0.06379744, 0.06396234]
pb207U = [0.46827978, 0.46972653, 0.4711747,  0.4726243,  0.47407533, 0.47552779,
 0.47698168, 0.478437, 0.47989376, 0.48135195, 0.48281158]
c.concordia(ti=0, tf=800000000.0, incre=100000000, labelincre=2)
plt.plot(pb207U, pb206U, "*r", label="Zircon Z1 à Z11")
plt.legend()
<matplotlib.legend.Legend at 0x6d262250>

png

On observe que les mesures des rapports pb/u effectuées sur les zircons (ici en rouge) s'alignent sur la concordia. Il n'y a pas eu de réouverture du système. L'âge de la roche est l'âge du plus vieux zircon lisible sur le graphique, soit 400 Ma. Il n'est donc pas nécessaire d'afficher la discordia.

Deuxième exemple:

On dispose des mesures suivantes effectuées sur 3 zircons d'une roche:

pb207U = [5, 3, 1]

pb206U = [0.3, 0.2, 0.1]

On affiche la concordia et on place nos mesures sur le graphique:

c.concordia(ti=0.4*10**9, tf=2.2*10**9, incre=100000000)
plt.plot(pb207U, pb206U, "*", label="Minéraux")
plt.legend()
<matplotlib.legend.Legend at 0x6d271530>

png

On observe que les mesures effectuées sur les minéraux ne sont pas alignées sur la concordia. Il y'a eu une ré-ouverture du système. Il faut donc afficher la discordia avec la méthode upb :

help(c.upb)
Help on method upb in module chrono:

upb(p207, p206, ti=100000000, tf=2500000000.0, delta=100000, incre=100000000, dotlabel=True, labelincre=3) method of chrono.Chrono instance
    Arguments:
    
            p206: liste valeurs 206Pb/238U mesurées (Y)
    
            p207: liste valeurs 207Pb/235U mesurées (X)
    
            ti: age présumé inférieur à réouverture
    
            tf: age présumé superieur à age roche
    
            delta: precision du calcul des intercepts en années
    
            incre: incrément de temps pour la construction de la concordia en années
    
            dotlabel: si True (defaut): indication age sur concordia
    
            labelincre: incrémentation des labels des points de la concordia pour ne pas surcharger trop
    
    Show: graphique concordia discordia             
    Return: {coeff correlation discordia, intersup, interlow, a, b}

dico = c.upb(pb207U, pb206U, ti=0.4*10**9, tf=2.2*10**9, incre=100000000, delta=1000)

png

dico
{'R': 1.0,
 'sup': 2033986000.0,
 'low': 514235000.0,
 'a': 0.05000000000000001,
 'b': 0.049999999999999954}

La méthode retourne un dictionnaire avec l'âge de la ré-ouverture du système correspondant à l'intercept inférieur: low = 515093000 +/- delta années; l'âge de la roche correspondant à l'intercept supérieur: sup = 2033278000 +/- delta années. La précision des âges obtenus dépend du paramètre "delta" passé en argument de la méthode upb.

Attention un delta faible augmente la précision, mais augmente aussi le temps des calculs.

La datation par Rb/Sr

On utilise la méthode rbsr()

help(c.rbsr)
Help on method rbsr in module chrono:

rbsr(rbsr, srsr) method of chrono.Chrono instance
    Arguments:
    
            rbsr: liste rapport 87Rb/86Sr
    
            srsr: liste rapport 87Sr/86Sr
    
    Show graphique isochrone
    
    Return equation droite reg, age (en annees)

rbsr = [0.288, 0.31, 0.996, 0.787, 0.898, 0.945, 0.901, 5.84, 3.36]

srsr = [0.7165, 0.7171, 0.7556, 0.7413, 0.7471, 0.7521, 0.7495, 1.0133, 0.8807]
y , a = c.rbsr(rbsr, srsr)

png

# fonction de la droite isochrone

print(y)
y = 0.05354693x + 0.7

# age 

print(a)
3673415678.371164

La méthode creatRbsr(age) permet de créer des couples de valeurs rbsr et srsr fictifs en fonction de l'âge désiré.

help(c.creatRbsr)
Help on method creatRbsr in module chrono:

creatRbsr(age, n=5, b=0.7, out='') method of chrono.Chrono instance
    Création de données RbSr en fonction de l’age.
    
    Arguments:
    
            age: age désiré
            
            n: nombre de couples SrSr, RbSr
            
            b: ordonnée à l’origine = Sr/Sr à t0
    
            out: path en .xlsx ou en .json pour sauvegarder les données crées
    
    Show graphique isochrone
    
    Return RbSr (liste), SrSr (liste)

rbsr, srsr = c.creatRbsr(5000000)

png

rbsr
[2.282978862971912,
 4.988633011174326,
 3.3254901343712717,
 2.091505108366786,
 4.1554292953511585]

srsr
[0.7001620972536554,
 0.7003542055179406,
 0.7002361181816367,
 0.7001485021344575,
 0.7002950459539775]

On peut vérifier l'âge en ré utilisant la méthode rbsr():

y , a = c.rbsr(rbsr, srsr)

png

a
4999999.999987198
Raw
chrono.py
#import os
#import matplotlib
#if 'DISPLAY' not in os.environ:
#matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
class Chrono:
def __init__(self):
# constantes radioactives lambda
self.l = {
"u238": 1.55 * 10**-10,
"u235": 9.8485 * 10**-10,
"rb": 1.42 * 10**-11,
"c14": 1.21 * 10**-4,
"kar": 5.81 * 10**-11,
"kca": 4.962 * 10**-10
}
# Périodes (demie-vie)
self.T = {k: np.log(2) / self.l[k] for k in list(self.l.keys())}
def interc_upper(self, p206, p207, tf=10**9, delta=100000):
l = np.polyfit(p207, p206, 1)
a = l[0]
b = l[1]
cc = np.corrcoef(p207, p206)
r = cc[0, 1]
while 0 == 0:
cp206 = np.exp(self.l["u238"] * tf) - 1 # y = p206 / u238
cp207 = np.exp(self.l["u235"] * tf) - 1 # x = p207 / u235
if cp206 > a *cp207 + b:
return tf
else:
tf -= delta
def interc_lower(self, p206, p207, ti=100*10**6, delta=100000):
l = np.polyfit(p207, p206, 1)
a = l[0]
b = l[1]
cc = np.corrcoef(p207, p206)
r = cc[0, 1]
while 0 == 0:
cp206 = np.exp(self.l["u238"] * ti) - 1 # y = p206 / u238
cp207 = np.exp(self.l["u235"] * ti) - 1 # x = p207 / u235
if cp206 > a *cp207 + b:
return ti
else:
ti += delta
def intercept(self, p206, p207, ti=100*10**6, tf=10**9, delta=100000):
return self.interc_upper(p206, p207, tf=tf, delta=delta), self.interc_lower(p206, p207, ti=ti, delta=delta)
def concordia(self, ti=0, tf=10**9, incre=10**7, dotlabel=True, labelincre=3):
"""
Arguments:
ti: age présumé inférieur à réouverture
tf: age présumé superieur à age roche
incre: incrément de temps pour la construction de la concordia en années
dotlabel: si True (defaut): indication age sur concordia
labelincre: icrémentation des labels des points de la concordia pour ne pas surcharger trop
Return: array des points de la concordia.
"""
dt = np.arange(ti, tf, incre)
# creation de la liste des dotlabel
dt_str = [" " + str(round(i*10**-9, 3)) + " Ga" if tf >= 10**9 else " " + str(round(i*10**-6, 3)) + " Ma" for i in dt]
dt_str = [j if i%labelincre == 0 else "" for i, j in enumerate(dt_str)]
# calcul concordia sur dt
p206 = np.exp(self.l["u238"] * dt) - 1 # y = p206 / u238
p207 = np.exp(self.l["u235"] * dt) - 1 # x = p207 / u235
# representation graphique concordia
plt.close()
plt.plot(p207, p206, "--*k", markersize=6, label="Concordia")
plt.xlabel("207p / 235u", fontsize=8)
plt.ylabel("206p / 238u", fontsize=8)
if dotlabel == True:
k = 0
for xy in zip(p207, p206):
plt.annotate(dt_str[k], xy=xy, textcoords='data', fontsize=6, verticalalignment='top')
k += 1
return p206, p207
def discordia(self, p207, p206, ti=100*10**6, tf=2.5*10**9, delta=100000):
l = np.polyfit(p207, p206, 1)
a = l[0]
b = l[1]
cc = np.corrcoef(p207, p206)
r = cc[0, 1]
#calcule position intercept
inter = self.intercept(p206, p207, ti=ti, tf=tf, delta=delta)
p207.append(np.exp(self.l["u235"] * inter[1]) - 1)
p207.append(np.exp(self.l["u235"] * inter[0]) - 1)
p206.append(np.exp(self.l["u238"] * inter[1]) - 1)
p206.append(np.exp(self.l["u238"] * inter[0]) - 1)
plt.plot(p207, p206, '*b', label='', markersize=6)
plt.plot(p207, [a * i + b for i in p207], '--b', markersize=3, label='Discordia')
plt.legend(fontsize=6)
return {"R": r, "sup": inter[0], "low": inter[1], "a": a, "b": b}
def upb(self, p207, p206, ti=100*10**6, tf=2.5*10**9, delta=100000, incre=10**8, dotlabel=True, labelincre=3):
"""
Arguments:
p206: liste valeurs 206Pb/238U mesurées (Y)
p207: liste valeurs 207Pb/235U mesurées (X)
ti: age présumé inférieur à réouverture
tf: age présumé superieur à age roche
delta: precision du calcul des intercepts en années
incre: incrément de temps pour la construction de la concordia en années
dotlabel: si True (defaut): indication age sur concordia
labelincre: incrémentation des labels des points de la concordia pour ne pas surcharger trop
Show: graphique concordia discordia
Return: {coeff correlation discordia, intersup, interlow, a, b}
"""
self.concordia(ti=ti, tf=tf, incre=incre, dotlabel=True, labelincre=labelincre)
r = self.discordia(p207, p206, ti=ti, tf=tf, delta=delta)
return r
def rbsr(self, rbsr, srsr):
"""
Arguments:
rbsr: liste rapport 87Rb/86Sr
srsr: liste rapport 87Sr/86Sr
Show graphique isochrone
Return equation droite reg, age (en annees)
"""
l = np.polyfit(rbsr, srsr, 1)
a = l[0]
b = l[1]
cc = np.corrcoef(rbsr, srsr)
r = cc[0, 1]
f = 'y = ' + str(round(a, 8)) + 'x + ' + str(round(
b, 2)) # equation droit reg lin y = ax + b
t = (np.log(a + 1)) / self.l["rb"] # age = ln(a+1) / lambda
plt.close()
plt.plot(rbsr, srsr, '^k', label='', markersize=6)
plt.plot(rbsr, [a * i + b for i in rbsr], '--b', label='Isochrone RbSr')
plt.xlabel("87Rb / 86Sr", fontsize=8)
plt.ylabel("87Sr / 86Sr", fontsize=8)
plt.legend(fontsize=6)
return f, t
def creatRbsr(self, age, n=5, b=0.7, out=""):
"""
Création de données RbSr en fonction de l’age.
Arguments:
age: age désiré
n: nombre de couples SrSr, RbSr
b: ordonnée à l’origine = Sr/Sr à t0
out: path en .xlsx ou en .json pour sauvegarder les données crées
Show graphique isochrone
Return RbSr (liste), SrSr (liste)
"""
a = np.exp(self.l["rb"] * age) - 1
Rbsr = np.random.uniform(0.00200, 5.00000, n)
Srsr = [a * x + b for x in Rbsr]
if out[len(out) - 5:] == '.xlsx':
df = pd.DataFrame({"RbSr": Rbsr, "SrSr": Srsr})
df.to_excel(out)
elif out[len(out) - 5:] == '.json':
df = pd.DataFrame({"RbSr": Rbsr, "SrSr": Srsr})
df.to_json(out)
plt.close()
self.rbsr(Rbsr, Srsr)
return Rbsr.tolist(), Srsr
def c14(self, pt=0, p0=13.56, xmax=50000):
"""Méthode 14C
Arguments:
p0: concentration initial P0 (defaut: 13.56 dpm)
pt: concentration P à l’instant t (mesure)
xmax: age max présumé pour la courbe
Show graphique décroissance
Return age
"""
xmax = np.arange(0, xmax, 1)
p = p0 * np.exp(-self.l["c14"] * xmax)
plt.close()
plt.plot(xmax, p, "--k", markersize=6, label="14C dpm")
plt.xlabel("Temps (années)", fontsize=8)
plt.ylabel("14C dpm", fontsize=8)
if pt > 0:
age = np.log(p0 / pt) / self.l["c14"]
plt.plot([age], [pt], "*r", markersize=8)
return age
def kar(self, k, ar, ymax=0.2):
""" Méthode K Ar:
Arguments:
k: quantité de 40K mesurée dans l’échantillon
ar: quantité de 40Ar mesurée dans l’échantillon
ymax: limite sup de la courbe de référence
Show graphique Ar/k = f(t)
Return age (années)
"""
# clalcul age echantillon
lam = self.l["kar"] + self.l["kca"]
age = (1/lam) * np.log(1 + (ar/k) * (1 + self.l["kca"]/self.l["kar"]))
# construction courbe reférence
ark = np.arange(0.0, ymax, 0.01)
t = [(1/lam) * np.log(1 + i * (1 + self.l["kca"]/self.l["kar"])) for i in ark]
plt.close()
plt.plot(t, ark, "--k", markersize=6, label="Ar/K")
plt.plot(age, ar/k, "*b", markersize=6, label="")
plt.xlabel("Temps (années)", fontsize=8)
plt.ylabel("40Ar/40K", fontsize=8)
return age
Display the source blob
Display the rendered blob
Raw
jupyter_netbook_chrono.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment