Skip to content
Snippets Groups Projects
Commit 2dd3153a authored by Arnaud HEMMERLE's avatar Arnaud HEMMERLE
Browse files

Add first version nb

parent 4fd4d1d4
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
xrr_analysis_refnx.ipynb 0 → 100644
+ 470
− 0
View file @ 2dd3153a
%% Cell type:markdown id:cdcdc51f tags:
# Import data
%% Cell type:code id:98e73446 tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from refnx.dataset import Data1D
from refnx.util import refplot
data_init = np.loadtxt('SIRIUS_2024_09_20_4446-4530_XRR.dat')
qz_exp = data_init[:,0]/10. #A^-1
R_exp = data_init[:,1]
R_err_exp = data_init[:,2]
data = Data1D(data=(qz_exp, R_exp, R_err_exp))
#refplot(data)
```
%% Cell type:code id:b7dbf4ea tags:
``` python
plt.errorbar(qz_exp, R_exp, R_err_exp, marker = 'x', linestyle = '')
plt.yscale('log')
plt.xlabel(r'$q_z$ (A$^{—1}$)')
plt.ylabel(r'R')
plt.show()
```
%% Output
%% Cell type:markdown id:e3a440cd tags:
# Creating model
%% Cell type:markdown id:d081153a tags:
## Define SLD
%% Cell type:code id:bca2c9cf tags:
``` python
from refnx.reflect import SLD, Slab, ReflectModel
# SLD = r_e*rho_el
r_el = 2.81794e-5 # in A
# SLD in refnx should be in 1e-6 A^-2
SLD_water = 0.334*r_el*1e6
SLD_helium = 0.
water = SLD(SLD_water, name = 'water')
helium = SLD(SLD_helium, name = 'helium')
```
%% Cell type:markdown id:f25844f4 tags:
## Define slabs
%% Cell type:code id:85f3f29a tags:
``` python
# first number is thickness, second number is roughness (with respect to top layer)
# lengths in A
# Semi-infinite bulk : thickness = 0
water_bulk = water(0,1)
helium_atm = helium(0,0)
# Constructed from top to bottom
structure = helium_atm | water_bulk
fig, ax1 = plt.subplots()
# First y-axis (SLD)
ax1.plot(*structure.sld_profile())
ax1.set_xlabel('Distance / $\AA$')
ax1.set_ylabel('SLD / $10^{-6} \AA^{-2}$', color='b')
ax1.tick_params(axis='y', labelcolor='b')
# Define transformation functions
def sld_to_ed(sld): # sld in 10^-6 Å⁻²
return (sld * 1e-6) / r_el # e⁻/ų
def ed_to_sld(ed): # e⁻/ų
return (ed * r_el) * 1e6 # back to 10^-6 Å⁻²
# Add secondary y-axis
ax2 = ax1.secondary_yaxis('right', functions=(sld_to_ed, ed_to_sld))
ax2.set_ylabel('Electron Density / e$^- \ \AA^{-3}$', color='r')
ax2.tick_params(axis='y', labelcolor='r')
plt.title('SLD and Electron Density Profile')
plt.tight_layout()
plt.show()
```
%% Output
%% Cell type:markdown id:5a2c4843 tags:
## Plot initial guess
%% Cell type:code id:3ad6f334 tags:
``` python
model = ReflectModel(structure, bkg=0)
plt.errorbar(qz_exp, R_exp, R_err_exp, marker = 'x', linestyle = '')
plt.plot(qz_exp, model(qz_exp), 'k-')
plt.xlabel(r'$q_z$ (A$^{—1}$)')
plt.ylabel(r'R')
plt.yscale('log')
plt.show()
```
%% Output
%% Cell type:markdown id:30d14073 tags:
# Fit
%% Cell type:markdown id:4914e6dc tags:
## Define fit parameters
%% Cell type:code id:d973e8ea tags:
``` python
water_bulk.rough.setp(bounds=(0.5, 7), vary=True)
```
%% Cell type:markdown id:f59599a3 tags:
## Fit the data
%% Cell type:code id:f0138907 tags:
``` python
from refnx.analysis import Transform, CurveFitter, Objective, Model, Parameter
objective = Objective(model, data, transform=Transform("logY"))
fitter = CurveFitter(objective)
fitter.fit("differential_evolution")
print(objective)
```
%% Output
56813.855182359344: : 4it [00:00, 5.03it/s]
________________________________________________________________________________
Objective - 140218423146320
Dataset = <None>, 85 points
datapoints = 85
chi2 = 114299.55496834047
Weighted = True
Transform = Transform('logY')
________________________________________________________________________________
Parameters: ''
________________________________________________________________________________
Parameters: 'instrument parameters'
<Parameter: 'scale' , value=1 (fixed) , bounds=[-inf, inf]>
<Parameter: 'bkg' , value=0 (fixed) , bounds=[-inf, inf]>
<Parameter:'dq - resolution', value=5 (fixed) , bounds=[-inf, inf]>
<Parameter: 'q_offset' , value=0 (fixed) , bounds=[-inf, inf]>
________________________________________________________________________________
Parameters: 'Structure - '
________________________________________________________________________________
Parameters: 'helium'
<Parameter:'helium - thick', value=0 (fixed) , bounds=[-inf, inf]>
________________________________________________________________________________
Parameters: 'helium'
<Parameter:'helium - sld' , value=0 (fixed) , bounds=[-inf, inf]>
<Parameter:'helium - isld', value=0 (fixed) , bounds=[-inf, inf]>
<Parameter:'helium - rough', value=0 (fixed) , bounds=[-inf, inf]>
<Parameter:'helium - volfrac solvent', value=0 (fixed) , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters: 'water'
<Parameter:'water - thick', value=0 (fixed) , bounds=[-inf, inf]>
________________________________________________________________________________
Parameters: 'water'
<Parameter: 'water - sld' , value=9.41192 (fixed) , bounds=[-inf, inf]>
<Parameter:'water - isld' , value=0 (fixed) , bounds=[-inf, inf]>
<Parameter:'water - rough', value=3.39189 +/- 0.0124, bounds=[0.5, 7.0]>
<Parameter:'water - volfrac solvent', value=0 (fixed) , bounds=[0.0, 1.0]>
%% Cell type:markdown id:44f8d495 tags:
## Plot the results
%% Cell type:code id:b340cc4f tags:
``` python
objective.plot()
plt.xlabel(r'$q_z$ ($\AA^{—1}$)')
plt.ylabel("logR")
plt.show()
fig, ax1 = plt.subplots()
# First y-axis (SLD)
ax1.plot(*structure.sld_profile())
ax1.set_xlabel('Distance / $\AA$')
ax1.set_ylabel('SLD / $10^{-6} \AA^{-2}$', color='b')
ax1.tick_params(axis='y', labelcolor='b')
# Define transformation functions
def sld_to_ed(sld): # sld in 10^-6 Å⁻²
return (sld * 1e-6) / r_el # e⁻/ų
def ed_to_sld(ed): # e⁻/ų
return (ed * r_el) * 1e6 # back to 10^-6 Å⁻²
# Add secondary y-axis
ax2 = ax1.secondary_yaxis('right', functions=(sld_to_ed, ed_to_sld))
ax2.set_ylabel('Electron Density / e$^- \ \AA^{-3}$', color='r')
ax2.tick_params(axis='y', labelcolor='r')
plt.title('SLD and Electron Density Profile')
plt.tight_layout()
plt.show()
```
%% Output
%% Cell type:markdown id:bdfd1460 tags:
# Statistics (for later use)
%% Cell type:code id:c6b0fc85 tags:
``` python
fitter.sample(100, pool=-1)
```
%% Output
100%|██████████| 100/100 [00:23<00:00, 4.24it/s]
[MCMCResult(name='water - rough', param=Parameter(value=3.3920854514375325, name='water - rough', vary=True, bounds=Interval(lb=0.5, ub=7.0), constraint=None), stderr=0.01276232987084347, chain=array([[3.39184015, 3.39438269, 3.41649244, ..., 3.39849706, 3.41172303,
3.39824873],
[3.39176784, 3.39438269, 3.40740042, ..., 3.39664135, 3.41172303,
3.39816917],
[3.38919273, 3.38919182, 3.40740042, ..., 3.39664135, 3.40834449,
3.40034635],
...,
[3.39544814, 3.38676058, 3.39825959, ..., 3.39790449, 3.39701153,
3.39476713],
[3.39365428, 3.38676058, 3.39731078, ..., 3.39797642, 3.39647369,
3.39699511],
[3.39364197, 3.39110629, 3.39731078, ..., 3.39819214, 3.40205574,
3.39701509]]), median=3.3920854514375325)]
%% Cell type:code id:a7fa301e tags:
``` python
objective.plot(samples=100);
```
%% Output
%% Cell type:code id:89b2d45c tags:
``` python
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

The place to be for your source code