Skip to content
Snippets Groups Projects
Commit f9eb5400 authored by Roger Wolf's avatar Roger Wolf
Browse files

cleanup of tools

parent 0709dcdc
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Resonanz/tools/resonance.ipynb
+ 2
− 49
View file @ f9eb5400
%% Cell type:markdown id:001d51ab-bdf8-40b5-9aad-54e33181e05c tags:
# Anpassung der Resonanzkurve und Bestimmung von $\Delta\Omega$
%% Cell type:markdown id:fffee2ad-f784-4614-bddf-491ce9de9f05 tags:
Mit dem folgenden Code-Fragment zeigen wir Ihnen:
* Wie man ein **geeignetes Modell** an die prozessierten Daten anpasst.
* Wie man die **angepasste Funktion** gewinnt.
* Wie man in mit Hilfe der *scypi*-Funktionen [fmin](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin.html) und [fsolve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html) den **$x$-Wert für das Maximum/Minimum oder bestimmte Schnittpunkte** bestimmt.
%% Cell type:code id:edc41a5d-2c39-4391-8f25-3a8a3a0e3271 tags:
``` python
import numpy as np
from kafe2 import XYContainer
# We
Omega = [ 1.0, 1.8, 2.6, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 4.0, 4.5, 5.0]
phi0 = [0.05, 0.10, 0.20, 0.30, 0.40, 0.75, 1.50, 1.75, 1.10, 0.47, 0.20, 0.10, 0.05]
# Fill XYContainer
xy_data = XYContainer(x_data=Omega, y_data=phi0)
# We add some arbitrary uncertaintes here
xy_data.add_error(axis="x", err_val=0.0125)
xy_data.add_error(axis="y", err_val=0.0500)
xy_data.label = "Data"
from kafe2 import Fit, Plot, ContoursProfiler
# Fit model, we start of with a damped oscillation, and potentially extend by
# a linear extra term to additional frictional effects of the wheel. The POIs
# (parameters of interest) are lambda0 and omega0.
def resonance(x, k=1.5, omega0=3.1, lambda0=0.2):
return k/np.sqrt((omega0**2-x**2)**2+(2*lambda0*x)**2)
# Declare and do the fit with a bit of cosmetic customization
fit = Fit(xy_data, model_function=resonance)
#fit.assign_model_function_latex_name(r"\Omega_{X}")
#fit.assign_parameter_latex_names(
# lambda0 = r"\lambda_{0}",
# x0 = r"x_{0}",
# omega0 = r"\omega_{0}",
# phi0 = r"\phi_{0}",
# a = r"a"
#)
# Set limits for parameters if needed.
#fit.limit_parameter("phi0", lower=-180., upper=180.)
#fit.limit_parameter("x0", lower=0.)
# Do the fit
results=fit.do_fit()
# Return a reasonably clean report
fit.report(show_data=False, show_model=False)
# Plot the result with a bit of customization
plot = Plot(fit_objects=fit)
plot.customize('model_line', 'label', [(0, r'Model Pohls Wheel ($\Omega_{X}$)')])
plot.customize('model_line', 'color', [(0, 'darkblue')])
plot.customize('model_line', 'linestyle', [(0, 'solid')])
plot.customize('model_error_band', 'label', [(0, r'$\pm1\sigma$')])
plot.customize('model_error_band', 'color', [(0, 'lightsteelblue')])
plot.x_label = r'$\Omega\,(2\pi\,\mathrm{Hz})$'
plot.y_label = r'$\varphi_{0}(\Omega)\,(\mathrm{rad})$'
plot.plot()
# Create contour plots, this step should only be used when using less than 300
# data points, since it may take really long otherwise
#cpf = ContoursProfiler(fit)
#cpf.plot_profiles_contours_matrix()
# Show plot in notebook
plot.show()
```
%% Output
###############
# Fit Results #
###############
Model Parameters
================
k = 0.873 +/- 0.048
omega0 = 3.3756 +/- 0.0063
lambda0 = -0.0639 +/- 0.0093
Model Parameter Correlations
============================
k omega0 lambda0
======= ======= =======
k 1.0 0.2624 -0.7683
omega0 0.2624 1.0 -0.4319
lambda0 -0.7683 -0.4319 1.0
Cost Function
=============
Cost function: chi-square (with covariance matrix)
chi2 / ndf = 13.08 / 10 = 1.308
chi2 probability = 0.220
OrderedDict({'k': np.float64(0.8729181747989906), 'omega0': np.float64(3.375592203784703), 'lambda0': np.float64(-0.06393978582236191)})
%% Cell type:code id:8561b826-d5f3-4ae4-bf68-1ba5cd32c955 tags:
``` python
from scipy.optimize import fsolve, fmin
# In case you want to evaluate the fitted function at any point per hand (this
# can be very useful to convince yourself that what is done below does what you
# expect it to do)
#fnc_values = fit.eval_model_function(x=np.array([1.0, 2.5, 3.4, 4.5]))
# One way to define the function with the fitted parameters
def fitted_fnc(x):
return resonance(x,
k = results["parameter_values"]["k"],
omega0 = results["parameter_values"]["omega0"],
lambda0 = results["parameter_values"]["lambda0"]
)
# Find the location of the maximum of fitted_fnc on the x-axis
x_max = fmin(lambda x: -fitted_fnc(x), 0)
# One way to define DeltaOmega=x_upper-x_lower. Be sure that the root finding
# by fslove is close to the roots that you expect. In this case we chose start
# values for the search 5% left and right of x_max
x_upper = fsolve(lambda x: fitted_fnc(x)-fitted_fnc(x_max)/np.sqrt(2),1.05*x_max)
x_lower = fsolve(lambda x: fitted_fnc(x)-fitted_fnc(x_max)/np.sqrt(2),0.95*x_max)
print("Results of the search for DeltaOmega:")
print("-------------------------------------")
print("Omega_max = ", x_max)
print("DeltaOmega = ", x_upper-x_lower)
print("Q = ", x_max/(x_upper-x_lower))
```
......
Resonanz/tools/swing.ipynb
+ 1
− 1
View file @ f9eb5400
%% Cell type:markdown id:001d51ab-bdf8-40b5-9aad-54e33181e05c tags:
# Werkzeuge zur Darstellung und Modellanpassungen für den Versuch Resonanz
%% Cell type:markdown id:a7fd56cb-ed0a-4516-86e6-9fe10a9dd806 tags:
Mit den folgenden Code-Fragmenten zeigen wir Ihnen:
* Wie man die Daten des CASSY-Messsystems mit PhyPraKit **einließt**.
* Wie man die eingelesenen Daten auf ein handhabbares Maß **reduziert**.
* Wie man die Daten mit Hilfe von [*Spline*](https://de.wikipedia.org/wiki/Spline)-Funktionen **interpoliert**.
* Wie man die angepassten *Spline*-Funktionen numerisch **differenziert**.
* Wie man die eingelesenen Daten mit Hilfe von *matplotlib* **darstellt**.
* Wie man an die eingelesenen Daten ein **geeignetes Modell anpasst**, um die Parameter $\Omega_{0}$ und $\lambda_{0}$ zu bestimmen.
Alle im Folgenden gezeigten Code-Fragmente lassen sich geeignet kombinieren.
%% Cell type:markdown id:d7c21a7c-b21d-474a-8692-715e51386cd1 tags:
## Einlesen der Daten des CASSY-Messsystems mit PhyPraKit
%% Cell type:markdown id:6b48a8c8-ffe3-484b-8434-0ffbc9006ece tags:
Mit dem folgenden Code-Fragment zeigen wir Ihnen:
* Wie man mit der PhyPrakit Funktion *labxParser* Daten aus einer *labx*-Datei des CASSY-Messsystems **einließt**.
* Bei der vom CASSY-Messsystem ausgegebenen *labx*-Datei handelt es sich um eine **gezippte *xml*-Datei**.
* Wenn Sie die Datei entpacken findent Sie die *xml*-Datei in einem entsprechend entpackten Verzeichnis.
* Diese Datei können Sie mit der Funktion *labxParser* wie in der folgenden Code-Zelle gezeigt einlesen.
* Wir verwenden hierzu die Beispiel-Datei `CassyExample.xml`, aus dem *tools*-Verzeichnis.
%% Cell type:code id:c5723f54-17d7-4ba2-ac75-56ce00cbac87 tags:
``` python
import numpy as np
from PhyPraKit.PhyPraKit.phyTools import labxParser
from PhyPraKit.phyTools import labxParser
# Read data from CASSY;
# column_names: hosts the column names;
# data: hosts the individual columns as lists.
# The argument prlevel steers how much information is printed to screen. We
# set it to 1 here, so that you can learn something about the structure of
# the labx file. If you want to shut the function up set prlevel=0. Check
# help(labxParser) to learn more about this function
column_names, data = labxParser("CassyExample.xml", prlevel=1)
# For this example we collect the oservables "Zeit" and "Winkel"
t=[]; phi=[]
# We cut off the first 400 and last 450 lines to reduce the amount of data
# points from ~1300 to <500. Downsampling is not really an option for this
# example, since the oscillation period is very small
LOWER_CUT= 400
UPPER_CUT=-450
for i, tag in enumerate(column_names):
# Column names are encoded in a more cryptic way than necessary by Leybold
# experts :-(; uncomment the following line to check this encoding
#print(tag.split(":"))
tcn = tag.split(":")[1]
if tcn == "Zeit":
# Copy full column in variable t
t = np.array(data[i][LOWER_CUT:UPPER_CUT])
if tcn == "Winkel":
# Copy full column in variable phi
phi = np.array(data[i][LOWER_CUT:UPPER_CUT])
#print(t, phi)
```
%% Cell type:markdown id:2d624120-49a5-4f05-aa56-3aa200a04726 tags:
## Dastellung der eingelesenen Daten
%% Cell type:markdown id:3dcacf72-00f4-4bff-8aa9-cba32b5425da tags:
Mit dem folgenden Code-Fragment zeigen wir Ihnen:
* Wie man sich mit Hilfe von *numpy* die (numerische) Ableitung $\partial_{t}\varphi$ verschafft.
* Wie man $\varphi(t)$, $\dot{\varphi}(t)$ und ein Phasenraumportrait $(\varphi, \dot{\varphi})(t)$ mit Hilfe von *matplotlib* darstellt.
%% Cell type:code id:f9886e0b-835c-4c91-8b3e-23d662322e3a tags:
``` python
# Choose difference between first an second entry in t as estimate for dt
dt = t[1]-t[0]
# Apply numerical derivative to phi
dphi = np.gradient(phi, dt)
```
%% Cell type:code id:3c27f647-1393-40a6-9c4f-fbe6b4cda016 tags:
``` python
import matplotlib.pyplot as plt
# Display phi and dphi
fig = plt.figure("Winkel und Winkelgeschwindigkeit", figsize=(12.0, 12.0))
ax1 = fig.add_subplot(2, 1, 1)
ax1.set_xlabel("t (s)")
ax1.set_ylabel(r"$\varphi(t)$ (rad)")
ax1.plot(t, phi , color="darkblue")
ax1.grid()
ax2 = fig.add_subplot(2, 1, 2)
ax2.set_xlabel("t (s)")
ax2.set_ylabel(r"$\dot{\varphi}(t)$ (rad/s)")
ax2.plot(t, dphi , color="darkcyan")
ax2.grid()
plt.show()
```
%% Cell type:code id:750e98b2-bea9-4bd0-912b-8abb19d1f8db tags:
``` python
from scipy import interpolate
# Do a nice cubic spline interpolation with help of the scipy function
# interpolate
cs_phi = interpolate.UnivariateSpline(t, phi, s=0)
# Get dphi as derivative of the spline approximation to phi
cs_dphi = cs_phi.derivative()
# Display Phasenraumprotrait
tplt = np.linspace(t[0], t[-1], 50000)
fig = plt.figure("Phasenraumportrait", figsize=(6.0, 6.0))
plt.xlabel(r"$\varphi(t)$ (rad)")
plt.ylabel(r"$\dot{\varphi}(t)$ (rad/s)")
plt.plot(cs_phi(tplt), cs_dphi(tplt) , color="darkblue")
plt.grid()
plt.show()
```
%% Cell type:markdown id:f41154f9-1c50-4c21-bc4d-3d270ab3c455 tags:
## Anpassung eines geeigneten Modells an die Daten
%% Cell type:markdown id:0122e806-5095-436e-9c0a-2224e42b8688 tags:
Mit dem folgenden Code-Fragment zeigen wir Ihnen:
* Wie man ein geeignetes Modell an die extrahierten Daten anpasst.
%% Cell type:code id:22588b29-966b-4adf-9756-c53eec7d7a17 tags:
``` python
from kafe2 import XYContainer
# Fill XYContainer
xy_data = XYContainer(x_data=t, y_data=phi)
# We add some typical uncertaintes here
xy_data.add_error(axis="x", err_val=0.0125)
xy_data.add_error(axis="y", err_val=0.0500)
xy_data.label = "Data"
from kafe2 import Fit, Plot, ContoursProfiler
# Fit model, we start of with a damped oscillation, and potentially extend by
# a linear extra term to additional frictional effects of the wheel. The POIs
# (parameters of interest) are lambda0 and omega0.
def harmonic_damped(x, lambda0=0.01, x0=3.5, omega0=3.1, phi0=0., a=-1.0):
return x0*np.exp(-x*lambda0)*np.cos(omega0*x+phi0)+a*x
# Declare and do the fit with a bit of cosmetic customization
fit = Fit(xy_data, model_function=harmonic_damped)
fit.assign_model_function_latex_name(r"\Omega_{X}")
fit.assign_parameter_latex_names(
lambda0 = r"\lambda_{0}",
x0 = r"x_{0}",
omega0 = r"\omega_{0}",
phi0 = r"\phi_{0}",
a = r"a"
)
# Set limits for parameters if needed.
fit.limit_parameter("phi0", lower=-180., upper=180.)
fit.limit_parameter("x0", lower=0.)
# Do the fit
fit.do_fit()
# Return a reasonably clean report
fit.report(show_data=False, show_model=False)
# Plot the result with a bit of customization
plot = Plot(fit_objects=fit)
plot.customize('model_line', 'label', [(0, r'Model Pohls Wheel ($\Omega_{X}$)')])
plot.customize('model_line', 'color', [(0, 'darkblue')])
plot.customize('model_line', 'linestyle', [(0, 'solid')])
plot.customize('model_error_band', 'label', [(0, r'$\pm1\sigma$')])
plot.customize('model_error_band', 'color', [(0, 'lightsteelblue')])
plot.x_label = r'$t\,(\mathrm{s})$'
plot.y_label = r'$\varphi\,(rad)$'
plot.plot()
# Create contour plots, this step should only be used when using less than 300
# data points, since it may take really long otherwise
#cpf = ContoursProfiler(fit)
#cpf.plot_profiles_contours_matrix()
# Show plot in notebook
plot.show()
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Help | Imprint | Privacy policy | Accessibility | Contact