Radio chronologie 14C 40K/40Ar U/Pb Rb/Sr

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 ¹⁴C

On mesure une concentration en ¹⁴C 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 ¹⁴C avec la position de la mesure.

age = c.c14(3)

age

12467.041271439173

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

c.c14()

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

Datation par K/Ar

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

age = c.kar(0.008, 0.0002)

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

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

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 ]))

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>

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>

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)

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)

# 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)

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)

a

4999999.999987198

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

jupyter_netbook_chrono.ipynb