xrs_utilities.py 161 KB
Newer Older
1
2
3
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
4
from six.moves import range
5
#!/usr/bin/python
christoph's avatar
christoph committed
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
# Filename: xrs_utilities.py

#/*##########################################################################
#
# The XRStools software package for XRS spectroscopy
#
# Copyright (c) 2013-2014 European Synchrotron Radiation Facility
#
# This file is part of the XRStools XRS spectroscopy package developed at
# the ESRF by the DEC and Software group and contains practical functions, 
# most of which are translated from Matlab functions from the University of
# Helsinki Electronic Structure Laboratory.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
#############################################################################*/
__author__ = "Christoph J. Sahle - ESRF"
__contact__ = "christoph.sahle@esrf.fr"
__license__ = "MIT"
__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
42
43
44

import os
import math
christoph's avatar
christoph committed
45
import copy
46
47
48
49

import numpy as np
import array as arr
import matplotlib.pyplot as plt
christoph's avatar
christoph committed
50
import pickle
51
52
import traceback
import sys
53

christoph's avatar
christoph committed
54
from matplotlib.widgets import Cursor
55
56
57
58
59
from itertools import groupby
from scipy.integrate import trapz
from scipy import interpolate, signal, integrate, constants, optimize
from re import findall
from scipy.ndimage import measurements
60
from scipy.optimize import leastsq, fmin, fsolve, minimize, nnls
61
from scipy.interpolate import Rbf, RectBivariateSpline
christoph's avatar
christoph committed
62
63
from scipy.integrate import odeint

christoph's avatar
christoph committed
64
# data_installation_dir = os.path.join( os.path.dirname(os.path.abspath(__file__)),"..","..","..","..","share","xrstools","data")
65
# data_installation_dir = os.path.abspath('.')
christoph's avatar
christoph committed
66
67
68

# whne you test the file from its source directory you, /data sits on level above. In this case you can
# work around it by creating a link in ./ to ../data
Alessandro Mirone's avatar
OK    
Alessandro Mirone committed
69
data_installation_dir = os.path.join( os.path.dirname(os.path.abspath(__file__)),"resources",  'data')
70

71
# os.path.join(getattr(install_cmd, 'install_lib'),"xrstools"+version,"..","..","..","..","share","xrstools","data")
christoph's avatar
christoph committed
72
73


74
75
76
77
78
def diode(current, energy, thickness=0.03):
    """ **diode**
    Calculates the number of photons incident for a Si PIPS diode.

    Args:
alex's avatar
alex committed
79
80
81
      * current (float): Diode current in [pA].
      * energy (float): Photon energy in [keV].
      * thickness (float): Thickness of Si active layer in [cm].
82
83

    Returns:
alex's avatar
alex committed
84
      * flux (float): Number of photons per second.
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

    Function adapted from Matlab function by S. Huotari.
    """
    t = thickness # thickness of diode in cm
    
    # Total cross-section of absorbed energy for Si (Storm & Israel)
    sx = np.array([2,3,4,5,6,8,10,15,20,30,40,50,60,80,100,150,200])
    s  = np.array([125000,44600,20600,11000,6600,2890,1510,448,186, \
                       53.3,21.9,11.2, 6.61,3.18,2.06,1.41,1.34])

    my = np.exp(np.interp(energy, sx , np.log(s)))
    my = (0.02144*2.32)*my

    n_ph= energy*(1.0-np.exp(-tauphoto(14,energy)*2.32*t))
    n_ph=n_ph + (energy)*(1.0-np.exp(-sigmainc(14,energy)*2.32*t))
    n_ph=n_ph/0.0036

    n_ph = energy*(1.0-np.exp(-my*t))/0.0036 # number of electrons/photon
    n_pA = current/1.6022e-7                     # number of electrons per pA
    print('The photon flux is: %E'%(n_pA/n_ph))
    return  n_pA/n_ph

def cshift(w1, th):
    """ **cshift**
    Calculates Compton peak position.

    Args:
alex's avatar
alex committed
112
113
      *  w1 (float, array): Incident energy in [keV].
      *  th (float): Scattering angle in [deg].
114
115

    Returns:
alex's avatar
alex committed
116
      *  w2 (foat, array): Energy of Compton peak in [keV].
117
118
119
120
121
122
123
124
125
126
127
128

    Funktion adapted from Keijo Hamalainen.

    """
    return w1/(1+w1/510.967*(1-np.cos(th/180*np.pi)))


def tauphoto(Z, energy, logtablefile=os.path.join(data_installation_dir,'logtable.dat')):
    """ **tauphoto**
    Calculates Photoelectric Cross Section in cm^2/g using Log-Log Fit.

    Args:
alex's avatar
alex committed
129
130
      * z (int or string): Element number or elements symbol.
      * energy (float or array): Energy (can be number or vector)
131
132

    Returns:
alex's avatar
alex committed
133
      * tau (float or array): Photoelectric cross section in [cm**2/g]
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166

    Adapted from original Matlab function of Keijo Hamalainen.
    """

    en = np.array([])
    en = np.append(en,energy)
    logtable = np.loadtxt(logtablefile)
    # find the right places in logtable
    if not isinstance(Z,int):
        Z = element(Z)
    try:
        ind = list(logtable[:,0]).index(Z)
    except:
        print( 'no such element in logtable.dat')

    c = np.array(logtable[ind:ind+5,:]) # 5 lines that corresponds to the element
    tau_i = np.zeros((4, len(en)))
    for ii in range(4):
        for jj in range(4):
            tau_i[ii,:] = tau_i[ii,:] + c[jj+1,ii]*np.log(en)**(jj)
    tau2 = np.zeros(len(en))
    tau2 = (en<c[0,1])*1*tau_i[0,:] + \
            np.logical_and(en>=c[0,1] , en<c[0,2])*1*tau_i[1,:] + \
            np.logical_and(en>=c[0,2] , en<c[0,3])*1*tau_i[2,:] + \
            (en>=c[0,3])*1*tau_i[3,:]
    return np.exp(tau2)*0.6022/c[0,4]

def sigmainc(Z, energy, logtablefile=os.path.join(data_installation_dir,'logtable.dat')):
    """ **sigmainc**
    Calculates the Incoherent Scattering Cross Section 
    in cm^2/g using Log-Log Fit.

    Args:
alex's avatar
alex committed
167
168
      * z (int or string): Element number or elements symbol.
      * energy (float or array): Energy (can be number or vector)
169
170

    Returns:
alex's avatar
alex committed
171
      * tau (float or array): Photoelectric cross section in [cm**2/g]
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194

    Adapted from original Matlab function of Keijo Hamalainen.

    """

    en = np.array([])
    en = np.append(en,energy)
    logtable = np.loadtxt(logtablefile)
    # find the right places in logtable
    if not isinstance(Z,int):
        Z = element(Z)
    try:
        ind = list(logtable[:,0]).index(Z)
    except:
        print( 'no such element in logtable.dat')

    c = np.array(logtable[ind:ind+5,:]) # 5 lines that corresponds to the element
    sigmai=0
    for jj in range(4):
        sigmai = sigmai + c[jj+1,5]*np.log(energy)**(jj)

    return np.exp(sigmai)*0.6022/c[0,4]

195
196
197
198
199
def Rx(chi, degrees=True):
    """ **Rx**
    Rotation matrix for vector rotations around the [1,0,0]-direction.

    Args:
alex's avatar
alex committed
200
201
      * chi   (float) : Angle of rotation.
      * degrees(bool) : Angle given in radians or degrees.
202
203

    Returns:
alex's avatar
alex committed
204
       * 3x3 rotation matrix.
205
206
207
208
209
210
211
212
213
214
    """
    if degrees:
        chi = np.radians(chi)
    return np.array([[1,0,0],[0, np.cos(chi), -np.sin(chi)], [0, np.sin(chi), np.cos(chi)]])

def Ry(phi, degrees=True):
    """ **Ry**
    Rotation matrix for vector rotations around the [0,1,0]-direction.

    Args:
alex's avatar
alex committed
215
216
      * phi   (float) : Angle of rotation.
      * degrees(bool) : Angle given in radians or degrees.
217
218

    Returns:
alex's avatar
alex committed
219
       * 3x3 rotation matrix.
220
221
222
223
224
225
226
227
228
229
    """
    if degrees:
        phi = np.radians(phi)
    return np.array([[np.cos(phi), 0, np.sin(phi)],[0, 1, 0],[-np.sin(phi), 0, np.cos(phi)]])

def Rz(omega, degrees=True):
    """ **Rz**
    Rotation matrix for vector rotations around the [0,0,1]-direction.

    Args:
alex's avatar
alex committed
230
231
      * omega (float) : Angle of rotation.
      * degrees(bool) : Angle given in radians or degrees.
232
233

    Returns:
alex's avatar
alex committed
234
      * 3x3 rotation matrix.
235
236
237
238
239
    """
    if degrees:
        omega = np.radians(omega)
    return np.array([[np.cos(omega), -np.sin(omega), 0],[np.sin(omega), np.cos(omega), 0],[0,0,1]])

christoph's avatar
christoph committed
240
241

def Phi(phi, degrees=True):
242
    """ rotation around (0,1,0), neg sense
christoph's avatar
christoph committed
243
244
245
    """
    if degrees:
        phi = np.radians(phi)
246
247
    return np.array([[np.cos(phi), 0, np.sin(phi)],[0, 1, 0],[-np.sin(phi), 0, np.cos(phi)]])
    
christoph's avatar
christoph committed
248
def Chi(chi, degrees=True):
249
    """ rotation around (1,0,0), pos sense
christoph's avatar
christoph committed
250
251
252
    """
    if degrees:
        chi = np.radians(chi)
253
254
    return np.array([[1,0,0],[0, np.cos(chi), -np.sin(chi)], [0, np.sin(chi), np.cos(chi)]])
    
christoph's avatar
christoph committed
255
def Omega(omega, degrees=True):
256
    """ rotation around (0,0,1), pos sense
christoph's avatar
christoph committed
257
258
259
    """
    if degrees:
        omega = np.radians(omega)
260
    return np.array([[np.cos(omega), -np.sin(omega), 0],[np.sin(omega), np.cos(omega), 0],[0,0,1]])
christoph's avatar
christoph committed
261
262

def get_UB_Q(tthv, tthh, phi, chi, omega, **kwargs):
263
264
265
266
267
268
    """ **get_UB_Q**
    Returns the momentum transfer and scattering vectors for given
    FOURC spectrometer and sample angles. U-, B-matrices and 
    incident/scattered wavelength are passed as keyword-arguments.

    Args:
alex's avatar
alex committed
269
270
271
272
273
274
275
276
277
278
279
280
      * tthv (float): Spectrometer vertical 2Theta angle.
      * tthh (float): Spectrometer horizontal 2Theta angle.
      * chi (float): Sample rotation around x-direction.
      * phi (float): Sample rotation around y-direction.
      * omega (float): Sample rotation around z-direction.

      * kwargs (dict): Dictionary with key-word arguments:
          * kwargs['U'] (array): 3x3 U-matrix Lab-to-sample transformation.
          * kwargs['B'] (array): 3x3 B-matrix reciprocal lattice 
            to absolute units transformation.
          * kwargs['lambdai'] (float): Incident x-ray wavelength in Angstrom.
          * kwargs['lambdao'] (float): Scattered x-ray wavelength in Angstrom.
281
282

    Returns:
alex's avatar
alex committed
283
284
285
      * Q_sample  (array): Momentum transfer in sample coordinates.
      * Ki_sample (array): Incident beam direction in sample coordinates.
      * Ko_sample (array): Scattered beam direction in sample coordinates.
286
287
288
289

    """
    U       = kwargs['U']
    B       = kwargs['B']
290
291
    Lab     = kwargs['Lab']
    beam_in = kwargs['beam_in']
292
293
294
    lambdai = kwargs['lambdai']
    lambdao = kwargs['lambdao']
    # scattering vectors in laboratory frame
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
    Ki_test = 2.0*np.pi/lambdai * beam_in/np.linalg.norm(beam_in)
    Ko_test = 2.0*np.pi/lambdao * Rz(tthh).dot(Ry(tthv)).dot(beam_in/np.linalg.norm(beam_in))
    # h_lab = Omega*Chi*Phi*U*B*h_cryst (Busing equ. 19)
    # invert
    # h_cryst = B_inv*U_inv*Phi_inv*Chi_inv*Omega_inv*h_lab
    Q_test    = Ko_test - Ki_test
    Phi_inv   = Phi(phi).T #np.linalg.inv(Phi(phi))
    print(Phi_inv)
    Chi_inv   = Chi(chi).T #np.linalg.inv(Chi(chi))
    Omega_inv = Omega(omega).T #np.linalg.inv(Omega(omega))
    U_inv     = np.linalg.inv(U)
    B_inv     = np.linalg.inv(B)
    Q_sample  = np.matmul(B_inv ,np.matmul(U_inv , np.matmul( Phi_inv , np.matmul(Chi_inv, np.matmul( Omega_inv, Q_test)))))
    Ki_sample = np.matmul(B_inv ,np.matmul(U_inv , np.matmul( Phi_inv , np.matmul(Chi_inv, np.matmul( Omega_inv, Ki_test)))))
    Ko_sample = np.matmul(B_inv ,np.matmul(U_inv , np.matmul( Phi_inv , np.matmul(Chi_inv, np.matmul( Omega_inv, Ko_test)))))
310
    return Q_sample, Ki_sample, Ko_sample
311
    
christoph's avatar
christoph committed
312
def find_diag_angles(q, x0, U, B, Lab, beam_in, lambdai, lambdao, tol=1e-8, method='BFGS'):
313
314
315
316
    """ **find_diag_angles**
    Finds the FOURC spectrometer and sample angles for a desired q.

    Args:
alex's avatar
alex committed
317
318
319
320
321
322
323
324
      * q (array): Desired momentum transfer in Lab coordinates.
      * x0 (list): Guesses for the angles (tthv, tthh, chi, phi, omega).
      * U (array): 3x3 U-matrix Lab-to-sample transformation.
      * B (array): 3x3 B-matrix reciprocal lattice to absolute units transformation.
      * lambdai (float): Incident x-ray wavelength in Angstrom.
      * lambdao (float): Scattered x-ray wavelength in Angstrom.
      * tol (float): Toleranz for minimization (see scipy.optimize.minimize)
      * method (str): Method for minimization (see scipy.optimize.minimize)
325
326

    Returns:
alex's avatar
alex committed
327
       * ans (array): tthv, tthh, phi, chi, omega
328
329
    """
    # put UB matrix and energies into keyword argument for minimization
330
    kwargs = {'U': U, 'B': B, 'Lab': Lab, 'beam_in':beam_in, 'lambdai': lambdai, 'lambdao': lambdao}
331
332
333
    # least square minimization between wanted and guessed q
    fitfctn = lambda x: np.sum(( q - get_UB_Q(x[0], x[1], x[2], x[3], x[4], \
                                                  **kwargs)[0]  )**2)
334
    ans=minimize(fitfctn, x0, bounds=((0.,0.),(-10.,110.),(-7.,7.),(-7.,7.),(None,None)), tol=tol, method=method)
335
336
337
338
339
340
341
342
    print( ans )
    return ans.x

def get_gnuplot_rgb( start=None, end=None, length=None ):
    """ **get_gnuplot_rgb**
    Prints out a progression of RGB hex-keys to use in Gnuplot.

    Args:
alex's avatar
alex committed
343
344
345
      * start (array): RGB code to start from (must be numbers out of [0,1]).
      * end   (array): RGB code to end at (must be numbers out of [0,1]).
      * length  (int): How many colors to print out.
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
    """
    if start==None and end==None and length==None:
        rgb =  [[0,      0,           1], [0,      0.2353, 0.7647], \
                [0,      0.4706, 0.5294], [0,      0.7059, 0.2941], \
                [0,      0.9412, 0.0588], [0.1765,  0.8235, 0    ], \
                [0.4118,  0.5882, 0    ], [0.6471,  0.3529, 0    ], \
                [1,      0,      0     ] ]
    else:
        rgb = np.zeros((length,3))
        rgb[:,0] = np.linspace(start[0], end[0], length)
        rgb[:,1] = np.linspace(start[1], end[1], length)
        rgb[:,2] = np.linspace(start[2], end[2], length)
    #
    for ii in range(len(rgb)):
        print( 'set style line %d lt -1 lc rgb \'#%02x%02x%02x\' lw 1.5'%(ii+1, rgb[ii][0]*255., rgb[ii][1]*255., rgb[ii][2]*255.) )

christoph's avatar
christoph committed
362
363
364
def hex2rgb( hex_val ):
    return tuple( int(hex_val.lstrip('#')[i:i+2], 16) for i in (0, 2, 4) )

christoph's avatar
christoph committed
365
class maxipix_det:
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
    """
    Class to store some useful values from the detectors used. To be used for arranging the ROIs.
    """
    def __init__(self,name,spot_arrangement):
        self.name = name
        assert spot_arrangement in ['3x4','vertical'], 'unknown ROI arrangement, select \'3x4\' or \'vertical\'.'
        self.spot_arrangement = spot_arrangement

        self.roi_indices   = []
        self.roi_x_indices = []
        self.roi_y_indices = []
        self.roi_x_means   = []
        self.roi_y_means   = []
        self.pixel_size    = [256,256]
        self.PIXEL_RANGE   = {'VD': [0,256,0,256],  'VU': [0,256,256,512],  'VB': [0,256,512,768],
                            'HR': [256,512,0,256],'HL': [256,512,256,512],'HB': [256,512,512,768]}
    def get_pixel_range(self):
        return self.PIXEL_RANGE[self.name]

    def get_det_name(self):
        return self.name
christoph's avatar
christoph committed
387

388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490

class bragg_refl:
    """
    Dynamical theory of diffraction.
    """
    def __init__(self, crystal, hkl, alpha=0.0 ):

        # constants
        self.hc = constants.h*constants.c
        self.C  = 1.0
        self.r0 = constants.e**2/ \
          (4.0*np.pi*constants.epsilon_0* \
            constants.m_e*constants.c**2)* \
            1.0e10 # classical electron radius expressed in Angstrom
        self.P = 1 # polarization factor
        self.alpha = alpha
            
        # params
        self.hkl     = hkl
        self.crystal = crystal
        self.dspace  = dspace( self.hkl, self.crystal )
        self.ff_energy, self.f1_energy, self.f2_energy = \
          self.get_nff()

    def get_reflectivity_bent(self, energy, delta_theta, R):
        #
        refl,e,dev,e0 = taupgen( energy, self.hkl, crystals=self.crystal, R=R,
                    dev=delta_theta, alpha=self.alpha )
        return refl, theta_B*180/np.pi
        
    def get_reflectivity(self, energy, delta_theta, case='sigma'):
        energy = energy*1e3
        wavelength = self.hc/(energy*constants.e)*1e10
        print(wavelength)
        theta_B = np.arcsin(wavelength/(2.0*self.dspace)) # Bragg angle
        theta_inc = theta_B + np.radians(self.alpha) # incidence angle
        theta_ref = theta_B - np.radians(self.alpha) # reflection angle
        b = -np.sin(theta_ref)/np.sin(theta_inc) # asymmetry factor
        # polarization factor
        P = self.get_polarization_factor(2.*theta_B, case=case)
        # chi
        chi_0, chi_h, chi_hbar = self.get_chi(energy, self.crystal, self.hkl)
        #    index of refraction
        n = 1 + chi_0/2
        n_delta = 1 - np.real(n)
        n_beta = -np.imag(n)
        # linear absorption coefficient
        mu = -2*np.pi*np.imag(chi_0)/wavelength 
        # Bragg angle correction
        theta_B_correction = -chi_0*(1 - b)/(2*np.sin(2*theta_B)) 
        # width of the total reflection domain
        delta = np.abs(P)*np.sqrt(np.abs(b)*chi_h*chi_hbar)/np.sin(2*theta_B) 
        Darwin_width = 2*np.real(delta)
        lambda_B = wavelength*np.abs(np.sin(theta_inc))/ \
          (2*np.pi*np.real(delta)*np.sin(2*theta_B))
        delta_theta = delta_theta/1e6 # delta_theta is input in microrad
        eta = (delta_theta - theta_B_correction)/delta # deviation parameter
        # reflectivity curve
        reflectivity_curve = np.abs(eta - np.sign(np.real(eta))* \
                                        np.sqrt(eta**2 - 1))**2
        return reflectivity_curve, theta_B*180/np.pi, theta_B_correction
        
    def get_nff(self,
        nff_path = os.path.join(data_installation_dir,'atomic_form_factors')):
        fname    = os.path.join(nff_path, self.crystal.lower()+'.nff')
        table = np.loadtxt(fname, unpack = True, skiprows = 1)
        table = np.transpose(table)
        ff_energy = table[:,0]
        f1_energy = table[:,1]
        f2_energy = table[:,2]
        return ff_energy, f1_energy, f2_energy

    def get_chi(self, energy, crystal=None, hkl=None):
        path       =  os.path.join(data_installation_dir,'chitable_')
        hkl_string = str(int(hkl[0])) + str(int(hkl[1])) + str(int(hkl[2]))
        filestring = path + crystal.lower() + hkl_string + '.dat'
        self.chi        = np.loadtxt(filestring)
        self.chi_0 = complex(np.interp(energy, self.chi[:,0], self.chi[:,1]), \
                    np.interp(energy, self.chi[:,0],self.chi[:,2]))
        self.chi_h = complex(np.interp(energy, self.chi[:,0], self.chi[:,3]), \
                    np.interp(energy, self.chi[:,0], self.chi[:,4]))
        self.chi_hbar = np.conj(self.chi_h)
        return self.chi_0, self.chi_h, self.chi_hbar

    def get_polarization_factor(self, tth, case='sigma'):
        """
        Calculate polarization factor.
        """
        if case == 'sigma':
            P = 1.0
            self.P = P
            return P
        elif case == 'pi':
            P = np.cos(tth)**2
            self.P = P
            return P
        elif case == None:
            P = (1 + np.cos(tth)**2)/2.0
            self.P = P
            return P

        
        
christoph's avatar
christoph committed
491
492
493
class dtxrd:
    """ class to hold all things dynamic theory of diffraction.
    """
494
495
496
497
498
499
500
    def __init__(self, hkl, energy, crystal='Si', asym_angle=0.0, angular_range=[-0.5e-3, 0.5e-3] , angular_step=1e-8 ):

        # constants
        self.hc = 12.3984191
        self.C  = 1.0

        # params
christoph's avatar
christoph committed
501
        self.hkl     = hkl
502
503
504
505
        self.energy  = energy
        self.lam     = self.hc/self.energy
        self.crystal = crystal
        self.dspace  = dspace( self.hkl, self.crystal )
christoph's avatar
christoph committed
506
507
508
509
510
511

        # load Chi from tables:
        path       =  os.path.join(data_installation_dir,'chitable_')
        hkl_string = str(int(hkl[0])) + str(int(hkl[1])) + str(int(hkl[2]))
        filestring = path + crystal.lower() + hkl_string + '.dat'
        self.chi        = np.loadtxt(filestring)
512
513
514
515
516
517
518
519
520
521
522
523
524
        self.chi0 = complex(np.interp(self.energy, self.chi[:,0], self.chi[:,1]), \
                    np.interp(self.energy, self.chi[:,0],self.chi[:,2]))
        self.chih = complex(np.interp(self.energy, self.chi[:,0], self.chi[:,3]), \
                    np.interp(self.energy, self.chi[:,0], self.chi[:,4]))
        self.chihbar = np.conj(self.chih)

        # set Bragg angle:
        self.thetab  = float( bragg( hkl, energy , xtal=crystal ) )
        self.thetabd = float( braggd( hkl, energy , xtal=crystal ) )

        # set asymmetry:
        self.set_asymmetry(asym_angle)

525
526
527
        # set reduced deviation parameter
        self.angular_range = angular_range
        self.angular_step  = angular_step
528
529
        self.get_eta(angular_range, angular_step)

christoph's avatar
christoph committed
530
531
532
533
534
        self.lam_ext = None
        self.mu0     = None
        self.mus     = None
        self.R       = None
        self.omega_h = None
535
        self.omega_0 = None
christoph's avatar
christoph committed
536
537

    def set_asymmetry(self, alpha):
538
539
540
        """
        negative alpha -> more grazing incidence
        """
541
542
        self.alpha  = alpha
        self.gammah = -np.sin(self.thetab + alpha) # Krisch et al. convention 
christoph's avatar
christoph committed
543
        self.gamma0 =  np.sin(self.thetab - alpha) # Krisch et al. convention 
544
545
        self.gamma  =  self.gamma0/np.abs(self.gammah)                             
        self.beta   =  self.gammah/self.gamma0
christoph's avatar
christoph committed
546
547
548
549

    def set_hkl(self, hkl):
        self.hkl = hkl

550
551
    def set_energy(self, energy):
        self.energy = energy
christoph's avatar
christoph committed
552

553
554
555
556
557
558
    def get_reflectivity(self, angular_range=None, angular_step=None):
        if not angular_range:
            angular_range = self.angular_range
        if not angular_step:
            angular_step = self.angular_step

559
560
561
562
563
564
565
566
        self.get_eta(angular_range, angular_step)
        
        pre_factor = (np.sqrt(self.chih*self.chihbar))/(self.chihbar) * \
          1.0/np.sqrt(np.abs(self.gamma)) * self.gammah/np.abs(self.gammah) * self.C/np.abs(self.C)
          
        self.R = np.abs(pre_factor * np.piecewise(self.eta, self.eta.real>0, \
                                          [lambda x: x - np.sqrt( x**2 + self.gammah/np.abs(self.gammah)  )  , \
                                            lambda x: x + np.sqrt( x**2 + self.gammah/np.abs(self.gammah)  ) ]))
christoph's avatar
christoph committed
567
            
568
569
570
571
572
    def get_eta(self, angular_range, angular_step=1e-8):
        self.theta = np.arange( self.thetab+angular_range[0], self.thetab+angular_range[1], angular_step )
        self.eta = ( (self.theta-self.thetab)*np.sin(2.*self.thetab) + self.chi0/2.0*(1-self.gamma)  ) / \
          (np.abs(self.C)*np.sqrt(np.abs(self.gamma))* np.sqrt(self.chih*self.chihbar))

573
574
575
576
    def get_anomalous_absorption(self, energy=None):
        if not energy:
            energy = self.energy

577
        # photoelectric absorption
578
579
580
581
582
583
584
585
        mu0 = myprho( energy, self.crystal )[0][0][0]*0.602252/ \
          myprho( energy, self.crystal )[2]

        omega = np.arctan( (1.0 -  np.abs(self.R)**2) / (1.0 + np.abs(self.R)**2) * np.tan(self.thetab) )

        mus = mu0 * np.cos(omega)/np.cos(self.thetab) * ( 1.0 + 2.0*np.abs(self.C) * \
                    self.chih.imag/self.chi0.imag * self.R.real/(1.0 + np.abs(self.R)**2)  )
        self.mus = mus
586

587
    def get_extinction_length(self, energy=None):
588

589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
        if energy:
            lam = energy/self.hc
        else:
            lam = self.lam

        kappa = -np.abs( (self.chih.imag)/(self.chih.real) )
        
        self.lam_ext = (lam*np.sqrt(self.gamma0 * np.abs(self.gammah))* np.cos(kappa) ) / \
          (2.0*np.pi*np.abs(self.C) * np.sqrt( np.abs(self.chih.real*self.chihbar.real)  ))

    def get_reflection_width(self):
        if not self.lam_ext:
            self.get_extinction_length()
        
        self.omega_h = self.lam*np.abs(self.gammah)/( np.pi*self.lam_ext*np.sin(2.0*self.thetab) )
        self.omega_0 = self.lam*self.gamma0/( np.pi*self.lam_ext*np.sin(2.0*self.thetab) )
605

christoph's avatar
christoph committed
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
    
def dtxrd_reflectivity( energy, hkl, alpha=0.0, crystal='Si', angular_range=np.arange(-0.5e-3, 0.5e-3) ):

    # constants
    hc, lam, dsp, thetab, alpha, C = \
      get_dtxrd_constants( energy, hkl, crystal=crystal, alpha=alpha )   

    # theta scale
    theta  = np.arange( thetab+angular_range[0], \
                        thetab+angular_range[1], 1e-8 )

    # load chi
    chi0, chih, chihbar = get_dtxrd_chi( energy, hkl, crystal=crystal )

    # asymmetry parameter
    gamma, gamma0, gammah, beta = get_dtxrd_assymmetry_params(thetab, alpha)

    # calculate eta
    eta = get_dtxrd_eta(theta, thetab, chi0, chih, chihbar)

    # calculate reflectivity
    pre_factor = (np.sqrt(chih*chihbar))/(chihbar) * \
        1.0/np.sqrt(np.abs(gamma)) * gammah/np.abs(gammah) * C/np.abs(C)
    r = pre_factor * np.piecewise(eta, eta.real>0, \
        [lambda x: x - np.sqrt( x**2 + gammah/np.abs(gammah)  )  , \
         lambda x: x + np.sqrt( x**2 + gammah/np.abs(gammah)  ) ])

    return theta, r

def dtxrd_anomalous_absorption( energy, hkl, alpha=0.0, crystal='Si', angular_range=np.arange(-0.5e-3, 0.5e-3) ):

    # constants
    hc     = 12.3984191
    lam    = hc/energy * 1e-10
    dsp    = dspace(hkl, crystal)
    thetab = float( bragg( hkl, energy , xtal=crystal ) )
    alpha  = np.radians(alpha)
    C      = 1.0

    # theta scale
    theta  = np.arange( thetab+angular_range[0], \
                        thetab+angular_range[1], 1e-8 )

    # load chi
    path       = os.path.join(data_installation_dir,'chitable_')
    hkl_string = str(int(hkl[0])) + str(int(hkl[1])) + str(int(hkl[2]))
    filestring = path + crystal.lower() + hkl_string + '.dat'
    chi        = np.loadtxt(filestring)
    chi0 = complex(np.interp(energy,chi[:,0],chi[:,1]),
                   np.interp(energy,chi[:,0],chi[:,2]))
    chih = complex(np.interp(energy,chi[:,0],chi[:,3]),
                   np.interp(energy,chi[:,0],chi[:,4]))
    chihbar = np.conj(chih)
    
    # asymmetry parameter
    gammah = -np.sin(thetab + alpha) # Krisch et al. convention
    gamma0 = np.sin(thetab - alpha) # Krisch et al. convention
    beta   = gamma0/np.abs(gammah)
    gamma  = gammah/gamma0

    # calculate eta
    eta = ( (theta-thetab)*np.sin(2.*thetab) + chi0/2.0*(1-gamma)  ) / \
    (np.abs(C)*np.sqrt(np.abs(gamma))* np.sqrt(chih*chihbar))

    # calculate reflectivity
    pre_factor = (np.sqrt(chih*chihbar))/(chihbar) * \
        1.0/np.sqrt(np.abs(gamma)) * gammah/np.abs(gammah) * C/np.abs(C)
    r = pre_factor * np.piecewise(eta, eta.real>0, \
        [lambda x: x - np.sqrt( x**2 + gammah/np.abs(gammah)  )  , \
         lambda x: x + np.sqrt( x**2 + gammah/np.abs(gammah)  ) ])

    # photoelectric absorption
    mu0 = myprho( energy, crystal )[0][0][0]*0.602252/ \
        myprho( energy, crystal )[2]
    mus = mu0 * np.cos(omega)/np.cos(thetab) * ( 1.0 + 2.0*np.abs(C) * chih.imag/chi0.imag * r.real/(1.0 + np.abs(r)**2)  )

    return theta, mus

def dtxrd_extinction_length( energy, hkl, alpha=0.0, crystal='Si' ):
    pass

687
def delE_dicedAnalyzerIntrinsic(E, Dw, Theta):
688
689
    """Calculates the intrinsic energy resolution of a diced crystal
    analyzer.
690

691
692
693
694
    Args:
        E     (float): Working energy in [eV].
        Dw    (float): Darwin width of the used reflection [microRad].
        Theta (float): Analyzer Bragg angle [degree].
695

696
697
698
699
700
    Returns:
        Intrinsic energy resolution of a perfect analyzer crystal.
    """
    Dw = Dw/1000000.0 # conversion to radians
    return E * (Dw)/(np.tan(np.radians(Theta)))
701
702

def delE_JohannAberration(E, A, R, Theta):
703
    """Calculates the Johann aberration of a spherical analyzer crystal.
704

705
706
707
708
709
    Args:
        E     (float): Working energy in [eV].
        A     (float): Analyzer aperture [mm].
        R     (float): Radius of the Rowland circle [mm].
        Theta (float): Analyzer Bragg angle [degree].
710

711
712
713
714
    Returns:
        Johann abberation in [eV].
    """
    return E/2.0 * ((A)/(2.0*R*np.tan(np.radians(Theta))))**2
715
716

def delE_pixelSize(E, p, R, Theta):
717
718
    """Calculates the pixel size contribution to the resolution function
    of a diced analyzer crystal.
719

720
721
722
723
724
    Args:
        E     (float): Working energy in [eV].
        p     (float): Pixel size in [mm].
        R     (float): Radius of the Rowland circle [mm].
        Theta (float): Analyzer Bragg angle [degree].
725

726
727
728
729
730
    Returns:
        Pixel size contribution in [eV] to the energy resolution for a diced analyzer
        crystal.
    """
    return E * (p/(2.0*R*np.sin(np.radians(Theta))))/(np.tan(np.radians(Theta)))
731
732

def delE_sourceSize(E, s, R, Theta):
733
    """Calculates the source size contribution to the resolution function.
734

735
736
737
738
739
    Args:
        E     (float): Working energy in [eV].
        s     (float): Source size in [mm].
        R     (float): Radius of the Rowland circle [mm].
        Theta (float): Analyzer Bragg angle [degree].
740

741
742
743
744
    Returns:
        Source size contribution in [eV] to the energy resolution.
    """
    return E * (s/(R*np.sin(np.radians(Theta))))/(np.tan(np.radians(Theta)))
745
746

def delE_offRowland(E, z, A, R, Theta):
747
    """Calculates the off-Rowland contribution of a spherical analyzer crystal.
748

749
750
751
752
753
754
    Args:
        E     (float): Working energy in [eV].
        z     (float): Off-Rowland distance [mm].
        A     (float): Analyzer aperture [mm].
        R     (float): Radius of the Rowland circle [mm].
        Theta (float): Analyzer Bragg angle [degree].
755

756
757
    Returns:
        Off-Rowland contribution in [eV] to the energy resolution.
758

759
760
    """
    return E * (z*A)/( (R*np.sin(np.radians(Theta)) + z)**2 ) * (1.0)/(np.tan(np.radians(Theta)))
761
762

def delE_stressedCrystal(E, t, v, R, Theta):
763
764
    """Calculates the stress induced contribution to the resulution function
    of a spherically bent crystal analyzer.
765

766
767
768
769
770
771
    Args:
        E     (float): Working energy in [eV].
        t     (float): Absorption length in the analyzer material [mm].
        v     (float): Poisson ratio of the analyzer material.
        R     (float): Radius of the Rowland circle [mm].
        Theta (float): Analyzer Bragg angle [degree].
772

773
774
    Returns:
        Stress-induced contribution in [eV] to the energy resolution.
775

776
777
    """
    return E * t/R * np.absolute((1.0)/(np.tan(np.radians(Theta))**2) - 2.0*v)
778
779

def get_num_of_MD_steps(time_ps,time_step):
780
781
    """Calculates the number of steps in an MD simulation for a desired 
    time (in ps) and given step size (in a.u.)
782

783
784
785
    Args:
        time_ps   (float): Desired time span (ps).
        time_step (float): Chosen time step (a.u.).
786

787
788
789
790
    Returns:
        The number of steps required to span the desired time span.
    """
    return time_ps / time_step /1.0e12 / 2.418884326505e-17
791

christoph's avatar
christoph committed
792
def nonzeroavg(y=None):
793
794
795
796
797
798
799
800
801
802
803
804
    dim1 = y.shape[0]
    yavg = np.zeros(dim1,float)
    for ii in range(dim1): 
        length = 0
        rowsum = 0.
        for ij in range(y.shape[1]):        
            if (not np.isnan(y[ii, ij])):                
                length += 1
                rowsum += y[ii, ij]            
        yavg[ii] = rowsum / float(length)
    yavg = yavg * y.shape[1]
    return(yavg)
christoph's avatar
christoph committed
805

christoph's avatar
christoph committed
806
def fermi(rs):
807
808
809
810
811
    """ **fermi**
    Calculates the plasmon energy (in eV), Fermi energy (in eV), Fermi 
    momentum (in a.u.), and critical plasmon cut-off vector (in a.u.).

    Args:
alex's avatar
alex committed
812
     * rs (float): electron separation parameter
813
814

    Returns:
alex's avatar
alex committed
815
816
817
818
       * wp (float): plasmon energy (in eV)
       * ef (float): Fermi energy (in eV)
       * kf (float): Fermi momentum (in a.u.)
       * kc (float): critical plasmon cut-off vector (in a.u.)
819
820
821
822
823
824
825
826
827
828
829

    Based on Matlab function from A. Soininen.
    """
    au   = 27.212
    alfa = (9.0*np.pi/4.0)**(1.0/3.0)
    kf = alfa/rs
    ef = kf*kf/2.0
    wp = np.sqrt(3.0/rs/rs/rs)
    kc = kf * (np.sqrt(1.0+wp/ef)-1.0)
    wp = wp*au
    ef = ef*au
830
    return wp, ef, kf, kcem
christoph's avatar
christoph committed
831
832

def lindhard_pol(q,w,rs=3.93,use_corr=False, lifetime=0.28):
833
834
835
836
837
    """ **lindhard_pol**
    Calculates the Lindhard polarizability function (RPA) for 
    certain q (a.u.), w (a.u.) and rs (a.u.).

    Args:
alex's avatar
alex committed
838
839
840
841
842
      * q (float): momentum transfer (in a.u.)
      * w (float): energy (in a.u.)
      * rs (float): electron parameter
      * use_corr (boolean): if True, uses Bernardo's calculation for n(k) instead of the Fermi function.
      * lifetime (float): life time (default is 0.28 eV for Na).
843
844
845

    Based on Matlab function by S. Huotari.
    """
846
847
    if type(w) in [float, int]:
        w = np.array([w])
christoph's avatar
christoph committed
848
    wp, ef, kf, kc = fermi(rs)
849
850
851
852
    ef     = ef/27.212
    gammal = lifetime/27.212  # lifetime  (0.28 eV for Na)
    th     = np.arange( 0.0, np.pi, np.pi/700.0 )
    k      = np.arange( 0.0, 2.0*kf, kf/1000.0 )
853
    [K,TH] = np.meshgrid(k,th)
854
855
856

    ek     = K**2/2.0
    ekq    = ( K**2+q**2+2*q*K*np.cos(TH) )/2.0
857
858
859
860
861
862
863
    if not use_corr:
        fek = np.zeros(np.shape(ek)) 
        fek[ek<=ef]=1.0
        fekq=np.zeros(np.shape(ekq)) 
        fekq[ekq<=ef] = 1.0
    if use_corr:
        print('Not implemented yet!')
864
    x = np.zeros_like(w, dtype='complex')
865
    for ii in range(len(w)):
866
867
868
        y=np.sin(TH)*(fek-fekq)/(w[ii]+ek-ekq+np.complex(0,1)*gammal)
        y=np.trapz( y, th, axis=0 )
        y=np.trapz( k**2.0*y, k, axis=0 )
869
870
871
872
        x[ii]=y
    x = 4.0*np.pi*x
    x = x/(2.0*np.pi)**3
    return x 
christoph's avatar
christoph committed
873

christoph's avatar
christoph committed
874
def energy(d,ba):
875
876
877
878
879
880
881
882
883
884
885
    """
    % ENERGY  Calculates energy corrresponing to Bragg angle for given d-spacing
    %         function e=energy(dspace,bragg_angle)
    %
    %      dspace for reflection
    %      bragg_angle in DEG
    %
    %         KH 28.09.93
    """
    hc = 12.3984191 # CODATA 2002 physics.nist.gov/constants
    return (2.0*d*np.sin(ba/180.0*np.pi)/hc)**(-1)
christoph's avatar
christoph committed
886
887

def dspace(hkl=[6,6,0],xtal='Si'):
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
    """
    % DSPACE Gives d-spacing for given xtal
    %     d=dspace(hkl,xtal)
    %     hkl can be a matrix i.e. hkl=[1,0,0 ; 1,1,1];
    %     xtal='Si','Ge','LiF','InSb','C','Dia','Li' (case insensitive)
    %     if xtal is number this is user as a d0
    %
    %     KH 28.09.93 
    %        SH 2005
    %
    """
    # create a database of lattice constants (could be a shelf)
    xtable = {}
    xtable['SI'] = 5.43102088
    xtable['GE'] = 5.657
    xtable['SIXOP'] = 5.430919
    xtable['SIKOH'] = 5.430707
    xtable['LIF'] = 4.027
    xtable['INSB'] = 6.4784
    xtable['C'] = 6.708
    xtable['DIA'] = 3.57
    xtable['LI'] = 3.41
    xtable['TCNE'] = 9.736
    xtable['CU'] = 3.61
    xtable['PB'] = 4.95
    xtable['NA'] = 4.2906
    xtable['AL'] = 4.0495

    if isinstance(xtal,str):
        try:
            a0 = xtable[xtal.upper()]
        except KeyError:
            print( 'Lattice constant is not in database')
            return
    else: 
        a0 = xtal # if number is provided, it's taken as lattice constant

    return a0/np.sqrt(np.sum(np.array(hkl)**2.0))
christoph's avatar
christoph committed
926
927

def bragg(hkl,e,xtal='Si'):
928
929
930
931
932
933
934
935
936
937
938
939
    """
    % BRAGG  Calculates Bragg angle for given reflection in RAD
    %      output=bangle(hkl,e,xtal)
    %        hkl can be a matrix i.e. hkl=[1,0,0 ; 1,1,1];
    %      e=energy in keV
    %      xtal='Si', 'Ge', etc. (check dspace.m) or d0 (Si default)
    %
    %      KH 28.09.93
    %
    """
    hc = 12.3984191 # CODATA 2002 recommended value, physics.nist.gov/constants
    return np.real(np.arcsin((2.0*dspace(hkl,xtal)*e/hc)**(-1.0)))
christoph's avatar
christoph committed
940
941

def braggd(hkl,e,xtal='Si'):
942
943
944
945
946
947
948
949
950
951
952
    """
    # BRAGGD  Calculates Bragg angle for given reflection in deg
    #      Call BRAGG.M
    #      output=bangle(hkl,e,xtal)
    #        hkl can be a matrix i.e. hkl=[1,0,0 ; 1,1,1];
    #      e=energy in keV
    #      xtal='Si', 'Ge', etc. (check dspace.m) or d0 (Si default)
    #
    #      KH 28.09.93
    """
    return bragg(hkl,e,xtal)/np.pi*180.0
christoph's avatar
christoph committed
953
954

def addch(xold,yold,n,n0=0,errors=None):
955
956
957
958
959
960
961
962
963
964
965
966
967
968
    """
    # ADDCH     Adds contents of given adjacent channels together
    #
    #           [x2,y2] = addch(x,y,n,n0)
    #           x  = original x-scale  (row or column vector)
    #           y  = original y-values (row or column vector)
    #           n  = number of channels to be summed up
    #            n0 = offset for adding, default is 0
    #           x2 = new x-scale 
    #           y2 = new y-values
    #
    #           KH 17.09.1990
    #        Modified 29.05.1995 to include offset
    """
969
    n0=int(n0-np.fix(n0/n)*n)
970
971
972
973
    if n0<0:
         n0 = (n + n0)
    datalen = np.floor( (len(xold) - n0) / n)

974
975
976
977
    xnew = np.zeros(int(np.min([datalen,len(xold)])))
    ynew = np.zeros(int(np.min([datalen,len(xold)])))
    errnew = np.zeros(int(np.min([datalen,len(xold)])))

978
979
980
981
982
983
984
    for i in range(int(datalen)):
        xnew[i] = np.sum(xold[i*n+n0:i*n+n+n0])/n
        ynew[i] = np.sum(yold[i*n+n0:i*n+n+n0])/n
        if np.any(errors):
            errnew[i] = np.sqrt(np.sum(errors[i*n+n0:i*n+n+n0]**2.0))
            return xnew, ynew, errnew
    return xnew, ynew
985
986

def fwhm(x,y):
987
988
989
990
991
992
993
994
995
996
997
    """
    finds full width at half maximum of the curve y vs. x
    returns 
    f  = FWHM
    x0 = position of the maximum
    """
    if x[-1] < x[0]:
        x = np.flipud(x)
        y = np.flipud(y)

    y0 = np.amax(y)
998
999
1000
    i0 = np.where(y == y0)[0]
    
    if len(i0)>1: