source: trunk/GSASIIpwd.py @ 2717

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