source: trunk/GSASIIElem.py @ 264

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

begin implementation of pdf generation - work in progress
modify image to have azimuth=0 as vertical "up"
add textctrls for max and min image intensity
fix to some lattice parameter defaults for some Laue groups

  • Property svn:keywords set to Date Author Revision URL Id
File size: 14.1 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-03-31 21:41:13 +0000 (Thu, 31 Mar 2011) $
6# $Author: vondreele $
7# $Revision: 264 $
8# $URL: trunk/GSASIIElem.py $
9# $Id: GSASIIElem.py 264 2011-03-31 21:41:13Z 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 ComptonFac(ComptonCoeff,SQ):
228    """compute Compton scattering factor
229    @param ComptonCoeff: list [Z, a1:a5, b1:b5]
230    @param SQ: (sin-theta/lambda)**2
231    @return: comp: compton scattering factor
232    """   
233    ca = np.array(ComptonCoeff[1:6])
234    cb = np.array(ComptonCoeff[6:11])
235    t = -cb*SQ
236    return ComptonCoeff[0]-np.sum(ca*np.exp(t)) 
237           
238def FPcalc(Orbs, KEv):
239    """Compute real & imaginary resonant X-ray scattering factors
240    @param Orbs: list of orbital dictionaries as defined in GetXsectionCoeff
241    @param KEv: x-ray energy in keV
242    @return: C: (f',f",mu): real, imaginary parts of resonant scattering & atomic absorption coeff.
243    """
244    def Aitken(Orb, LKev):
245        Nval = Orb['Nval']
246        j = Nval-1
247        LEner = Orb['LEner']
248        for i in range(Nval):
249            if LEner[i] <= LKev: j = i
250        if j > Nval-3: j= Nval-3
251        T = [0,0,0,0,0,0]
252        LXSect = Orb['LXSect']
253        for i in range(3):
254           T[i] = LXSect[i+j]
255           T[i+3] = LEner[i+j]-LKev
256        T[1] = (T[0]*T[4]-T[1]*T[3])/(LEner[j+1]-LEner[j])
257        T[2] = (T[0]*T[5]-T[2]*T[3])/(LEner[j+2]-LEner[j])
258        T[2] = (T[1]*T[5]-T[2]*T[4])/(LEner[j+2]-LEner[j+1])
259        C = T[2]
260        return C
261   
262    def DGauss(Orb,CX,RX,ISig):
263        ALG = (0.11846344252810,0.23931433524968,0.284444444444,
264        0.23931433524968,0.11846344252810)
265        XLG = (0.04691007703067,0.23076534494716,0.5,
266        0.76923465505284,0.95308992296933)
267       
268        D = 0.0
269        B2 = Orb['BB']**2
270        R2 = RX**2
271        XSecIP = Orb['XSecIP']
272        for i in range(5):
273            X = XLG[i]
274            X2 = X**2
275            XS = XSecIP[i]
276            if ISig == 0:
277                S = BB*(XS*(B2/X2)-CX*R2)/(R2*X2-B2)
278            elif ISig == 1:
279                S = 0.5*BB*B2*XS/(math.sqrt(X)*(R2*X2-X*B2))
280            elif ISig == 2:
281                T = X*X2*R2-B2/X
282                S = 2.0*BB*(XS*B2/(T*X2**2)-(CX*R2/T))
283            else:
284                S = BB*B2*(XS-Orb['SEdge']*X2)/(R2*X2**2-X2*B2)
285            A = ALG[i]
286            D += A*S
287        return D
288   
289    AU = 2.80022e+7
290    C1 = 0.02721
291    C = 137.0367
292    FP = 0.0
293    FPP = 0.0
294    Mu = 0.0
295    LKev = math.log(KEv)
296    RX = KEv/C1
297    if Orbs:
298        for Orb in Orbs:
299            CX = 0.0
300            BB = Orb['BB']
301            BindEn = Orb['BindEn']
302            if Orb['IfBe'] != 0: ElEterm = Orb['ElEterm']
303            if BindEn <= KEv:
304                CX = math.exp(Aitken(Orb,LKev))
305                Mu += CX
306                CX /= AU
307            Corr = 0.0
308            if Orb['IfBe'] == 0 and BindEn >= KEv:
309                CX = 0.0
310                FPI = DGauss(Orb,CX,RX,3)
311                Corr = 0.5*Orb['SEdge']*BB**2*math.log((RX-BB)/(-RX-BB))/RX
312            else:
313                FPI = DGauss(Orb,CX,RX,Orb['IfBe'])
314                if CX != 0.0: Corr = -0.5*CX*RX*math.log((RX+BB)/(RX-BB))
315            FPI = (FPI+Corr)*C/(2.0*math.pi**2)
316            FPPI = C*CX*RX/(4.0*math.pi)
317            FP += FPI
318            FPP += FPPI
319        FP -= ElEterm
320   
321    return (FP, FPP, Mu)
322   
323
324class PickElement(wx.Dialog):
325    "Makes periodic table widget for picking element - caller maintains element list"
326    Elem=None
327    def _init_ctrls(self, prnt,oneOnly):
328        wx.Dialog.__init__(self, id=-1, name='PickElement',
329              parent=prnt, pos=wx.DefaultPosition, 
330              style=wx.DEFAULT_DIALOG_STYLE, title='Pick Element')
331        import ElementTable as ET
332        self.SetClientSize(wx.Size(770, 250))
333       
334        i=0
335        for E in ET.ElTable:
336            if oneOnly:
337                color=E[4]
338            else:
339                color=E[6]
340            PickElement.ElButton(self,name=E[0],
341               pos=wx.Point(E[1]*40+25,E[2]*24+24),tip=E[3],color=color,oneOnly=oneOnly)
342            i+=1
343
344    def __init__(self, parent,oneOnly=False):
345        self._init_ctrls(parent,oneOnly)
346       
347    def ElButton(self, name, pos, tip, color, oneOnly):
348        Black = wx.Colour(0,0,0)
349        if oneOnly:
350            El = wscs.ColourSelect(label=name[0], parent=self,colour=color,
351                pos=pos, size=wx.Size(40,23), style=wx.RAISED_BORDER)
352#            El.SetLabel(name)
353            El.Bind(wx.EVT_BUTTON, self.OnElButton)
354        else:
355            El = wx.ComboBox(choices=name, parent=self, pos=pos, size=wx.Size(40,23),
356                style=wx.CB_READONLY, value=name[0])
357            El.Bind(wx.EVT_COMBOBOX,self.OnElButton)
358       
359        El.SetBackgroundColour(color)
360        El.SetToolTipString(tip)
361
362    def OnElButton(self, event):
363        El = event.GetEventObject().GetLabel()
364        self.Elem = El
365        self.EndModal(wx.ID_OK)       
366       
367class DeleteElement(wx.Dialog):
368    "Delete element from selected set widget"
369    def _init_ctrls(self, parent):
370        l = len(DeleteElement.Elems)-1
371        wx.Dialog.__init__(self, id=-1, name='Delete', parent=parent,
372              pos=wx.DefaultPosition, size=wx.Size(max(128,64+l*24), 87),
373              style=wx.DEFAULT_DIALOG_STYLE, title='Delete Element')
374        self.Show(True)
375        self.SetAutoLayout(True)
376        self.SetHelpText('Select element to delete')
377        self.SetWindowVariant(wx.WINDOW_VARIANT_SMALL)
378
379        i = 0
380        Elem = []
381        for Elem in DeleteElement.Elems:
382            name = Elem[0].lower().capitalize()
383            self.ElButton(id=-1,name=name,pos=wx.Point(16+i*24, 16))
384            i+=1
385             
386    def __init__(self, parent):
387        DeleteElement.Elems = parent.Elems
388        DeleteElement.El = ' '
389        self._init_ctrls(parent)
390
391    def ElButton(self, id, name, pos):
392        White = wx.Colour(255, 255, 255)
393        El = wscs.ColourSelect(label=name, parent=self, colour = White,
394            pos=pos, size=wx.Size(24, 23), style=wx.RAISED_BORDER)
395        El.Bind(wx.EVT_BUTTON, self.OnDeleteButton)
396   
397    def OnDeleteButton(self, event):
398        DeleteElement.El=event.GetEventObject().GetLabel()
399        self.EndModal(wx.ID_OK)
400       
401    def GetDeleteElement(self):
402        return DeleteElement.El
403       
404
Note: See TracBrowser for help on using the repository browser.