source: trunk/GSASIIpwd.py @ 2183

Last change on this file since 2183 was 2183, checked in by vondreele, 6 years ago

fix General - atomic wt. now updated on isotope change
add new Stacking fault option to make sequence of DIFFaX runs varying one variable - under development
fix problem in G2image & G2imageGUI where change in FlatBkg?, Background or Dark didn't change thresholds or image plot limits - these are now set to new image limits.
Some of these used tricky changes to sizer values

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