source: trunk/GSASIIElem.py @ 151

Last change on this file since 151 was 151, checked in by vondreel, 11 years ago

removed ElTable? - now in own ElementTable?.py file, removed GetAtomColors? & GetElemColor? - not to be used (ever) - now new GetAtomInfo? fixed to work with valence choices, etc.

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