# -*- coding: utf-8 -*-
"""
basic variables
"""
from __future__ import print_function, division
from scipy.optimize import fsolve
from scipy.optimize import leastsq, curve_fit
from scipy import constants as con
# from matplotlib.pylab import *
# from io import StringIO
# import fileinput
import numpy as np
# ----------------------------------------------------------------------------
# --------------------- Konstanten und Flags ---------------------------------
# ----------------------------------------------------------------------------
A = con.angstrom #: Angstroem in meters
m_He3 = con.m_u * 3.0160293 # Mass of a helium atoms (modified He3)
m_He4 = con.m_u * 4.0026
m_Te = con.m_u * 127.6
a_Bi2Te3_111 = 4.38 * A
td_Bi2Te3_111 = 80 # K maybe wrong
hb = con.hbar # reduced Planck's constant
[docs]def h2m(m):
"""Energy-Momentum-factor"""
return hb**2/(2*m)
meV = mel = con.eV/1E3 # milli-electronvolt
kb = 1.3806504E-23
[docs]def gauss(p, sectx):
return p[1]*np.exp(-(sectx-p[0])**2/p[2]**2)
[docs]def gigagauss(p, sectx):
"""
Funktion to model a more complex signal using a threshold,
a peak at dE=0 and several gauss functions given in p
"""
A1 = p[0]+gauss(np.array([0, abs(p[1]), p[2]]), sectx)
for ind in np.arange((p.size-3)/3)+1:
ing = ind*3
A1 = A1 + gauss(np.array([p[ing], abs(p[ing+1]), p[ing+2]]), sectx)
return A1
[docs]def gigagauss_generic(p, sectx):
"""
Funktion to model a more complex signal using only
several gauss functions given in p
"""
A1 = np.zeros(sectx.size)
for ind in np.arange((p.size)/3):
ing = ind*3
A1 = A1 + gauss(np.array([p[ing], abs(p[ing+1]), p[ing+2]]), sectx)
return A1
[docs]def fitgigagauss_generic(sectx, secty, p):
# Definition of a Fitfunction
def errf(p, x, y):
return gigagauss_generic(p, x)-y
# Optimieren
p, success = leastsq(errf, p, args=(sectx, secty), maxfev=10000)
return p, success
[docs]def fitgauss(sectx, secty, x0, a=1, s=0.9):
# Setzen der initialen Parameterliste
p = np.array([x0, a, s])
# Definieren der Fitfunktion und der dazugehoerigen errorfunct
def errf(p, x, y):
return gauss(p, x)-y
# Optimieren
p, success = leastsq(errf, p, args=(sectx, secty), maxfev=10000)
return [p[0], p[1], p[2], success]
[docs]def fitgauss_cf(sectx, secty, x0, a=1, s=0.9):
# Setzen der initialen Parameterliste
p = np.array([x0, a, s])
# Definieren der Fitfunktion und der dazugehoerigen errorfunct
def errf(p, x):
return gauss(p, x)
# Optimieren
p, _ = curve_fit(errf, sectx, secty, p0=p, maxfev=10000)
return p
[docs]def fitgigagauss(sectx, secty, p):
# Definition of a Fitfunction
def errf(p, x, y):
return gigagauss(p, x)-y
# Optimieren
p, success = leastsq(errf, p, args=(sectx, secty), maxfev=10000)
return p, success
[docs]def find_inters(func1, func2, x0=0.0):
# Schnittpunkt-Auffindefunktion
def f(x):
return func1(x)-func2(x)
y, dict, du, my = fsolve(f, x0, full_output=1, maxfev=1000)
return [y, du]
[docs]def myunique(arra, thr_pos=3, col=0):
thresh = int(10**thr_pos)
if arra.ndim > 1:
index = [int(x) for x in arra[:, col]*thresh]
uniq = np.unique(index, return_index=True)
return arra[uniq[1], :]
else:
return arra
[docs]def tempsweep(Tofobj, Spec1, Spec2):
"""
tempsweep
Function to calculate the Nozzle temperature of a TOFmeasurement
using two specular measurements (one before, one after)
"""
TEMP1 = Spec1.TN
TIME1 = Spec1.mes_time
TEMP2 = Spec2.TN
TIME2 = Spec2.mes_time
# Calculatin linear parameters
k = (TEMP2-TEMP1)/(TIME2-TIME1)
d = TEMP2 - k*TIME2
# Returning correspondant Temperature value
return k*Tofobj.mes_time+d
[docs]def findQ(TOFobj, AppObj, E, BZ):
kb = 1.3806504E-23 # Bolzmanns Constant
hb = 1.054571628E-34 # H-Bar (Planck's constant)
mel = 1.602E-22 # Millielectronvolt - Energy in Joule
m = 1.660538782E-27*4.002602 # Mass of a helium atom in kg
Ang = (AppObj.TSD/2.0-(TOFobj.theta_app-AppObj.specpos))*np.pi/180.0
Tem = TOFobj.TN
TSD = AppObj.TSD*np.pi/180.0
Ei = 5*kb*Tem/2.0 # Helium Atom Energy in J
ki = np.sqrt(2*m*Ei)/hb # Helium momentum m**-1
Q = ki*np.sin(Ang) * \
(np.sqrt((E*mel/Ei+1)*((np.sin(TSD-Ang)**2)/(np.sin(Ang)**2)))-1)
# remapping
Q = abs(Q)
hugo = True
while hugo:
hugo = False
if Q > BZ:
Q = 2*BZ-Q
Q = abs(Q)
hugo = True
return Q
# Define scanc function
[docs]def sancurvec(appa, TN, theta, Gbound):
kb = 1.3806504E-23 # Bolzmanns Constant
hb = 1.054571628E-34 # H-Bar (Planck's constant)
mel = 1.602E-22 # Millielectronvolt - Energy in Joule
m = 1.660538782E-27*4.002602 # Mass of a helium atom in kg
TSD = appa.TSD*np.pi/180.0
Ang = (theta)*np.pi/180.0
Tem = TN
Ei = 5*kb*Tem/2.0 # Helium Atom Energy in J
ki = np.sqrt(2*m*Ei)/hb # Helium momentum m**-1
Q = np.arange(-4E10, 4E10, 1E8)
# fig = figure()
# ax = fig.add_subplot(111)
dE = (hb**2/(2*m)) * \
((ki*np.sin(Ang)+Q)**2/(np.sin(TSD-Ang)**2)) - \
(hb**2/(2*m))*ki**2
dE = dE/mel
# remap
Q = abs(Q)
dE = abs(dE)
GG = Gbound
hugo = True
while hugo:
hugo = False
for run in range(Q.size):
if Q[run] > GG:
Q[run] = 2*GG-Q[run]
Q = abs(Q)
hugo = True
return (Q*1E-10, dE)
[docs]def scanc(TN, theta, GG=14583570019.449413, TSD=91.75, m=m_He4):
"""
calculate the scancurve
@param TN: Nozzle Temp in K
@param theta: incident angle in degrees
@param GG: Brillouin zone boundary in 1/m
@param TSD: Source detector angle in degrees
@param m: beam mass im kg
@return: tuple of Q and E in 1e9 1/m and meV
"""
TSD = TSD*np.pi/180.0
Ang = theta*np.pi/180.0
Ei = 5*kb*TN/2.0 # Helium Atom Energy in J
ki = np.sqrt(2*m*Ei)/hb # Helium momentum m**-1
Q = np.linspace(-2*ki, 2*ki, 4001)
dE = (hb**2/(2*m))*((ki*np.sin(Ang)+Q)**2/(np.sin(TSD-Ang)**2)-ki**2)
dE = dE/mel
# from matlab
ind_p = dE.argmin()
dE = dE[ind_p:]
Q = Q[ind_p:]
# remap
# dE = abs(dE)
Q = abs(Q)
hugo = True
while hugo:
hugo = False
for run in range(Q.size):
if Q[run] > GG:
Q[run] = 2*GG-Q[run]
Q = abs(Q)
hugo = True
return (Q*1E-9, dE)
[docs]def prettyprint_elas(CONTAINER):
"""
prettyprint(CONTAINER)
CONTAINER ... output of sim_elastic measurement
Function to provide easy-to-read access to the relevant elastic positions
Patrick Kraus, 11.09.2012
"""
LABCONT = {
-1: "Elastic Scattering Peaks",
0: "Threshold Resonance Positions",
1: "Resonance with n=0",
2: "Resonance with n=1",
3: "Resonance with n=2"}
for ind in LABCONT.keys():
dummy = np.find(CONTAINER[:, 1] == ind)
try:
dummy = dummy[0]
print(LABCONT[ind]+":")
POSS = ""
VALS = ""
for jnd in np.where(CONTAINER[:, 1] == ind)[0]:
POSS = POSS + "("+str(int(CONTAINER[jnd, 2])) + \
","+str(int(CONTAINER[jnd, 3]))+")"+"\t"
for jnd in np.where(CONTAINER[:, 1] == ind)[0]:
VALS = VALS + str(CONTAINER[jnd, 0])[0:5]+"\t"
print(POSS)
print(VALS)
print("\n")
except:
pass