Source code for haslib.science_basics

# -*- 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