source: trunk/GSASIIElem.py @ 409

Last change on this file since 409 was 409, checked in by vondreele, 10 years ago

make it do CW neutrons - 1st pass at it

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