source: trunk/GSASIIElem.py @ 271

Last change on this file since 271 was 265, checked in by vondreele, 14 years ago

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