| 1 | # -*- coding: utf-8 -*- |
|---|
| 2 | """ |
|---|
| 3 | *GSASIIElem: functions for element types* |
|---|
| 4 | ----------------------------------------- |
|---|
| 5 | |
|---|
| 6 | """ |
|---|
| 7 | # Copyright: 2008, Robert B. Von Dreele & Brian H. Toby (Argonne National Laboratory) |
|---|
| 8 | ########### SVN repository information ################### |
|---|
| 9 | # $Date: 2018-09-17 15:03:54 +0000 (Mon, 17 Sep 2018) $ |
|---|
| 10 | # $Author: vondreele $ |
|---|
| 11 | # $Revision: 3605 $ |
|---|
| 12 | # $URL: trunk/GSASIIElem.py $ |
|---|
| 13 | # $Id: GSASIIElem.py 3605 2018-09-17 15:03:54Z vondreele $ |
|---|
| 14 | ########### SVN repository information ################### |
|---|
| 15 | |
|---|
| 16 | import math |
|---|
| 17 | import sys |
|---|
| 18 | import os.path |
|---|
| 19 | import GSASIIpath |
|---|
| 20 | GSASIIpath.SetVersionNumber("$Revision: 3605 $") |
|---|
| 21 | import numpy as np |
|---|
| 22 | import atmdata |
|---|
| 23 | import GSASIImath as G2mth |
|---|
| 24 | |
|---|
| 25 | getElSym = lambda sym: sym.split('+')[0].split('-')[0].capitalize() |
|---|
| 26 | def GetFormFactorCoeff(El): |
|---|
| 27 | """Read X-ray form factor coefficients from `atomdata.py` file |
|---|
| 28 | |
|---|
| 29 | :param str El: element 1-2 character symbol, case irrevelant |
|---|
| 30 | :return: `FormFactors`: list of form factor dictionaries |
|---|
| 31 | |
|---|
| 32 | Each X-ray form factor dictionary is: |
|---|
| 33 | |
|---|
| 34 | * `Symbol`: 4 character element symbol with valence (e.g. 'NI+2') |
|---|
| 35 | * `Z`: atomic number |
|---|
| 36 | * `fa`: 4 A coefficients |
|---|
| 37 | * `fb`: 4 B coefficients |
|---|
| 38 | * `fc`: C coefficient |
|---|
| 39 | |
|---|
| 40 | """ |
|---|
| 41 | |
|---|
| 42 | Els = El.capitalize().strip() |
|---|
| 43 | valences = [ky for ky in atmdata.XrayFF.keys() if Els == getElSym(ky)] |
|---|
| 44 | FormFactors = [atmdata.XrayFF[val] for val in valences] |
|---|
| 45 | for Sy,FF in zip(valences,FormFactors): |
|---|
| 46 | FF.update({'Symbol':Sy.upper()}) |
|---|
| 47 | return FormFactors |
|---|
| 48 | |
|---|
| 49 | def GetFFtable(atomTypes): |
|---|
| 50 | ''' returns a dictionary of form factor data for atom types found in atomTypes |
|---|
| 51 | |
|---|
| 52 | :param list atomTypes: list of atom types |
|---|
| 53 | :return: FFtable, dictionary of form factor data; key is atom type |
|---|
| 54 | |
|---|
| 55 | ''' |
|---|
| 56 | FFtable = {} |
|---|
| 57 | for El in atomTypes: |
|---|
| 58 | FFs = GetFormFactorCoeff(getElSym(El)) |
|---|
| 59 | for item in FFs: |
|---|
| 60 | if item['Symbol'] == El.upper(): |
|---|
| 61 | FFtable[El] = item |
|---|
| 62 | return FFtable |
|---|
| 63 | |
|---|
| 64 | def GetMFtable(atomTypes,Landeg): |
|---|
| 65 | ''' returns a dictionary of magnetic form factor data for atom types found in atomTypes |
|---|
| 66 | |
|---|
| 67 | :param list atomTypes: list of atom types |
|---|
| 68 | :param list Landeg: Lande g factors for atomTypes |
|---|
| 69 | :return: FFtable, dictionary of form factor data; key is atom type |
|---|
| 70 | |
|---|
| 71 | ''' |
|---|
| 72 | MFtable = {} |
|---|
| 73 | for El,gfac in zip(atomTypes,Landeg): |
|---|
| 74 | MFs = GetMagFormFacCoeff(getElSym(El)) |
|---|
| 75 | for item in MFs: |
|---|
| 76 | if item['Symbol'] == El.upper(): |
|---|
| 77 | item['gfac'] = gfac |
|---|
| 78 | MFtable[El] = item |
|---|
| 79 | return MFtable |
|---|
| 80 | |
|---|
| 81 | def GetBLtable(General): |
|---|
| 82 | ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General |
|---|
| 83 | |
|---|
| 84 | :param dict General: dictionary of phase info.; includes AtomTypes & Isotopes |
|---|
| 85 | :return: BLtable, dictionary of scattering length data; key is atom type |
|---|
| 86 | ''' |
|---|
| 87 | atomTypes = General['AtomTypes'] |
|---|
| 88 | BLtable = {} |
|---|
| 89 | isotope = General['Isotope'] |
|---|
| 90 | for El in atomTypes: |
|---|
| 91 | ElS = getElSym(El) |
|---|
| 92 | if 'Nat' in isotope[El]: |
|---|
| 93 | BLtable[El] = [isotope[El],atmdata.AtmBlens[ElS+'_']] |
|---|
| 94 | else: |
|---|
| 95 | BLtable[El] = [isotope[El],atmdata.AtmBlens[ElS+'_'+isotope[El]]] |
|---|
| 96 | return BLtable |
|---|
| 97 | |
|---|
| 98 | def getFFvalues(FFtables,SQ,ifList=False): |
|---|
| 99 | 'Needs a doc string' |
|---|
| 100 | if ifList: |
|---|
| 101 | FFvals = [] |
|---|
| 102 | for El in FFtables: |
|---|
| 103 | FFvals.append(ScatFac(FFtables[El],SQ)[0]) |
|---|
| 104 | else: |
|---|
| 105 | FFvals = {} |
|---|
| 106 | for El in FFtables: |
|---|
| 107 | FFvals[El] = ScatFac(FFtables[El],SQ)[0] |
|---|
| 108 | return FFvals |
|---|
| 109 | |
|---|
| 110 | def getBLvalues(BLtables,ifList=False): |
|---|
| 111 | 'Needs a doc string' |
|---|
| 112 | if ifList: |
|---|
| 113 | BLvals = [] |
|---|
| 114 | for El in BLtables: |
|---|
| 115 | if 'BW-LS' in El: |
|---|
| 116 | BLvals.append(BLtables[El][1]['BW-LS'][0]) |
|---|
| 117 | else: |
|---|
| 118 | BLvals.append(BLtables[El][1]['SL'][0]) |
|---|
| 119 | else: |
|---|
| 120 | BLvals = {} |
|---|
| 121 | for El in BLtables: |
|---|
| 122 | if 'BW-LS' in El: |
|---|
| 123 | BLvals[El] = BLtables[El][1]['BW-LS'][0] |
|---|
| 124 | else: |
|---|
| 125 | BLvals[El] = BLtables[El][1]['SL'][0] |
|---|
| 126 | return BLvals |
|---|
| 127 | |
|---|
| 128 | def getMFvalues(MFtables,SQ,ifList=False): |
|---|
| 129 | 'Needs a doc string' |
|---|
| 130 | if ifList: |
|---|
| 131 | MFvals = [] |
|---|
| 132 | for El in MFtables: |
|---|
| 133 | MFvals.append(MagScatFac(MFtables[El],SQ)[0]) |
|---|
| 134 | else: |
|---|
| 135 | MFvals = {} |
|---|
| 136 | for El in MFtables: |
|---|
| 137 | MFvals[El] = MagScatFac(MFtables[El],SQ)[0] |
|---|
| 138 | return MFvals |
|---|
| 139 | |
|---|
| 140 | def GetFFC5(ElSym): |
|---|
| 141 | '''Get 5 term form factor and Compton scattering data |
|---|
| 142 | |
|---|
| 143 | :param ElSym: str(1-2 character element symbol with proper case); |
|---|
| 144 | :return El: dictionary with 5 term form factor & compton coefficients |
|---|
| 145 | ''' |
|---|
| 146 | import FormFactors as FF |
|---|
| 147 | El = {} |
|---|
| 148 | FF5 = FF.FFac5term[ElSym] |
|---|
| 149 | El['fa'] = FF5[:5] |
|---|
| 150 | El['fc'] = FF5[5] |
|---|
| 151 | El['fb'] = FF5[6:] |
|---|
| 152 | Cmp5 = FF.Compton[ElSym] |
|---|
| 153 | El['cmpz'] = Cmp5[0] |
|---|
| 154 | El['cmpa'] = Cmp5[1:6] |
|---|
| 155 | El['cmpb'] = Cmp5[6:] |
|---|
| 156 | return El |
|---|
| 157 | |
|---|
| 158 | def CheckElement(El): |
|---|
| 159 | '''Check if element El is in the periodic table |
|---|
| 160 | |
|---|
| 161 | :param str El: One or two letter element symbol, capitaliztion ignored |
|---|
| 162 | :returns: True if the element is found |
|---|
| 163 | |
|---|
| 164 | ''' |
|---|
| 165 | import ElementTable as ET |
|---|
| 166 | Elements = [] |
|---|
| 167 | for elem in ET.ElTable: |
|---|
| 168 | Elements.append(elem[0][0]) |
|---|
| 169 | if El.capitalize() in Elements: |
|---|
| 170 | return True |
|---|
| 171 | else: |
|---|
| 172 | return False |
|---|
| 173 | |
|---|
| 174 | def FixValence(El): |
|---|
| 175 | 'Returns the element symbol, even when a valence is present' |
|---|
| 176 | if '+' in El[-1]: #converts An+/- to A+/-n |
|---|
| 177 | num = El[-2] |
|---|
| 178 | El = El.split(num)[0]+'+'+num |
|---|
| 179 | if '+0' in El: |
|---|
| 180 | El = El.split('+0')[0] |
|---|
| 181 | if '-' in El[-1]: |
|---|
| 182 | num = El[-2] |
|---|
| 183 | El = El.split(num)[0]+'-'+num |
|---|
| 184 | if '-0' in El: |
|---|
| 185 | El = El.split('-0')[0] |
|---|
| 186 | return El |
|---|
| 187 | |
|---|
| 188 | def GetAtomInfo(El,ifMag=False): |
|---|
| 189 | 'reads element information from atmdata.py' |
|---|
| 190 | import ElementTable as ET |
|---|
| 191 | Elem = ET.ElTable |
|---|
| 192 | if ifMag: |
|---|
| 193 | Elem = ET.MagElTable |
|---|
| 194 | Elements = [elem[0][0] for elem in Elem] |
|---|
| 195 | AtomInfo = {} |
|---|
| 196 | ElS = getElSym(El) |
|---|
| 197 | if El not in atmdata.XrayFF and El not in atmdata.MagFF: |
|---|
| 198 | if ElS not in atmdata.XrayFF: |
|---|
| 199 | print('Atom type '+El+' not found, using H') |
|---|
| 200 | ElS = 'H' |
|---|
| 201 | # return # not sure what this element should be! |
|---|
| 202 | print('Atom type '+El+' not found, using '+ElS) |
|---|
| 203 | El = ElS |
|---|
| 204 | AtomInfo.update(dict(zip(['Drad','Arad','Vdrad','Hbrad'],atmdata.AtmSize[ElS]))) |
|---|
| 205 | AtomInfo['Symbol'] = El |
|---|
| 206 | AtomInfo['Color'] = ET.ElTable[Elements.index(ElS)][6] |
|---|
| 207 | AtomInfo['Z'] = atmdata.XrayFF[ElS]['Z'] |
|---|
| 208 | isotopes = [ky for ky in atmdata.AtmBlens.keys() if ElS == ky.split('_')[0]] |
|---|
| 209 | isotopes.sort() |
|---|
| 210 | AtomInfo['Mass'] = atmdata.AtmBlens[isotopes[0]]['Mass'] #default to nat. abund. |
|---|
| 211 | AtomInfo['Isotopes'] = {} |
|---|
| 212 | for isotope in isotopes: |
|---|
| 213 | data = atmdata.AtmBlens[isotope] |
|---|
| 214 | if isotope == ElS+'_': |
|---|
| 215 | AtomInfo['Isotopes']['Nat. Abund.'] = data |
|---|
| 216 | else: |
|---|
| 217 | AtomInfo['Isotopes'][isotope.split('_')[1]] = data |
|---|
| 218 | AtomInfo['Lande g'] = 2.0 |
|---|
| 219 | return AtomInfo |
|---|
| 220 | |
|---|
| 221 | def GetElInfo(El,inst): |
|---|
| 222 | ElemSym = El.strip().capitalize() |
|---|
| 223 | if 'X' in inst['Type'][0]: |
|---|
| 224 | keV = 12.397639/G2mth.getWave(inst) |
|---|
| 225 | FpMu = FPcalc(GetXsectionCoeff(ElemSym), keV) |
|---|
| 226 | ElData = GetFormFactorCoeff(ElemSym)[0] |
|---|
| 227 | ElData['FormulaNo'] = 0.0 |
|---|
| 228 | ElData.update(GetAtomInfo(ElemSym)) |
|---|
| 229 | ElData.update(dict(zip(['fp','fpp','mu'],FpMu))) |
|---|
| 230 | ElData.update(GetFFC5(El)) |
|---|
| 231 | else: #'N'eutron |
|---|
| 232 | ElData = {} |
|---|
| 233 | ElData.update(GetAtomInfo(ElemSym)) |
|---|
| 234 | ElData['FormulaNo'] = 0.0 |
|---|
| 235 | ElData.update({'mu':0.0,'fp':0.0,'fpp':0.0}) |
|---|
| 236 | return ElData |
|---|
| 237 | |
|---|
| 238 | def GetXsectionCoeff(El): |
|---|
| 239 | """Read atom orbital scattering cross sections for fprime calculations via Cromer-Lieberman algorithm |
|---|
| 240 | |
|---|
| 241 | :param El: 2 character element symbol |
|---|
| 242 | :return: Orbs: list of orbitals each a dictionary with detailed orbital information used by FPcalc |
|---|
| 243 | |
|---|
| 244 | each dictionary is: |
|---|
| 245 | |
|---|
| 246 | * 'OrbName': Orbital name read from file |
|---|
| 247 | * 'IfBe' 0/2 depending on orbital |
|---|
| 248 | * 'BindEn': binding energy |
|---|
| 249 | * 'BB': BindEn/0.02721 |
|---|
| 250 | * 'XSectIP': 5 cross section inflection points |
|---|
| 251 | * 'ElEterm': energy correction term |
|---|
| 252 | * 'SEdge': absorption edge for orbital |
|---|
| 253 | * 'Nval': 10/11 depending on IfBe |
|---|
| 254 | * 'LEner': 10/11 values of log(energy) |
|---|
| 255 | * 'LXSect': 10/11 values of log(cross section) |
|---|
| 256 | |
|---|
| 257 | """ |
|---|
| 258 | AU = 2.80022e+7 |
|---|
| 259 | C1 = 0.02721 |
|---|
| 260 | ElS = El.upper() |
|---|
| 261 | ElS = ElS.ljust(2) |
|---|
| 262 | filename = os.path.join(os.path.split(__file__)[0],'Xsect.dat') |
|---|
| 263 | try: |
|---|
| 264 | xsec = open(filename,'Ur') |
|---|
| 265 | except: |
|---|
| 266 | print ('**** ERROR - File Xsect.dat not found in directory %s'%os.path.split(filename)[0]) |
|---|
| 267 | sys.exit() |
|---|
| 268 | S = '1' |
|---|
| 269 | Orbs = [] |
|---|
| 270 | while S: |
|---|
| 271 | S = xsec.readline() |
|---|
| 272 | if S[:2] == ElS: |
|---|
| 273 | S = S[:-1]+xsec.readline()[:-1]+xsec.readline() |
|---|
| 274 | OrbName = S[9:14] |
|---|
| 275 | S = S[14:] |
|---|
| 276 | IfBe = int(S[0]) |
|---|
| 277 | S = S[1:] |
|---|
| 278 | val = S.split() |
|---|
| 279 | BindEn = float(val[0]) |
|---|
| 280 | BB = BindEn/C1 |
|---|
| 281 | Orb = {'OrbName':OrbName,'IfBe':IfBe,'BindEn':BindEn,'BB':BB} |
|---|
| 282 | Energy = [] |
|---|
| 283 | XSect = [] |
|---|
| 284 | for i in range(11): |
|---|
| 285 | Energy.append(float(val[2*i+1])) |
|---|
| 286 | XSect.append(float(val[2*i+2])) |
|---|
| 287 | XSecIP = [] |
|---|
| 288 | for i in range(5): XSecIP.append(XSect[i+5]/AU) |
|---|
| 289 | Orb['XSecIP'] = XSecIP |
|---|
| 290 | if IfBe == 0: |
|---|
| 291 | Orb['SEdge'] = XSect[10]/AU |
|---|
| 292 | Nval = 11 |
|---|
| 293 | else: |
|---|
| 294 | Orb['ElEterm'] = XSect[10] |
|---|
| 295 | del Energy[10] |
|---|
| 296 | del XSect[10] |
|---|
| 297 | Nval = 10 |
|---|
| 298 | Orb['SEdge'] = 0.0 |
|---|
| 299 | Orb['Nval'] = Nval |
|---|
| 300 | D = dict(zip(Energy,XSect)) |
|---|
| 301 | Energy.sort() |
|---|
| 302 | X = [] |
|---|
| 303 | for key in Energy: |
|---|
| 304 | X.append(D[key]) |
|---|
| 305 | XSect = X |
|---|
| 306 | LEner = [] |
|---|
| 307 | LXSect = [] |
|---|
| 308 | for i in range(Nval): |
|---|
| 309 | LEner.append(math.log(Energy[i])) |
|---|
| 310 | if XSect[i] > 0.0: |
|---|
| 311 | LXSect.append(math.log(XSect[i])) |
|---|
| 312 | else: |
|---|
| 313 | LXSect.append(0.0) |
|---|
| 314 | Orb['LEner'] = LEner |
|---|
| 315 | Orb['LXSect'] = LXSect |
|---|
| 316 | Orbs.append(Orb) |
|---|
| 317 | xsec.close() |
|---|
| 318 | return Orbs |
|---|
| 319 | |
|---|
| 320 | def GetMagFormFacCoeff(El): |
|---|
| 321 | """Read magnetic form factor data from atmdata.py |
|---|
| 322 | |
|---|
| 323 | :param El: 2 character element symbol |
|---|
| 324 | :return: MagFormFactors: list of all magnetic form factors dictionaries for element El. |
|---|
| 325 | |
|---|
| 326 | each dictionary contains: |
|---|
| 327 | |
|---|
| 328 | * 'Symbol':Symbol |
|---|
| 329 | * 'Z':Z |
|---|
| 330 | * 'mfa': 4 MA coefficients |
|---|
| 331 | * 'nfa': 4 NA coefficients |
|---|
| 332 | * 'mfb': 4 MB coefficients |
|---|
| 333 | * 'nfb': 4 NB coefficients |
|---|
| 334 | * 'mfc': MC coefficient |
|---|
| 335 | * 'nfc': NC coefficient |
|---|
| 336 | |
|---|
| 337 | """ |
|---|
| 338 | Els = El.capitalize().strip() |
|---|
| 339 | MagFormFactors = [] |
|---|
| 340 | mags = [ky for ky in atmdata.MagFF.keys() if Els == getElSym(ky)] |
|---|
| 341 | for mag in mags: |
|---|
| 342 | magData = {} |
|---|
| 343 | data = atmdata.MagFF[mag] |
|---|
| 344 | magData['Symbol'] = mag.upper() |
|---|
| 345 | magData['Z'] = atmdata.XrayFF[getElSym(mag)]['Z'] |
|---|
| 346 | magData['mfa'] = [data['M'][i] for i in [0,2,4,6]] |
|---|
| 347 | magData['mfb'] = [data['M'][i] for i in [1,3,5,7]] |
|---|
| 348 | magData['mfc'] = data['M'][8] |
|---|
| 349 | magData['nfa'] = [data['N'][i] for i in [0,2,4,6]] |
|---|
| 350 | magData['nfb'] = [data['N'][i] for i in [1,3,5,7]] |
|---|
| 351 | magData['nfc'] = data['N'][8] |
|---|
| 352 | MagFormFactors.append(magData) |
|---|
| 353 | return MagFormFactors |
|---|
| 354 | |
|---|
| 355 | def ScatFac(El, SQ): |
|---|
| 356 | """compute value of form factor |
|---|
| 357 | |
|---|
| 358 | :param El: element dictionary defined in GetFormFactorCoeff |
|---|
| 359 | :param SQ: (sin-theta/lambda)**2 |
|---|
| 360 | :return: real part of form factor |
|---|
| 361 | """ |
|---|
| 362 | fa = np.array(El['fa']) |
|---|
| 363 | fb = np.array(El['fb']) |
|---|
| 364 | t = -fb[:,np.newaxis]*SQ |
|---|
| 365 | return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc'] |
|---|
| 366 | |
|---|
| 367 | def MagScatFac(El, SQ): |
|---|
| 368 | """compute value of form factor |
|---|
| 369 | |
|---|
| 370 | :param El: element dictionary defined in GetFormFactorCoeff |
|---|
| 371 | :param SQ: (sin-theta/lambda)**2 |
|---|
| 372 | :param gfac: Lande g factor (normally = 2.0) |
|---|
| 373 | :return: real part of form factor |
|---|
| 374 | """ |
|---|
| 375 | mfa = np.array(El['mfa']) |
|---|
| 376 | mfb = np.array(El['mfb']) |
|---|
| 377 | nfa = np.array(El['nfa']) |
|---|
| 378 | nfb = np.array(El['nfb']) |
|---|
| 379 | mt = -mfb[:,np.newaxis]*SQ |
|---|
| 380 | nt = -nfb[:,np.newaxis]*SQ |
|---|
| 381 | MMF = np.sum(mfa[:,np.newaxis]*np.exp(mt)[:],axis=0)+El['mfc'] |
|---|
| 382 | MMF0 = np.sum(mfa)+El['mfc'] |
|---|
| 383 | NMF = np.sum(nfa[:,np.newaxis]*np.exp(nt)[:],axis=0)+El['nfc'] |
|---|
| 384 | NMF0 = np.sum(nfa)+El['nfc'] |
|---|
| 385 | MF0 = MMF0+(2.0/El['gfac']-1.0)*NMF0 |
|---|
| 386 | return (MMF+(2.0/El['gfac']-1.0)*NMF)/MF0 |
|---|
| 387 | |
|---|
| 388 | def BlenResCW(Els,BLtables,wave): |
|---|
| 389 | FP = np.zeros(len(Els)) |
|---|
| 390 | FPP = np.zeros(len(Els)) |
|---|
| 391 | for i,El in enumerate(Els): |
|---|
| 392 | BL = BLtables[El][1] |
|---|
| 393 | if 'BW-LS' in BL: |
|---|
| 394 | Re,Im,E0,gam,A,E1,B,E2 = BL['BW-LS'][1:] |
|---|
| 395 | Emev = 81.80703/wave**2 |
|---|
| 396 | T0 = Emev-E0 |
|---|
| 397 | T1 = Emev-E1 |
|---|
| 398 | T2 = Emev-E2 |
|---|
| 399 | D0 = T0**2+gam**2 |
|---|
| 400 | D1 = T1**2+gam**2 |
|---|
| 401 | D2 = T2**2+gam**2 |
|---|
| 402 | FP[i] = Re*(T0/D0+A*T1/D1+B*T2/D2)+BL['BW-LS'][0] |
|---|
| 403 | FPP[i] = -Im*(1/D0+A/D1+B/D2) |
|---|
| 404 | else: |
|---|
| 405 | FPP[i] = BL['SL'][1] #for Li, B, etc. |
|---|
| 406 | return FP,FPP |
|---|
| 407 | |
|---|
| 408 | def BlenResTOF(Els,BLtables,wave): |
|---|
| 409 | FP = np.zeros((len(Els),len(wave))) |
|---|
| 410 | FPP = np.zeros((len(Els),len(wave))) |
|---|
| 411 | BL = [BLtables[el][1] for el in Els] |
|---|
| 412 | for i,El in enumerate(Els): |
|---|
| 413 | if 'BW-LS' in BL[i]: |
|---|
| 414 | Re,Im,E0,gam,A,E1,B,E2 = BL[i]['BW-LS'][1:] |
|---|
| 415 | Emev = 81.80703/wave**2 |
|---|
| 416 | T0 = Emev-E0 |
|---|
| 417 | T1 = Emev-E1 |
|---|
| 418 | T2 = Emev-E2 |
|---|
| 419 | D0 = T0**2+gam**2 |
|---|
| 420 | D1 = T1**2+gam**2 |
|---|
| 421 | D2 = T2**2+gam**2 |
|---|
| 422 | FP[i] = Re*(T0/D0+A*T1/D1+B*T2/D2)+BL[i]['BW-LS'][0] |
|---|
| 423 | FPP[i] = -Im*(1/D0+A/D1+B/D2) |
|---|
| 424 | else: |
|---|
| 425 | FPP[i] = np.ones(len(wave))*BL[i]['SL'][1] #for Li, B, etc. |
|---|
| 426 | return FP,FPP |
|---|
| 427 | |
|---|
| 428 | def ComptonFac(El,SQ): |
|---|
| 429 | """compute Compton scattering factor |
|---|
| 430 | |
|---|
| 431 | :param El: element dictionary |
|---|
| 432 | :param SQ: (sin-theta/lambda)**2 |
|---|
| 433 | :return: compton scattering factor |
|---|
| 434 | """ |
|---|
| 435 | ca = np.array(El['cmpa']) |
|---|
| 436 | cb = np.array(El['cmpb']) |
|---|
| 437 | t = -cb[:,np.newaxis]*SQ |
|---|
| 438 | return El['cmpz']-np.sum(ca[:,np.newaxis]*np.exp(t),axis=0) |
|---|
| 439 | |
|---|
| 440 | def FPcalc(Orbs, KEv): |
|---|
| 441 | """Compute real & imaginary resonant X-ray scattering factors |
|---|
| 442 | |
|---|
| 443 | :param Orbs: list of orbital dictionaries as defined in GetXsectionCoeff |
|---|
| 444 | :param KEv: x-ray energy in keV |
|---|
| 445 | :return: C: (f',f",mu): real, imaginary parts of resonant scattering & atomic absorption coeff. |
|---|
| 446 | """ |
|---|
| 447 | def Aitken(Orb, LKev): |
|---|
| 448 | Nval = Orb['Nval'] |
|---|
| 449 | j = Nval-1 |
|---|
| 450 | LEner = Orb['LEner'] |
|---|
| 451 | for i in range(Nval): |
|---|
| 452 | if LEner[i] <= LKev: j = i |
|---|
| 453 | if j > Nval-3: j= Nval-3 |
|---|
| 454 | T = [0,0,0,0,0,0] |
|---|
| 455 | LXSect = Orb['LXSect'] |
|---|
| 456 | for i in range(3): |
|---|
| 457 | T[i] = LXSect[i+j] |
|---|
| 458 | T[i+3] = LEner[i+j]-LKev |
|---|
| 459 | T[1] = (T[0]*T[4]-T[1]*T[3])/(LEner[j+1]-LEner[j]) |
|---|
| 460 | T[2] = (T[0]*T[5]-T[2]*T[3])/(LEner[j+2]-LEner[j]) |
|---|
| 461 | T[2] = (T[1]*T[5]-T[2]*T[4])/(LEner[j+2]-LEner[j+1]) |
|---|
| 462 | C = T[2] |
|---|
| 463 | return C |
|---|
| 464 | |
|---|
| 465 | def DGauss(Orb,CX,RX,ISig): |
|---|
| 466 | ALG = (0.11846344252810,0.23931433524968,0.284444444444, |
|---|
| 467 | 0.23931433524968,0.11846344252810) |
|---|
| 468 | XLG = (0.04691007703067,0.23076534494716,0.5, |
|---|
| 469 | 0.76923465505284,0.95308992296933) |
|---|
| 470 | |
|---|
| 471 | D = 0.0 |
|---|
| 472 | B2 = Orb['BB']**2 |
|---|
| 473 | R2 = RX**2 |
|---|
| 474 | XSecIP = Orb['XSecIP'] |
|---|
| 475 | for i in range(5): |
|---|
| 476 | X = XLG[i] |
|---|
| 477 | X2 = X**2 |
|---|
| 478 | XS = XSecIP[i] |
|---|
| 479 | if ISig == 0: |
|---|
| 480 | S = BB*(XS*(B2/X2)-CX*R2)/(R2*X2-B2) |
|---|
| 481 | elif ISig == 1: |
|---|
| 482 | S = 0.5*BB*B2*XS/(math.sqrt(X)*(R2*X2-X*B2)) |
|---|
| 483 | elif ISig == 2: |
|---|
| 484 | T = X*X2*R2-B2/X |
|---|
| 485 | S = 2.0*BB*(XS*B2/(T*X2**2)-(CX*R2/T)) |
|---|
| 486 | else: |
|---|
| 487 | S = BB*B2*(XS-Orb['SEdge']*X2)/(R2*X2**2-X2*B2) |
|---|
| 488 | A = ALG[i] |
|---|
| 489 | D += A*S |
|---|
| 490 | return D |
|---|
| 491 | |
|---|
| 492 | AU = 2.80022e+7 |
|---|
| 493 | C1 = 0.02721 |
|---|
| 494 | C = 137.0367 |
|---|
| 495 | FP = 0.0 |
|---|
| 496 | FPP = 0.0 |
|---|
| 497 | Mu = 0.0 |
|---|
| 498 | LKev = math.log(KEv) |
|---|
| 499 | RX = KEv/C1 |
|---|
| 500 | if Orbs: |
|---|
| 501 | for Orb in Orbs: |
|---|
| 502 | CX = 0.0 |
|---|
| 503 | BB = Orb['BB'] |
|---|
| 504 | BindEn = Orb['BindEn'] |
|---|
| 505 | if Orb['IfBe'] != 0: ElEterm = Orb['ElEterm'] |
|---|
| 506 | if BindEn <= KEv: |
|---|
| 507 | CX = math.exp(Aitken(Orb,LKev)) |
|---|
| 508 | Mu += CX |
|---|
| 509 | CX /= AU |
|---|
| 510 | Corr = 0.0 |
|---|
| 511 | if Orb['IfBe'] == 0 and BindEn >= KEv: |
|---|
| 512 | CX = 0.0 |
|---|
| 513 | FPI = DGauss(Orb,CX,RX,3) |
|---|
| 514 | Corr = 0.5*Orb['SEdge']*BB**2*math.log((RX-BB)/(-RX-BB))/RX |
|---|
| 515 | else: |
|---|
| 516 | FPI = DGauss(Orb,CX,RX,Orb['IfBe']) |
|---|
| 517 | if CX != 0.0: Corr = -0.5*CX*RX*math.log((RX+BB)/(RX-BB)) |
|---|
| 518 | FPI = (FPI+Corr)*C/(2.0*math.pi**2) |
|---|
| 519 | FPPI = C*CX*RX/(4.0*math.pi) |
|---|
| 520 | FP += FPI |
|---|
| 521 | FPP += FPPI |
|---|
| 522 | FP -= ElEterm |
|---|
| 523 | |
|---|
| 524 | return (FP, FPP, Mu) |
|---|
| 525 | |
|---|
| 526 | |
|---|
