diff --git a/.gitignore b/.gitignore
index aee26506690023b02e63b3f2b626e96aeff42820..e23ad529d75c8d41c1b90aab4a5ad50ff027b132 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,7 +16,7 @@ MATLAB/*.asv
MATLAB/**/*.asv
# Python generated files
-py/__pycache__
+py/**/__pycache__
# Automatic generated file on Windows
desktop.ini
diff --git a/py/datafit.py b/py/datafit.py
deleted file mode 100644
index 8e48066cffd6305d1fd4206f45dd02b2edb994fb..0000000000000000000000000000000000000000
--- a/py/datafit.py
+++ /dev/null
@@ -1,154 +0,0 @@
-"""
-In this file, the four measurements from Julian are analysed and temperature
-dependent models fitted onto these results.
-"""
-
-import numpy as np
-import matplotlib.pyplot as plt
-import scipy.optimize
-from scripts import linreg, chisq_ndf
-from measurementclasses import Experiment
-
-# As their layer was really thin, assume 2D theories
-dimensions = 2;
-
-# Literature values for melting temperature, latent heat and interfacial energy
-class Literature:
- Tm = 817
- L = 878.9e6;
- IE = 5.5e-2 *3
-
-# Some physial constants
-kB = 1.380649e-23
-eV = 1.6021766208e-19
-zeroC = 273.15
-
-# Function for bulk free energy
-dGv = lambda T,L,Tm: L*(Tm-T)/Tm*2*T/(Tm+T)
-
-
-# The temperatures and geometries the data was taken at
-Ts = np.array([115,120,125,130]) + zeroC
-height = 30e-9
-
-# Read all the data
-A = Experiment("../data/experiment/Sb2Te_asd_20W_115C.csv", height)
-B = Experiment("../data/experiment/Sb2Te_asd_20W_120C.csv", height)
-C = Experiment("../data/experiment/Sb2Te_asd_20W_125C.csv", height)
-D = Experiment("../data/experiment/Sb2Te_asd_20W_130C.csv", height)
-
-
-# Get all JMAK fits and growth velocities and calculate nucleation rate from
-# them; get all nucleation rates and calculate nucleation rates from them too
-Js, vs, uvs = np.ones(Ts.shape), np.ones(Ts.shape), np.ones(Ts.shape)
-for i,M in enumerate([A,B,C,D]):
- v, uv = M.all_rfits(show=False)
- vs[i] = v
- uvs[i] = uv
- [(dNdt,_)] = M.nucleationratefit(show=False)
- if dimensions == 2:
- J = dNdt / M.A
- elif dimensions == 3:
- J = dNdt / (M.A * height)
- else:
- print("Dimensions need to be 2 or 3.")
- Js[i] = J
-
-
-def velocityfit(show=True):
- """ Plot velocities over temperature and fit arrhenius behaviour.
- Argument show can be 2 to show the arrhenius plot. """
- fitfunc = lambda x,a,b: a*np.exp(-b/x)
- m0, _, c0, _, _, _ = linreg(1/Ts, np.log(vs), uvs)
- popt, pcov = scipy.optimize.curve_fit(fitfunc, Ts, vs, p0=[np.exp(c0),-m0],
- sigma=uvs, absolute_sigma=True)
- if show==2:
- plt.errorbar(1/Ts, np.log(vs), uvs/vs, linestyle="none", marker="_",
- label="data")
- m, c = -popt[1], np.log(popt[0])
- plt.plot(1/Ts, m/Ts+c, label="fit")
- plt.xlabel("$1/(T \\ [{}^{\circ}C])$")
- plt.ylabel("$\\log(v \\ [m/s])$")
- elif show:
- plt.errorbar(Ts, vs, uvs, linestyle="none", marker="_", label="data")
- pxs = np.linspace(np.min(Ts)-1,np.max(Ts)+1,100)
- plt.plot(pxs,fitfunc(pxs,*popt), label="fit")
- plt.xlabel("$T \\ [{}^{\circ}C]$")
- plt.ylabel("$v \\ [m/s]$")
- if show:
- print("Fitting v = a*exp(-b/T)")
- print("a:",popt[0],"+-",np.sqrt(pcov[0,0]))
- print("b:",popt[1],"+-",np.sqrt(pcov[1,1]))
- Ea, uEa = [x/eV*kB for x in [popt[1], np.sqrt(pcov[1,1])]]
- print(" = (",Ea,"+-",uEa,")eV/kB")
- print("chisq/ndf:", chisq_ndf(Ts,vs,uvs,fitfunc,popt))
- plt.legend()
- plt.tight_layout()
- plt.show()
- return popt, pcov
-
-def nucleationratefit(plotrange=None, show=True):
- """
- Fit classical nucleation theory formula to J(T).
- Use either JMAK theory fit results(second argument False) or counted
- nuclei (second argument True).
- A range of temperatures to plot the calculated rate over can be
- given manually as a 2 tuple.
- """
-
- # Formulas for critical energy and nucleation rate
- dGc = lambda T,IE,L,Tm: 16*np.pi*IE**3/(3*dGv(T,L,Tm)**2)
- J = lambda T,b,IE,L,Tm : (b*dGv(T,L,Tm)**2/(8*np.pi) *
- 1/np.sqrt(IE**3*kB*T) *
- np.exp(-dGc(T,IE,L,Tm)/kB/T))
- # Fit with a free prefactor and interfacial energy, other should be known
- fitfunc = lambda T,b,IE : J(T,b,IE,Literature.L,Literature.Tm)
-# fitfunc = lambda T,b,IE,L,Tm : J(T,b,IE,L,Tm)
-
- # Guess for IE is literature value, choose initial prefactor to fit some
- # arbitary point. Do the fit.
- b0 = Js[3]/fitfunc(Ts[3],1,Literature.IE)
- popt, pcov = scipy.optimize.curve_fit(fitfunc, Ts, Js, maxfev=60000,
- p0=[b0,Literature.IE])
-# p0=[b0,Literature.IE,Literature.L,Literature.Tm]
- if show:
- print("Fitting J = b*dGv**2/(8*pi)/sqrt(IE**3*kB*T) * exp(-dGc/kB/T)")
- print(" dGc = 16*pi*IE^3/(3*dGv)")
- print(" dGv = L*(Tm-T)/Tm*2*T/(Tm+T)")
- print("b: ",popt[0],"+-",np.sqrt(pcov[0,0]))
- print("IE:",popt[1],"+-",np.sqrt(pcov[1,1]))
-# print("L: ",popt[2]," +- ",np.sqrt(pcov[2,2]))
-# print("Tm:",popt[3]," +- ",np.sqrt(pcov[3,3]))
- plt.plot(Ts, Js, linestyle="none", marker="x", label="mean rates")
- if plotrange is None:
- pxs = np.linspace(np.min(Ts)-3,np.max(Ts)+3,50)
- else:
- pxs = np.linspace(plotrange[0],plotrange[-1],200)
- plt.plot(pxs,fitfunc(pxs,*popt), label="fit")
- plt.xlabel("$T \\ [{}^{\circ}C]$")
- plt.ylabel("$J \\ [1/m^{}/s]$".format(dimensions))
- plt.tight_layout()
- plt.legend()
- plt.show()
- return popt, pcov
-
-def mobilityfit(show=True):
- """ Calculate mobility (as defined in the AIST paper), plot over
- temperature and fit arrhenius behaviour """
-
- M = v/dGv(Ts,Literature.L,Literature.Tm)
-
- plt.plot(Ts, M, linestyle="none", marker="x", label="data")
- expfunc = lambda x,a,b: a*np.exp(b/x)
- m0, _, c0, _, _, _ = linreg(1/Ts, np.log(M), uvs)
- popt, pcov = scipy.optimize.curve_fit(expfunc, Ts, M, p0=[np.exp(c0),m0],
- sigma=uvs)
- pxs = np.linspace(np.min(Ts)-3,np.max(Ts)+3,50)
- plt.plot(pxs,expfunc(pxs,np.exp(c0),m0),label="fit")
- plt.xlabel("T [C]")
- plt.ylabel("M [m/s / (J/m^3)]")
-
-#A.nucleationratefit(splitat=[1.2e3,10e3,19e3,28e3])
-#B.nucleationratefit(splitat=[1000,8888])
-#C.nucleationratefit(splitat=[2000,3500])
-#D.nucleationratefit(splitat=[250,1100,1900])
diff --git a/py/experimentalfits.py b/py/experimentalfits.py
new file mode 100644
index 0000000000000000000000000000000000000000..c3ed03beef8ea5fb07dda76224ccb56050ce48a4
--- /dev/null
+++ b/py/experimentalfits.py
@@ -0,0 +1,115 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import fits
+
+# Set font size and stuff for images in the tex document
+plt.rcParams.update({'font.size': 16,
+ 'lines.linewidth': 2,
+ 'lines.markersize': 10,
+ 'lines.markeredgewidth': 2})
+
+# Read the experiments
+height = 30e-9
+Ex115 = fits.Experiment("../data/experiment/Sb2Te_asd_20W_115C.csv", height)
+Ex120 = fits.Experiment("../data/experiment/Sb2Te_asd_20W_120C.csv", height)
+Ex125 = fits.Experiment("../data/experiment/Sb2Te_asd_20W_125C.csv", height)
+Ex130 = fits.Experiment("../data/experiment/Sb2Te_asd_20W_130C.csv", height)
+
+# Methods performing fits on this data, giving images and information ready to
+# use in tex
+
+def show_incubationtime():
+ axes = Ex120.plot_N()
+ fig = axes.figure
+ line = axes.lines[0]
+ line.set_linestyle(":")
+ axes.set_xlim(-120, 3500)
+ axes.set_ylim(-5, 120)
+ axes.set_title("Sb${}_2$Te, 120°C")
+ fig.canvas.draw()
+ fig.tight_layout()
+
+def fit_v(Ex=Ex130, i=2):
+ vfit = Ex.velocityFit()
+ vfit.fitall()
+ vfit.info()
+ axes = vfit.plot(i)
+ Ex.plot_r(i, axes)
+ [fitline, expline] = axes.lines
+ expline.set_color("red")
+ fitline.set_color("black")
+ axes.legend()
+ fig = axes.figure
+ fig.tight_layout()
+
+def fit_J_115():
+ jfit = Ex115.nucleationRateFit([1.7e3,9.7e3,20e3]).fit()
+ jfit.info()
+ axes = jfit.plotJ()
+ axes.set_xlim(-500,31800)
+ axes.set_ylim(-0.001, 0.048)
+ fig = axes.figure
+ fig.tight_layout()
+
+def fit_Nc_115():
+ jfit = Ex115.nucleationRateFit([1.7e3,9.7e3,20e3]).fit()
+ jfit.info()
+ axes = jfit.plotN()
+ axes.figure.tight_layout()
+
+def comp_120():
+ axes = Ex120.plot_fc()
+ axes.lines[0].set_markersize(6)
+ fit = Ex120.jmakFit().fit(2)
+ rec = Ex120.jmakBuild(2)
+ fit.plot(axes)
+ rec.plot(axes)
+ axes.legend()
+ axes.figure.tight_layout()
+
+def comp_125():
+ axes = Ex125.plot_fc()
+ axes.lines[0].set_markersize(6)
+ fit = Ex125.jmakFit().fit(2)
+ rec = Ex125.jmakBuild(2)
+ fit.plot(axes)
+ rec.plot(axes)
+ axes.legend()
+ axes.figure.tight_layout()
+
+# Combine the measurements and perform fits to their relations
+
+S = fits.SeriesOfMeasurements(np.array([115,120,125,130]) + 273.15,
+ [Ex115,Ex120,Ex125,Ex130], 2)
+S.load()
+V = fits.ArrheniusFit(S).fit()
+C = fits.SimpleCntFit(S, 296.4e6, 837, 9e-2).fit()
+
+def fit_cnt():
+ C.info()
+ fig, axes = plt.subplots(1,1)
+ axes.set_xlim(384.7,405.3)
+ axes.set_ylim(-1e12,4e13)
+ C.plot(axes)
+ axes.legend()
+ fig.tight_layout()
+
+def fit_arrhenius():
+ V.info()
+ axes = V.plot()
+ axes.legend()
+ axes.figure.tight_layout()
+
+def fit_arrhenius_lin():
+ V.info()
+ axes = V.plot_lin()
+ axes.set_xticks(np.linspace(2.48e-3, 2.58e-3, 3))
+ axes.legend()
+ axes.figure.tight_layout()
+
+def fit_mobility(L=296.4e6, Tm=837):
+ mf = S.mobilityFit(lambda T: L*(Tm-T)/(2*Tm)*2*T/(Tm+T) ).fit()
+ mf.info()
+ axes = mf.plot()
+ axes.legend()
+ axes.figure.tight_layout()
diff --git a/py/fits/__init__.py b/py/fits/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..269a1989138b65b1e992df994243b93b52297e25
--- /dev/null
+++ b/py/fits/__init__.py
@@ -0,0 +1,15 @@
+"""
+Analysis of experiments on phase change materials.
+"""
+
+# Import classes for single measurements and fits
+from .measurements import *
+from .measurements import __all__
+
+# Import classes for combining several measurements
+from .series import *
+from .series import __all__ as __all2__
+
+# Correct list of available classes
+__all__.extend(__all2__)
+del __all2__
diff --git a/py/fits/measurements.py b/py/fits/measurements.py
new file mode 100644
index 0000000000000000000000000000000000000000..c7d6ec352e51b17ae563dbe88c75da344b1cae0a
--- /dev/null
+++ b/py/fits/measurements.py
@@ -0,0 +1,551 @@
+"""
+This file provides classes to analyse experimental or simulated data on the
+crystallisation of a phase change material.
+The first type of classes represent measurements:
+ Measurement - Parent class for all types of measurements
+ Experiment - Class for experimental results as I got from Julian
+ Simulation - Class for simulation results as returned from the export
+ methods in my simulation code.
+The second type represents fits to several quantities:
+ VelocityFit - Fit to the grain sizes
+ NucleationRateFit - Fit to the number of counted nuclei
+ JmakFit - Fit of a JMAK model to the crystallised fraction
+ JmakBuild - Conterpart to JmakFit: Get a crystallised fraction
+ from JMAK model and fit results for v and J
+"""
+
+import csv
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.optimize
+from .scripts import linreg, weightedmean, chisq_ndf
+
+NaN = float('nan')
+
+def lininterp(a):
+ """
+ Linear inperpolation of an array.
+ The length of the returned array is one smaller as the one of the given.
+ Each element is the mean of two neighbouring elements of a.
+ """
+ return (a[1:] + a[:-1]) / 2
+
+class Measurement:
+ """
+ A class for a measurement for one temperature as I got from Julian or as
+ returned from one simulation run.
+ Properties:
+ times - times of measurements [s]
+ nuclei - number of seen nuclei
+ radii - 2D array of radii [m]
+ afract - amorphous fraction [1]
+ cfract - crystallized fraction [1]
+ """
+
+ def plot_fc(self, axes=None):
+ """ Plot the crystallised fraction over time. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ axes.plot(self.times, self.cfract, linestyle="none", marker="x",
+ label="data")
+ axes.set_xlabel("$t \\ [s]$")
+ axes.set_ylabel("$f_c$")
+ return axes
+
+ def plot_N(self, axes=None):
+ """ Plot the counted nuclei over time. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ axes.plot(self.times, self.nuclei, linestyle="none", marker="x",
+ label="data")
+ axes.set_xlabel("$t \\ [s]$")
+ axes.set_ylabel("$N$")
+ return axes
+
+ def plot_r(self, i, axes=None):
+ """ Plot the size of the i-th nucleus over time. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ axes.plot(self.times, self.radii[i]*1e6, linestyle="none", marker="x",
+ label="$r_{{{}}}$".format(i))
+ axes.set_xlabel("$t \\ [s]$")
+ axes.set_ylabel("$r \\ [\\mu m]$")
+ return axes
+
+ def jmakFit(self):
+ """ Return an object for a JMAK fit to measurement. """
+ return JmakFit(self)
+
+ def velocityFit(self):
+ """ Return an object for a velocity fit to measurement. """
+ return VelocityFit(self)
+
+ def nucleationRateFit(self, splitat=None):
+ """ Return an object for a nucleation rate fit to measurement. """
+ return NucleationRateFit(self,splitat)
+
+ def jmakBuild(self, dimensions):
+ """
+ Return an object giving a function for the crystallised fraction
+ calculated in JMAK theory with velocity and nucleation rate from fits.
+ """
+ return JmakBuild(self.velocityFit().fitall(),
+ self.nucleationRateFit(splitat=[]).fit(),
+ dimensions)
+
+
+class VelocityFit:
+ """
+ A class to fit a linear growht velocity to one measurement, to get
+ information about the fit and to visualise it.
+ Properties:
+ vs - Array of velocities for the individual grains [m/s]
+ uvs - Array of uncertainties to vs [m/s]
+ v - Mean velocity [m/s]
+ uv - Uncertainty to v [m/s]
+ """
+
+ def __init__(self, M):
+ """ Build an object for a measurement M. """
+ self.M = M
+ # Extract the necessary data
+ self.rs = M.radii
+ self.ts = M.times
+ # Initialise results
+ self.vs = np.zeros(len(M.radii)) * NaN
+ self.uvs = np.zeros(len(M.radii)) * NaN
+ self._cs = np.zeros(len(M.radii)) * NaN
+ self.v, self.uv = NaN, NaN
+
+ def fit(self, i, maxtime=None):
+ """
+ Fit a linear groth to the i-th grain.
+ If the growth is not linear anymore after a time maxtime, the fit only
+ takes earlier points into account.
+ """
+ if maxtime is None:
+ maxtime = float('inf')
+ # Extract points where t < tmax
+ cond = np.logical_and(np.logical_not( np.isnan(self.rs[i]) ),
+ self.ts < maxtime)
+ ts = np.compress(cond, self.ts)
+ rs = np.compress(cond, self.rs[i])
+ # Linear regression
+ v, uv, c, _, _, chiq = linreg(ts, rs, np.ones(rs.shape))
+ # Scale to make chisquared/ndf equals 1
+ chiq = chiq/(len(ts)-2)
+ uv = uv * np.sqrt(chiq)
+ # Save results
+ self.vs[i], self.uvs[i], self._cs[i] = v, uv, c
+ return self
+
+ def fitall(self):
+ """ Fit all grains and calculate mean growht velocity """
+ for i,_ in enumerate(self.vs):
+ self.fit(i)
+ self.v, self.uv = np.mean(self.vs), np.std(self.vs)
+ return self
+
+ def info(self):
+ """ Print some information about this fit. """
+ print("Velocity fit:")
+ for i, (v, uv) in enumerate(zip(self.vs, self.uvs)):
+ print(" {}-th grain: v={}+-{}".format(i,v,uv))
+ print("mean uncertainty:", np.mean(self.uvs))
+ print("overall result: v =", self.v, "+-", self.uv)
+ return self
+
+ def plot(self, i, axes=None):
+ """ Plot the fit to the i-th grain """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ ts, rs, v, c = self.ts, self.rs[i], self.vs[i], self._cs[i]
+ ts = np.compress(np.logical_not(np.isnan(rs)), ts)
+ axes.plot(ts, (v*ts+c)*1e6, label="$v_{{{}}}\\cdot t$".format(i),
+ alpha=0.7)
+ axes.set_xlabel("$t \\ [s]$")
+ axes.set_ylabel("$r \\ [\\mu m]$")
+ return axes
+
+class NucleationRateFit:
+ """
+ Fit one or more nucleation rates on the data.
+ Properties:
+ Js - List of (J,uJ) tupels where J is a nucleation rate [1/s] and uJ
+ its uncertainty
+ domainboundaries - List of (tmin, tmax) tupels defining the borders of
+ the domains a constant rate was fittet to
+ """
+
+ def __init__(self, M, splitat):
+ """ Build an object for a measurement M. """
+ self.M = M
+
+ # If no points to split at are given, create an empty list
+ if splitat is None:
+ splitat = []
+ # Calculate domain boundarys: Insert first time at front, cut out every
+ # time thats out of available times and insert last time at end
+ splitat = np.insert(splitat, 0, M.times[0])
+ splitat = np.compress(splitat<M.times[-1], splitat)
+ splitat = np.append(splitat,M.times[-1])
+ self.domainboundaries = list(zip(splitat[0:-1], splitat[1:]))
+
+ # Extract necessary data
+ self._dN = np.diff(M.nuclei)
+ self._dt = np.diff(M.times)
+ self._tm = lininterp(M.times) # interpolated times
+ self._ta = M.times # all times
+ self._fa = lininterp(M.afract)
+ self._Js = self._dN/self._dt * 1/self._fa
+ self._uJs = np.sqrt(np.mean(self._Js)/(self._dt * self._fa))
+
+ self.Js = [(NaN,NaN)]
+
+ def fit(self):
+ """ Perform the fit. """
+ res = []
+ # For each of the domains
+ for tmin, tmax in self.domainboundaries:
+ # extract data in this domain
+ cond = np.logical_and(self._tm >= tmin, self._tm <= tmax)
+ Js_i = np.compress(cond, self._Js)
+ dN_i = np.compress(cond, self._dN)
+ dt_i = np.compress(cond, self._dt)
+ fa_i = np.compress(cond, self._fa)
+ # Assume poisson uncertainty with the mean rate on all points,
+ # override the original guessed uncertainty
+ uJ_i = np.sqrt(np.maximum(np.mean(dN_i),1))/dt_i * 1/fa_i
+ self._uJs[cond] = uJ_i
+ # Mean rate
+ J, uJ = weightedmean(Js_i, uJ_i)
+ res.append((J,uJ))
+ self.Js = res
+ return self
+
+ def info(self):
+ """ Print some information about this fit. """
+ for (tmin, tmax), (J,uJ) in zip(self.domainboundaries, self.Js):
+ print("from",tmin,"to",tmax,":")
+ print(" rate:",J,"+-",uJ)
+ return self
+
+ def plotN(self, axes=None):
+ """ Plot the integrated rates. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ Js, dt = self._Js, self._dt
+ # Integrate again
+ N2 = np.cumsum(Js * dt)
+ N2 = np.insert(N2,0,0)
+ # Assume poisson uncertainty
+ uN2 = np.maximum(np.sqrt(N2), 1)
+ axes.errorbar(self._ta, N2, uN2, linestyle=":", label="data")
+ y0, y1 = 0, 0
+ for (tmin, tmax), (J,_) in zip(self.domainboundaries, self.Js):
+ y0, y1 = y1, y1+J*(tmax-tmin)
+ plt.plot([tmin,tmax], [y0, y1], label="integrated mean rate")
+ axes.set_xlabel("$t \\ [s]$")
+ axes.set_ylabel("$N_{corrected}$")
+ return axes
+
+ def plotJ(self, axes=None):
+ """ Plot momentary and mean rates. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ t, Js, uJs = self._tm, self._Js, self._uJs
+ axes.errorbar(t, Js, uJs, linestyle="none", marker="_")
+ for (tmin, tmax), (J,_) in zip(self.domainboundaries, self.Js):
+ axes.hlines(J,tmin,tmax)
+ axes.set_xlabel("$t\\ [s]$")
+ axes.set_ylabel("$\\frac{1}{f_a}\\frac{\\Delta N}{\\Delta t}$")
+ return axes
+
+
+class JmakFit:
+ """
+ A class to fit a JMAK model to the crystallised fraction of a measurement.
+ Properties:
+ n - JMAK exponent [1]
+ k - JMAK coefficient [s^(-n)]
+ un - Uncertainty on n [1]
+ uk - Uncertainty on k [s^(-n)]
+ corr - Correlation coefficient of n and k
+ chisq - Chisquared of the fit
+ """
+
+ def __init__(self, M):
+ """ Build an object for a measurement M. """
+ self.M = M
+ # Initialise used values
+ self.n, self.un = NaN, NaN
+ self.k, self.uk = NaN, NaN
+ self.corr, self.chisq = NaN, NaN
+
+ def linearise(self):
+ """
+ Get the linearised vaiables:
+ xs = log( t )
+ ys = log( - log( 1 - f_c ))
+ sigmas = const / ( log(f_a) * f_a )
+ """
+ # Extract points where log() is well defined
+ M = self.M
+ cond = np.logical_and( np.logical_and( M.afract > 0,
+ M.afract < 1),
+ M.times > 0)
+ # Linearise
+ xs = np.log( np.compress(cond, M.times))
+ fs = np.compress(cond, M.afract)
+ ys = np.log( - np.log( fs ))
+ # Assume equal uncertainty (=1) on area values, transform to y
+ sigmas = np.abs(np.ones(xs.shape)/(np.log(fs)*fs))
+ return xs, ys, sigmas
+
+ def linfit(self):
+ """
+ Perform a fit to the linearised problem and save the result.
+ """
+ xs, ys, sigmas = self.linearise()
+ n, un, c, uc, corr, chisq = linreg(xs, ys, sigmas)
+ chisqndf = chisq / (len(xs)-2)
+ # correct uncertainties, calculate k
+ un, uc = un*np.sqrt(chisqndf), uc*np.sqrt(chisqndf)
+ k, uk = np.exp(c), np.exp(c)*uc
+ self.n, self.un = n, un
+ self.k, self.uk = k, uk
+ self.corr, self.chisq = corr, 1
+ return self
+
+ def fit(self, dimensions=None):
+ """
+ Perform a fit on the original data, considering also values where the
+ linear problem is not defined
+ """
+ # Initial values from linear problem
+ if self.k is NaN:
+ self.linfit()
+ n0, k0 = self.n, self.k
+ # Get data
+ times, fa = self.M.times, self.M.afract
+ # Define function to fit onto
+ fitfunc = lambda t,n,k: np.exp(-k*t**n)
+ if dimensions is None:
+ popt, pcov = scipy.optimize.curve_fit(fitfunc, times, fa,
+ p0=[n0,k0])
+ corr = pcov[0,1]/np.sqrt(pcov[0,0]*pcov[1,1])
+ else:
+ # If a dimension is given, calculate exponent and put it into
+ # the fitfunction
+ n = dimensions + 1
+ fitfunc2 = lambda t,k: fitfunc(t,n,k)
+ popt, pcov = scipy.optimize.curve_fit(fitfunc2, times, fa, p0=[k0])
+ # Fill popt and pcov to the same size as from above
+ popt = np.insert(popt, 0, n)
+ pcov = np.array([[0, 0], [0, pcov[0,0]]])
+ corr = 0
+ # Get chisquared and correct uncertainties, save
+ chisq = chisq_ndf(times, fa, np.ones(fa.shape), fitfunc, popt)
+ pcov = pcov*chisq
+ self.n, self.k = popt
+ self.un, self.uk = np.sqrt(pcov[0,0]), np.sqrt(pcov[1,1])
+ self.corr, self.chisq = corr, 1
+ return self
+
+ def info(self):
+ """ Print some information about the fit. """
+ print("f_a = exp(-k*t^n)")
+ print("<=> ln(-ln(f_a)) = n * ln(t) + ln(k)")
+ print("Result:")
+ print(" n:", self.n, "+-", self.un)
+ print(" k:", self.k, "+-", self.uk)
+ print(" correlation:", self.corr)
+ return self
+
+ def resfun(self, t):
+ """
+ Return the crystalline fraction to a time t as calculated from the
+ model and the results from the fit.
+ """
+ return 1-np.exp(-self.k*t**self.n)
+
+ def linplot(self, axes=None):
+ """ Plot the linearised problem and the fit. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ xs, ys, _ = self.linearise()
+ n, k = self.n, self.k
+ c = np.log(k)
+ axes.plot(xs, ys, linestyle="none", marker="x", label="data")
+ axes.plot(xs, n*xs+c, label="fit")
+ axes.set_xlabel('$\ln(t \\ [s])$')
+ axes.set_ylabel('$\ln( - \ln( 1 - f_c ))$')
+ return axes
+
+ def plot(self, axes=None):
+ """ Plot the fit. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ plottimes = np.linspace(np.min(self.M.times), np.max(self.M.times), 100)
+ axes.plot(plottimes, self.resfun(plottimes), label="JMAK fit")
+ axes.set_xlabel("$t \\ [s]$")
+ axes.set_ylabel("$f_c$")
+ return axes
+
+
+class JmakBuild:
+ """
+ This class gives a counterpart to JmakFit coming from a defined velocity
+ and nucleation rate.
+ Properties:
+ M - Measurement
+ n - JMAK exponent [1]
+ k - JMAK coefficient [s^(-n)]
+ un - Uncertainty on n [1]
+ uk - Uncertainty on k [s^(-n)]
+ """
+
+ def __init__(self, vFit, jFit, dimensions):
+ """
+ Initialise. Needs a velocity and a nucleation rate fit object
+ and a number of dimensions to use.
+ """
+ # Save the fits
+ self.M = jFit.M
+ self.vFit = vFit
+ self.jFit = jFit
+ # Extract velocity with uncertainty
+ v, uv = vFit.v, vFit.uv
+ if v is NaN:
+ raise AttributeError("Velocity not yet fitted.")
+ # Try to get a single value for J
+ try:
+ [(dNdt, udNdt)] = jFit.Js
+ except ValueError as err:
+ print("Probably more than one one nucleation rate was fitted.")
+ raise err
+ # Try to get an area and a height
+ try:
+ A, h = self.M.A, self.M.h
+ except AttributeError as err:
+ print("Data regarding the volume missing in the used measurement.")
+ raise err
+ # Calculate J
+ J = dNdt / (A*h)
+ uJ = udNdt / (A*h)
+ # Calculate JMAK formula according to choosen dimension
+ if dimensions == 1:
+ self.k = J*A*v
+ self.uk = np.sqrt( (uJ*A*v)**2 + (J*A*uv)**2 )
+ elif dimensions == 2:
+ self.k = np.pi/3 * J * h * v**2
+ self.uk = np.pi/3 * np.sqrt( (uJ*h*v**2)**2 + (2*J*h*v*uv)**2 )
+ elif dimensions == 3:
+ self.k = np.pi/3 * J * v**3
+ self.uk = np.pi/3 * np.sqrt( (uJ*v**3)**2 + (3*J*v**2*uv)**2 )
+ else:
+ raise ValueError("Wrong number of dimensions")
+ # Exponent with uncertainty
+ self.n, self.un = dimensions + 1, 0
+
+ def resfun(self, t):
+ """
+ Time dependtent equation for the crystallised fraction from the fit
+ results.
+ """
+ return 1-np.exp(-self.k*t**self.n)
+
+ def plot(self, axes=None):
+ """ Plot the fit. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ plottimes = np.linspace(np.min(self.M.times), np.max(self.M.times), 100)
+ axes.plot(plottimes, self.resfun(plottimes),
+ label="JMAK $\\left(v,\\ \\frac{dN}{dt} \\right)$")
+ axes.set_xlabel("$t \\ [s]$")
+ axes.set_ylabel("$f_c$")
+ return axes
+
+
+# Calculate a length from the values in pixels, according to the mail I got.
+# First method is the correction I got in the second mail
+correctedpx = lambda px: 300/72*px
+px2len = lambda p: (0.3876 * correctedpx(p) - 0.3802) * 1e-6
+
+class Experiment(Measurement):
+ """
+ Class for experimental data as I got from Julian.
+ Additional properties:
+ columntitles - column titles from the csv
+ aarea - amorphous area [m^2]
+ A - overall observed area [m^2]
+ h - thickness of the material [m]
+ """
+
+ def __init__(self, filename, height):
+ """ Read all data from the file and build variables of interest. """
+
+ # Open file and build csv reader
+ file = open(filename, newline="")
+ reader = csv.reader(file, delimiter=",")
+ # First line has titles to the columns
+ self.columntitles = list(map(lambda x: x.strip(), next(reader)))
+ # Rest is data structured in columns
+ columns = [ [] for _ in self.columntitles]
+ for row in reader:
+ for column, value in zip(columns, row):
+ # Replace empty cells with nan to parse them to float
+ if value == "":
+ value = "nan"
+ column.append(value.replace(",","."))
+ # First column are times in minutes, convert to seconds
+ self.times = np.array(columns[0], dtype=np.float) * 60
+ # Second column is amporphous area. The square root is a length and can
+ # be converted according to the mail, then square the result again
+ self.aarea = (px2len(np.sqrt(np.array(columns[1], dtype=np.float))))**2
+ # Third column has number of grains, here no value means probably zero
+ self.nuclei = np.array(columns[2], dtype=np.float)
+ self.nuclei[np.isnan(self.nuclei)] = 0
+ # Rest of the file are radii of grains
+ self.radii = px2len( np.array(columns[3:], dtype=np.float) )
+ # Overall observed area
+ self.A = max(self.aarea)
+ self.h = height
+ # Fractions of amorphous and crystaline area
+ self.afract = self.aarea/self.A
+ self.cfract = 1-self.afract
+
+
+class Simulation(Measurement):
+ """
+ Class for simulation results as received from my code.
+ Additional properties:
+ columntitles - column titles from the csv
+ """
+
+ def __init__(self, filename):
+ """ Read all data from the file and build variables of interest. """
+
+ # Open file and build csv reader
+ file = open(filename, newline="")
+ reader = csv.reader(file, delimiter=",")
+ # First line has titles to the columns
+ self.columntitles = list(map(lambda x: x.strip(), next(reader)))
+ # Rest is data structured in columns
+ lines = []
+ for row in reader:
+ lines.append(row)
+ data = np.array(lines, dtype=np.float)
+ # First column are times
+ self.times = data[:,0]
+ # Second column has number of grains
+ self.nuclei = data[:,1]
+ # Third column is crystalline fraction
+ self.cfract = data[:,2]
+ self.afract = 1 - self.cfract
+ # Rest of the file are sizes of grains
+ self.radii = np.sqrt(data[:,3:]/np.pi).transpose()
+
+
+__all__ = ["Measurement", "Experiment", "Simulation", "VelocityFit",
+ "NucleationRateFit", "JmakFit", "JmakBuild"]
diff --git a/py/scripts.py b/py/fits/scripts.py
similarity index 70%
rename from py/scripts.py
rename to py/fits/scripts.py
index 44783ba2219207009638cd52f85dde2884375386..219bbb510e6374a196c70b7a7d502499ded93c32 100644
--- a/py/scripts.py
+++ b/py/fits/scripts.py
@@ -1,4 +1,10 @@
-# Some math stuff
+"""
+This file provides some general mathematical methods:
+ linreg - Linear regression
+ weightedmean - Weighted mean value of values with uncertainties
+ chisq_ndf - Chisquared per degree of freedom for given data, function and
+ parameters
+"""
import numpy as np
@@ -41,4 +47,13 @@ def chisq_ndf(x,y,uncert_y,f,popt):
Returns the chisquared per degree of freedom for data x, y with
uncertainty uncert_y and a function f(x,*popt) with parameters popt.
"""
- return np.sum(((y-f(x,*popt))/uncert_y)**2) / (len(x) - len(popt))
\ No newline at end of file
+ return np.sum(((y-f(x,*popt))/uncert_y)**2) / (len(x) - len(popt))
+
+def enlargedLinspace(left,right,num=50,factor=1):
+ """
+ Get evenly spaced numbers over the interval from left to right, enlarged by
+ a factor.
+ Default factor is one, default amount of points is 50.
+ """
+ shift = (right-left)*factor
+ return np.linspace(left-shift, right+shift, num)
\ No newline at end of file
diff --git a/py/fits/series.py b/py/fits/series.py
new file mode 100644
index 0000000000000000000000000000000000000000..f6eb2dea942916bc6035ce831442e23e5deb53ff
--- /dev/null
+++ b/py/fits/series.py
@@ -0,0 +1,298 @@
+"""
+This file provides classes to fit temperature dependent models to data obtained
+from several measurements at different temperatures.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.optimize
+from .scripts import linreg, chisq_ndf
+
+# Some physial constants
+kB = 1.380649e-23
+eV = 1.6021766208e-19
+
+class SeriesOfMeasurements:
+ """
+ A class for a series of several measurements at different temperatures.
+ """
+
+ def __init__(self, temperatures, measurements, dimensions):
+ """
+ Initialise a series of measurements (given as a list) at different
+ temperatures (given as a list).
+ Analyse assuming a given number of dimensions.
+ """
+ self.Ts = np.array(temperatures)
+ self.Ms = measurements
+ self.dimensions = dimensions
+
+ def load(self):
+ """
+ Do all the fits for all the measurements.
+ """
+ self.vs = np.zeros(len(self.Ms))
+ self.uvs = np.zeros(len(self.Ms))
+ self.Js = np.zeros(len(self.Ms))
+ self.uJs = np.zeros(len(self.Ms))
+ for i,M in enumerate(self.Ms):
+ vFit = M.velocityFit().fitall()
+ self.vs[i], self.uvs[i] = vFit.v, vFit.uv
+ volume = M.A * M.h
+ jFit = M.nucleationRateFit().fit()
+ [(J,uJ)] = jFit.Js
+ self.Js[i], self.uJs[i] = J/volume, uJ/volume
+ return self
+
+ def arrheniusFit(self):
+ """
+ Return an arrhenius behaviour fit for the growth velocities of this
+ series of measurements.
+ """
+ return ArrheniusFit(self)
+
+ def simpleCntFit(self, L, Tm, IE):
+ """
+ Return an object to fit classical nucleation theory to the rates of
+ this series, using given values for the latent heat L, the melting
+ temperature Tm and an initial guess for the interfacial energy IE.
+ """
+ return SimpleCntFit(self, L, Tm, IE)
+
+ def mobilityFit(self, dGv):
+ """
+ Return an arrhenius behaviour fit for the mobility v/dGv of this
+ series of measurements, using a given temperature dependent function
+ for the free energy.
+ """
+ return MobilityFit(self, dGv)
+
+
+class ArrheniusFit:
+ """ Class to fit an Arrhenius behaviour to the growth velocity. """
+
+ def __init__(self, series):
+ """ Build an object of a given series of measurements. """
+ self.Ts = series.Ts
+ self.vs = series.vs
+ self.uvs = series.uvs
+
+ def linfit(self):
+ """ Fit in the linearised problem. """
+ m, um, c, uc, _, _ = linreg(1/self.Ts, np.log(self.vs), self.uvs)
+ self.vinf, self.uvinf = np.exp(c), np.exp(c)*uc
+ self.Ea, self.uEa = -m, um
+ return self
+
+ def fit(self):
+ """ Perform a fit. """
+ self.linfit()
+ fitfunc = lambda x,a,b: a*np.exp(-b/x)
+ Ts, vs, uvs = self.Ts, self.vs, self.uvs
+ p0 = [self.vinf, self.Ea]
+ popt, pcov = scipy.optimize.curve_fit(fitfunc, Ts, vs, p0=p0,
+ sigma=uvs, absolute_sigma=True)
+ self.vinf, self.Ea = popt[0], popt[1]
+ self.uvinf, self.uEa = np.sqrt(pcov[0,0]), np.sqrt(pcov[1,1])
+ return self
+
+ def v(self,T):
+ """
+ Velocity as a function of temperature according to the fit result.
+ """
+ return self.vinf * np.exp(-self.Ea/T)
+
+ def info(self):
+ """ Print some information. """
+ print("Fitting v = vinf*exp(-Ea/T)")
+ print("vinf = {} +- {}".format(self.vinf, self.uvinf))
+ print("Ea = {} +- {}".format(self.Ea, self.uEa))
+ Ea, uEa = [x/eV*kB for x in [self.Ea, self.uEa]]
+ print(" = (",Ea,"+-",uEa,")eV/kB")
+ print("chisq/ndf:", chisq_ndf(self.Ts, self.vs, self.uvs,
+ (lambda x,a,b: a*np.exp(-b/x)),
+ [self.vinf, self.Ea]))
+ return self
+
+ def plot_lin(self, axes=None):
+ """ Plot the linearised problem. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ invTs = 1/self.Ts
+ m, c = -self.Ea, np.log(self.vinf)
+ axes.errorbar(invTs, np.log(self.vs), self.uvs/self.vs,
+ linestyle="none", marker="_", label="data")
+ axes.plot(invTs, m*invTs+c, label="fit")
+ axes.set_xlabel("$1/(T \\ [K])$")
+ axes.set_ylabel("$\\log(v \\ [m/s])$")
+ return axes
+
+ def plot(self, axes=None):
+ """ Show the data and fit. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ axes.errorbar(self.Ts, self.vs, self.uvs, linestyle="none", marker="_",
+ label="data")
+ pxs = np.linspace(np.min(self.Ts)-1,np.max(self.Ts)+1,100)
+ axes.plot(pxs,self.v(pxs), label="fit")
+ axes.set_xlabel("$T \\ [K]$")
+ axes.set_ylabel("$v \\ [m/s]$")
+ return axes
+
+
+# Function for bulk free energy and critical energy
+def bulkFreeEnergy(T,L,Tm):
+ return L*(Tm-T)/Tm*2*T/(Tm+T)
+
+def criticalFreeEnergy(T,IE,L,Tm):
+ return 16*np.pi*IE**3/(3*bulkFreeEnergy(T,L,Tm)**2)
+
+class SimpleCntFit:
+ """
+ A simple fit of classical nucleation theory onto the nucleation rates,
+ using given values for latent heat L and melting temperature Tm and leaving
+ the prefactor b and interfacial energy IE as fit parameters.
+ """
+
+ def __init__(self, series, L, Tm, IE):
+ """
+ Build an object of a given series of measurements, using given values
+ for the latent heat L, melting temperature Tm and an initial guess for
+ the interfacial energy IE.
+ """
+ self.Ts = series.Ts
+ self.Js = series.Js
+ self.uJs = series.uJs
+ self.L, self.Tm, self.IE = L, Tm, IE
+
+ def fit(self):
+ """ Perform the fit. """
+ Ts, Js, uJs = self.Ts, self.Js, self.uJs
+ # Functions for nucleation rate from CNT
+ J = lambda T,b,IE,L,Tm: (b*bulkFreeEnergy(T,L,Tm)**2/(8*np.pi) *
+ 1/np.sqrt(IE**3*kB*T) *
+ np.exp(-criticalFreeEnergy(T,IE,L,Tm)/kB/T))
+ # Fitfunction with L and Tm from literature
+ fitfunc = lambda T,b,IE: J(T,b,IE,self.L,self.Tm)
+ # Guess for IE is literature value, choose initial prefactor to fit some
+ # arbitary point. Do the fit.
+ b0 = Js[len(Js)//2]/fitfunc(Ts[len(Js)//2],1,self.IE)
+ popt, pcov = scipy.optimize.curve_fit(fitfunc, Ts, Js, maxfev=60000,
+ p0=[b0,self.IE])
+ self.b, self.IE = popt
+ self.ub, self.uIE = np.sqrt(pcov[0,0]), np.sqrt(pcov[1,1])
+# self.b, self.IE = b0, self.IE
+# self.ub, self.uIE = 1, 1
+ return self
+
+ def J(self,T):
+ """
+ Nucleation rate as a temperature dependent function from this fit.
+ """
+ b, IE, L, Tm = self.b, self.IE, self.L, self.Tm
+ return (b*bulkFreeEnergy(T,L,Tm)**2/(8*np.pi) *
+ 1/np.sqrt(IE**3*kB*T) *
+ np.exp(-criticalFreeEnergy(T,IE,L,Tm)/kB/T))
+
+ def info(self):
+ """ Print some information. """
+ print("Fitting J = b*dGv**2/(8*pi)/sqrt(IE**3*kB*T) * exp(-dGc/kB/T)")
+ print(" dGc = 16*pi*IE^3/(3*dGv)")
+ print(" dGv = L*(Tm-T)/Tm*2*T/(Tm+T)")
+ print("b: ",self.b,"+-",self.ub)
+ print("IE:",self.IE,"+-",self.uIE)
+ return self
+
+ def plot(self, axes=None):
+ """ Plot the problem and fit. """
+ Ts, Js = self.Ts, self.Js
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ plotrange = np.min(Ts), np.max(Ts)
+ else:
+ plotrange = axes.get_xlim()
+ axes.plot(Ts, Js, linestyle="none", marker=".", label="mean rates")
+ pxs = np.linspace(plotrange[0], plotrange[1], 200)
+ axes.plot(pxs, self.J(pxs), label="fit")
+ axes.set_xlabel("$T \\ [K]$")
+ axes.set_ylabel("$J \\ [1/m^3/s]$")
+ return axes
+
+class MobilityFit:
+ """ Fitting an Arrhenius behaviour to the mobility M = v / dGv """
+
+ def __init__(self, series, dGv):
+ """
+ Build an object of a given series of measurements and a given model for
+ the free energy dGv as a function fo temperature.
+ """
+ self.Ts = series.Ts
+ self.vs = series.vs
+ self.uvs = series.uvs
+ self.Ms = self.vs / dGv(self.Ts)
+ self.uMs = self.uvs / dGv(self.Ts)
+
+ def linfit(self):
+ """ Fit the linear problem. """
+ m, um, c, uc, _, _ = linreg(1/self.Ts, np.log(self.Ms), self.uMs)
+ self.Minf, self.uMinf = np.exp(c), np.exp(c)*uc
+ self.Ea, self.uEa = -m, um
+ return self
+
+ def fit(self):
+ """ Perform a fit. """
+ self.linfit()
+ fitfunc = lambda x,a,b: a*np.exp(-b/x)
+ Ts, Ms, uMs = self.Ts, self.Ms, self.uMs
+ p0 = [self.Minf, self.Ea]
+ popt, pcov = scipy.optimize.curve_fit(fitfunc, Ts, Ms, p0=p0,
+ sigma=uMs, absolute_sigma=True)
+ self.Minf, self.Ea = popt[0], popt[1]
+ self.uMinf, self.uEa = np.sqrt(pcov[0,0]), np.sqrt(pcov[1,1])
+ return self
+
+ def M(self,T):
+ """
+ Mobility as temperature dependent function according to this fit results.
+ """
+ return self.Minf * np.exp(-self.Ea/T)
+
+ def info(self):
+ """ Print some information. """
+ print("Fitting M = Minf*exp(-Ea/T)")
+ print("Minf = {} +- {}".format(self.Minf, self.uMinf))
+ print("Ea = {} +- {}".format(self.Ea, self.uEa))
+ Ea, uEa = [x/eV*kB for x in [self.Ea, self.uEa]]
+ print(" = (",Ea,"+-",uEa,")eV/kB")
+ print("chisq/ndf:", chisq_ndf(self.Ts, self.Ms, self.uMs,
+ (lambda x,a,b: a*np.exp(-b/x)),
+ [self.Minf, self.Ea]))
+ return self
+
+ def plot_lin(self, axes=None):
+ """ Plot the linearised problem. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ invTs = 1/self.Ts
+ m, c = -self.Ea, np.log(self.Minf)
+ axes.errorbar(invTs, np.log(self.Ms), self.uMs/self.Ms,
+ linestyle="none", marker="_", label="data")
+ axes.plot(invTs, m*invTs+c, label="fit")
+ axes.set_xlabel("$1/(T \\ [K])$")
+ axes.set_ylabel("$\\log(M \\ \\left[\\frac{m^3}{J s}\\right])$")
+ return axes
+
+ def plot(self, axes=None):
+ """ Show a plot. """
+ if axes is None:
+ fig, axes = plt.subplots(1,1)
+ axes.errorbar(self.Ts, self.Ms, self.uMs, linestyle="none", marker="_",
+ label="data")
+ pxs = np.linspace(np.min(self.Ts)-1,np.max(self.Ts)+1,100)
+ axes.plot(pxs,self.M(pxs), label="fit")
+ axes.set_xlabel("$T \\ [K]$")
+ axes.set_ylabel("$M \\ \\left[\\frac{m^3}{J s}\\right]$")
+ return axes
+
+__all__ = ["SeriesOfMeasurements", "ArrheniusFit", "SimpleCntFit",
+ "MobilityFit"]
diff --git a/py/measurementclasses.py b/py/measurementclasses.py
deleted file mode 100644
index 89513e90a4ba2915f2814fe19e48b5ff055b0260..0000000000000000000000000000000000000000
--- a/py/measurementclasses.py
+++ /dev/null
@@ -1,351 +0,0 @@
-"""
-This file provides classes to analyse experimental or simulated data on the
-crystallisation of a phase change material.
-"""
-
-import csv
-import numpy as np
-import matplotlib.pyplot as plt
-import scipy.optimize
-from scripts import linreg, weightedmean, chisq_ndf
-
-
-class Measurement:
- """
- A class for a measurement for one temperature as I got from Julian or as
- returned from one simulation run.
- Properties:
- times - times of measurements [s]
- nuclei - number of seen nuclei
- radii - 2D array of radii [m]
- afract - amorphous fraction
- cfract - crystallized fraction
- Functions:
- plot_fc
- plot_N
- plot_r
- JMAK_linfit
- JMAK_fit
- rfit
- all_rfits
- nucleationratefit
- """
-
- def plot_fc(self):
- """ Plot the crystallised fraction over time. """
- plt.plot(self.times, self.cfract, linestyle="none", marker="x")
- plt.xlabel("$t \\ [s]$")
- plt.ylabel("$f_c$")
- plt.tight_layout()
-
- def plot_N(self):
- """ Plot the counted nuclei over time. """
- plt.plot(self.times, self.nuclei, linestyle="none", marker="x")
- plt.xlabel("$t \\ [s]$")
- plt.ylabel("$N$")
- plt.tight_layout()
-
- def plot_r(self, i):
- """ Plot the size of the i-th nucleus over time. """
- plt.plot(self.times, self.radii[i]*1e6, linestyle="none", marker="x")
- plt.xlabel("$t \\ [s]$")
- plt.ylabel("$r \\ [\\mu m]$")
- plt.tight_layout()
-
- def JMAK_linfit(self, show=True):
- """
- Linearise the fit to a JMAK model and use linear regression.
- Uncertainties are given so that chisquared per degree of freedom is 1.
- Returns (n, k) for f_a = exp(-k*t^n)
- """
-
- # Extract points where log() is well defined
- cond = np.logical_and( np.logical_and( self.afract > 0,
- self.afract < 1),
- self.times > 0)
- xs = np.log( np.compress(cond, self.times))
- fs = np.compress(cond, self.afract)
- ys = np.log( - np.log( fs ))
- # Assume equal uncertainty on area values, transform to y
- sigmas = np.abs(np.ones(xs.shape)/np.log(fs)/fs)
- n, un, c, uc, corr, chisq = linreg(xs, ys, sigmas)
- chisqndf = chisq / (len(xs)-2)
- # correct uncertainties, calculate k
- un, uc = un*np.sqrt(chisqndf), uc*np.sqrt(chisqndf)
- k, uk = np.exp(c), np.exp(c)*uc
- if show:
- print("f_a = exp(-k*t^n), linearise:")
- print("ln(-ln(f_a)) = n * ln(t) + ln(k)")
- print("n:", n, "+-", un)
- print("k:", k, "+-", uk)
- plt.plot(xs, ys, linestyle="none", marker="x", label="data")
- plt.plot(xs, n*xs+c, label="fit")
- plt.xlabel('$\ln(t \\ [s])$')
- plt.ylabel('$\ln( - \ln( 1 - f_c ))$')
- plt.legend()
- plt.tight_layout()
- return n, k
-
-
- def JMAK_fit(self, dimensions=None, show=True):
- """
- Fit to a JMAK model.
- Either with a fixed number of dimensions or using this as another free
- parameter.
- Uncertainties are given so that chisquared per degree of freedom is 1.
- Returns ([n,k],pcov) where pcov is the covariance matix.
- """
-
- # Initial values from linear problem, function to fit onto
- n0, k0 = self.JMAK_linfit(show=False)
- fitfunc = lambda t,n,k: np.exp(-k*t**n)
- if dimensions is None:
- popt, pcov = scipy.optimize.curve_fit(fitfunc, self.times,
- self.afract, p0=[n0,k0])
- else:
- # If dimension is given, calculate exponent and put into fitfunction
- n = dimensions + 1
- fitfunc2 = lambda t,k: fitfunc(t,n,k)
- popt, pcov = scipy.optimize.curve_fit(fitfunc2, self.times,
- self.afract, p0=[k0])
- # Fill popt, pcov to the same size as from above
- popt = np.insert(popt, 0, n)
- pcov = np.array([[0, 0], [0, pcov[0,0]]])
- # Get chisquared and correct uncertainties
- chisq = chisq_ndf(self.times, self.afract, np.ones(self.afract.shape),
- fitfunc, popt)
- pcov = pcov*chisq
- if show:
- print("f_a = exp(-k*t^n)")
- print("n:", popt[0], "+-", np.sqrt(pcov[0,0]))
- print("k:", popt[1], "+-", np.sqrt(pcov[1,1]))
- print("covariance: ", pcov[0,1])
- plt.plot(self.times, self.afract, linestyle="none", marker="x",
- label="data")
- plottimes = np.linspace(np.min(self.times), np.max(self.times), 100)
- plt.plot(plottimes, fitfunc(plottimes, *popt), label="fit")
- plt.xlabel("$t \\ [s]$")
- plt.ylabel("$1-f_c$")
- plt.legend()
- plt.tight_layout()
- return popt, pcov
-
- def rfit(self, i, show=True):
- """
- Fit a linear growth to the i-th grain.
- Uncertainties are given so that chisquared per degree of freedom is 1.
- Returns (v, sigma_v).
- """
-
- cond = np.logical_not( np.isnan(self.radii[i]) )
- ts = np.compress(cond, self.times)
- rs = np.compress(cond, self.radii[i])
- v, uv, c, uc, corr, chiq = linreg(ts, rs, np.ones(rs.shape))
- # Scale to make chisquared/ndf equals 1
- chiq = chiq/(len(ts)-2)
- uv, uc = [x * np.sqrt(chiq) for x in [uv, uc]]
- if show:
- print("v:", v, "+-", uv)
- plt.plot(ts, rs*1e6, linestyle="none", marker="x", label="data")
- plt.plot(ts, (v*ts+c)*1e6, label="fit", color="black", alpha=0.7)
- plt.xlabel("$t \\ [s]$")
- plt.ylabel("$r \\ [\\mu m]$")
- if show != 2:
- plt.legend()
- plt.tight_layout()
- plt.show()
- return v, uv
-
- def all_rfits(self, show=True):
- """ Fit all grains and calculate mean growht velocity """
-
- # Initialise results array and do first fit
- vs = np.zeros(len(self.radii))
- vs[0] = self.rfit(0, show=show)[0]
- # For all other fits, do not show a legend anymore
- for i in range(1,len(self.radii)):
- vs[i] = self.rfit(i, show=show*2)[0]
- v, uv = np.mean(vs), np.std(vs)
- return v, uv
-
- def nucleationratefit(self, splitat=None, show=True):
- """
- Fit a nucleation rate to the number of grains. Either to the whole
- domain, or give points to split at. Will return a list of 2-tuples
- each containing a rate and an uncertainty.
- """
-
- # Calculate the new nuclei per time step
- dN = np.diff(self.nuclei)
- dt = np.diff(self.times)
- meantimes = (self.times[1:]+self.times[:-1])/2
- # Scale with fraction of material where nucleation can occour
- fa = (self.afract[1:] + self.afract[:-1]) / 2
- Js = dN/dt * 1/fa
- # Integrate again
- N2 = np.cumsum(Js * dt)
- N2 = np.insert(N2,0,0)
- # Assume poisson uncertainty
- uN2 = np.maximum(np.sqrt(N2), 1)
-
- # If no points to split at are given, create an empty list
- if splitat is None:
- splitat = []
- # Calculate domain boundarys: Insert first time at front, cut out every
- # time thats out of available times and insert last time at end
- splitat = np.insert(splitat, 0, self.times[0])
- splitat = np.compress(splitat<self.times[-1], splitat)
- splitat = np.append(splitat,self.times[-1])
-
- # Show the data
- if show:
- plt.errorbar(self.times, N2, uN2, linestyle=":", label="data")
- plt.xlabel("$t \\ [s]$")
- plt.ylabel("$N_{corrected}$")
- plt.tight_layout()
-
- # For each of the domains
- res = []
- y0, y1 = 0, 0
- for tmin, tmax in zip(splitat[0:-1],splitat[1:]):
- # extract data in this domain
- cond = np.logical_and(meantimes >= tmin,
- meantimes <= tmax)
- Js_i = np.compress(cond, Js)
- dN_i = np.compress(cond, dN)
- dt_i = np.compress(cond, dt)
- fa_i = np.compress(cond,fa)
-
- # Assume poisson uncertainty with the mean rate on all points
- uJ_i = np.sqrt(np.maximum(np.mean(dN_i),1))/dt_i * 1/fa_i
-
- # Mean rate
- J, uJ = weightedmean(Js_i, uJ_i)
-
- if show:
- print("from",tmin,"to",tmax,":")
- print(" rate:",J,"+-",uJ)
- y0, y1 = y1, y1+J*(tmax-tmin)
- plt.plot([tmin,tmax], [y0, y1], label="integrated mean rate")
-
- res.append((J,uJ))
-
- return res
-
-
-# Calculate a length from the values in pixels, according to the mail I got.
-# First method is the correction I got in the second mail
-correctedpx = lambda px: 300/72*px
-px2len = lambda p: (0.3876 * correctedpx(p) - 0.3802) * 1e-6
-
-class Experiment(Measurement):
- """
- Class for experimental data as I got from Julian.
- Additional properties:
- columntitles - column titles from the csv
- aarea - amorphous area [m^2]
- A - overall observed area [m^2]
- h - thickness of the material [m]
- Additional methods:
- compare_results
- """
-
- def __init__(self, filename, height):
- """ Read all data from the file and build variables of interest. """
-
- # Open file and build csv reader
- file = open(filename, newline="")
- reader = csv.reader(file, delimiter=",")
- # First line has titles to the columns
- self.columntitles = list(map(lambda x: x.strip(), next(reader)))
- # Rest is data structured in columns
- columns = [ [] for _ in self.columntitles]
- for row in reader:
- for column, value in zip(columns, row):
- # Replace empty cells with nan to parse them to float
- if value == "":
- value = "nan"
- column.append(value.replace(",","."))
- # First column are times in minutes, convert to seconds
- self.times = np.array(columns[0], dtype=np.float) * 60
- # Second column is amporphous area. The square root is a length and can
- # be converted according to the mail, then square the result again
- self.aarea = (px2len(np.sqrt(np.array(columns[1], dtype=np.float))))**2
- # Third column has number of grains, here no value means probably zero
- self.nuclei = np.array(columns[2], dtype=np.float)
- self.nuclei[np.isnan(self.nuclei)] = 0
- # Rest of the file are radii of grains
- self.radii = px2len( np.array(columns[3:], dtype=np.float) )
- # Overall observed area
- self.A = max(self.aarea)
- self.h = height
- # Fractions of amorphous and crystaline area
- self.afract = self.aarea/self.A
- self.cfract = 1-self.afract
-
- def compare_results(self, dimensions):
- """
- Compare the results from the JMAK fit with the ones obtained from
- fitting on grain count and sizes.
- """
- v, uv = self.all_rfits(show=False)
- [(dNdt, udNdt)] = self.nucleationratefit(show=False)
- J = dNdt / (self.A*self.h)
- uJ = udNdt / (self.A*self.h)
- [n,k1], pcov = self.JMAK_fit(dimensions,show=False)
-
- if dimensions == 1:
- func = lambda t: np.exp(-J*self.h**2*v*t**2)
- k2 = J*self.h**2*v
- uk2 = k2 * np.sqrt(uJ**2/J**2 + uv**2/v**2)
- elif dimensions == 2:
- func = lambda t: np.exp(-np.pi/3*J*self.h*v**2*t**3)
- k2 = np.pi/3*J*self.h*v**2
- uk2 = k2 * np.sqrt(uJ**2/J**2 + 4*uv**2/v**2)
- elif dimensions == 3:
- func = lambda t: np.exp(-np.pi/3*J*v**3*t**4)
- k2 = np.pi/3*J*v**3
- uk2 = k2 * np.sqrt(uJ**2/J**2 + 9*uv**2/v**2)
-
- print("JMAK fit: k = {} +- {}".format(k1,np.sqrt(pcov[0,0])))
- print("v & J fits: k = {} +- {}".format(k2,uk2))
-
- plottimes = np.linspace(np.min(self.times), np.max(self.times), 100)
- plt.plot(self.times, self.afract, linestyle="none", marker="x",
- label="data")
- plt.plot(plottimes, np.exp(-k1*plottimes**n),
- label="fit on $f_c$")
- plt.plot(plottimes, func(plottimes),
- label="fits on $\\frac{dN}{dt}$ and $v$")
- plt.legend()
-
-
-class Simulation(Measurement):
- """
- Class for simulation results as received from my code.
- Additional properties:
- columntitles - column titles from the csv
- """
-
- def __init__(self, filename):
- """ Read all data from the file and build variables of interest. """
-
- # Open file and build csv reader
- file = open(filename, newline="")
- reader = csv.reader(file, delimiter=",")
- # First line has titles to the columns
- self.columntitles = list(map(lambda x: x.strip(), next(reader)))
- # Rest is data structured in columns
- lines = []
- for row in reader:
- lines.append(row)
- data = np.array(lines, dtype=np.float)
- # First column are times
- self.times = data[:,0]
- # Second column has number of grains
- self.nuclei = data[:,1]
- # Third column is crystalline fraction
- self.cfract = data[:,2]
- self.afract = 1 - self.cfract
- # Rest of the file are sizes of grains
- self.radii = np.sqrt(data[:,3:]/np.pi).transpose()