source: trunk/GSASIIpwd.py @ 289

Last change on this file since 289 was 289, checked in by vondreele, 11 years ago
File size: 32.2 KB
Line 
1#GSASII powder calculation module
2########### SVN repository information ###################
3# $Date: 2011-04-20 13:09:53 -0500 (Wed, 20 Apr 2011) $
4# $Author: vondreele $
5# $Revision: 267 $
6# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIIpwd.py $
7# $Id: GSASIIpwd.py 267 2011-04-20 18:09:53Z vondreele $
8########### SVN repository information ###################
9import sys
10import math
11import wx
12import time
13import numpy as np
14import scipy as sp
15import numpy.linalg as nl
16import scipy.interpolate as si
17import scipy.stats as st
18import GSASIIpath
19import pypowder as pyp              #assumes path has been amended to include correctr bin directory
20import GSASIIplot as G2plt
21import GSASIIlattice as G2lat
22import GSASIIElem as G2elem
23import GSASIIgrid as G2gd
24
25# trig functions in degrees
26sind = lambda x: math.sin(x*math.pi/180.)
27asind = lambda x: 180.*math.asin(x)/math.pi
28tand = lambda x: math.tan(x*math.pi/180.)
29atand = lambda x: 180.*math.atan(x)/math.pi
30atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
31cosd = lambda x: math.cos(x*math.pi/180.)
32acosd = lambda x: 180.*math.acos(x)/math.pi
33rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
34#numpy versions
35npsind = lambda x: np.sin(x*np.pi/180.)
36npasind = lambda x: 180.*np.arcsin(x)/math.pi
37npcosd = lambda x: np.cos(x*math.pi/180.)
38npacosd = lambda x: 180.*np.arccos(x)/math.pi
39nptand = lambda x: np.tan(x*math.pi/180.)
40npatand = lambda x: 180.*np.arctan(x)/np.pi
41npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
42npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
43npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
44   
45np.seterr(divide='ignore')      #this is for the FCJ functions
46
47#Peak shape definitions
48# Finger-Cox_Jephcoat D(2phi,2th) function for S/L = H/L
49
50class fcjde_gen(st.rv_continuous):
51   
52    """
53    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
54    Ref: J. Appl. Cryst. (1994) 27, 892-900.
55    Parameters
56    -----------------------------------------
57    x: array like 2-theta positions
58    t: 2-theta position of peak
59    s: sum(S/L,H/L); S: sample height, H: detector opening,
60        L: sample to detector opening distance
61    Result for fcj.pdf
62    -----------------------------------------
63    if x < t & s = S/L+H/L:
64        fcj.pdf = [1/sqrt({cos(x)**2/cos(t)**2}-1) - 1/s]/cos(x)   
65    if x >= t:
66        fcj.pdf = 0   
67    """
68    def _pdf(self,x,t,s):
69        ax = npcosd(x)**2
70        bx = npcosd(t)**2
71        bx = np.where(ax>bx,bx,ax)
72        fx = np.where(ax>bx,(1./np.sqrt((ax/bx)-1.)-1./s)/npcosd(x),0.0)
73        return np.where(((x < t) & (fx > 0)),fx,0.0)
74#    def _cdf(self, x):
75#    def _ppf(self, q):
76#    def _sf(self, x):
77#    def _isf(self, q):
78#    def _stats(self):
79#    def _entropy(self):
80       
81fcjde = fcjde_gen(name='fcjde')
82               
83       
84# Finger-Cox_Jephcoat D(2phi,2th) function for S/L != H/L
85
86class fcjd_gen(st.rv_continuous):
87    """
88    Finger-Cox-Jephcoat D(2phi,2th) function for S/L != H/L
89    Ref: J. Appl. Cryst. (1994) 27, 892-900.
90    Parameters
91    -----------------------------------------
92    x: array like 2-theta positions
93    t: 2-theta position of peak
94    h: detector opening height/sample to detector opening distance
95    s: sample height/sample to detector opening distance
96    Result for fcj2.pdf
97    -----------------------------------------
98    infl = acos(cos(t)*sqrt((h-s)**2+1))
99    if x < infl:   
100        fcj.pdf = [1/sqrt({cos(x)**2/cos(t)**2}-1) - 1/shl]/cos(2phi)
101   
102    for 2phi < 2tth & shl = S/L+H/L
103   
104    fcj.pdf(x,tth,shl) = 0
105   
106    for 2phi >= 2th
107    """
108    def _pdf(self,x,t,h,s):
109        a = npcosd(t)*(np.sqrt((h-s)**2+1.))
110        infl = np.where((a <= 1.),npacosd(a),t)
111        ax = npcosd(x)**2
112        bx = npcosd(t)**2
113        bx = np.where(ax>bx,bx,ax)
114        H = np.where(ax>bx,np.sqrt((ax/bx)-1.),0.0)
115        W1 = h+s-H
116        W2 = np.where ((h > s),2.*s,2.*h)
117        fx = 2.*h*np.sqrt((ax/bx)-1.)*npcosd(x)
118        fx = np.where(fx>0.0,1./fx,0.0)
119        fx = np.where((x < infl),fx*W1,fx*W2)
120        return np.where((fx > 0.),fx,0.0)
121#    def _cdf(self, x):
122#    def _ppf(self, q):
123#    def _sf(self, x):
124#    def _isf(self, q):
125#    def _stats(self):
126#    def _entropy(self):
127       
128fcjd = fcjd_gen(name='fcjd')
129               
130# Finger-Cox_Jephcoat D(2phi,2th) function for S/L != H/L using sum & difference
131
132class fcjdsd_gen(st.rv_continuous):
133    """
134    Finger-Cox-Jephcoat D(2phi,2th) function for S/L != H/L using sum & difference
135   
136    fcj.pdf(x,tth,shl) = [1/sqrt({cos(2phi)**2/cos(2th)**2}-1) - 1/shl]/cos(2phi)
137   
138    for 2phi < 2tth & shl = S/L+H/L
139   
140    fcj.pdf(x,tth,shl) = 0
141   
142    for 2phi >= 2th
143    """
144    def _argcheck(self,t,s,d):
145        return (t > 0)&(s > 0)&(abs(d) < s)
146    def _pdf(self,x,t,s,d):
147        a = npcosd(t)*np.sqrt(d**2+1.)
148        infl = np.where((a < 1.),npacosd(a),t)
149        ax = npcosd(x)**2
150        bx = npcosd(t)**2
151        bx = np.where(ax>bx,bx,ax)
152        H = np.where(ax>bx,np.sqrt((ax/bx)-1.),0.0)
153        W1 = s-H
154        W2 = np.where ((d > 0),s-d,s+d)
155        fx = np.where(ax>bx,1./((s+d)*np.sqrt((ax/bx)-1.)*npcosd(x)),0.0)
156        fx = np.where((x < infl),fx*W1,fx*W2)
157        return np.where((fx > 0.),fx,0.0)
158#    def _cdf(self, x):
159#    def _ppf(self, q):
160#    def _sf(self, x):
161#    def _isf(self, q):
162#    def _stats(self):
163#    def _entropy(self):
164       
165fcjdsd = fcjdsd_gen(name='fcjdsd')
166               
167def factorize(num):
168    ''' Provide prime number factors for integer num
169    Returns dictionary of prime factors (keys) & power for each (data)
170    '''
171    factors = {}
172    orig = num
173
174    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
175    i, sqi = 2, 4
176    while sqi <= num:
177        while not num%i:
178            num /= i
179            factors[i] = factors.get(i, 0) + 1
180
181        sqi += 2*i + 1
182        i += 1
183
184    if num != 1 and num != orig:
185        factors[num] = factors.get(num, 0) + 1
186
187    if factors:
188        return factors
189    else:
190        return {num:1}          #a prime number!
191           
192def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
193    ''' Provide list of optimal data sizes for FFT calculations
194    Input:
195        nmin: minimum data size >= 1
196        nmax: maximum data size > nmin
197        thresh: maximum prime factor allowed
198    Returns:
199        list of data sizes where the maximum prime factor is < thresh
200    ''' 
201    plist = []
202    nmin = max(1,nmin)
203    nmax = max(nmin+1,nmax)
204    for p in range(nmin,nmax):
205        if max(factorize(p).keys()) < thresh:
206            plist.append(p)
207    return plist
208
209def Transmission(Geometry,Abs,Diam):
210#Calculate sample transmission
211#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
212#   Abs: absorption coeff in cm-1
213#   Diam: sample thickness/diameter in mm
214    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
215        MuR = Abs*Diam/20.0
216        if MuR <= 3.0:
217            T0 = 16/(3.*math.pi)
218            T1 = -0.045780
219            T2 = -0.02489
220            T3 = 0.003045
221            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
222            if T < -20.:
223                return 2.06e-9
224            else:
225                return math.exp(T)
226        else:
227            T1 = 1.433902
228            T2 = 0.013869+0.337894
229            T3 = 1.933433+1.163198
230            T4 = 0.044365-0.04259
231            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
232            return T/100.
233    elif 'plate' in Geometry:
234        MuR = Abs*Diam/10.
235        return math.exp(-MuR)
236    elif 'Bragg' in Geometry:
237        return 0.0
238
239def Absorb(Geometry,Abs,Diam,Tth,Phi=0,Psi=0):
240#Calculate sample absorption
241#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
242#   Abs: absorption coeff in cm-1
243#   Diam: sample thickness/diameter in mm
244#   Tth: 2-theta scattering angle - can be numpy array
245#   Phi: flat plate tilt angle - future
246#   Psi: flat plate tilt axis - future
247    Sth2 = npsind(Tth/2.0)**2
248    Cth2 = 1.-Sth2
249    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
250        MuR = Abs*Diam/20.0
251        if MuR < 3.0:
252            T0 = 16.0/(3*np.pi)
253            T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
254                0.109561*np.sqrt(Sth2)-26.04556
255            T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
256                1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
257            T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
258            Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
259            return np.exp(Trns)
260        else:
261            T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
262                10.02088*Sth2**3-3.36778*Sth2**4
263            T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
264                (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
265            T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
266                0.13576*np.sqrt(Sth2)+1.163198
267            T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
268            Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
269            return Trns/100.
270    elif 'Bragg' in Geometry:
271        return 1.0
272    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
273        # and only defined for 2theta < 90
274        MuR = Abs*Diam/10.0
275        T1 = np.exp(-MuR)
276        T2 = np.exp(-MuR/npcosd(Tth))
277        Tb = MuR-MuR/npcosd(Tth)
278        return (T2-T1)/Tb
279    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
280        MuR = Abs*Diam/10.0
281        cth = npcosd(Tth/2.0)
282        return np.exp(-MuR/cth)/cth
283       
284def Polarization(Pola,Tth,Azm=0.0):
285#   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
286#   Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
287#   Azm: azimuthal angle e.g. 0.0 in plane of polarization
288#   Tth: 2-theta scattering angle - can be numpy array
289#       which (if either) of these is "right"?
290#    return (Pola*npcosd(Azm)**2+(1.-Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \
291#        Pola*npsind(Azm)**2+(1.-Pola)*npcosd(Azm)**2
292    return ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
293        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
294   
295def Oblique(ObCoeff,Tth):
296    if ObCoeff:
297        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
298    else:
299        return 1.0
300               
301def Ruland(RulCoff,wave,Q,Compton):
302    C = 2.9978e8
303    D = 1.5e-3
304    hmc = 0.024262734687
305    sinth2 = (Q*wave/(4.0*np.pi))**2
306    dlam = (wave**2)*Compton*Q/C
307    dlam_c = 2.0*hmc*sinth2-D*wave**2
308    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
309   
310def LorchWeight(Q):
311    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
312           
313def GetAsfMean(ElList,Sthl2):
314#   Calculate various scattering factor terms for PDF calcs
315#   ElList: element dictionary contains scattering factor coefficients, etc.
316#   Sthl2: numpy array of sin theta/lambda squared values
317#   returns: mean(f^2), mean(f)^2, mean(compton)
318    sumNoAtoms = 0.0
319    FF = np.zeros_like(Sthl2)
320    FF2 = np.zeros_like(Sthl2)
321    CF = np.zeros_like(Sthl2)
322    for El in ElList:
323        sumNoAtoms += ElList[El]['FormulaNo']
324    for El in ElList:
325        el = ElList[El]
326        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
327        cf = G2elem.ComptonFac(el,Sthl2)
328        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
329        FF2 += ff2*el['FormulaNo']/sumNoAtoms
330        CF += cf*el['FormulaNo']/sumNoAtoms
331    return FF2,FF**2,CF
332   
333def GetNumDensity(ElList,Vol):
334    sumNoAtoms = 0.0
335    for El in ElList:
336        sumNoAtoms += ElList[El]['FormulaNo']
337    return sumNoAtoms/Vol
338           
339def MultGetQ(Tth,MuT,Geometry,b=88.0,a=0.01):
340    NS = 500
341    Gama = np.linspace(0.,np.pi/2.,NS,False)[1:]
342    Cgama = np.cos(Gama)[:,np.newaxis]
343    Sgama = np.sin(Gama)[:,np.newaxis]
344    CSgama = 1.0/Sgama
345    Delt = Gama[1]-Gama[0]
346    emc = 7.94e-26
347    Navo = 6.023e23
348    Cth = npcosd(Tth/2.0)
349    CTth = npcosd(Tth)
350    Sth = npcosd(Tth/2.0)
351    STth = npsind(Tth)
352    CSth = 1.0/Sth
353    CSTth = 1.0/STth
354    SCth = 1.0/Cth
355    SCTth = 1.0/CTth
356    if 'Bragg' in Geometry:
357        G = 8.0*Delt*Navo*emc*Sth/((1.0-CTth**2)*(1.0-np.exp(-2.0*MuT*CSth)))
358        Y1 = np.pi
359        Y2 = np.pi/2.0
360        Y3 = 3.*np.pi/8. #3pi/4?
361        W = 1.0/(Sth+np.fabs(Sgama))
362        W += np.exp(-MuT*CSth)*(2.0*np.fabs(Sgama)*np.exp(-MuT*np.fabs(CSgama))-
363            (Sth+np.fabs(Sgama))*np.exp(-MuT*CSth))/(Sth**2-Sgama**2)
364        Fac0 = Sth**2*Sgama**2
365        X = Fac0+(Fac0+CTth)**2/2
366        Y = Cgama**2*Cth**2*(1.0-Fac0-CTth)
367        Z = Cgama**4*Cth**4/2.0
368        E = 2.0*(1.0-a)/(b*Cgama/Cth)
369        F1 = (2.0+b*(1.0+Sth*Sgama))/(b*Cth*Cgama) #trouble if < 1
370        F2 = (2.0+b*(1.0-Sth*Sgama))/(b*Cth*Cgama)
371        T1 = np.pi/np.sqrt(F1**2-1.0)
372        T2 = np.pi/np.sqrt(F2**2-1.0)
373        Y4 = T1+T2
374        Y5 = F1**2*T1+F2**2*T2-np.pi*(F1+F2)
375        Y6 = F1**4*T1+F2**4*T2-np.pi*(F1+F2)/2.0-np.pi*(F1**3+F2**3)
376        Y7 = (T2-T1)/(F1-F2)
377        YT = F2**2*T2-F1**2*T1
378        Y8 = Y1+YT/(F1-F2)
379        Y9 = Y2+(F2**4*T2-F1**4*T1)/(F1-F2)+Y1*((F1+F2)**2-F1*F2)
380        M = (a**2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E**2*(X*Y7+Y*Y8+Z*Y9))*Cgama
381       
382        Q = np.sum(W*M,axis=0)
383        return Q*G*NS/(NS-1.)
384#
385#      cos2th=2.0d*costh^2 - 1.0d
386#      G= delta * 8.0d * Na * emc * sinth/(1.0d + cos2th^2)/(1.0d - exp(-2.0d*mut*cscth))
387#      Y1=3.1415926d
388#      Y2=Y1*0.5d
389#      Y3=Y2*0.75d
390#      for i=1,num_steps-1 do begin
391#         cosgama=double(cos(gama[i]))
392#         singama=double(sin(gama[i]))
393#         cscgama=1.0d / singama
394#
395#         W=1.0d/(sinth+abs(singama))
396#         W=W+exp(-1.0*mut*cscth)*(2.0d*abs(singama)*exp(-1.0d*mut*abs(cscgama))- $
397#                 (sinth+abs(singama))*exp(-1.0d*mut*cscth))/(sinth^2-singama^2)
398#
399#         factor0=sinth^2*singama^2
400#         X=factor0+(factor0+cos2th)^2/2.0d
401#         Y=cosgama^2*(1.0d - factor0-cos2th)*costh^2
402#         Z=cosgama^4/2.0d*costh^4
403#         E=2.0d*(1.0-a)/b/cosgama/costh
404#
405#         F1=1.0d/b/cosgama*(2.0d + b*(1.0+sinth*singama))/costh
406#         F2=1.0d/b/cosgama*(2.0d + b*(1.0-sinth*singama))/costh
407#
408#         T1=3.14159/sqrt(F1^2-1.0d)
409#         T2=3.14159/sqrt(F2^2-1.0d)
410#         Y4=T1+T2
411#         Y5=F1^2*T1+F2^2*T2-3.14159*(F1+F2)
412#         Y6=F1^4*T1+F2^4*T2-3.14159*(F1+F2)/2.0-3.14159*(F1^3+F2^3)
413#         Y7=(T2-T1)/(F1-F2)
414#         Y8=Y1+(F2^2*T2-F1^2*T1)/(F1-F2)
415#         Y9=Y2+(F2^4*T2-F1^4*T1)/(F1-F2)+Y1*((F1+F2)^2-F1*F2)
416#         M=(a^2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E^2* $
417#                      (X*Y7+Y*Y8+Z*Y9))*cosgama
418#
419#         Q=Q+W*M
420#
421#      endfor
422#      Q=double(num_steps)/Double(num_steps-1)*Q*G
423#      end
424    elif 'Cylinder' in Geometry:
425        Q = 0.
426    elif 'Fixed' in Geometry:   #Dwiggens & Park, Acta Cryst. A27, 264 (1971) with corrections
427        EMA = np.exp(-MuT*(1.0-SCTth))
428        Fac1 = (1.-EMA)/(1.0-SCTth)
429        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*Fac1)
430        Fac0 = Cgama/(1-Sgama*SCTth)
431        Wp = Fac0*(Fac1-(EMA-np.exp(-MuT*(CSgama-SCTth)))/(CSgama-1.0))
432        Fac0 = Cgama/(1.0+Sgama*SCTth)
433        Wm = Fac0*(Fac1+(np.exp(-MuT*(1.0+CSgama))-1.0)/(CSgama+1.0))
434        X = (Sgama**2+CTth**2*(1.0-Sgama**2+Sgama**4))/2.0
435        Y = Sgama**3*Cgama*STth*CTth
436        Z = Cgama**2*(1.0+Sgama**2)*STth**2/2.0
437        Fac2 = 1.0/(b*Cgama*STth)
438        U = 2.0*(1.0-a)*Fac2
439        V = (2.0+b*(1.0-CTth*Sgama))*Fac2
440        Mp = 2.0*np.pi*(a+2.0*(1.0-a)/(2.0+b*(1.0-Sgama)))*(a*X+a*Z/2.0-U*Y+U*(X+Y*V+Z*V**2)/np.sqrt(V**2-1.0)-U*Z*V)
441        V = (2.0+b*(1.0+CTth*Sgama))*Fac2
442        Y = -Y
443        Mm = 2.0*np.pi*(a+2.0*(1.0-a)/(2.0+b*(1.0+Sgama)))*(a*X+a*Z/2.0-U*Y+U*(X+Y*V+Z*V**2)/np.sqrt(V**2-1.0)-U*Z*V)
444        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
445        return Q*G*NS/(NS-1.)
446    elif 'Tilting' in Geometry:
447        EMA = np.exp(-MuT*SCth)
448        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*EMA)
449#          Wplus = expmutsecth/(1.0d - singama*secth) + singama/mut/(1.0 -singama*secth)/(1.0-singama*secth)* $
450#                                                       (Exp(-1.0*mut*cscgama) - expmutsecth)
451#          Wminus = expmutsecth/(1.0d + singama*secth) - singama/mut/(1.0 +singama*secth)/(1.0+singama*secth)* $
452#                                                        expmutsecth*(1.0d - expmutsecth*Exp(-1.0*mut*cscgama))
453        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)/(1.0-Sgama*SCth)*(np.exp(-MuT*CSgama)-EMA)
454#        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)**2*(np.exp(-MuT*CSgama)-EMA)
455        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)/(1.0+Sgama*SCth)*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
456#        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)**2*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
457        X = 0.5*(Cth**2*(Cth**2*Sgama**4-4.0*Sth**2*Cgama**2)+1.0)
458        Y = Cgama**2*(1.0+Cgama**2)*Cth**2*Sth**2
459        Z = 0.5*Cgama**4*Sth**4
460#          X = 0.5*(costh*costh*(costh*costh*singama*singama*singama*singama - $
461#                           4.0*sinth*sinth*cosgama*cosgama) +1.0d)
462#
463#          Y = cosgama*cosgama*(1.0 + cosgama*cosgama)*costh*costh*sinth*sinth
464#
465#          Z= 0.5 *cosgama*cosgama*cosgama*cosgama* (sinth^4)
466#
467        AlP = 2.0+b*(1.0-Cth*Sgama)
468        AlM = 2.0+b*(1.0+Cth*Sgama)
469#          alphaplus = 2.0 + b*(1.0 - costh*singama)
470#          alphaminus = 2.0 + b*(1.0 + costh*singama)
471        BeP = np.sqrt(np.fabs(AlP**2-(b*Cgama*Sth)**2))
472        BeM = np.sqrt(np.fabs(AlM**2-(b*Cgama*Sth)**2))
473#          betaplus = Sqrt(Abs(alphaplus*alphaplus - b*b*cosgama*cosgama*sinth*sinth))
474#          betaminus = Sqrt(Abs(alphaminus*alphaminus - b*b*cosgama*cosgama*sinth*sinth))
475        Mp = Cgama*(np.pi*a**2*(2.0*X+Y+0.75*Z)+(2.0*np.pi*(1.0-a))*(1.0-a+a*AlP)* \
476            (4.0*X/AlP/BeP+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlP/BeP-1.0)+
477            2.0/b**4*AlP/BeP*AlP**2-2.0/b**4*AlP**2-Cgama**2/b/b*Sth*2))
478#          Mplus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
479#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaplus)* $
480#                   (4.0*x/alphaplus/betaplus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaplus/betaplus -1.0) + $
481#                   2.0/(b^4)*alphaplus/betaplus*alphaplus*alphaplus - 2.0/(b^4)*alphaplus*alphaplus - $
482#                   cosgama*cosgama/b/b*sinth*sinth))
483        Mm =Cgama*(np.pi*a**2*(2.0*X+Y+0.75*Z)+(2.0*np.pi*(1.0-a))*(1.0-a+a*AlM)* \
484            (4.0*X/AlM/BeM+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlM/BeM-1.0)+
485            2.0/b**4*AlM/BeM*AlM**2-2.0/b**4*AlM**2-Cgama**2/b/b*Sth*2))
486#          Mminus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
487#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaminus)* $
488#                   (4.0*x/alphaminus/betaminus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaminus/betaminus -1.0) + $
489#                   2.0/(b^4)*alphaminus/betaminus*alphaminus*alphaminus - 2.0/(b^4)*alphaminus*alphaminus - $
490#                   cosgama*cosgama/b/b*sinth*sinth))
491        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
492        return Q*G*NS/(NS-1.)
493#       expmutsecth = Exp(-1.0*mut*secth)
494#       G= delta * 2.0 * Na * emc /(1.0+costth^2)/expmutsecth
495#       for i=1, num_steps-1 do begin
496#          cosgama=double(cos(gama[i]))
497#          singama=double(sin(gama[i]))
498#          cscgama=1.0d/singama
499#
500#
501#     ; print, "W", min(wplus), max(wplus), min(wminus), max(wminus)
502#
503#
504#
505#
506#    ;               print, a,b
507#  ; print, "M", min(mplus), max(mplus), min(mminus), max(mminus)
508#          Q=Q+ Wplus*Mplus + Wminus*Mminus
509#      endfor
510#      Q=double(num_steps)/double(num_steps-1)*Q*G
511#   ;   print, min(q), max(q)
512#      end
513
514def MultiScattering(Geometry,ElList,Tth):
515    BN = BD = 0.0
516    Amu = 0.0
517    for El in ElList:
518        el = ElList[El]
519        BN += el['Z']*el['FormulaNo']
520        BD += el['FormulaNo']
521        Amu += el['FormulaNo']*el['mu']
522       
523
524def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
525    # returns value(esd) string; nTZ=True for no trailing zeros
526    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
527    #get the 2 significant digits in the esd
528    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
529    #get the number of digits to represent them
530    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
531   
532    mdec = lambda esd: -int(math.log10(abs(esd)))
533    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
534    if esd > 0:
535        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
536        print fmt,ndec(esd),esd*10**(mdec(esd)+1)
537        return fmt%(value,int(esd*10**(mdec(esd)+1)))
538    elif esd < 0:
539         return str(round(value,mdec(esd)))
540    else:
541        text = "%F"%(value)
542        if nTZ:
543            return text.rstrip('0')
544        else:
545            return text
546
547       
548#GSASII peak fitting routine: Thompson, Cox & Hastings; Finger, Cox & Jephcoat model       
549
550def DoPeakFit(peaks,background,limits,inst,data):
551   
552    def backgroundPrint(background,sigback):
553        if background[1]:
554            print 'Background coefficients for',background[0],'function'
555            ptfmt = "%12.5f"
556            ptstr =  'values:'
557            sigstr = 'esds  :'
558            for i,back in enumerate(background[3:]):
559                ptstr += ptfmt % (back)
560                sigstr += ptfmt % (sigback[i+3])
561            print ptstr
562            print sigstr
563        else:
564            print 'Background not refined'
565   
566    def instPrint(instVal,siginst,insLabels):
567        print 'Instrument Parameters:'
568        ptfmt = "%12.6f"
569        ptlbls = 'names :'
570        ptstr =  'values:'
571        sigstr = 'esds  :'
572        for i,value in enumerate(instVal):
573            ptlbls += "%s" % (insLabels[i].center(12))
574            ptstr += ptfmt % (value)
575            if siginst[i]:
576                sigstr += ptfmt % (siginst[i])
577            else:
578                sigstr += 12*' '
579        print ptlbls
580        print ptstr
581        print sigstr
582   
583    def peaksPrint(peaks,sigpeaks):
584        print 'Peak list='
585
586    begin = time.time()
587    Np = len(peaks[0])
588    DataType = inst[1][0]
589    instVal = inst[1][1:]
590    Insref = inst[2][1:]
591    insLabels = inst[3][1:]
592    Ka2 = False
593    Ioff = 3
594    if len(instVal) == 12:
595        lamratio = instVal[1]/instVal[0]
596        Ka2 = True
597        Ioff = 5
598    insref = Insref[len(Insref)-7:-1]               #just U,V,W,X,Y,SH/L
599    for peak in peaks:
600        dip = []
601        dip.append(tand(peak[0]/2.0)**2)
602        dip.append(tand(peak[0]/2.0))
603        dip.append(1.0/cosd(peak[0]/2.0))
604        dip.append(tand(peak[0]/2.0))
605        peak.append(dip)
606    B = background[2]
607    bcof = background[3:3+B]
608    Bv = 0
609    if background[1]:
610        Bv = B
611    x,y,w,yc,yb,yd = data               #these are numpy arrays!
612    V = []
613    A = []
614    swobs = 0.0
615    smin = 0.0
616    first = True
617    GoOn = True
618    Go = True
619    dlg = wx.ProgressDialog("Elapsed time","Fitting peaks to pattern",len(x), \
620        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
621    screenSize = wx.DisplaySize()
622    Size = dlg.GetSize()
623    dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0))
624    try:
625        i = 0
626        for xi in x :
627            Go = dlg.Update(i)[0]
628            if GoOn:
629                GoOn = Go
630            if limits[0] <= xi <= limits[1]:
631                yb[i] = 0.0
632                dp = []
633                for j in range(B):
634                    t = (xi-limits[0])**j
635                    yb[i] += t*bcof[j]
636                    if background[1]:
637                        dp.append(t)
638                yc[i] = yb[i]
639                Iv = 0
640                for j in range(6):
641                    if insref[j]:
642                        dp.append(0.0)
643                        Iv += 1
644                for peak in peaks:
645                    dip = peak[-1]
646                    f = pyp.pypsvfcj(peak[2],xi-peak[0],peak[0],peak[4],peak[6],instVal[-2],0.0)
647                    yc[i] += f[0]*peak[2]
648                    if f[0] > 0.0:
649                        j = 0
650                        if insref[0]:              #U
651                            dp[Bv+j] += f[3]*dip[0]
652                            j += 1
653                        if insref[1]:              #V
654                            dp[Bv+j] += f[3]*dip[1]
655                            j += 1
656                        if insref[2]:              #W
657                            dp[Bv+j] += f[3]
658                            j += 1
659                        if insref[3]:              #X
660                            dp[Bv+j] += f[4]*dip[2]
661                            j += 1
662                        if insref[4]:              #Y
663                            dp[Bv+j] += f[4]*dip[3]
664                            j += 1
665                        if insref[5]:              #SH/L
666                            dp[Bv+j] += f[5]
667                    if Ka2:
668                       pos2 = 2.0*asind(lamratio*sind(peak[0]/2.0))
669                       f2 = pyp.pypsvfcj(peak[2],xi-pos2,peak[0],peak[4],peak[6],instVal[-2],0.0)
670                       yc[i] += f2[0]*peak[2]*instVal[3]
671                       if f[0] > 0.0:
672                           j = 0
673                           if insref[0]:              #U
674                               dp[Bv+j] += f2[3]*dip[0]*instVal[3]
675                               j += 1
676                           if insref[1]:              #V
677                               dp[Bv+j] += f2[3]*dip[1]*instVal[3]
678                               j += 1
679                           if insref[2]:              #W
680                               dp[Bv+j] += f2[3]*instVal[3]
681                               j += 1
682                           if insref[3]:              #X
683                               dp[Bv+j] += f2[4]*dip[2]*instVal[3]
684                               j += 1
685                           if insref[4]:              #Y
686                               dp[Bv+j] += f2[4]*dip[3]*instVal[3]
687                               j += 1
688                           if insref[5]:              #SH/L
689                               dp[Bv+j] += f2[5]*instVal[3]                       
690                    for j in range(0,Np,2):
691                        if peak[j+1]: dp.append(f[j/2+1])
692                yd[i] = y[i]-yc[i]
693                swobs += w[i]*y[i]**2
694                t2 = w[i]*yd[i]
695                smin += t2*yd[i]
696                if first:
697                    first = False
698                    M = len(dp)
699                    A = np.zeros(shape=(M,M))
700                    V = np.zeros(shape=(M))
701                A,V = pyp.buildmv(t2,w[i],M,dp,A,V)
702            i += 1
703    finally:
704        dlg.Destroy()
705    Rwp = smin/swobs
706    Rwp = math.sqrt(Rwp)*100.0
707    norm = np.diag(A)
708    for i,elm in enumerate(norm):
709        if elm <= 0.0:
710            print norm
711            return False,0,0,0,False
712    for i in xrange(len(V)):
713        norm[i] = 1.0/math.sqrt(norm[i])
714        V[i] *= norm[i]
715        a = A[i]
716        for j in xrange(len(V)):
717            a[j] *= norm[i]
718    A = np.transpose(A)
719    for i in xrange(len(V)):
720        a = A[i]
721        for j in xrange(len(V)):
722            a[j] *= norm[i]
723    b = nl.solve(A,V)
724    A = nl.inv(A)
725    sig = np.diag(A)
726    for i in xrange(len(V)):
727        b[i] *= norm[i]
728        sig[i] *= norm[i]*norm[i]
729        sig[i] = math.sqrt(abs(sig[i]))
730    sigback = [0,0,0]
731    for j in range(Bv):
732        background[j+3] += b[j]
733        sigback.append(sig[j])
734    backgroundPrint(background,sigback)
735    k = 0
736    delt = []
737    if Ka2:
738        siginst = [0,0,0,0,0]
739    else:
740        siginst = [0,0,0]
741    for j in range(6):
742        if insref[j]:
743            instVal[j+Ioff] += b[Bv+k]
744            siginst.append(sig[Bv+k])
745            delt.append(b[Bv+k])
746            k += 1
747        else:
748            delt.append(0.0)
749            siginst.append(0.0)
750    delt.append(0.0)                    #dummies for azm
751    siginst.append(0.0)
752    instPrint(instVal,siginst,insLabels)
753    inst[1] = [DataType,]
754    for val in instVal:
755        inst[1].append(val)
756    B = Bv+Iv
757    for peak in peaks:
758        del peak[-1]                        # remove dip from end
759        delsig = delt[0]*tand(peak[0]/2.0)**2+delt[1]*tand(peak[0]/2.0)+delt[2]
760        delgam = delt[3]/cosd(peak[0]/2.0)+delt[4]*tand(peak[0]/2.0)
761        for j in range(0,len(peak[:-1]),2):
762            if peak[j+1]: 
763                peak[j] += b[B]
764                B += 1
765        peak[4] += delsig
766        if peak[4] < 0.0:
767            print 'ERROR - negative sigma'
768            return False,0,0,0,False           
769        peak[6] += delgam
770        if peak[6] < 0.0:
771            print 'ERROR - negative gamma'
772            return False,0,0,0,False
773    runtime = time.time()-begin   
774    data = [x,y,w,yc,yb,yd]
775    return True,smin,Rwp,runtime,GoOn
776
777def CalcPDF(data,inst,xydata):
778    auxPlot = []
779    import copy
780    import scipy.fftpack as ft
781    #subtract backgrounds - if any
782    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
783    if data['Sample Bkg.']['Name']:
784        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
785            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
786    if data['Container']['Name']:
787        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
788        if data['Container Bkg.']['Name']:
789            xycontainer += (xydata['Container Bkg.'][1][1]+
790                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
791        xydata['IofQ'][1][1] += xycontainer
792    #get element data & absorption coeff.
793    ElList = data['ElList']
794    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
795    #Apply angle dependent corrections
796    Tth = xydata['Sample'][1][0]
797    dt = (Tth[1]-Tth[0])
798    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth)
799    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])
800    if data['DetType'] == 'Image plate':
801        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
802    XY = xydata['IofQ'][1]   
803    #convert to Q
804    hc = 12.397639
805    if 'Lam' in inst:
806        wave = inst['Lam']
807    else:
808        wave = inst['Lam1']
809    keV = hc/wave
810    minQ = npT2q(Tth[0],wave)
811    maxQ = npT2q(Tth[-1],wave)   
812    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
813    dq = Qpoints[1]-Qpoints[0]
814    XY[0] = npT2q(XY[0],wave)   
815#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
816    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
817    Qdata = T(Qpoints)
818   
819    qLimits = data['QScaleLim']
820    minQ = np.searchsorted(Qpoints,qLimits[0])
821    maxQ = np.searchsorted(Qpoints,qLimits[1])
822    newdata = []
823    xydata['IofQ'][1][0] = Qpoints
824    xydata['IofQ'][1][1] = Qdata
825    for item in xydata['IofQ'][1]:
826        newdata.append(item[:maxQ])
827    xydata['IofQ'][1] = newdata
828   
829
830    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
831    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
832    Q = xydata['SofQ'][1][0]
833    ruland = Ruland(data['Ruland'],wave,Q,CF)
834#    auxPlot.append([Q,ruland,'Ruland'])     
835    CF *= ruland
836#    auxPlot.append([Q,CF,'CF-Corr'])
837    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
838    xydata['SofQ'][1][1] *= scale
839    xydata['SofQ'][1][1] -= CF
840    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
841    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
842    xydata['SofQ'][1][1] *= scale
843   
844    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
845    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
846    if data['Lorch']:
847        xydata['FofQ'][1][1] *= LorchWeight(Q)
848   
849    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
850    nR = len(xydata['GofR'][1][1])
851    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
852    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
853   
854       
855    return auxPlot
856       
Note: See TracBrowser for help on using the repository browser.