source: trunk/GSASIIpwd.py @ 2614

Last change on this file since 2614 was 2614, checked in by toby, 5 years ago

add PDF correction optimizer; fix vdWaals bug after adding element in structure plot

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 88.9 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASII powder calculation module*
5==================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2017-01-06 19:04:03 +0000 (Fri, 06 Jan 2017) $
10# $Author: toby $
11# $Revision: 2614 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 2614 2017-01-06 19:04:03Z toby $
14########### SVN repository information ###################
15import sys
16import math
17import time
18import os
19import subprocess as subp
20
21import numpy as np
22import numpy.linalg as nl
23import numpy.ma as ma
24import random as rand
25from numpy.fft import ifft, fft, fftshift
26import scipy.interpolate as si
27import scipy.stats as st
28import scipy.optimize as so
29
30import GSASIIpath
31GSASIIpath.SetVersionNumber("$Revision: 2614 $")
32import GSASIIlattice as G2lat
33import GSASIIspc as G2spc
34import GSASIIElem as G2elem
35import GSASIImath as G2mth
36import pypowder as pyd
37try:
38    import pydiffax as pyx
39except ImportError:
40    print 'pydiffax is not available for this platform - under develpment'
41
42   
43# trig functions in degrees
44tand = lambda x: math.tan(x*math.pi/180.)
45atand = lambda x: 180.*math.atan(x)/math.pi
46atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
47cosd = lambda x: math.cos(x*math.pi/180.)
48acosd = lambda x: 180.*math.acos(x)/math.pi
49rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
50#numpy versions
51npsind = lambda x: np.sin(x*np.pi/180.)
52npasind = lambda x: 180.*np.arcsin(x)/math.pi
53npcosd = lambda x: np.cos(x*math.pi/180.)
54npacosd = lambda x: 180.*np.arccos(x)/math.pi
55nptand = lambda x: np.tan(x*math.pi/180.)
56npatand = lambda x: 180.*np.arctan(x)/np.pi
57npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
58npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave    #=d*
59npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)    #=2pi*d*
60ateln2 = 8.0*math.log(2.0)
61
62################################################################################
63#### Powder utilities
64################################################################################
65
66def PhaseWtSum(G2frame,histo):
67    '''
68    Calculate sum of phase mass*shase fraction for PWDR data (exclude magnetic phases)
69   
70    :param G2frame: GSASII main frame structure
71    :param str histo: histogram name
72    :returns sum(scale*mass) for phases in histo
73    '''
74    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
75    wtSum = 0.0
76    for phase in Phases:
77        if Phases[phase]['General']['Type'] != 'magnetic':
78            if histo in Phases[phase]['Histograms']:
79                mass = Phases[phase]['General']['Mass']
80                phFr = Phases[phase]['Histograms'][histo]['Scale'][0]
81                wtSum += mass*phFr
82    return wtSum
83   
84################################################################################
85#GSASII pdf calculation routines
86################################################################################
87       
88def Transmission(Geometry,Abs,Diam):
89    '''
90    Calculate sample transmission
91
92    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
93    :param float Abs: absorption coeff in cm-1
94    :param float Diam: sample thickness/diameter in mm
95    '''
96    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
97        MuR = Abs*Diam/20.0
98        if MuR <= 3.0:
99            T0 = 16/(3.*math.pi)
100            T1 = -0.045780
101            T2 = -0.02489
102            T3 = 0.003045
103            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
104            if T < -20.:
105                return 2.06e-9
106            else:
107                return math.exp(T)
108        else:
109            T1 = 1.433902
110            T2 = 0.013869+0.337894
111            T3 = 1.933433+1.163198
112            T4 = 0.044365-0.04259
113            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
114            return T/100.
115    elif 'plate' in Geometry:
116        MuR = Abs*Diam/10.
117        return math.exp(-MuR)
118    elif 'Bragg' in Geometry:
119        return 0.0
120       
121def SurfaceRough(SRA,SRB,Tth):
122    ''' Suortti (J. Appl. Cryst, 5,325-331, 1972) surface roughness correction
123    :param float SRA: Suortti surface roughness parameter
124    :param float SRB: Suortti surface roughness parameter
125    :param float Tth: 2-theta(deg) - can be numpy array
126   
127    '''
128    sth = npsind(Tth/2.)
129    T1 = np.exp(-SRB/sth)
130    T2 = SRA+(1.-SRA)*np.exp(-SRB)
131    return (SRA+(1.-SRA)*T1)/T2
132   
133def SurfaceRoughDerv(SRA,SRB,Tth):
134    ''' Suortti surface roughness correction derivatives
135    :param float SRA: Suortti surface roughness parameter (dimensionless)
136    :param float SRB: Suortti surface roughness parameter (dimensionless)
137    :param float Tth: 2-theta(deg) - can be numpy array
138    :return list: [dydSRA,dydSRB] derivatives to be used for intensity derivative
139    '''
140    sth = npsind(Tth/2.)
141    T1 = np.exp(-SRB/sth)
142    T2 = SRA+(1.-SRA)*np.exp(-SRB)
143    Trans = (SRA+(1.-SRA)*T1)/T2
144    dydSRA = ((1.-T1)*T2-(1.-np.exp(-SRB))*Trans)/T2**2
145    dydSRB = ((SRA-1.)*T1*T2/sth-Trans*(SRA-T2))/T2**2
146    return [dydSRA,dydSRB]
147
148def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
149    '''Calculate sample absorption
150    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
151    :param float MuR: absorption coeff * sample thickness/2 or radius
152    :param Tth: 2-theta scattering angle - can be numpy array
153    :param float Phi: flat plate tilt angle - future
154    :param float Psi: flat plate tilt axis - future
155    '''
156   
157    def muRunder3(MuR,Sth2):
158        T0 = 16.0/(3.*np.pi)
159        T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
160            0.109561*np.sqrt(Sth2)-26.04556
161        T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
162            1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
163        T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
164        Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
165        return np.exp(Trns)
166   
167    def muRover3(MuR,Sth2):
168        T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
169            10.02088*Sth2**3-3.36778*Sth2**4
170        T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
171            (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
172        T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
173            0.13576*np.sqrt(Sth2)+1.163198
174        T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
175        Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
176        return Trns/100.
177       
178    Sth2 = npsind(Tth/2.0)**2
179    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
180        if 'array' in str(type(MuR)):
181            MuRSTh2 = np.concatenate((MuR,Sth2))
182            AbsCr = np.where(MuRSTh2[0]<=3.0,muRunder3(MuRSTh2[0],MuRSTh2[1]),muRover3(MuRSTh2[0],MuRSTh2[1]))
183            return AbsCr
184        else:
185            if MuR <= 3.0:
186                return muRunder3(MuR,Sth2)
187            else:
188                return muRover3(MuR,Sth2)
189    elif 'Bragg' in Geometry:
190        return 1.0
191    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
192        # and only defined for 2theta < 90
193        MuT = 2.*MuR
194        T1 = np.exp(-MuT)
195        T2 = np.exp(-MuT/npcosd(Tth))
196        Tb = MuT-MuT/npcosd(Tth)
197        return (T2-T1)/Tb
198    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
199        MuT = 2.*MuR
200        cth = npcosd(Tth/2.0)
201        return np.exp(-MuT/cth)/cth
202       
203def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
204    'needs a doc string'
205    dA = 0.001
206    AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
207    if MuR:
208        AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
209        return (AbsP-AbsM)/(2.0*dA)
210    else:
211        return (AbsP-1.)/dA
212       
213def Polarization(Pola,Tth,Azm=0.0):
214    """   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
215
216    :param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
217    :param Azm: azimuthal angle e.g. 0.0 in plane of polarization
218    :param Tth: 2-theta scattering angle - can be numpy array
219      which (if either) of these is "right"?
220    :return: (pola, dpdPola)
221      * pola = ((1-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \
222        (1-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
223      * dpdPola: derivative needed for least squares
224
225    """
226    pola = ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
227        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
228    dpdPola = -npsind(Tth)**2*(npsind(Azm)**2-npcosd(Azm)**2)
229    return pola,dpdPola
230   
231def Oblique(ObCoeff,Tth):
232    'currently assumes detector is normal to beam'
233    if ObCoeff:
234        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
235    else:
236        return 1.0
237               
238def Ruland(RulCoff,wave,Q,Compton):
239    'needs a doc string'
240    C = 2.9978e8
241    D = 1.5e-3
242    hmc = 0.024262734687
243    sinth2 = (Q*wave/(4.0*np.pi))**2
244    dlam = (wave**2)*Compton*Q/C
245    dlam_c = 2.0*hmc*sinth2-D*wave**2
246    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
247   
248def LorchWeight(Q):
249    'needs a doc string'
250    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
251           
252def GetAsfMean(ElList,Sthl2):
253    '''Calculate various scattering factor terms for PDF calcs
254
255    :param dict ElList: element dictionary contains scattering factor coefficients, etc.
256    :param np.array Sthl2: numpy array of sin theta/lambda squared values
257    :returns: mean(f^2), mean(f)^2, mean(compton)
258    '''
259    sumNoAtoms = 0.0
260    FF = np.zeros_like(Sthl2)
261    FF2 = np.zeros_like(Sthl2)
262    CF = np.zeros_like(Sthl2)
263    for El in ElList:
264        sumNoAtoms += ElList[El]['FormulaNo']
265    for El in ElList:
266        el = ElList[El]
267        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
268        cf = G2elem.ComptonFac(el,Sthl2)
269        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
270        FF2 += ff2*el['FormulaNo']/sumNoAtoms
271        CF += cf*el['FormulaNo']/sumNoAtoms
272    return FF2,FF**2,CF
273   
274def GetNumDensity(ElList,Vol):
275    'needs a doc string'
276    sumNoAtoms = 0.0
277    for El in ElList:
278        sumNoAtoms += ElList[El]['FormulaNo']
279    return sumNoAtoms/Vol
280           
281def CalcPDF(data,inst,limits,xydata):
282    '''Computes I(Q), S(Q) & G(r) from Sample, Bkg, etc. diffraction patterns loaded into
283    dict xydata; results are placed in xydata.
284    Calculation parameters are found in dicts data and inst and list limits.
285    The return value is at present an empty list.
286    '''
287    auxPlot = []
288    import copy
289    import scipy.fftpack as ft
290    Ibeg = np.searchsorted(xydata['Sample'][1][0],limits[0])
291    Ifin = np.searchsorted(xydata['Sample'][1][0],limits[1])+1
292    #subtract backgrounds - if any & use PWDR limits
293#    GSASIIpath.IPyBreak()
294    IofQ = copy.deepcopy(xydata['Sample'])
295    IofQ[1] = np.array(IofQ[1])[:,Ibeg:Ifin]
296    if data['Sample Bkg.']['Name']:
297        IofQ[1][1] += (xydata['Sample Bkg.'][1][1][Ibeg:Ifin]+
298            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
299    if data['Container']['Name']:
300        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
301        if data['Container Bkg.']['Name']:
302            xycontainer += (xydata['Container Bkg.'][1][1][Ibeg:Ifin]+
303                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
304        IofQ[1][1] += xycontainer[Ibeg:Ifin]
305    data['IofQmin'] = IofQ[1][1][Ifin-1]
306    IofQ[1][1] -= data.get('Flat Bkg',0.)
307    #get element data & absorption coeff.
308    ElList = data['ElList']
309    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
310    #Apply angle dependent corrections
311    Tth = IofQ[1][0]
312    MuR = Abs*data['Diam']/20.0
313    IofQ[1][1] /= Absorb(data['Geometry'],MuR,Tth)
314    if 'X' in inst['Type'][0]:
315        IofQ[1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
316    if data['DetType'] == 'Image plate':
317        IofQ[1][1] *= Oblique(data['ObliqCoeff'],Tth)
318    XY = IofQ[1]   
319    #convert to Q
320#    nQpoints = len(XY[0])     #points for Q interpolation
321    nQpoints = 5000
322    if 'C' in inst['Type'][0]:
323        wave = G2mth.getWave(inst)
324        minQ = npT2q(Tth[0],wave)
325        maxQ = npT2q(Tth[-1],wave)   
326        Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
327        dq = Qpoints[1]-Qpoints[0]
328        XY[0] = npT2q(XY[0],wave)
329    elif 'T' in inst['Type'][0]:
330        difC = inst['difC'][1]
331        minQ = 2.*np.pi*difC/Tth[-1]
332        maxQ = 2.*np.pi*difC/Tth[0]
333        Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
334        dq = Qpoints[1]-Qpoints[0]
335        XY[0] = 2.*np.pi*difC/XY[0]
336    Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0])    #interpolate I(Q)
337    Qdata -= np.min(Qdata)*data['BackRatio']
338   
339    qLimits = data['QScaleLim']
340    minQ = np.searchsorted(Qpoints,qLimits[0])
341    maxQ = np.searchsorted(Qpoints,qLimits[1])
342    newdata = []
343    xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],IofQ[2]]
344    for item in xydata['IofQ'][1]:
345        newdata.append(item[:maxQ])
346    xydata['IofQ'][1] = newdata
347   
348    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
349    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
350    Q = xydata['SofQ'][1][0]
351#    auxPlot.append([Q,np.copy(CF),'CF-unCorr'])
352    ruland = Ruland(data['Ruland'],wave,Q,CF)
353#    auxPlot.append([Q,ruland,'Ruland'])     
354    CF *= ruland
355#    auxPlot.append([Q,CF,'CF-Corr'])
356    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
357    xydata['SofQ'][1][1] *= scale
358    xydata['SofQ'][1][1] -= CF
359    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
360    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
361    xydata['SofQ'][1][1] *= scale
362    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
363    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
364    if data['Lorch']:
365        xydata['FofQ'][1][1] *= LorchWeight(Q)   
366    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
367    nR = len(xydata['GofR'][1][1])
368    mul = int(round(2.*np.pi*nR/(data.get('Rmax',100.)*qLimits[1])))
369    xydata['GofR'][1][0] = 2.*np.pi*np.linspace(0,nR,nR)/(mul*qLimits[1])
370    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],mul*nR)[:nR])
371    if data.get('noRing',True):
372        xydata['GofR'][1][1] = np.where(xydata['GofR'][1][0]<0.5,0.,xydata['GofR'][1][1])
373    return auxPlot
374   
375def MakeRDF(RDFcontrols,background,inst,pwddata):
376    import scipy.fftpack as ft
377    import scipy.signal as signal
378    auxPlot = []
379    if 'C' in inst['Type'][0]:
380        Tth = pwddata[0]
381        wave = G2mth.getWave(inst)
382        minQ = npT2q(Tth[0],wave)
383        maxQ = npT2q(Tth[-1],wave)
384        powQ = npT2q(Tth,wave) 
385    elif 'T' in inst['Type'][0]:
386        TOF = pwddata[0]
387        difC = inst['difC'][1]
388        minQ = 2.*np.pi*difC/TOF[-1]
389        maxQ = 2.*np.pi*difC/TOF[0]
390        powQ = 2.*np.pi*difC/TOF
391    piDQ = np.pi/(maxQ-minQ)
392    Qpoints = np.linspace(minQ,maxQ,len(pwddata[0]),endpoint=True)
393    if RDFcontrols['UseObsCalc']:
394        Qdata = si.griddata(powQ,pwddata[1]-pwddata[3],Qpoints,method=RDFcontrols['Smooth'],fill_value=0.)
395    else:
396        Qdata = si.griddata(powQ,pwddata[1]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
397    Qdata *= np.sin((Qpoints-minQ)*piDQ)/piDQ
398    Qdata *= 0.5*np.sqrt(Qpoints)       #Qbin normalization
399#    GSASIIpath.IPyBreak()
400    dq = Qpoints[1]-Qpoints[0]
401    nR = len(Qdata)
402    R = 0.5*np.pi*np.linspace(0,nR,nR)/(4.*maxQ)
403    iFin = np.searchsorted(R,RDFcontrols['maxR'])
404    bBut,aBut = signal.butter(4,0.01)
405    Qsmooth = signal.filtfilt(bBut,aBut,Qdata)
406#    auxPlot.append([Qpoints,Qdata,'interpolate:'+RDFcontrols['Smooth']])
407#    auxPlot.append([Qpoints,Qsmooth,'interpolate:'+RDFcontrols['Smooth']])
408    DofR = dq*np.imag(ft.fft(Qsmooth,16*nR)[:nR])
409    auxPlot.append([R[:iFin],DofR[:iFin],'D(R)'])   
410    return auxPlot
411
412################################################################################       
413#GSASII peak fitting routines: Finger, Cox & Jephcoat model       
414################################################################################
415
416def factorize(num):
417    ''' Provide prime number factors for integer num
418    :returns: dictionary of prime factors (keys) & power for each (data)
419    '''
420    factors = {}
421    orig = num
422
423    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
424    i, sqi = 2, 4
425    while sqi <= num:
426        while not num%i:
427            num /= i
428            factors[i] = factors.get(i, 0) + 1
429
430        sqi += 2*i + 1
431        i += 1
432
433    if num != 1 and num != orig:
434        factors[num] = factors.get(num, 0) + 1
435
436    if factors:
437        return factors
438    else:
439        return {num:1}          #a prime number!
440           
441def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
442    ''' Provide list of optimal data sizes for FFT calculations
443
444    :param int nmin: minimum data size >= 1
445    :param int nmax: maximum data size > nmin
446    :param int thresh: maximum prime factor allowed
447    :Returns: list of data sizes where the maximum prime factor is < thresh
448    ''' 
449    plist = []
450    nmin = max(1,nmin)
451    nmax = max(nmin+1,nmax)
452    for p in range(nmin,nmax):
453        if max(factorize(p).keys()) < thresh:
454            plist.append(p)
455    return plist
456
457np.seterr(divide='ignore')
458
459# Normal distribution
460
461# loc = mu, scale = std
462_norm_pdf_C = 1./math.sqrt(2*math.pi)
463class norm_gen(st.rv_continuous):
464    'needs a doc string'
465     
466    def pdf(self,x,*args,**kwds):
467        loc,scale=kwds['loc'],kwds['scale']
468        x = (x-loc)/scale
469        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
470       
471norm = norm_gen(name='norm',longname='A normal',extradoc="""
472
473Normal distribution
474
475The location (loc) keyword specifies the mean.
476The scale (scale) keyword specifies the standard deviation.
477
478normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
479""")
480
481## Cauchy
482
483# median = loc
484
485class cauchy_gen(st.rv_continuous):
486    'needs a doc string'
487
488    def pdf(self,x,*args,**kwds):
489        loc,scale=kwds['loc'],kwds['scale']
490        x = (x-loc)/scale
491        return 1.0/np.pi/(1.0+x*x) / scale
492       
493cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
494
495Cauchy distribution
496
497cauchy.pdf(x) = 1/(pi*(1+x**2))
498
499This is the t distribution with one degree of freedom.
500""")
501   
502   
503#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
504
505
506class fcjde_gen(st.rv_continuous):
507    """
508    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
509    Ref: J. Appl. Cryst. (1994) 27, 892-900.
510
511    :param x: array -1 to 1
512    :param t: 2-theta position of peak
513    :param s: sum(S/L,H/L); S: sample height, H: detector opening,
514      L: sample to detector opening distance
515    :param dx: 2-theta step size in deg
516
517    :returns: for fcj.pdf
518
519     * T = x*dx+t
520     * s = S/L+H/L
521     * if x < 0::
522
523        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|
524
525     * if x >= 0: fcj.pdf = 0   
526    """
527    def _pdf(self,x,t,s,dx):
528        T = dx*x+t
529        ax2 = abs(npcosd(T))
530        ax = ax2**2
531        bx = npcosd(t)**2
532        bx = np.where(ax>bx,bx,ax)
533        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
534        fx = np.where(fx > 0.,fx,0.0)
535        return fx
536             
537    def pdf(self,x,*args,**kwds):
538        loc=kwds['loc']
539        return self._pdf(x-loc,*args)
540       
541fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
542               
543def getWidthsCW(pos,sig,gam,shl):
544    '''Compute the peak widths used for computing the range of a peak
545    for constant wavelength data. On low-angle side, 50 FWHM are used,
546    on high-angle side 75 are used, low angle side extended for axial divergence
547    (for peaks above 90 deg, these are reversed.)
548    '''
549    widths = [np.sqrt(sig)/100.,gam/100.]
550    fwhm = 2.355*widths[0]+widths[1]
551    fmin = 50.*(fwhm+shl*abs(npcosd(pos)))
552    fmax = 75.0*fwhm
553    if pos > 90:
554        fmin,fmax = [fmax,fmin]         
555    return widths,fmin,fmax
556   
557def getWidthsTOF(pos,alp,bet,sig,gam):
558    '''Compute the peak widths used for computing the range of a peak
559    for constant wavelength data. 50 FWHM are used on both sides each
560    extended by exponential coeff.
561    '''
562    widths = [np.sqrt(sig),gam]
563    fwhm = 2.355*widths[0]+2.*widths[1]
564    fmin = 50.*fwhm*(1.+1./alp)   
565    fmax = 50.*fwhm*(1.+1./bet)
566    return widths,fmin,fmax
567   
568def getFWHM(pos,Inst):
569    '''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
570    via getgamFW(g,s).
571   
572    :param pos: float peak position in deg 2-theta or tof in musec
573    :param Inst: dict instrument parameters
574   
575    :returns float: total FWHM of pseudoVoigt in deg or musec
576    ''' 
577   
578    sig = lambda Th,U,V,W: np.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W))
579    sigTOF = lambda dsp,S0,S1,S2,Sq:  S0+S1*dsp**2+S2*dsp**4+Sq/dsp**2
580    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))
581    gamTOF = lambda dsp,X,Y: X*dsp+Y*dsp**2
582    if 'C' in Inst['Type'][0]:
583        s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
584        g = gam(pos/2.,Inst['X'][1],Inst['Y'][1])
585        return getgamFW(g,s)/100.  #returns FWHM in deg
586    else:
587        dsp = pos/Inst['difC'][0]
588        s = sigTOF(dsp,Inst['sig-0'][1],Inst['sig-1'][1],Inst['sig-2'][1],Inst['sig-q'][1])
589        g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1])
590        return getgamFW(g,s)
591   
592def getgamFW(g,s):
593    '''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
594    lambda fxn needs FWHM for both Gaussian & Lorentzian components
595   
596    :param g: float Lorentzian gamma = FWHM(L)
597    :param s: float Gaussian sig
598   
599    :returns float: total FWHM of pseudoVoigt
600    ''' 
601    gamFW = lambda s,g: np.exp(np.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
602    return gamFW(2.35482*s,g)   #sqrt(8ln2)*sig = FWHM(G)
603               
604def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
605    'needs a doc string'
606    DX = xdata[1]-xdata[0]
607    widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
608    x = np.linspace(pos-fmin,pos+fmin,256)
609    dx = x[1]-x[0]
610    Norm = norm.pdf(x,loc=pos,scale=widths[0])
611    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
612    arg = [pos,shl/57.2958,dx,]
613    FCJ = fcjde.pdf(x,*arg,loc=pos)
614    if len(np.nonzero(FCJ)[0])>5:
615        z = np.column_stack([Norm,Cauchy,FCJ]).T
616        Z = fft(z)
617        Df = ifft(Z.prod(axis=0)).real
618    else:
619        z = np.column_stack([Norm,Cauchy]).T
620        Z = fft(z)
621        Df = fftshift(ifft(Z.prod(axis=0))).real
622    Df /= np.sum(Df)
623    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
624    return intens*Df(xdata)*DX/dx
625
626def getBackground(pfx,parmDict,bakType,dataType,xdata):
627    'needs a doc string'
628    if 'T' in dataType:
629        q = 2.*np.pi*parmDict[pfx+'difC']/xdata
630    elif 'C' in dataType:
631        wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
632        q = npT2q(xdata,wave)
633    yb = np.zeros_like(xdata)
634    nBak = 0
635    cw = np.diff(xdata)
636    cw = np.append(cw,cw[-1])
637    sumBk = [0.,0.,0]
638    while True:
639        key = pfx+'Back;'+str(nBak)
640        if key in parmDict:
641            nBak += 1
642        else:
643            break
644#empirical functions
645    if bakType in ['chebyschev','cosine']:
646        dt = xdata[-1]-xdata[0]   
647        for iBak in range(nBak):
648            key = pfx+'Back;'+str(iBak)
649            if bakType == 'chebyschev':
650                ybi = parmDict[key]*(-1.+2.*(xdata-xdata[0])/dt)**iBak
651            elif bakType == 'cosine':
652                ybi = parmDict[key]*npcosd(180.*xdata*iBak/xdata[-1])
653            yb += ybi
654        sumBk[0] = np.sum(yb)
655    elif bakType in ['Q^2 power series','Q^-2 power series']:
656        QT = 1.
657        yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
658        for iBak in range(nBak-1):
659            key = pfx+'Back;'+str(iBak+1)
660            if '-2' in bakType:
661                QT *= (iBak+1)*q**-2
662            else:
663                QT *= q**2/(iBak+1)
664            yb += QT*parmDict[key]
665        sumBk[0] = np.sum(yb)
666    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
667        if nBak == 1:
668            yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
669        elif nBak == 2:
670            dX = xdata[-1]-xdata[0]
671            T2 = (xdata-xdata[0])/dX
672            T1 = 1.0-T2
673            yb = parmDict[pfx+'Back;0']*T1+parmDict[pfx+'Back;1']*T2
674        else:
675            if bakType == 'lin interpolate':
676                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
677            elif bakType == 'inv interpolate':
678                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
679            elif bakType == 'log interpolate':
680                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
681            bakPos[0] = xdata[0]
682            bakPos[-1] = xdata[-1]
683            bakVals = np.zeros(nBak)
684            for i in range(nBak):
685                bakVals[i] = parmDict[pfx+'Back;'+str(i)]
686            bakInt = si.interp1d(bakPos,bakVals,'linear')
687            yb = bakInt(ma.getdata(xdata))
688        sumBk[0] = np.sum(yb)
689#Debye function       
690    if pfx+'difC' in parmDict:
691        ff = 1.
692    else:       
693        try:
694            wave = parmDict[pfx+'Lam']
695        except KeyError:
696            wave = parmDict[pfx+'Lam1']
697        SQ = (q/(4.*np.pi))**2
698        FF = G2elem.GetFormFactorCoeff('Si')[0]
699        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
700    iD = 0       
701    while True:
702        try:
703            dbA = parmDict[pfx+'DebyeA;'+str(iD)]
704            dbR = parmDict[pfx+'DebyeR;'+str(iD)]
705            dbU = parmDict[pfx+'DebyeU;'+str(iD)]
706            ybi = ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
707            yb += ybi
708            sumBk[1] += np.sum(ybi)
709            iD += 1       
710        except KeyError:
711            break
712#peaks
713    iD = 0
714    while True:
715        try:
716            pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
717            pkI = parmDict[pfx+'BkPkint;'+str(iD)]
718            pkS = parmDict[pfx+'BkPksig;'+str(iD)]
719            pkG = parmDict[pfx+'BkPkgam;'+str(iD)]
720            if 'C' in dataType:
721                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
722            else: #'T'OF
723                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
724            iBeg = np.searchsorted(xdata,pkP-fmin)
725            iFin = np.searchsorted(xdata,pkP+fmax)
726            lenX = len(xdata)
727            if not iBeg:
728                iFin = np.searchsorted(xdata,pkP+fmax)
729            elif iBeg == lenX:
730                iFin = iBeg
731            else:
732                iFin = np.searchsorted(xdata,pkP+fmax)
733            if 'C' in dataType:
734                ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])
735                yb[iBeg:iFin] += ybi
736            else:   #'T'OF
737                ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
738                yb[iBeg:iFin] += ybi
739            sumBk[2] += np.sum(ybi)
740            iD += 1       
741        except KeyError:
742            break
743        except ValueError:
744            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
745            break       
746    return yb,sumBk
747   
748def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata):
749    'needs a doc string'
750    if 'T' in dataType:
751        q = 2.*np.pi*parmDict[hfx+'difC']/xdata
752    elif 'C' in dataType:
753        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
754        q = 2.*np.pi*npsind(xdata/2.)/wave
755    nBak = 0
756    while True:
757        key = hfx+'Back;'+str(nBak)
758        if key in parmDict:
759            nBak += 1
760        else:
761            break
762    dydb = np.zeros(shape=(nBak,len(xdata)))
763    dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
764    dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
765    cw = np.diff(xdata)
766    cw = np.append(cw,cw[-1])
767
768    if bakType in ['chebyschev','cosine']:
769        dt = xdata[-1]-xdata[0]   
770        for iBak in range(nBak):   
771            if bakType == 'chebyschev':
772                dydb[iBak] = (-1.+2.*(xdata-xdata[0])/dt)**iBak
773            elif bakType == 'cosine':
774                dydb[iBak] = npcosd(180.*xdata*iBak/xdata[-1])
775    elif bakType in ['Q^2 power series','Q^-2 power series']:
776        QT = 1.
777        dydb[0] = np.ones_like(xdata)
778        for iBak in range(nBak-1):
779            if '-2' in bakType:
780                QT *= (iBak+1)*q**-2
781            else:
782                QT *= q**2/(iBak+1)
783            dydb[iBak+1] = QT
784    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
785        if nBak == 1:
786            dydb[0] = np.ones_like(xdata)
787        elif nBak == 2:
788            dX = xdata[-1]-xdata[0]
789            T2 = (xdata-xdata[0])/dX
790            T1 = 1.0-T2
791            dydb = [T1,T2]
792        else:
793            if bakType == 'lin interpolate':
794                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
795            elif bakType == 'inv interpolate':
796                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
797            elif bakType == 'log interpolate':
798                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
799            bakPos[0] = xdata[0]
800            bakPos[-1] = xdata[-1]
801            for i,pos in enumerate(bakPos):
802                if i == 0:
803                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
804                elif i == len(bakPos)-1:
805                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
806                else:
807                    dydb[i] = np.where(xdata>bakPos[i],
808                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
809                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
810    if hfx+'difC' in parmDict:
811        ff = 1.
812    else:
813        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
814        q = npT2q(xdata,wave)
815        SQ = (q/(4*np.pi))**2
816        FF = G2elem.GetFormFactorCoeff('Si')[0]
817        ff = np.array(G2elem.ScatFac(FF,SQ)[0])*np.pi**2    #needs pi^2~10. for cw data (why?)
818    iD = 0       
819    while True:
820        try:
821            if hfx+'difC' in parmDict:
822                q = 2*np.pi*parmDict[hfx+'difC']/xdata
823            dbA = parmDict[hfx+'DebyeA;'+str(iD)]
824            dbR = parmDict[hfx+'DebyeR;'+str(iD)]
825            dbU = parmDict[hfx+'DebyeU;'+str(iD)]
826            sqr = np.sin(q*dbR)/(q*dbR)
827            cqr = np.cos(q*dbR)
828            temp = np.exp(-dbU*q**2)
829            dyddb[3*iD] = ff*sqr*temp
830            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(dbR)
831            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2
832            iD += 1
833        except KeyError:
834            break
835    iD = 0
836    while True:
837        try:
838            pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
839            pkI = parmDict[hfx+'BkPkint;'+str(iD)]
840            pkS = parmDict[hfx+'BkPksig;'+str(iD)]
841            pkG = parmDict[hfx+'BkPkgam;'+str(iD)]
842            if 'C' in dataType:
843                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
844            else: #'T'OF
845                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
846            iBeg = np.searchsorted(xdata,pkP-fmin)
847            iFin = np.searchsorted(xdata,pkP+fmax)
848            lenX = len(xdata)
849            if not iBeg:
850                iFin = np.searchsorted(xdata,pkP+fmax)
851            elif iBeg == lenX:
852                iFin = iBeg
853            else:
854                iFin = np.searchsorted(xdata,pkP+fmax)
855            if 'C' in dataType:
856                Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
857                dydpk[4*iD][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdp
858                dydpk[4*iD+1][iBeg:iFin] += 100.*cw[iBeg:iFin]*Df
859                dydpk[4*iD+2][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFds
860                dydpk[4*iD+3][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdg
861            else:   #'T'OF
862                Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
863                dydpk[4*iD][iBeg:iFin] += pkI*dFdp
864                dydpk[4*iD+1][iBeg:iFin] += Df
865                dydpk[4*iD+2][iBeg:iFin] += pkI*dFds
866                dydpk[4*iD+3][iBeg:iFin] += pkI*dFdg
867            iD += 1       
868        except KeyError:
869            break
870        except ValueError:
871            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
872            break       
873    return dydb,dyddb,dydpk
874
875#use old fortran routine
876def getFCJVoigt3(pos,sig,gam,shl,xdata):
877    'needs a doc string'
878   
879    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
880#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
881    Df /= np.sum(Df)
882    return Df
883
884def getdFCJVoigt3(pos,sig,gam,shl,xdata):
885    'needs a doc string'
886   
887    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
888#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
889    return Df,dFdp,dFds,dFdg,dFdsh
890
891def getPsVoigt(pos,sig,gam,xdata):
892    'needs a doc string'
893   
894    Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
895    Df /= np.sum(Df)
896    return Df
897
898def getdPsVoigt(pos,sig,gam,xdata):
899    'needs a doc string'
900   
901    Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
902    return Df,dFdp,dFds,dFdg
903
904def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
905    'needs a doc string'
906    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
907    Df /= np.sum(Df)
908    return Df 
909   
910def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
911    'needs a doc string'
912    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
913    return Df,dFdp,dFda,dFdb,dFds,dFdg   
914
915def ellipseSize(H,Sij,GB):
916    'needs a doc string'
917    HX = np.inner(H.T,GB)
918    lenHX = np.sqrt(np.sum(HX**2))
919    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
920    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
921    lenR = np.sqrt(np.sum(R**2))
922    return lenR
923
924def ellipseSizeDerv(H,Sij,GB):
925    'needs a doc string'
926    lenR = ellipseSize(H,Sij,GB)
927    delt = 0.001
928    dRdS = np.zeros(6)
929    for i in range(6):
930        Sij[i] -= delt
931        lenM = ellipseSize(H,Sij,GB)
932        Sij[i] += 2.*delt
933        lenP = ellipseSize(H,Sij,GB)
934        Sij[i] -= delt
935        dRdS[i] = (lenP-lenM)/(2.*delt)
936    return lenR,dRdS
937
938def getHKLpeak(dmin,SGData,A,Inst=None):
939    'needs a doc string'
940    HKL = G2lat.GenHLaue(dmin,SGData,A)       
941    HKLs = []
942    for h,k,l,d in HKL:
943        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
944        if not ext:
945            if Inst == None:
946                HKLs.append([h,k,l,d,0,-1])
947            else:
948                HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
949    return HKLs
950
951def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
952    'needs a doc string'
953    HKLs = []
954    vec = np.array(Vec)
955    vstar = np.sqrt(G2lat.calc_rDsq(vec,A))     #find extra needed for -n SS reflections
956    dvec = 1./(maxH*vstar+1./dmin)
957    HKL = G2lat.GenHLaue(dvec,SGData,A)       
958    SSdH = [vec*h for h in range(-maxH,maxH+1)]
959    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
960    for h,k,l,d in HKL:
961        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
962        if not ext and d >= dmin:
963            HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
964        for dH in SSdH:
965            if dH:
966                DH = SSdH[dH]
967                H = [h+DH[0],k+DH[1],l+DH[2]]
968                d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
969                if d >= dmin:
970                    HKLM = np.array([h,k,l,dH])
971                    if G2spc.checkSSextc(HKLM,SSGData):
972                        HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])   
973#    GSASIIpath.IPyBreak()
974    return G2lat.sortHKLd(HKLs,True,True,True)
975
976def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
977    'needs a doc string'
978   
979    yb = getBackground('',parmDict,bakType,dataType,xdata)[0]
980    yc = np.zeros_like(yb)
981    cw = np.diff(xdata)
982    cw = np.append(cw,cw[-1])
983    if 'C' in dataType:
984        shl = max(parmDict['SH/L'],0.002)
985        Ka2 = False
986        if 'Lam1' in parmDict.keys():
987            Ka2 = True
988            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
989            kRatio = parmDict['I(L2)/I(L1)']
990        iPeak = 0
991        while True:
992            try:
993                pos = parmDict['pos'+str(iPeak)]
994                tth = (pos-parmDict['Zero'])
995                intens = parmDict['int'+str(iPeak)]
996                sigName = 'sig'+str(iPeak)
997                if sigName in varyList:
998                    sig = parmDict[sigName]
999                else:
1000                    sig = G2mth.getCWsig(parmDict,tth)
1001                sig = max(sig,0.001)          #avoid neg sigma^2
1002                gamName = 'gam'+str(iPeak)
1003                if gamName in varyList:
1004                    gam = parmDict[gamName]
1005                else:
1006                    gam = G2mth.getCWgam(parmDict,tth)
1007                gam = max(gam,0.001)             #avoid neg gamma
1008                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1009                iBeg = np.searchsorted(xdata,pos-fmin)
1010                iFin = np.searchsorted(xdata,pos+fmin)
1011                if not iBeg+iFin:       #peak below low limit
1012                    iPeak += 1
1013                    continue
1014                elif not iBeg-iFin:     #peak above high limit
1015                    return yb+yc
1016                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1017                if Ka2:
1018                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1019                    iBeg = np.searchsorted(xdata,pos2-fmin)
1020                    iFin = np.searchsorted(xdata,pos2+fmin)
1021                    if iBeg-iFin:
1022                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1023                iPeak += 1
1024            except KeyError:        #no more peaks to process
1025                return yb+yc
1026    else:
1027        Pdabc = parmDict['Pdabc']
1028        difC = parmDict['difC']
1029        iPeak = 0
1030        while True:
1031            try:
1032                pos = parmDict['pos'+str(iPeak)]               
1033                tof = pos-parmDict['Zero']
1034                dsp = tof/difC
1035                intens = parmDict['int'+str(iPeak)]
1036                alpName = 'alp'+str(iPeak)
1037                if alpName in varyList:
1038                    alp = parmDict[alpName]
1039                else:
1040                    if len(Pdabc):
1041                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1042                    else:
1043                        alp = G2mth.getTOFalpha(parmDict,dsp)
1044                alp = max(0.0001,alp)
1045                betName = 'bet'+str(iPeak)
1046                if betName in varyList:
1047                    bet = parmDict[betName]
1048                else:
1049                    if len(Pdabc):
1050                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1051                    else:
1052                        bet = G2mth.getTOFbeta(parmDict,dsp)
1053                bet = max(0.0001,bet)
1054                sigName = 'sig'+str(iPeak)
1055                if sigName in varyList:
1056                    sig = parmDict[sigName]
1057                else:
1058                    sig = G2mth.getTOFsig(parmDict,dsp)
1059                gamName = 'gam'+str(iPeak)
1060                if gamName in varyList:
1061                    gam = parmDict[gamName]
1062                else:
1063                    gam = G2mth.getTOFgamma(parmDict,dsp)
1064                gam = max(gam,0.001)             #avoid neg gamma
1065                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1066                iBeg = np.searchsorted(xdata,pos-fmin)
1067                iFin = np.searchsorted(xdata,pos+fmax)
1068                lenX = len(xdata)
1069                if not iBeg:
1070                    iFin = np.searchsorted(xdata,pos+fmax)
1071                elif iBeg == lenX:
1072                    iFin = iBeg
1073                else:
1074                    iFin = np.searchsorted(xdata,pos+fmax)
1075                if not iBeg+iFin:       #peak below low limit
1076                    iPeak += 1
1077                    continue
1078                elif not iBeg-iFin:     #peak above high limit
1079                    return yb+yc
1080                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1081                iPeak += 1
1082            except KeyError:        #no more peaks to process
1083                return yb+yc
1084           
1085def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
1086    'needs a doc string'
1087# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1088    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1089    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,dataType,xdata)
1090    if 'Back;0' in varyList:            #background derivs are in front if present
1091        dMdv[0:len(dMdb)] = dMdb
1092    names = ['DebyeA','DebyeR','DebyeU']
1093    for name in varyList:
1094        if 'Debye' in name:
1095            parm,id = name.split(';')
1096            ip = names.index(parm)
1097            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
1098    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1099    for name in varyList:
1100        if 'BkPk' in name:
1101            parm,id = name.split(';')
1102            ip = names.index(parm)
1103            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
1104    cw = np.diff(xdata)
1105    cw = np.append(cw,cw[-1])
1106    if 'C' in dataType:
1107        shl = max(parmDict['SH/L'],0.002)
1108        Ka2 = False
1109        if 'Lam1' in parmDict.keys():
1110            Ka2 = True
1111            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1112            kRatio = parmDict['I(L2)/I(L1)']
1113        iPeak = 0
1114        while True:
1115            try:
1116                pos = parmDict['pos'+str(iPeak)]
1117                tth = (pos-parmDict['Zero'])
1118                intens = parmDict['int'+str(iPeak)]
1119                sigName = 'sig'+str(iPeak)
1120                if sigName in varyList:
1121                    sig = parmDict[sigName]
1122                    dsdU = dsdV = dsdW = 0
1123                else:
1124                    sig = G2mth.getCWsig(parmDict,tth)
1125                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1126                sig = max(sig,0.001)          #avoid neg sigma
1127                gamName = 'gam'+str(iPeak)
1128                if gamName in varyList:
1129                    gam = parmDict[gamName]
1130                    dgdX = dgdY = 0
1131                else:
1132                    gam = G2mth.getCWgam(parmDict,tth)
1133                    dgdX,dgdY = G2mth.getCWgamDeriv(tth)
1134                gam = max(gam,0.001)             #avoid neg gamma
1135                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1136                iBeg = np.searchsorted(xdata,pos-fmin)
1137                iFin = np.searchsorted(xdata,pos+fmin)
1138                if not iBeg+iFin:       #peak below low limit
1139                    iPeak += 1
1140                    continue
1141                elif not iBeg-iFin:     #peak above high limit
1142                    break
1143                dMdpk = np.zeros(shape=(6,len(xdata)))
1144                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1145                for i in range(1,5):
1146                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
1147                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
1148                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1149                if Ka2:
1150                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1151                    iBeg = np.searchsorted(xdata,pos2-fmin)
1152                    iFin = np.searchsorted(xdata,pos2+fmin)
1153                    if iBeg-iFin:
1154                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1155                        for i in range(1,5):
1156                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
1157                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
1158                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
1159                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1160                for parmName in ['pos','int','sig','gam']:
1161                    try:
1162                        idx = varyList.index(parmName+str(iPeak))
1163                        dMdv[idx] = dervDict[parmName]
1164                    except ValueError:
1165                        pass
1166                if 'U' in varyList:
1167                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1168                if 'V' in varyList:
1169                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1170                if 'W' in varyList:
1171                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1172                if 'X' in varyList:
1173                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1174                if 'Y' in varyList:
1175                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1176                if 'SH/L' in varyList:
1177                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
1178                if 'I(L2)/I(L1)' in varyList:
1179                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1180                iPeak += 1
1181            except KeyError:        #no more peaks to process
1182                break
1183    else:
1184        Pdabc = parmDict['Pdabc']
1185        difC = parmDict['difC']
1186        iPeak = 0
1187        while True:
1188            try:
1189                pos = parmDict['pos'+str(iPeak)]               
1190                tof = pos-parmDict['Zero']
1191                dsp = tof/difC
1192                intens = parmDict['int'+str(iPeak)]
1193                alpName = 'alp'+str(iPeak)
1194                if alpName in varyList:
1195                    alp = parmDict[alpName]
1196                else:
1197                    if len(Pdabc):
1198                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1199                        dada0 = 0
1200                    else:
1201                        alp = G2mth.getTOFalpha(parmDict,dsp)
1202                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1203                betName = 'bet'+str(iPeak)
1204                if betName in varyList:
1205                    bet = parmDict[betName]
1206                else:
1207                    if len(Pdabc):
1208                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1209                        dbdb0 = dbdb1 = dbdb2 = 0
1210                    else:
1211                        bet = G2mth.getTOFbeta(parmDict,dsp)
1212                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1213                sigName = 'sig'+str(iPeak)
1214                if sigName in varyList:
1215                    sig = parmDict[sigName]
1216                    dsds0 = dsds1 = dsds2 = dsds3 = 0
1217                else:
1218                    sig = G2mth.getTOFsig(parmDict,dsp)
1219                    dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1220                gamName = 'gam'+str(iPeak)
1221                if gamName in varyList:
1222                    gam = parmDict[gamName]
1223                    dsdX = dsdY = 0
1224                else:
1225                    gam = G2mth.getTOFgamma(parmDict,dsp)
1226                    dsdX,dsdY = G2mth.getTOFgammaDeriv(dsp)
1227                gam = max(gam,0.001)             #avoid neg gamma
1228                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1229                iBeg = np.searchsorted(xdata,pos-fmin)
1230                lenX = len(xdata)
1231                if not iBeg:
1232                    iFin = np.searchsorted(xdata,pos+fmax)
1233                elif iBeg == lenX:
1234                    iFin = iBeg
1235                else:
1236                    iFin = np.searchsorted(xdata,pos+fmax)
1237                if not iBeg+iFin:       #peak below low limit
1238                    iPeak += 1
1239                    continue
1240                elif not iBeg-iFin:     #peak above high limit
1241                    break
1242                dMdpk = np.zeros(shape=(7,len(xdata)))
1243                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1244                for i in range(1,6):
1245                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1246                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1247                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1248                for parmName in ['pos','int','alp','bet','sig','gam']:
1249                    try:
1250                        idx = varyList.index(parmName+str(iPeak))
1251                        dMdv[idx] = dervDict[parmName]
1252                    except ValueError:
1253                        pass
1254                if 'alpha' in varyList:
1255                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1256                if 'beta-0' in varyList:
1257                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1258                if 'beta-1' in varyList:
1259                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1260                if 'beta-q' in varyList:
1261                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1262                if 'sig-0' in varyList:
1263                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1264                if 'sig-1' in varyList:
1265                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1266                if 'sig-2' in varyList:
1267                    dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1268                if 'sig-q' in varyList:
1269                    dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1270                if 'X' in varyList:
1271                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1272                if 'Y' in varyList:
1273                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
1274                iPeak += 1
1275            except KeyError:        #no more peaks to process
1276                break
1277    return dMdv
1278       
1279def Dict2Values(parmdict, varylist):
1280    '''Use before call to leastsq to setup list of values for the parameters
1281    in parmdict, as selected by key in varylist'''
1282    return [parmdict[key] for key in varylist] 
1283   
1284def Values2Dict(parmdict, varylist, values):
1285    ''' Use after call to leastsq to update the parameter dictionary with
1286    values corresponding to keys in varylist'''
1287    parmdict.update(zip(varylist,values))
1288   
1289def SetBackgroundParms(Background):
1290    'needs a doc string'
1291    if len(Background) == 1:            # fix up old backgrounds
1292        Background.append({'nDebye':0,'debyeTerms':[]})
1293    bakType,bakFlag = Background[0][:2]
1294    backVals = Background[0][3:]
1295    backNames = ['Back;'+str(i) for i in range(len(backVals))]
1296    Debye = Background[1]           #also has background peaks stuff
1297    backDict = dict(zip(backNames,backVals))
1298    backVary = []
1299    if bakFlag:
1300        backVary = backNames
1301
1302    backDict['nDebye'] = Debye['nDebye']
1303    debyeDict = {}
1304    debyeList = []
1305    for i in range(Debye['nDebye']):
1306        debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1307        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1308        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1309    debyeVary = []
1310    for item in debyeList:
1311        if item[1]:
1312            debyeVary.append(item[0])
1313    backDict.update(debyeDict)
1314    backVary += debyeVary
1315
1316    backDict['nPeaks'] = Debye['nPeaks']
1317    peaksDict = {}
1318    peaksList = []
1319    for i in range(Debye['nPeaks']):
1320        peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1321        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1322        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1323    peaksVary = []
1324    for item in peaksList:
1325        if item[1]:
1326            peaksVary.append(item[0])
1327    backDict.update(peaksDict)
1328    backVary += peaksVary   
1329    return bakType,backDict,backVary
1330   
1331def DoCalibInst(IndexPeaks,Inst):
1332   
1333    def SetInstParms():
1334        dataType = Inst['Type'][0]
1335        insVary = []
1336        insNames = []
1337        insVals = []
1338        for parm in Inst:
1339            insNames.append(parm)
1340            insVals.append(Inst[parm][1])
1341            if parm in ['Lam','difC','difA','difB','Zero',]:
1342                if Inst[parm][2]:
1343                    insVary.append(parm)
1344        instDict = dict(zip(insNames,insVals))
1345        return dataType,instDict,insVary
1346       
1347    def GetInstParms(parmDict,Inst,varyList):
1348        for name in Inst:
1349            Inst[name][1] = parmDict[name]
1350       
1351    def InstPrint(Inst,sigDict):
1352        print 'Instrument Parameters:'
1353        if 'C' in Inst['Type'][0]:
1354            ptfmt = "%12.6f"
1355        else:
1356            ptfmt = "%12.3f"
1357        ptlbls = 'names :'
1358        ptstr =  'values:'
1359        sigstr = 'esds  :'
1360        for parm in Inst:
1361            if parm in  ['Lam','difC','difA','difB','Zero',]:
1362                ptlbls += "%s" % (parm.center(12))
1363                ptstr += ptfmt % (Inst[parm][1])
1364                if parm in sigDict:
1365                    sigstr += ptfmt % (sigDict[parm])
1366                else:
1367                    sigstr += 12*' '
1368        print ptlbls
1369        print ptstr
1370        print sigstr
1371       
1372    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1373        parmDict.update(zip(varyList,values))
1374        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1375
1376    peakPos = []
1377    peakDsp = []
1378    peakWt = []
1379    for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1380        if peak[2] and peak[3] and sig > 0.:
1381            peakPos.append(peak[0])
1382            peakDsp.append(peak[-1])    #d-calc
1383            peakWt.append(peak[-1]**2/sig**2)   #weight by d**2
1384    peakPos = np.array(peakPos)
1385    peakDsp = np.array(peakDsp)
1386    peakWt = np.array(peakWt)
1387    dataType,insDict,insVary = SetInstParms()
1388    parmDict = {}
1389    parmDict.update(insDict)
1390    varyList = insVary
1391    if not len(varyList):
1392        print '**** ERROR - nothing to refine! ****'
1393        return False
1394    while True:
1395        begin = time.time()
1396        values =  np.array(Dict2Values(parmDict, varyList))
1397        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1398            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1399        ncyc = int(result[2]['nfev']/2)
1400        runtime = time.time()-begin   
1401        chisq = np.sum(result[2]['fvec']**2)
1402        Values2Dict(parmDict, varyList, result[0])
1403        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
1404        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(peakPos),' Number of parameters: ',len(varyList)
1405        print 'calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1406        print 'chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)
1407        try:
1408            sig = np.sqrt(np.diag(result[1])*GOF)
1409            if np.any(np.isnan(sig)):
1410                print '*** Least squares aborted - some invalid esds possible ***'
1411            break                   #refinement succeeded - finish up!
1412        except ValueError:          #result[1] is None on singular matrix
1413            print '**** Refinement failed - singular matrix ****'
1414       
1415    sigDict = dict(zip(varyList,sig))
1416    GetInstParms(parmDict,Inst,varyList)
1417    InstPrint(Inst,sigDict)
1418    return True
1419           
1420def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
1421    '''Called to perform a peak fit, refining the selected items in the peak
1422    table as well as selected items in the background.
1423
1424    :param str FitPgm: type of fit to perform. At present this is ignored.
1425    :param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1426      four values followed by a refine flag where the values are: position, intensity,
1427      sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1428      tree entry, dict item "peaks"
1429    :param list Background: describes the background. List with two items.
1430      Item 0 specifies a background model and coefficients. Item 1 is a dict.
1431      From the Histogram/Background tree entry.
1432    :param list Limits: min and max x-value to use
1433    :param dict Inst: Instrument parameters
1434    :param dict Inst2: more Instrument parameters
1435    :param numpy.array data: a 5xn array. data[0] is the x-values,
1436      data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1437      calc, background and difference intensities, respectively.
1438    :param list prevVaryList: Used in sequential refinements to override the
1439      variable list. Defaults as an empty list.
1440    :param bool oneCycle: True if only one cycle of fitting should be performed
1441    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1442      and derivType = controls['deriv type']. If None default values are used.
1443    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1444      Defaults to None, which means no updates are done.
1445    '''
1446    def GetBackgroundParms(parmList,Background):
1447        iBak = 0
1448        while True:
1449            try:
1450                bakName = 'Back;'+str(iBak)
1451                Background[0][iBak+3] = parmList[bakName]
1452                iBak += 1
1453            except KeyError:
1454                break
1455        iDb = 0
1456        while True:
1457            names = ['DebyeA;','DebyeR;','DebyeU;']
1458            try:
1459                for i,name in enumerate(names):
1460                    val = parmList[name+str(iDb)]
1461                    Background[1]['debyeTerms'][iDb][2*i] = val
1462                iDb += 1
1463            except KeyError:
1464                break
1465        iDb = 0
1466        while True:
1467            names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1468            try:
1469                for i,name in enumerate(names):
1470                    val = parmList[name+str(iDb)]
1471                    Background[1]['peaksList'][iDb][2*i] = val
1472                iDb += 1
1473            except KeyError:
1474                break
1475               
1476    def BackgroundPrint(Background,sigDict):
1477        print 'Background coefficients for',Background[0][0],'function'
1478        ptfmt = "%12.5f"
1479        ptstr =  'value: '
1480        sigstr = 'esd  : '
1481        for i,back in enumerate(Background[0][3:]):
1482            ptstr += ptfmt % (back)
1483            if Background[0][1]:
1484                prm = 'Back;'+str(i)
1485                if prm in sigDict:
1486                    sigstr += ptfmt % (sigDict[prm])
1487                else:
1488                    sigstr += " "*12
1489            if len(ptstr) > 75:
1490                print ptstr
1491                if Background[0][1]: print sigstr
1492                ptstr =  'value: '
1493                sigstr = 'esd  : '
1494        if len(ptstr) > 8:
1495            print ptstr
1496            if Background[0][1]: print sigstr
1497
1498        if Background[1]['nDebye']:
1499            parms = ['DebyeA;','DebyeR;','DebyeU;']
1500            print 'Debye diffuse scattering coefficients'
1501            ptfmt = "%12.5f"
1502            print ' term       DebyeA       esd        DebyeR       esd        DebyeU        esd'
1503            for term in range(Background[1]['nDebye']):
1504                line = ' term %d'%(term)
1505                for ip,name in enumerate(parms):
1506                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1507                    if name+str(term) in sigDict:
1508                        line += ptfmt%(sigDict[name+str(term)])
1509                    else:
1510                        line += " "*12
1511                print line
1512        if Background[1]['nPeaks']:
1513            print 'Coefficients for Background Peaks'
1514            ptfmt = "%15.3f"
1515            for j,pl in enumerate(Background[1]['peaksList']):
1516                names =  'peak %3d:'%(j+1)
1517                ptstr =  'values  :'
1518                sigstr = 'esds    :'
1519                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1520                    val = pl[2*i]
1521                    prm = lbl+";"+str(j)
1522                    names += '%15s'%(prm)
1523                    ptstr += ptfmt%(val)
1524                    if prm in sigDict:
1525                        sigstr += ptfmt%(sigDict[prm])
1526                    else:
1527                        sigstr += " "*15
1528                print names
1529                print ptstr
1530                print sigstr
1531                           
1532    def SetInstParms(Inst):
1533        dataType = Inst['Type'][0]
1534        insVary = []
1535        insNames = []
1536        insVals = []
1537        for parm in Inst:
1538            insNames.append(parm)
1539            insVals.append(Inst[parm][1])
1540            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1541                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1542                    insVary.append(parm)
1543        instDict = dict(zip(insNames,insVals))
1544        instDict['X'] = max(instDict['X'],0.01)
1545        instDict['Y'] = max(instDict['Y'],0.01)
1546        if 'SH/L' in instDict:
1547            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1548        return dataType,instDict,insVary
1549       
1550    def GetInstParms(parmDict,Inst,varyList,Peaks):
1551        for name in Inst:
1552            Inst[name][1] = parmDict[name]
1553        iPeak = 0
1554        while True:
1555            try:
1556                sigName = 'sig'+str(iPeak)
1557                pos = parmDict['pos'+str(iPeak)]
1558                if sigName not in varyList:
1559                    if 'C' in Inst['Type'][0]:
1560                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1561                    else:
1562                        dsp = G2lat.Pos2dsp(Inst,pos)
1563                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1564                gamName = 'gam'+str(iPeak)
1565                if gamName not in varyList:
1566                    if 'C' in Inst['Type'][0]:
1567                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1568                    else:
1569                        dsp = G2lat.Pos2dsp(Inst,pos)
1570                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1571                iPeak += 1
1572            except KeyError:
1573                break
1574       
1575    def InstPrint(Inst,sigDict):
1576        print 'Instrument Parameters:'
1577        ptfmt = "%12.6f"
1578        ptlbls = 'names :'
1579        ptstr =  'values:'
1580        sigstr = 'esds  :'
1581        for parm in Inst:
1582            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1583                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1584                ptlbls += "%s" % (parm.center(12))
1585                ptstr += ptfmt % (Inst[parm][1])
1586                if parm in sigDict:
1587                    sigstr += ptfmt % (sigDict[parm])
1588                else:
1589                    sigstr += 12*' '
1590        print ptlbls
1591        print ptstr
1592        print sigstr
1593
1594    def SetPeaksParms(dataType,Peaks):
1595        peakNames = []
1596        peakVary = []
1597        peakVals = []
1598        if 'C' in dataType:
1599            names = ['pos','int','sig','gam']
1600        else:
1601            names = ['pos','int','alp','bet','sig','gam']
1602        for i,peak in enumerate(Peaks):
1603            for j,name in enumerate(names):
1604                peakVals.append(peak[2*j])
1605                parName = name+str(i)
1606                peakNames.append(parName)
1607                if peak[2*j+1]:
1608                    peakVary.append(parName)
1609        return dict(zip(peakNames,peakVals)),peakVary
1610               
1611    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1612        if 'C' in Inst['Type'][0]:
1613            names = ['pos','int','sig','gam']
1614        else:   #'T'
1615            names = ['pos','int','alp','bet','sig','gam']
1616        for i,peak in enumerate(Peaks):
1617            pos = parmDict['pos'+str(i)]
1618            if 'difC' in Inst:
1619                dsp = pos/Inst['difC'][1]
1620            for j in range(len(names)):
1621                parName = names[j]+str(i)
1622                if parName in varyList:
1623                    peak[2*j] = parmDict[parName]
1624                elif 'alpha' in parName:
1625                    peak[2*j] = parmDict['alpha']/dsp
1626                elif 'beta' in parName:
1627                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1628                elif 'sig' in parName:
1629                    if 'C' in Inst['Type'][0]:
1630                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1631                    else:
1632                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1633                elif 'gam' in parName:
1634                    if 'C' in Inst['Type'][0]:
1635                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1636                    else:
1637                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1638                       
1639    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1640        print 'Peak coefficients:'
1641        if 'C' in dataType:
1642            names = ['pos','int','sig','gam']
1643        else:   #'T'
1644            names = ['pos','int','alp','bet','sig','gam']           
1645        head = 13*' '
1646        for name in names:
1647            if name in ['alp','bet']:
1648                head += name.center(8)+'esd'.center(8)
1649            else:
1650                head += name.center(10)+'esd'.center(10)
1651        print head
1652        if 'C' in dataType:
1653            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1654        else:
1655            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1656        for i,peak in enumerate(Peaks):
1657            ptstr =  ':'
1658            for j in range(len(names)):
1659                name = names[j]
1660                parName = name+str(i)
1661                ptstr += ptfmt[name] % (parmDict[parName])
1662                if parName in varyList:
1663                    ptstr += ptfmt[name] % (sigDict[parName])
1664                else:
1665                    if name in ['alp','bet']:
1666                        ptstr += 8*' '
1667                    else:
1668                        ptstr += 10*' '
1669            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1670               
1671    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1672        parmdict.update(zip(varylist,values))
1673        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1674           
1675    def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):       
1676        parmdict.update(zip(varylist,values))
1677        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1678        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1679        if dlg:
1680            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1681            if not GoOn:
1682                return -M           #abort!!
1683        return M
1684       
1685    if controls:
1686        Ftol = controls['min dM/M']
1687    else:
1688        Ftol = 0.0001
1689    if oneCycle:
1690        Ftol = 1.0
1691    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1692    yc *= 0.                            #set calcd ones to zero
1693    yb *= 0.
1694    yd *= 0.
1695    xBeg = np.searchsorted(x,Limits[0])
1696    xFin = np.searchsorted(x,Limits[1])
1697    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1698    dataType,insDict,insVary = SetInstParms(Inst)
1699    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1700    parmDict = {}
1701    parmDict.update(bakDict)
1702    parmDict.update(insDict)
1703    parmDict.update(peakDict)
1704    parmDict['Pdabc'] = []      #dummy Pdabc
1705    parmDict.update(Inst2)      #put in real one if there
1706    if prevVaryList:
1707        varyList = prevVaryList[:]
1708    else:
1709        varyList = bakVary+insVary+peakVary
1710    fullvaryList = varyList[:]
1711    while True:
1712        begin = time.time()
1713        values =  np.array(Dict2Values(parmDict, varyList))
1714        Rvals = {}
1715        badVary = []
1716        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1717               args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1718        ncyc = int(result[2]['nfev']/2)
1719        runtime = time.time()-begin   
1720        chisq = np.sum(result[2]['fvec']**2)
1721        Values2Dict(parmDict, varyList, result[0])
1722        Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1723        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1724        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1725        if ncyc:
1726            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1727        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1728        sig = [0]*len(varyList)
1729        if len(varyList) == 0: break  # if nothing was refined
1730        try:
1731            sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1732            if np.any(np.isnan(sig)):
1733                print '*** Least squares aborted - some invalid esds possible ***'
1734            break                   #refinement succeeded - finish up!
1735        except ValueError:          #result[1] is None on singular matrix
1736            print '**** Refinement failed - singular matrix ****'
1737            Ipvt = result[2]['ipvt']
1738            for i,ipvt in enumerate(Ipvt):
1739                if not np.sum(result[2]['fjac'],axis=1)[i]:
1740                    print 'Removing parameter: ',varyList[ipvt-1]
1741                    badVary.append(varyList[ipvt-1])
1742                    del(varyList[ipvt-1])
1743                    break
1744            else: # nothing removed
1745                break
1746    if dlg: dlg.Destroy()
1747    sigDict = dict(zip(varyList,sig))
1748    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]
1749    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1750    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1751    GetBackgroundParms(parmDict,Background)
1752    if bakVary: BackgroundPrint(Background,sigDict)
1753    GetInstParms(parmDict,Inst,varyList,Peaks)
1754    if insVary: InstPrint(Inst,sigDict)
1755    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1756    if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList)
1757    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
1758
1759def calcIncident(Iparm,xdata):
1760    'needs a doc string'
1761
1762    def IfunAdv(Iparm,xdata):
1763        Itype = Iparm['Itype']
1764        Icoef = Iparm['Icoeff']
1765        DYI = np.ones((12,xdata.shape[0]))
1766        YI = np.ones_like(xdata)*Icoef[0]
1767       
1768        x = xdata/1000.                 #expressions are in ms
1769        if Itype == 'Exponential':
1770            for i in [1,3,5,7,9]:
1771                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1772                YI += Icoef[i]*Eterm
1773                DYI[i] *= Eterm
1774                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
1775        elif 'Maxwell'in Itype:
1776            Eterm = np.exp(-Icoef[2]/x**2)
1777            DYI[1] = Eterm/x**5
1778            DYI[2] = -Icoef[1]*DYI[1]/x**2
1779            YI += (Icoef[1]*Eterm/x**5)
1780            if 'Exponential' in Itype:
1781                for i in range(3,11,2):
1782                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1783                    YI += Icoef[i]*Eterm
1784                    DYI[i] *= Eterm
1785                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
1786            else:   #Chebyschev
1787                T = (2./x)-1.
1788                Ccof = np.ones((12,xdata.shape[0]))
1789                Ccof[1] = T
1790                for i in range(2,12):
1791                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1792                for i in range(1,10):
1793                    YI += Ccof[i]*Icoef[i+2]
1794                    DYI[i+2] =Ccof[i]
1795        return YI,DYI
1796       
1797    Iesd = np.array(Iparm['Iesd'])
1798    Icovar = Iparm['Icovar']
1799    YI,DYI = IfunAdv(Iparm,xdata)
1800    YI = np.where(YI>0,YI,1.)
1801    WYI = np.zeros_like(xdata)
1802    vcov = np.zeros((12,12))
1803    k = 0
1804    for i in range(12):
1805        for j in range(i,12):
1806            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1807            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1808            k += 1
1809    M = np.inner(vcov,DYI.T)
1810    WYI = np.sum(M*DYI,axis=0)
1811    WYI = np.where(WYI>0.,WYI,0.)
1812    return YI,WYI
1813   
1814################################################################################
1815# Stacking fault simulation codes
1816################################################################################
1817
1818def GetStackParms(Layers):
1819   
1820    Parms = []
1821#cell parms
1822    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
1823        Parms.append('cellA')
1824        Parms.append('cellC')
1825    else:
1826        Parms.append('cellA')
1827        Parms.append('cellB')
1828        Parms.append('cellC')
1829        if Layers['Laue'] != 'mmm':
1830            Parms.append('cellG')
1831#Transition parms
1832    for iY in range(len(Layers['Layers'])):
1833        for iX in range(len(Layers['Layers'])):
1834            Parms.append('TransP;%d;%d'%(iY,iX))
1835            Parms.append('TransX;%d;%d'%(iY,iX))
1836            Parms.append('TransY;%d;%d'%(iY,iX))
1837            Parms.append('TransZ;%d;%d'%(iY,iX))
1838    return Parms
1839
1840def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
1841    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
1842   
1843    param: Layers dict: 'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
1844                        'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
1845                        'Layers':[],'Stacking':[],'Transitions':[]}
1846    param: ctrls string: controls string to be written on DIFFaX controls.dif file
1847    param: scale float: scale factor
1848    param: background dict: background parameters
1849    param: limits list: min/max 2-theta to be calculated
1850    param: inst dict: instrument parameters dictionary
1851    param: profile list: powder pattern data
1852   
1853    all updated in place   
1854    '''
1855    import atmdata
1856    path = sys.path
1857    for name in path:
1858        if 'bin' in name:
1859            DIFFaX = name+'/DIFFaX.exe'
1860            print ' Execute ',DIFFaX
1861            break
1862    # make form factor file that DIFFaX wants - atom types are GSASII style
1863    sf = open('data.sfc','w')
1864    sf.write('GSASII special form factor file for DIFFaX\n\n')
1865    atTypes = Layers['AtInfo'].keys()
1866    if 'H' not in atTypes:
1867        atTypes.insert(0,'H')
1868    for atType in atTypes:
1869        if atType == 'H': 
1870            blen = -.3741
1871        else:
1872            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
1873        Adat = atmdata.XrayFF[atType]
1874        text = '%4s'%(atType.ljust(4))
1875        for i in range(4):
1876            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
1877        text += '%11.6f%11.6f'%(Adat['fc'],blen)
1878        text += '%3d\n'%(Adat['Z'])
1879        sf.write(text)
1880    sf.close()
1881    #make DIFFaX control.dif file - future use GUI to set some of these flags
1882    cf = open('control.dif','w')
1883    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
1884        x0 = profile[0]
1885        iBeg = np.searchsorted(x0,limits[0])
1886        iFin = np.searchsorted(x0,limits[1])
1887        if iFin-iBeg > 20000:
1888            iFin = iBeg+20000
1889        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
1890        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
1891        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
1892    else:
1893        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
1894        inst = {'Type':['XSC','XSC',]}
1895    cf.close()
1896    #make DIFFaX data file
1897    df = open('GSASII-DIFFaX.dat','w')
1898    df.write('INSTRUMENTAL\n')
1899    if 'X' in inst['Type'][0]:
1900        df.write('X-RAY\n')
1901    elif 'N' in inst['Type'][0]:
1902        df.write('NEUTRON\n')
1903    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
1904        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
1905        U = ateln2*inst['U'][1]/10000.
1906        V = ateln2*inst['V'][1]/10000.
1907        W = ateln2*inst['W'][1]/10000.
1908        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
1909        HW = np.sqrt(np.mean(HWHM))
1910    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
1911        if 'Mean' in Layers['selInst']:
1912            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
1913        elif 'Gaussian' in Layers['selInst']:
1914            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
1915        else:
1916            df.write('None\n')
1917    else:
1918        df.write('0.10\nNone\n')
1919    df.write('STRUCTURAL\n')
1920    a,b,c = Layers['Cell'][1:4]
1921    gam = Layers['Cell'][6]
1922    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
1923    laue = Layers['Laue']
1924    if laue == '2/m(ab)':
1925        laue = '2/m(1)'
1926    elif laue == '2/m(c)':
1927        laue = '2/m(2)'
1928    if 'unknown' in Layers['Laue']:
1929        df.write('%s %.3f\n'%(laue,Layers['Toler']))
1930    else:   
1931        df.write('%s\n'%(laue))
1932    df.write('%d\n'%(len(Layers['Layers'])))
1933    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
1934        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
1935    layerNames = []
1936    for layer in Layers['Layers']:
1937        layerNames.append(layer['Name'])
1938    for il,layer in enumerate(Layers['Layers']):
1939        if layer['SameAs']:
1940            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
1941            continue
1942        df.write('LAYER %d\n'%(il+1))
1943        if '-1' in layer['Symm']:
1944            df.write('CENTROSYMMETRIC\n')
1945        else:
1946            df.write('NONE\n')
1947        for ia,atom in enumerate(layer['Atoms']):
1948            [name,atype,x,y,z,frac,Uiso] = atom
1949            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
1950                frac /= 2.
1951            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
1952    df.write('STACKING\n')
1953    df.write('%s\n'%(Layers['Stacking'][0]))
1954    if 'recursive' in Layers['Stacking'][0]:
1955        df.write('%s\n'%Layers['Stacking'][1])
1956    else:
1957        if 'list' in Layers['Stacking'][1]:
1958            Slen = len(Layers['Stacking'][2])
1959            iB = 0
1960            iF = 0
1961            while True:
1962                iF += 68
1963                if iF >= Slen:
1964                    break
1965                iF = min(iF,Slen)
1966                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
1967                iB = iF
1968        else:
1969            df.write('%s\n'%Layers['Stacking'][1])   
1970    df.write('TRANSITIONS\n')
1971    for iY in range(len(Layers['Layers'])):
1972        sumPx = 0.
1973        for iX in range(len(Layers['Layers'])):
1974            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
1975            p = round(p,3)
1976            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
1977            sumPx += p
1978        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
1979            print 'ERROR - Layer probabilities sum to ',sumPx,' DIFFaX will insist it = 1.0'
1980            df.close()
1981            os.remove('data.sfc')
1982            os.remove('control.dif')
1983            os.remove('GSASII-DIFFaX.dat')
1984            return       
1985    df.close()   
1986    time0 = time.time()
1987    try:
1988        subp.call(DIFFaX)
1989    except OSError:
1990        print ' DIFFax.exe is not available for this platform - under development'
1991    print ' DIFFaX time = %.2fs'%(time.time()-time0)
1992    if os.path.exists('GSASII-DIFFaX.spc'):
1993        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
1994        iFin = iBeg+Xpat.shape[1]
1995        bakType,backDict,backVary = SetBackgroundParms(background)
1996        backDict['Lam1'] = G2mth.getWave(inst)
1997        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
1998        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
1999        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2000            rv = st.poisson(profile[3][iBeg:iFin])
2001            profile[1][iBeg:iFin] = rv.rvs()
2002            Z = np.ones_like(profile[3][iBeg:iFin])
2003            Z[1::2] *= -1
2004            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2005            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2006        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2007    #cleanup files..
2008        os.remove('GSASII-DIFFaX.spc')
2009    elif os.path.exists('GSASII-DIFFaX.sadp'):
2010        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
2011        Sadp = np.reshape(Sadp,(256,-1))
2012        Layers['Sadp']['Img'] = Sadp
2013        os.remove('GSASII-DIFFaX.sadp')
2014    os.remove('data.sfc')
2015    os.remove('control.dif')
2016    os.remove('GSASII-DIFFaX.dat')
2017   
2018def SetPWDRscan(inst,limits,profile):
2019   
2020    wave = G2mth.getMeanWave(inst)
2021    x0 = profile[0]
2022    iBeg = np.searchsorted(x0,limits[0])
2023    iFin = np.searchsorted(x0,limits[1])
2024    if iFin-iBeg > 20000:
2025        iFin = iBeg+20000
2026    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2027    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
2028    return iFin-iBeg
2029       
2030def SetStackingSF(Layers,debug):
2031# Load scattering factors into DIFFaX arrays
2032    import atmdata
2033    atTypes = Layers['AtInfo'].keys()
2034    aTypes = []
2035    for atype in atTypes:
2036        aTypes.append('%4s'%(atype.ljust(4)))
2037    SFdat = []
2038    for atType in atTypes:
2039        Adat = atmdata.XrayFF[atType]
2040        SF = np.zeros(9)
2041        SF[:8:2] = Adat['fa']
2042        SF[1:8:2] = Adat['fb']
2043        SF[8] = Adat['fc']
2044        SFdat.append(SF)
2045    SFdat = np.array(SFdat)
2046    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
2047   
2048def SetStackingClay(Layers,Type):
2049# Controls
2050    rand.seed()
2051    ranSeed = rand.randint(1,2**16-1)
2052    try:   
2053        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
2054            '6/m','6/mmm'].index(Layers['Laue'])+1
2055    except ValueError:  #for 'unknown'
2056        laueId = -1
2057    if 'SADP' in Type:
2058        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
2059        lmax = int(Layers['Sadp']['Lmax'])
2060    else:
2061        planeId = 0
2062        lmax = 0
2063# Sequences
2064    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
2065    try:
2066        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
2067    except ValueError:
2068        StkParm = -1
2069    if StkParm == 2:    #list
2070        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
2071        Nstk = len(StkSeq)
2072    else:
2073        Nstk = 1
2074        StkSeq = [0,]
2075    if StkParm == -1:
2076        StkParm = int(Layers['Stacking'][1])
2077    Wdth = Layers['Width'][0]
2078    mult = 1
2079    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
2080    LaueSym = Layers['Laue'].ljust(12)
2081    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
2082    return laueId,controls
2083   
2084def SetCellAtoms(Layers):
2085    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
2086# atoms in layers
2087    atTypes = Layers['AtInfo'].keys()
2088    AtomXOU = []
2089    AtomTp = []
2090    LayerSymm = []
2091    LayerNum = []
2092    layerNames = []
2093    Natm = 0
2094    Nuniq = 0
2095    for layer in Layers['Layers']:
2096        layerNames.append(layer['Name'])
2097    for il,layer in enumerate(Layers['Layers']):
2098        if layer['SameAs']:
2099            LayerNum.append(layerNames.index(layer['SameAs'])+1)
2100            continue
2101        else:
2102            LayerNum.append(il+1)
2103            Nuniq += 1
2104        if '-1' in layer['Symm']:
2105            LayerSymm.append(1)
2106        else:
2107            LayerSymm.append(0)
2108        for ia,atom in enumerate(layer['Atoms']):
2109            [name,atype,x,y,z,frac,Uiso] = atom
2110            Natm += 1
2111            AtomTp.append('%4s'%(atype.ljust(4)))
2112            Ta = atTypes.index(atype)+1
2113            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
2114    AtomXOU = np.array(AtomXOU)
2115    Nlayers = len(layerNames)
2116    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
2117    return Nlayers
2118   
2119def SetStackingTrans(Layers,Nlayers):
2120# Transitions
2121    TransX = []
2122    TransP = []
2123    for Ytrans in Layers['Transitions']:
2124        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
2125        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
2126    TransP = np.array(TransP,dtype='float').T
2127    TransX = np.array(TransX,dtype='float')
2128#    GSASIIpath.IPyBreak()
2129    pyx.pygettrans(Nlayers,TransP,TransX)
2130   
2131def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
2132# Scattering factors
2133    SetStackingSF(Layers,debug)
2134# Controls & sequences
2135    laueId,controls = SetStackingClay(Layers,'PWDR')
2136# cell & atoms
2137    Nlayers = SetCellAtoms(Layers)
2138    Volume = Layers['Cell'][7]   
2139# Transitions
2140    SetStackingTrans(Layers,Nlayers)
2141# PWDR scan
2142    Nsteps = SetPWDRscan(inst,limits,profile)
2143# result as Spec
2144    x0 = profile[0]
2145    profile[3] = np.zeros(len(profile[0]))
2146    profile[4] = np.zeros(len(profile[0]))
2147    profile[5] = np.zeros(len(profile[0]))
2148    iBeg = np.searchsorted(x0,limits[0])
2149    iFin = np.searchsorted(x0,limits[1])
2150    if iFin-iBeg > 20000:
2151        iFin = iBeg+20000
2152    Nspec = 20001       
2153    spec = np.zeros(Nspec,dtype='double')   
2154    time0 = time.time()
2155    pyx.pygetspc(controls,Nspec,spec)
2156    print ' GETSPC time = %.2fs'%(time.time()-time0)
2157    time0 = time.time()
2158    U = ateln2*inst['U'][1]/10000.
2159    V = ateln2*inst['V'][1]/10000.
2160    W = ateln2*inst['W'][1]/10000.
2161    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2162    HW = np.sqrt(np.mean(HWHM))
2163    BrdSpec = np.zeros(Nsteps)
2164    if 'Mean' in Layers['selInst']:
2165        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
2166    elif 'Gaussian' in Layers['selInst']:
2167        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
2168    else:
2169        BrdSpec = spec[:Nsteps]
2170    BrdSpec /= Volume
2171    iFin = iBeg+Nsteps
2172    bakType,backDict,backVary = SetBackgroundParms(background)
2173    backDict['Lam1'] = G2mth.getWave(inst)
2174    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2175    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
2176    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2177        rv = st.poisson(profile[3][iBeg:iFin])
2178        profile[1][iBeg:iFin] = rv.rvs()
2179        Z = np.ones_like(profile[3][iBeg:iFin])
2180        Z[1::2] *= -1
2181        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2182        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2183    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2184    print ' Broadening time = %.2fs'%(time.time()-time0)
2185   
2186def CalcStackingSADP(Layers,debug):
2187   
2188# Scattering factors
2189    SetStackingSF(Layers,debug)
2190# Controls & sequences
2191    laueId,controls = SetStackingClay(Layers,'SADP')
2192# cell & atoms
2193    Nlayers = SetCellAtoms(Layers)   
2194# Transitions
2195    SetStackingTrans(Layers,Nlayers)
2196# result as Sadp
2197    Nspec = 20001       
2198    spec = np.zeros(Nspec,dtype='double')   
2199    time0 = time.time()
2200    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
2201    Sapd = np.zeros((256,256))
2202    iB = 0
2203    for i in range(hkLim):
2204        iF = iB+Nblk
2205        p1 = 127+int(i*Incr)
2206        p2 = 128-int(i*Incr)
2207        if Nblk == 128:
2208            if i:
2209                Sapd[128:,p1] = spec[iB:iF]
2210                Sapd[:128,p1] = spec[iF:iB:-1]
2211            Sapd[128:,p2] = spec[iB:iF]
2212            Sapd[:128,p2] = spec[iF:iB:-1]
2213        else:
2214            if i:
2215                Sapd[:,p1] = spec[iB:iF]
2216            Sapd[:,p2] = spec[iB:iF]
2217        iB += Nblk
2218    Layers['Sadp']['Img'] = Sapd
2219    print ' GETSAD time = %.2fs'%(time.time()-time0)
2220#    GSASIIpath.IPyBreak()
2221   
2222#testing data
2223NeedTestData = True
2224def TestData():
2225    'needs a doc string'
2226#    global NeedTestData
2227    global bakType
2228    bakType = 'chebyschev'
2229    global xdata
2230    xdata = np.linspace(4.0,40.0,36000)
2231    global parmDict0
2232    parmDict0 = {
2233        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
2234        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
2235        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
2236        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
2237        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
2238        'Back0':5.384,'Back1':-0.015,'Back2':.004,
2239        }
2240    global parmDict1
2241    parmDict1 = {
2242        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
2243        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
2244        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
2245        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
2246        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
2247        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
2248        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
2249        'Back0':36.897,'Back1':-0.508,'Back2':.006,
2250        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2251        }
2252    global parmDict2
2253    parmDict2 = {
2254        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
2255        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
2256        'Back0':5.,'Back1':-0.02,'Back2':.004,
2257#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2258        }
2259    global varyList
2260    varyList = []
2261
2262def test0():
2263    if NeedTestData: TestData()
2264    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
2265    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
2266    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
2267    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
2268    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
2269    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
2270   
2271def test1():
2272    if NeedTestData: TestData()
2273    time0 = time.time()
2274    for i in range(100):
2275        getPeakProfile(parmDict1,xdata,varyList,bakType)
2276    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
2277   
2278def test2(name,delt):
2279    if NeedTestData: TestData()
2280    varyList = [name,]
2281    xdata = np.linspace(5.6,5.8,400)
2282    hplot = plotter.add('derivatives test for '+name).gca()
2283    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
2284    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2285    parmDict2[name] += delt
2286    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2287    hplot.plot(xdata,(y1-y0)/delt,'r+')
2288   
2289def test3(name,delt):
2290    if NeedTestData: TestData()
2291    names = ['pos','sig','gam','shl']
2292    idx = names.index(name)
2293    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
2294    xdata = np.linspace(5.6,5.8,800)
2295    dx = xdata[1]-xdata[0]
2296    hplot = plotter.add('derivatives test for '+name).gca()
2297    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
2298    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2299    myDict[name] += delt
2300    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2301    hplot.plot(xdata,(y1-y0)/delt,'r+')
2302
2303if __name__ == '__main__':
2304    import GSASIItestplot as plot
2305    global plotter
2306    plotter = plot.PlotNotebook()
2307#    test0()
2308#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
2309    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
2310        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
2311        test2(name,shft)
2312    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
2313        test3(name,shft)
2314    print "OK"
2315    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.