source: trunk/GSASIIpwd.py @ 2167

Last change on this file since 2167 was 2167, checked in by vondreele, 7 years ago

clean up of stacking fault stuff

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