External variable - Bloch example
[1]:
import nexus as nx
import numpy as np
import matplotlib.pyplot as plt
[2]:
vel_1000, int_1000 = nx.data.Load("spectrum_1000K.txt", x_index = 0, intensity_index = 1)
vel_800, int_800 = nx.data.Load("spectrum_800K.txt", x_index = 0, intensity_index = 1)
vel_600, int_600 = nx.data.Load("spectrum_600K.txt", x_index = 0, intensity_index = 1)
vel_400, int_400 = nx.data.Load("spectrum_400K.txt", x_index = 0, intensity_index = 1)
vel_200, int_200 = nx.data.Load("spectrum_200K.txt", x_index = 0, intensity_index = 1)
[3]:
plt.plot(vel_1000, int_1000, label = "1000 K")
plt.plot(vel_800, int_800 + 4e3, label = "800 K")
plt.plot(vel_600, int_600 + 8e3, label = "600 K")
plt.plot(vel_400, int_400 + 12e3, label = "400 K")
plt.plot(vel_200, int_200 + 16e3, label = "200 K")
plt.legend()
plt.show()
[4]:
# create external fit parameters
Tc = nx.Var(1040, min = 1000, max= 1100, fit = True, id = "Tc")
M_0 = nx.Var(40, min = 30, max= 50, fit = True, id ="M_0")
exponent = nx.Var(1.2, min = 1, max= 2, fit = True,id = "exponent")
# create a Bloch model function
# M(T) = M(0)* (1-(T/Tc)^3/2)
# we want to be able to fit all parameters
# Only the temperature T is known from the experiment
def Bloch(T):
# return a magnetic field value in Tesla
return M_0.value * (1 - (T/Tc.value)**(exponent.value))
#check the function call
print(Bloch(100))
37.59220236152488
[5]:
# create the experiment in a for loop
# we assume that the oterh hyperfine parameters do not change individually
# we also want to fit the isomer shift as a common variable to all measurements
isomer_shift = nx.Var(0.5, min = -1, max = 1, fit = True, id = "isomer all")
# It is also possible to create additional Vars
# either in the function to be fit individually
# we can use the same beam object in each experient
beam = nx.Beam()
beam.Unpolarized()
# because the hyperfine parameters change we have to create new samples for each temperature
# also we have to pass the parameters of the experiment T and velocity and intensity data
def setup_specturm(temperature, velocity_data, intensity_data):
site = nx.Hyperfine(weight=1,
isomer = isomer_shift,
# now assign the external function to the magnetic field of each experiment
# given by the Bloch function and the set temperature
magnetic_field = nx.Var(0, id="mag_at_"+str(temperature), equality = lambda: Bloch(temperature)),
isotropic = True)
mat_Fe = nx.Material.Template(nx.lib.material.Fe)
mat_Fe.hyperfine_sites = [site]
layer_Fe = nx.Layer(id = "Fe",
material = mat_Fe,
thickness = 3000)
sample = nx.Sample(layers = [layer_Fe])
exp = nx.Experiment(beam = beam,
objects = [sample],
isotope = nx.lib.moessbauer.Fe57)
spectrum = nx.MoessbauerSpectrum(experiment = exp,
velocity = velocity_data,
intensity_data = intensity_data)
return spectrum
[6]:
# create list to loop over
temps = [200, 400, 600, 800, 1000]
vels = [vel_200, vel_400, vel_600, vel_800, vel_1000]
dats = [int_200, int_400, int_600, int_800, int_1000]
measurements = []
# create all spectra
for t, v, d in zip(temps, vels, dats):
measurements.append(setup_specturm(t,v,d))
[7]:
fit = nx.Fit(measurements = measurements,
id = "Bloch fit",
external_fit_variables = [Tc, M_0, exponent])
[8]:
fit()
Run Fit instance with id: Bloch fit
Starting fit with 5 measurement data set(s) and 14 fit parameter(s):
no. | id | initial value | min | max
0 | ES scaling | 87159 | 1e-16 | 8.7159e+06
1 | ES backgr | 8602.84 | 0 | 860284
2 | isomer all | 0.5 | -1 | 1
3 | ES scaling | 87177 | 1e-16 | 8.7177e+06
4 | ES backgr | 8604.19 | 0 | 860419
5 | ES scaling | 87185 | 1e-16 | 8.7185e+06
6 | ES backgr | 8604.69 | 0 | 860469
7 | ES scaling | 87182 | 1e-16 | 8.7182e+06
8 | ES backgr | 8600.63 | 0 | 860063
9 | ES scaling | 87190 | 1e-16 | 8.719e+06
10 | ES backgr | 8606.13 | 0 | 860613
11 | Tc | 1040 | 1000 | 1100
12 | M_0 | 40 | 30 | 50
13 | exponent | 1.2 | 1 | 2
Using 5 equality constraint(s) on parameter(s):
no. | id | value
0 | mag_at_200 | 34.4683
1 | mag_at_400 | 27.2916
2 | mag_at_600 | 19.3271
3 | mag_at_800 | 10.8037
4 | mag_at_1000 | 1.83898
Using 0 inequality constraint(s).
Calling ceres solver with fit method LevMar
Ceres Solver Report: Iterations: 23, Initial cost: 4.429313e+04, Final cost: 3.090389e+02, Termination: CONVERGENCE
Gradient error analysis.
Fit performed with algorithm:
LevMar
Error analysis:
Gradient
Using 14 fit parameter(s):
no. | id | fit value | +/- std dev | initial value | min | max
0 | ES scaling | 99219.2 | 798.158 | 87159 | 1e-16 | 8.7159e+06
1 | ES backgr | 1156.36 | 682.626 | 8602.84 | 0 | 860284
2 | isomer all | 0.00049078 | 0.000492938 | 0.5 | -1 | 1
3 | ES scaling | 99312.6 | 797.726 | 87177 | 1e-16 | 8.7177e+06
4 | ES backgr | 1089.7 | 682.258 | 8604.19 | 0 | 860419
5 | ES scaling | 99458.7 | 791.2 | 87185 | 1e-16 | 8.7185e+06
6 | ES backgr | 970.047 | 676.672 | 8604.69 | 0 | 860469
7 | ES scaling | 99464.4 | 768.868 | 87182 | 1e-16 | 8.7182e+06
8 | ES backgr | 924.201 | 657.568 | 8600.63 | 0 | 860063
9 | ES scaling | 100155 | 524.807 | 87190 | 1e-16 | 8.719e+06
10 | ES backgr | 373.908 | 448.809 | 8606.13 | 0 | 860613
11 | Tc | 1042.85 | 0.165413 | 1040 | 1000 | 1100
12 | M_0 | 39.1217 | 0.0147747 | 40 | 30 | 50
13 | exponent | 1.49983 | 0.00149505 | 1.2 | 1 | 2
Values at boundaries:
Correlation matrix:
no. | 0 1 2 3 4 5 6 7 8 9 10 11 12 13
-----|------------------------------------------------------------------------------------------------------------------------------
0 | 1.000 -1.000 0.001 0.000 -0.000 -0.000 0.000 -0.000 0.000 0.000 -0.000 0.001 0.005 -0.003
1 | -1.000 1.000 -0.001 -0.000 0.000 0.000 -0.000 0.000 -0.000 -0.000 0.000 -0.001 -0.005 0.003
2 | 0.001 -0.001 1.000 -0.001 0.001 0.000 -0.000 -0.000 0.000 -0.000 0.000 -0.001 -0.000 0.000
3 | 0.000 -0.000 -0.001 1.000 -1.000 0.000 -0.000 0.000 -0.000 -0.000 0.000 -0.002 0.000 0.002
4 | -0.000 0.000 0.001 -1.000 1.000 -0.000 0.000 -0.000 0.000 0.000 -0.000 0.002 -0.000 -0.002
5 | -0.000 0.000 0.000 0.000 -0.000 1.000 -1.000 0.000 -0.000 -0.000 0.000 -0.003 -0.005 0.008
6 | 0.000 -0.000 -0.000 -0.000 0.000 -1.000 1.000 -0.000 0.000 0.000 -0.000 0.003 0.005 -0.008
7 | -0.000 0.000 -0.000 0.000 -0.000 0.000 -0.000 1.000 -1.000 0.001 -0.001 0.004 -0.005 0.005
8 | 0.000 -0.000 0.000 -0.000 0.000 -0.000 0.000 -1.000 1.000 -0.001 0.001 -0.004 0.005 -0.005
9 | 0.000 -0.000 -0.000 -0.000 0.000 -0.000 0.000 0.001 -0.001 1.000 -1.000 0.196 0.063 -0.104
10 | -0.000 0.000 0.000 0.000 -0.000 0.000 -0.000 -0.001 0.001 -1.000 1.000 -0.196 -0.063 0.104
11 | 0.001 -0.001 -0.001 -0.002 0.002 -0.003 0.003 0.004 -0.004 0.196 -0.196 1.000 0.414 -0.635
12 | 0.005 -0.005 -0.000 0.000 -0.000 -0.005 0.005 -0.005 0.005 0.063 -0.063 0.414 1.000 -0.902
13 | -0.003 0.003 0.000 0.002 -0.002 0.008 -0.008 0.005 -0.005 -0.104 0.104 -0.635 -0.902 1.000
Using 5 equality constraint(s) on parameter(s):
no. | id | value
0 | mag_at_200 | 35.835
1 | mag_at_400 | 29.827
2 | mag_at_600 | 22.047
3 | mag_at_800 | 12.835
4 | mag_at_1000 | 2.386
and 0 inequality constraint(s).
Total cost = 1.207e+00
cost for each FitMeasurement is:
no. | id | cost | %
0 | | 2.254e-01 | 18.670
1 | | 2.386e-01 | 19.761
2 | | 2.399e-01 | 19.875
3 | | 2.466e-01 | 20.428
4 | | 2.567e-01 | 21.266
Fit instance with id "Bloch fit" finished.
[9]:
offset = 0
for t, v, d, m in zip(temps, vels, dats, measurements):
int = m.Calculate()
offset += max(int) * 0.05
plt.plot(v, d + offset, label = str(t)+" K", marker='o', markersize = 2, linewidth=0)
plt.plot(v, int + offset, color='grey')
plt.legend()
plt.show()
