source: trunk/GSASIIpwd.py @ 1551

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

fixes to calibration

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