source: trunk/GSASIIpwd.py @ 2659

Last change on this file since 2659 was 2659, checked in by vondreele, 5 years ago

make lab data (2 x-ray wavelengths) instrument default 'Bragg-Brentano', all others 'Debye-Scherrer'
refactor PDF stuff to show PDF Controls & (new) PDF Peaks on G2 tree (removing I(Q)...).
Old gpx files with I(Q)... updated automatically to new scheme
Add new tree item for PDF Peaks - does nothing yet.
Fix FWHM calc for TOF so bins/FWHM on peak fitting make sense.

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