Update lattice strain fitting example

Author: morganjwilliams

docs/source/gallery/examples/geochem/mineral_lattice.py

@@ -53,7 +53,7 @@
 #
 fontsize = 8
 fig, ax = plt.subplots(1)
-
+ax.set(ylabel="$D_X$", xlabel="Radii ($\AA$)", yscale="log")
 site2labels = ["Na", "Ca", "Eu", "Sr"]
 # get the Shannon ionic radii for the elements in the 2+ site
 site2radii = [
@@ -103,7 +103,6 @@
     ax.annotate(
         l, xy=(r, d), xycoords="data", ha="right", va="bottom", fontsize=fontsize
     )
-ax.set(ylabel="$D_X$", xlabel="Radii ($\AA$)", yscale="log")
 fig
 ########################################################################################
 # As europium is commonly present as a mixture of both :math:`Eu^{2+}`
@@ -136,14 +135,17 @@
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # Given the lattice strain model and a partitioning profile (for e.g. REE data),
 # we can also fit a model to a given curve; here we fit to our REE data above,
-# for which we have some known parameters to compare to:
+# for which we have some known parameters to compare to. Note that this uses
+# the youngs modulus approximation behind the scenes where it isn't provided:
 #
-from pyrolite.mineral.lattice import fit_lattice_strain, _lattice_opt_function
+from pyrolite.mineral.lattice import fit_lattice_strain, youngs_modulus_approximation
 
-_ri, _tk, _D = fit_lattice_strain(np.array(site3radii), site3Ds, z=3)
+t0 = 273.15 + 700  # estimated temperature
+_ri, _tk, _D = fit_lattice_strain(site3radii, site3Ds, z=3, t0=t0)
 ########################################################################################
 # We can compare the results of this fit to our source parameters - the ionic radius of
-# La, 900°C and estimated :math:`D_{La}`:
+# La, 900°C and estimated :math:`D_{La}` - note that the temperature isn't being fit
+# well here, but the other parameters are recovered reasonably well:
 #
 import pandas as pd
 
@@ -154,27 +156,26 @@
 ).T
 
 ########################################################################################
-# We can also compare the curves visually:
+# From this point, we could We can also compare the curves visually:
 #
 from pyrolite.plot.spider import REE_v_radii
 
 ax = REE_v_radii(index="radii")
 ax.set(ylabel="$D_X$", xlabel="Radii ($\AA$)")
+
 ax.plot(site3radii, site3Ds, label="True", color="k")
-ax.plot(
-    site3radii,
-    _lattice_opt_function(site3radii, site3radii.mean(), 298 + 1000, 1, z=3),
-    label="Starting Guess",
-    ls="--",
-    color="0.5",
+
+r0 = site3radii.mean()
+starting_guess = strain_coefficient(
+    r0, site3radii, r0=r0, z=3, T=t0, E=youngs_modulus_approximation(3, r0)
 )
+ax.plot(site3radii, starting_guess, label="Starting Guess", ls="--", color="0.5")
 
-ax.plot(
-    site3radii,
-    _lattice_opt_function(site3radii, _ri, _tk, _D, z=3),
-    label="Fit",
-    color="r",
+fit_curve = _D * strain_coefficient(
+    _ri, site3radii, r0=_ri, T=_tk, z=3, E=youngs_modulus_approximation(3, _ri)
 )
+ax.plot(site3radii, fit_curve, color="r", ls="--", label="Fit")
+
 ax.legend()
 ax.figure
 ########################################################################################
@@ -202,4 +203,3 @@
 #             American Journal of Science 314, 1319–1372.
 #             doi: `10.2475/09.2014.04 <https://doi.org/10.2475/09.2014.04>`__
 #
-

pyrolite/mineral/lattice.py

@@ -227,6 +227,11 @@ def fit_lattice_strain(
     ri, tk, D : :class:`float`
         Radius, temperature and partition coefficeint describing the
         lattice strain fit.
+
+    Notes
+    -----
+    * Uses the youngs modulus approximation where `E` not provided.
+    * Passes keyword arguments to :func:`scipy.optimize_curve_fit`.
     """
 
     popt, pcov = curve_fit(