source: trunk/GSASIIElem.py @ 230

Last change on this file since 230 was 230, checked in by vondreele, 11 years ago

add svn keywords & add log intensity plotting

  • Property svn:keywords set to Date Author Revision URL Id
File size: 13.4 KB
Line 
1"""Element: functions for element types
2   Copyright: 2008, Robert B. Von Dreele (Argonne National Laboratory)
3"""
4########### SVN repository information ###################
5# $Date: 2011-01-07 19:27:24 +0000 (Fri, 07 Jan 2011) $
6# $Author: vondreele $
7# $Revision: 230 $
8# $URL: trunk/GSASIIElem.py $
9# $Id: GSASIIElem.py 230 2011-01-07 19:27:24Z vondreele $
10########### SVN repository information ###################
11
12import wx
13import math
14import sys
15import os.path
16import  wx.lib.colourselect as wscs
17import GSASIIpath
18import numpy as np
19
20def GetFormFactorCoeff(El):
21    """Read form factor coefficients from atomdata.asc file
22    @param El: element 1-2 character symbol case irrevelant
23    @return: FormFactors: list of form factor dictionaries
24    each dictionary is:
25    'Symbol':4 character element symbol with valence (e.g. 'NI+2')
26    'Z': atomic number
27    'fa': 4 A coefficients
28    'fb':4 B coefficients
29    'fc': C coefficient
30    """
31    ElS = El.upper()
32    ElS = ElS.rjust(2)
33    filename = os.path.join(sys.path[1],'atmdata.dat')
34    try:
35        FFdata = open(filename,'Ur')
36    except:
37        wx.MessageBox(message="File atmdata.dat not found in directory %s" % sys.path[0],
38            caption="No atmdata.dat file",style=wx.OK | wx.ICON_EXCLAMATION | wx.STAY_ON_TOP)
39        sys.exit()
40    S = '1'
41    FormFactors = []
42    while S:
43        S = FFdata.readline()
44        if S[3:5] == ElS:
45            if S[5:6] != '_':
46                Z=int(S[:2])
47                Symbol = S[3:7].strip()
48                S = S[12:]
49                fa = (float(S[:7]),float(S[14:21]),float(S[28:35]),float(S[42:49]))
50                fb = (float(S[7:14]),float(S[21:28]),float(S[35:42]),float(S[49:56]))
51                FormFac = {'Symbol':Symbol,'Z':Z,'fa':fa,'fb':fb,'fc':float(S[56:63])}
52                FormFactors.append(FormFac)               
53    FFdata.close()
54    return FormFactors
55   
56def GetAtomInfo(El):
57   
58    import ElementTable as ET
59    Elements = []
60    for elem in ET.ElTable:
61        Elements.append(elem[0][0])
62    if len(El) in [2,4]:
63        ElS = El.upper()[:2].rjust(2)
64    else:
65        ElS = El.upper()[:1].rjust(2)
66    filename = os.path.join(sys.path[1],'atmdata.dat')
67    try:
68        FFdata = open(filename,'Ur')
69    except:
70        wx.MessageBox(message="File atmdata.dat not found in directory %s" % sys.path[0],
71            caption="No atmdata.dat file",style=wx.OK | wx.ICON_EXCLAMATION | wx.STAY_ON_TOP)
72        sys.exit()
73    S = '1'
74    AtomInfo = {}
75    Mass = []
76    while S:
77        S = FFdata.readline()
78        if S[3:5] == ElS:
79            if S[5:6] == '_':
80                if not Mass:                                 #picks 1st one; natural abundance or 1st isotope
81                    Mass = float(S[10:19])
82                if S[5:9] == '_SIZ':
83                    Z=int(S[:2])
84                    Symbol = S[3:5].strip().lower().capitalize()
85                    Drad = float(S[12:22])
86                    Arad = float(S[22:32])
87                    Vdrad = float(S[32:38])
88                    Color = ET.ElTable[Elements.index(Symbol)][6]
89    FFdata.close()
90    AtomInfo={'Symbol':Symbol,'Mass':Mass,'Z':Z,'Drad':Drad,'Arad':Arad,'Vdrad':Vdrad,'Color':Color}   
91    return AtomInfo
92     
93def GetXsectionCoeff(El):
94    """Read atom orbital scattering cross sections for fprime calculations via Cromer-Lieberman algorithm
95    @param El: 2 character element symbol
96    @return: Orbs: list of orbitals each a dictionary with detailed orbital information used by FPcalc
97    each dictionary is:
98    'OrbName': Orbital name read from file
99    'IfBe' 0/2 depending on orbital
100    'BindEn': binding energy
101    'BB': BindEn/0.02721
102    'XSectIP': 5 cross section inflection points
103    'ElEterm': energy correction term
104    'SEdge': absorption edge for orbital
105    'Nval': 10/11 depending on IfBe
106    'LEner': 10/11 values of log(energy)
107    'LXSect': 10/11 values of log(cross section)
108    """
109    AU = 2.80022e+7
110    C1 = 0.02721
111    ElS = El.upper()
112    ElS = ElS.ljust(2)
113    filename = os.path.join(sys.path[1],'Xsect.dat')
114    try:
115        xsec = open(filename,'Ur')
116    except:
117        wx.MessageBox(message="File Xsect.dat not found in directory %s" % sys.path[0],
118            caption="No Xsect.dat file",style=wx.OK | wx.ICON_EXCLAMATION |wx.STAY_ON_TOP)
119        sys.exit()
120    S = '1'
121    Orbs = []
122    while S:
123        S = xsec.readline()
124        if S[:2] == ElS:
125            S = S[:-1]+xsec.readline()[:-1]+xsec.readline()
126            OrbName = S[9:14]
127            S = S[14:]
128            IfBe = int(S[0])
129            S = S[1:]
130            val = S.split()
131            BindEn = float(val[0])
132            BB = BindEn/C1
133            Orb = {'OrbName':OrbName,'IfBe':IfBe,'BindEn':BindEn,'BB':BB}
134            Energy = []
135            XSect = []
136            for i in range(11):
137                Energy.append(float(val[2*i+1]))
138                XSect.append(float(val[2*i+2]))
139            XSecIP = []
140            for i in range(5): XSecIP.append(XSect[i+5]/AU)
141            Orb['XSecIP'] = XSecIP
142            if IfBe == 0:
143                Orb['SEdge'] = XSect[10]/AU
144                Nval = 11
145            else:
146                Orb['ElEterm'] = XSect[10]
147                del Energy[10]
148                del XSect[10]
149                Nval = 10
150                Orb['SEdge'] = 0.0
151            Orb['Nval'] = Nval
152            D = dict(zip(Energy,XSect))
153            Energy.sort()
154            X = []
155            for key in Energy:
156                X.append(D[key])
157            XSect = X
158            LEner = []
159            LXSect = []
160            for i in range(Nval):
161                LEner.append(math.log(Energy[i]))
162                if XSect[i] > 0.0:
163                    LXSect.append(math.log(XSect[i]))
164                else:
165                    LXSect.append(0.0)
166            Orb['LEner'] = LEner
167            Orb['LXSect'] = LXSect
168            Orbs.append(Orb)
169    xsec.close()
170    return Orbs
171   
172def GetMagFormFacCoeff(El):
173    """Read magnetic form factor data from atomdata.asc file
174    @param El: 2 character element symbol
175    @return: MagFormFactors: list of all magnetic form factors dictionaries for element El.
176    each dictionary contains:
177    'Symbol':Symbol
178    'Z':Z
179    'mfa': 4 MA coefficients
180    'nfa': 4 NA coefficients
181    'mfb': 4 MB coefficients
182    'nfb': 4 NB coefficients
183    'mfc': MC coefficient
184    'nfc': NC coefficient
185    """
186    ElS = El.upper()
187    ElS = ElS.rjust(2)
188    filename = os.path.join(sys.path[1],'atmdata.dat')
189    try:
190        FFdata = open(filename,'Ur')
191    except:
192        wx.MessageBox(message="File atmdata.dat not found in directory %s" % sys.path[0],
193            caption="No atmdata.dat file",style=wx.OK | wx.ICON_EXCLAMATION |wx.STAY_ON_TOP)
194        sys.exit()
195    S = '1'
196    MagFormFactors = []
197    while S:
198        S = FFdata.readline()
199        if S[3:5] == ElS:
200            if S[8:9] == 'M':
201                SN = FFdata.readline()               #'N' is assumed to follow 'M' in Atomdata.asc
202                Z=int(S[:2])
203                Symbol = S[3:7]
204                S = S[12:]
205                SN = SN[12:]
206                mfa = (float(S[:7]),float(S[14:21]),float(S[28:35]),float(S[42:49]))
207                mfb = (float(S[7:14]),float(S[21:28]),float(S[35:42]),float(S[49:56]))
208                nfa = (float(SN[:7]),float(SN[14:21]),float(SN[28:35]),float(SN[42:49]))
209                nfb = (float(SN[7:14]),float(SN[21:28]),float(SN[35:42]),float(SN[49:56]))
210                FormFac = {'Symbol':Symbol,'Z':Z,'mfa':mfa,'nfa':nfa,'mfb':mfb,'nfb':nfb,
211                    'mfc':float(S[56:63]),'nfc':float(SN[56:63])}
212                MagFormFactors.append(FormFac)
213    FFdata.close()
214    return MagFormFactors
215
216def ScatFac(FormFac, SQ):
217    """compute value of form factor
218    @param FormFac: dictionary  defined in GetFormFactorCoeff
219    @param SQ: (sin-theta/lambda)**2
220    @return: f: real part of form factor
221    """
222    fa = np.array(FormFac['fa'])
223    fb = np.array(FormFac['fb'])
224    t = -fb*SQ
225    return np.sum(fa*np.exp(t))+FormFac['fc']
226           
227def FPcalc(Orbs, KEv):
228    """Compute real & imaginary resonant X-ray scattering factors
229    @param Orbs: list of orbital dictionaries as defined in GetXsectionCoeff
230    @param KEv: x-ray energy in keV
231    @return: C: (f',f",mu): real, imaginary parts of resonant scattering & atomic absorption coeff.
232    """
233    def Aitken(Orb, LKev):
234        Nval = Orb['Nval']
235        j = Nval-1
236        LEner = Orb['LEner']
237        for i in range(Nval):
238            if LEner[i] <= LKev: j = i
239        if j > Nval-3: j= Nval-3
240        T = [0,0,0,0,0,0]
241        LXSect = Orb['LXSect']
242        for i in range(3):
243           T[i] = LXSect[i+j]
244           T[i+3] = LEner[i+j]-LKev
245        T[1] = (T[0]*T[4]-T[1]*T[3])/(LEner[j+1]-LEner[j])
246        T[2] = (T[0]*T[5]-T[2]*T[3])/(LEner[j+2]-LEner[j])
247        T[2] = (T[1]*T[5]-T[2]*T[4])/(LEner[j+2]-LEner[j+1])
248        C = T[2]
249        return C
250   
251    def DGauss(Orb,CX,RX,ISig):
252        ALG = (0.11846344252810,0.23931433524968,0.284444444444,
253        0.23931433524968,0.11846344252810)
254        XLG = (0.04691007703067,0.23076534494716,0.5,
255        0.76923465505284,0.95308992296933)
256       
257        D = 0.0
258        B2 = Orb['BB']**2
259        R2 = RX**2
260        XSecIP = Orb['XSecIP']
261        for i in range(5):
262            X = XLG[i]
263            X2 = X**2
264            XS = XSecIP[i]
265            if ISig == 0:
266                S = BB*(XS*(B2/X2)-CX*R2)/(R2*X2-B2)
267            elif ISig == 1:
268                S = 0.5*BB*B2*XS/(math.sqrt(X)*(R2*X2-X*B2))
269            elif ISig == 2:
270                T = X*X2*R2-B2/X
271                S = 2.0*BB*(XS*B2/(T*X2**2)-(CX*R2/T))
272            else:
273                S = BB*B2*(XS-Orb['SEdge']*X2)/(R2*X2**2-X2*B2)
274            A = ALG[i]
275            D += A*S
276        return D
277   
278    AU = 2.80022e+7
279    C1 = 0.02721
280    C = 137.0367
281    FP = 0.0
282    FPP = 0.0
283    Mu = 0.0
284    LKev = math.log(KEv)
285    RX = KEv/C1
286    if Orbs:
287        for Orb in Orbs:
288            CX = 0.0
289            BB = Orb['BB']
290            BindEn = Orb['BindEn']
291            if Orb['IfBe'] != 0: ElEterm = Orb['ElEterm']
292            if BindEn <= KEv:
293                CX = math.exp(Aitken(Orb,LKev))
294                Mu += CX
295                CX /= AU
296            Corr = 0.0
297            if Orb['IfBe'] == 0 and BindEn >= KEv:
298                CX = 0.0
299                FPI = DGauss(Orb,CX,RX,3)
300                Corr = 0.5*Orb['SEdge']*BB**2*math.log((RX-BB)/(-RX-BB))/RX
301            else:
302                FPI = DGauss(Orb,CX,RX,Orb['IfBe'])
303                if CX != 0.0: Corr = -0.5*CX*RX*math.log((RX+BB)/(RX-BB))
304            FPI = (FPI+Corr)*C/(2.0*math.pi**2)
305            FPPI = C*CX*RX/(4.0*math.pi)
306            FP += FPI
307            FPP += FPPI
308        FP -= ElEterm
309   
310    return (FP, FPP, Mu)
311   
312
313class PickElement(wx.Dialog):
314    "Makes periodic table widget for picking element - caller maintains element list"
315    Elem=None
316    def _init_ctrls(self, prnt):
317        wx.Dialog.__init__(self, id=-1, name='PickElement',
318              parent=prnt, pos=wx.DefaultPosition, 
319              style=wx.DEFAULT_DIALOG_STYLE, title='Pick Element')
320        import ElementTable as ET
321        self.SetClientSize(wx.Size(770, 250))
322       
323        i=0
324        for E in ET.ElTable:
325            PickElement.ElButton(self,name=E[0],
326#            pos=wx.Point(E[1]*40+25,E[2]*24+24),tip=E[3],color=E[4])
327            pos=wx.Point(E[1]*40+25,E[2]*24+24),tip=E[3],color=E[6])
328            i+=1
329
330    def __init__(self, parent):
331        self._init_ctrls(parent)
332       
333    def ElButton(self, name, pos, tip, color):
334        Black = wx.Colour(0,0,0)
335        El = wx.ComboBox(choices=name, parent=self, pos=pos, size=wx.Size(40,23),
336            style=wx.CB_READONLY, value=name[0])
337        El.SetBackgroundColour(color)
338        El.SetToolTipString(tip)
339        El.Bind(wx.EVT_COMBOBOX, self.OnElButton)
340
341    def OnElButton(self, event):
342        El = event.GetEventObject().GetLabel()
343        self.Elem = (El)
344        self.EndModal(wx.ID_OK)       
345       
346class DeleteElement(wx.Dialog):
347    "Delete element from selected set widget"
348    def _init_ctrls(self, parent):
349        l = len(DeleteElement.Elems)-1
350        wx.Dialog.__init__(self, id=-1, name='Delete', parent=parent,
351              pos=wx.DefaultPosition, size=wx.Size(max(128,64+l*24), 87),
352              style=wx.DEFAULT_DIALOG_STYLE, title='Delete Element')
353        self.Show(True)
354        self.SetAutoLayout(True)
355        self.SetHelpText('Select element to delete')
356        self.SetWindowVariant(wx.WINDOW_VARIANT_SMALL)
357
358        i = 0
359        Elem = []
360        for Elem in DeleteElement.Elems:
361            name = Elem[0].lower().capitalize()
362            self.ElButton(id=-1,name=name,pos=wx.Point(16+i*24, 16))
363            i+=1
364             
365    def __init__(self, parent):
366        DeleteElement.Elems = parent.Elems
367        DeleteElement.El = ' '
368        self._init_ctrls(parent)
369
370    def ElButton(self, id, name, pos):
371        White = wx.Colour(255, 255, 255)
372        El = wscs.ColourSelect(label=name, parent=self, colour = White,
373            pos=pos, size=wx.Size(24, 23), style=wx.RAISED_BORDER)
374        El.Bind(wx.EVT_BUTTON, self.OnDeleteButton)
375   
376    def OnDeleteButton(self, event):
377        DeleteElement.El=event.GetEventObject().GetLabel()
378        self.EndModal(wx.ID_OK)
379       
380    def GetDeleteElement(self):
381        return DeleteElement.El
382       
383
Note: See TracBrowser for help on using the repository browser.