source: trunk/GSASIIpwd.py @ 1443

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

add calibration of lam, difC, etc. from index peak positions
new plotCalib routine

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 62.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-07-28 21:53:56 +0000 (Mon, 28 Jul 2014) $
10# $Author: vondreele $
11# $Revision: 1443 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 1443 2014-07-28 21:53:56Z 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: 1443 $")
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,Sq:  S0+S1*dsp**2+Sq*dsp
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-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,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 '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            shl = 0.002
592            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
593            iBeg = np.searchsorted(xdata,pkP-fmin)
594            iFin = np.searchsorted(xdata,pkP+fmax)
595            yb[iBeg:iFin] += pkI*getFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
596            iD += 1       
597        except KeyError:
598            break
599        except ValueError:
600            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
601            break       
602    return yb
603   
604def getBackgroundDerv(pfx,parmDict,bakType,xdata):
605    'needs a doc string'
606    nBak = 0
607    while True:
608        key = pfx+'Back:'+str(nBak)
609        if key in parmDict:
610            nBak += 1
611        else:
612            break
613    dydb = np.zeros(shape=(nBak,len(xdata)))
614    dyddb = np.zeros(shape=(3*parmDict[pfx+'nDebye'],len(xdata)))
615    dydpk = np.zeros(shape=(4*parmDict[pfx+'nPeaks'],len(xdata)))
616    cw = np.diff(xdata)
617    cw = np.append(cw,cw[-1])
618
619    if bakType in ['chebyschev','cosine']:
620        dt = xdata[-1]-xdata[0]   
621        for iBak in range(nBak):   
622            if bakType == 'chebyschev':
623                dydb[iBak] = (2.*(xdata-xdata[0])/dt-1.)**iBak
624            elif bakType == 'cosine':
625                dydb[iBak] = npcosd(xdata*iBak)
626    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
627        if nBak == 1:
628            dydb[0] = np.ones_like(xdata)
629        elif nBak == 2:
630            dX = xdata[-1]-xdata[0]
631            T2 = (xdata-xdata[0])/dX
632            T1 = 1.0-T2
633            dydb = [T1,T2]
634        else:
635            if bakType == 'lin interpolate':
636                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
637            elif bakType == 'inv interpolate':
638                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
639            elif bakType == 'log interpolate':
640                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
641            bakPos[0] = xdata[0]
642            bakPos[-1] = xdata[-1]
643            for i,pos in enumerate(bakPos):
644                if i == 0:
645                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
646                elif i == len(bakPos)-1:
647                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
648                else:
649                    dydb[i] = np.where(xdata>bakPos[i],
650                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
651                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
652    if 'difC' in parmDict:
653        ff = 1.
654    else:
655        try:
656            wave = parmDict[pfx+'Lam']
657        except KeyError:
658            wave = parmDict[pfx+'Lam1']
659        q = 4.0*np.pi*npsind(xdata/2.0)/wave
660        SQ = (q/(4*np.pi))**2
661        FF = G2elem.GetFormFactorCoeff('Si')[0]
662        ff = np.array(G2elem.ScatFac(FF,SQ)[0])
663    iD = 0       
664    while True:
665        try:
666            dbA = parmDict[pfx+'DebyeA:'+str(iD)]
667            dbR = parmDict[pfx+'DebyeR:'+str(iD)]
668            dbU = parmDict[pfx+'DebyeU:'+str(iD)]
669            sqr = np.sin(q*dbR)/(q*dbR)
670            cqr = np.cos(q*dbR)
671            temp = np.exp(-dbU*q**2)
672            dyddb[3*iD] = ff*sqr*temp/(np.pi*cw)
673            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(np.pi*dbR*cw)
674            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2/(np.pi*cw)
675            iD += 1       #ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
676        except KeyError:
677            break
678    iD = 0
679    while True:
680        try:
681            pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
682            pkI = parmDict[pfx+'BkPkint;'+str(iD)]
683            pkS = parmDict[pfx+'BkPksig;'+str(iD)]
684            pkG = parmDict[pfx+'BkPkgam;'+str(iD)]
685            shl = 0.002
686            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
687            iBeg = np.searchsorted(xdata,pkP-fmin)
688            iFin = np.searchsorted(xdata,pkP+fmax)
689            Df,dFdp,dFds,dFdg,dFdsh = getdFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
690            dydpk[4*iD][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdp
691            dydpk[4*iD+1][iBeg:iFin] += 100.*cw[iBeg:iFin]*Df
692            dydpk[4*iD+2][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFds
693            dydpk[4*iD+3][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdg
694            iD += 1       
695        except KeyError:
696            break
697        except ValueError:
698            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
699            break       
700    return dydb,dyddb,dydpk
701
702#use old fortran routine
703def getFCJVoigt3(pos,sig,gam,shl,xdata):
704    'needs a doc string'
705   
706    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
707#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
708    Df /= np.sum(Df)
709    return Df
710
711def getdFCJVoigt3(pos,sig,gam,shl,xdata):
712    'needs a doc string'
713   
714    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
715#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
716    sumDf = np.sum(Df)
717    return Df,dFdp,dFds,dFdg,dFdsh
718
719def getPsVoigt(pos,sig,gam,xdata):
720    'needs a doc string'
721   
722    Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
723    Df /= np.sum(Df)
724    return Df
725
726def getdPsVoigt(pos,sig,gam,xdata):
727    'needs a doc string'
728   
729    Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
730    sumDf = np.sum(Df)
731    return Df,dFdp,dFds,dFdg
732
733def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
734    'needs a doc string'
735    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
736    Df /= np.sum(Df)
737    return Df 
738   
739def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
740    'needs a doc string'
741    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
742    sumDf = np.sum(Df)
743    return Df,dFdp,dFda,dFdb,dFds,dFdg   
744
745def ellipseSize(H,Sij,GB):
746    'needs a doc string'
747    HX = np.inner(H.T,GB)
748    lenHX = np.sqrt(np.sum(HX**2))
749    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
750    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
751    lenR = np.sqrt(np.sum(R**2))
752    return lenR
753
754def ellipseSizeDerv(H,Sij,GB):
755    'needs a doc string'
756    lenR = ellipseSize(H,Sij,GB)
757    delt = 0.001
758    dRdS = np.zeros(6)
759    for i in range(6):
760        dSij = Sij[:]
761        dSij[i] += delt
762        dRdS[i] = (ellipseSize(H,dSij,GB)-lenR)/delt
763    return lenR,dRdS
764
765def getHKLpeak(dmin,SGData,A):
766    'needs a doc string'
767    HKL = G2lat.GenHLaue(dmin,SGData,A)       
768    HKLs = []
769    for h,k,l,d in HKL:
770        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
771        if not ext:
772            HKLs.append([h,k,l,d,-1])
773    return HKLs
774
775def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
776    'needs a doc string'
777   
778    yb = getBackground('',parmDict,bakType,xdata)
779    yc = np.zeros_like(yb)
780    cw = np.diff(xdata)
781    cw = np.append(cw,cw[-1])
782    if 'C' in dataType:
783        shl = max(parmDict['SH/L'],0.002)
784        Ka2 = False
785        if 'Lam1' in parmDict.keys():
786            Ka2 = True
787            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
788            kRatio = parmDict['I(L2)/I(L1)']
789        iPeak = 0
790        while True:
791            try:
792                pos = parmDict['pos'+str(iPeak)]
793                tth = (pos-parmDict['Zero'])
794                intens = parmDict['int'+str(iPeak)]
795                sigName = 'sig'+str(iPeak)
796                if sigName in varyList:
797                    sig = parmDict[sigName]
798                else:
799                    sig = G2mth.getCWsig(parmDict,tth)
800                sig = max(sig,0.001)          #avoid neg sigma
801                gamName = 'gam'+str(iPeak)
802                if gamName in varyList:
803                    gam = parmDict[gamName]
804                else:
805                    gam = G2mth.getCWgam(parmDict,tth)
806                gam = max(gam,0.001)             #avoid neg gamma
807                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
808                iBeg = np.searchsorted(xdata,pos-fmin)
809                iFin = np.searchsorted(xdata,pos+fmin)
810                if not iBeg+iFin:       #peak below low limit
811                    iPeak += 1
812                    continue
813                elif not iBeg-iFin:     #peak above high limit
814                    return yb+yc
815                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
816                if Ka2:
817                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
818                    iBeg = np.searchsorted(xdata,pos2-fmin)
819                    iFin = np.searchsorted(xdata,pos2+fmin)
820                    if iBeg-iFin:
821                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
822                iPeak += 1
823            except KeyError:        #no more peaks to process
824                return yb+yc
825    else:
826        Pdabc = parmDict['Pdabc']
827        difC = parmDict['difC']
828        iPeak = 0
829        while True:
830            try:
831                pos = parmDict['pos'+str(iPeak)]               
832                tof = pos-parmDict['Zero']
833                dsp = tof/difC
834                intens = parmDict['int'+str(iPeak)]
835                alpName = 'alp'+str(iPeak)
836                if alpName in varyList:
837                    alp = parmDict[alpName]
838                else:
839                    if len(Pdabc):
840                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
841                    else:
842                        alp = G2mth.getTOFalpha(parmDict,dsp)
843                betName = 'bet'+str(iPeak)
844                if betName in varyList:
845                    bet = parmDict[betName]
846                else:
847                    if len(Pdabc):
848                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
849                    else:
850                        bet = G2mth.getTOFbeta(parmDict,dsp)
851                sigName = 'sig'+str(iPeak)
852                if sigName in varyList:
853                    sig = parmDict[sigName]
854                else:
855                    sig = G2mth.getTOFsig(parmDict,dsp)
856                gamName = 'gam'+str(iPeak)
857                if gamName in varyList:
858                    gam = parmDict[gamName]
859                else:
860                    gam = G2mth.getTOFgamma(parmDict,dsp)
861                gam = max(gam,0.001)             #avoid neg gamma
862                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
863                iBeg = np.searchsorted(xdata,pos-fmin)
864                iFin = np.searchsorted(xdata,pos+fmax)
865                lenX = len(xdata)
866                if not iBeg:
867                    iFin = np.searchsorted(xdata,pos+fmax)
868                elif iBeg == lenX:
869                    iFin = iBeg
870                else:
871                    iFin = np.searchsorted(xdata,pos+fmax)
872                if not iBeg+iFin:       #peak below low limit
873                    iPeak += 1
874                    continue
875                elif not iBeg-iFin:     #peak above high limit
876                    return yb+yc
877                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
878                iPeak += 1
879            except KeyError:        #no more peaks to process
880                return yb+yc
881           
882def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
883    'needs a doc string'
884# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
885    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
886    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,xdata)
887    if 'Back:0' in varyList:            #background derivs are in front if present
888        dMdv[0:len(dMdb)] = dMdb
889    names = ['DebyeA','DebyeR','DebyeU']
890    for name in varyList:
891        if 'Debye' in name:
892            parm,id = name.split(':')
893            ip = names.index(parm)
894            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
895    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
896    for name in varyList:
897        if 'BkPk' in name:
898            parm,id = name.split(':')
899            ip = names.index(parm)
900            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
901    cw = np.diff(xdata)
902    cw = np.append(cw,cw[-1])
903    if 'C' in dataType:
904        shl = max(parmDict['SH/L'],0.002)
905        Ka2 = False
906        if 'Lam1' in parmDict.keys():
907            Ka2 = True
908            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
909            kRatio = parmDict['I(L2)/I(L1)']
910        iPeak = 0
911        while True:
912            try:
913                pos = parmDict['pos'+str(iPeak)]
914                tth = (pos-parmDict['Zero'])
915                intens = parmDict['int'+str(iPeak)]
916                sigName = 'sig'+str(iPeak)
917                if sigName in varyList:
918                    sig = parmDict[sigName]
919                    dsdU = dsdV = dsdW = 0
920                else:
921                    sig = G2mth.getCWsig(parmDict,tth)
922                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
923                sig = max(sig,0.001)          #avoid neg sigma
924                gamName = 'gam'+str(iPeak)
925                if gamName in varyList:
926                    gam = parmDict[gamName]
927                    dgdX = dgdY = 0
928                else:
929                    gam = G2mth.getCWgam(parmDict,tth)
930                    dgdX,dgdY = G2mth.getCWgamDeriv(tth)
931                gam = max(gam,0.001)             #avoid neg gamma
932                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
933                iBeg = np.searchsorted(xdata,pos-fmin)
934                iFin = np.searchsorted(xdata,pos+fmin)
935                if not iBeg+iFin:       #peak below low limit
936                    iPeak += 1
937                    continue
938                elif not iBeg-iFin:     #peak above high limit
939                    break
940                dMdpk = np.zeros(shape=(6,len(xdata)))
941                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
942                for i in range(1,5):
943                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
944                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
945                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
946                if Ka2:
947                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
948                    iBeg = np.searchsorted(xdata,pos2-fmin)
949                    iFin = np.searchsorted(xdata,pos2+fmin)
950                    if iBeg-iFin:
951                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
952                        for i in range(1,5):
953                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
954                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
955                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
956                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
957                for parmName in ['pos','int','sig','gam']:
958                    try:
959                        idx = varyList.index(parmName+str(iPeak))
960                        dMdv[idx] = dervDict[parmName]
961                    except ValueError:
962                        pass
963                if 'U' in varyList:
964                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
965                if 'V' in varyList:
966                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
967                if 'W' in varyList:
968                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
969                if 'X' in varyList:
970                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
971                if 'Y' in varyList:
972                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
973                if 'SH/L' in varyList:
974                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
975                if 'I(L2)/I(L1)' in varyList:
976                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
977                iPeak += 1
978            except KeyError:        #no more peaks to process
979                break
980    else:
981        Pdabc = parmDict['Pdabc']
982        difC = parmDict['difC']
983        iPeak = 0
984        while True:
985            try:
986                pos = parmDict['pos'+str(iPeak)]               
987                tof = pos-parmDict['Zero']
988                dsp = tof/difC
989                intens = parmDict['int'+str(iPeak)]
990                alpName = 'alp'+str(iPeak)
991                if alpName in varyList:
992                    alp = parmDict[alpName]
993                else:
994                    if len(Pdabc):
995                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
996                        dad0 = 0
997                    else:
998                        alp = G2mth.getTOFalpha(parmDict,dsp)
999                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1000                betName = 'bet'+str(iPeak)
1001                if betName in varyList:
1002                    bet = parmDict[betName]
1003                else:
1004                    if len(Pdabc):
1005                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1006                        dbdb0 = dbdb1 = dbdb2 = 0
1007                    else:
1008                        bet = G2mth.getTOFbeta(parmDict,dsp)
1009                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1010                sigName = 'sig'+str(iPeak)
1011                if sigName in varyList:
1012                    sig = parmDict[sigName]
1013                    dsds0 = dsds1 = dsds2 = 0
1014                else:
1015                    sig = G2mth.getTOFsig(parmDict,dsp)
1016                    dsds0,dsds1,dsds2 = G2mth.getTOFsigDeriv(dsp)
1017                gamName = 'gam'+str(iPeak)
1018                if gamName in varyList:
1019                    gam = parmDict[gamName]
1020                    dsdX = dsdY = 0
1021                else:
1022                    gam = G2mth.getTOFgamma(parmDict,dsp)
1023                    dsdX,dsdY = G2mth.getTOFgammaDeriv(dsp)
1024                gam = max(gam,0.001)             #avoid neg gamma
1025                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1026                iBeg = np.searchsorted(xdata,pos-fmin)
1027                lenX = len(xdata)
1028                if not iBeg:
1029                    iFin = np.searchsorted(xdata,pos+fmax)
1030                elif iBeg == lenX:
1031                    iFin = iBeg
1032                else:
1033                    iFin = np.searchsorted(xdata,pos+fmax)
1034                if not iBeg+iFin:       #peak below low limit
1035                    iPeak += 1
1036                    continue
1037                elif not iBeg-iFin:     #peak above high limit
1038                    break
1039                dMdpk = np.zeros(shape=(7,len(xdata)))
1040                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1041                for i in range(1,6):
1042                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1043                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1044                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1045                for parmName in ['pos','int','alp','bet','sig','gam']:
1046                    try:
1047                        idx = varyList.index(parmName+str(iPeak))
1048                        dMdv[idx] = dervDict[parmName]
1049                    except ValueError:
1050                        pass
1051                if 'alpha' in varyList:
1052                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1053                if 'beta-0' in varyList:
1054                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1055                if 'beta-1' in varyList:
1056                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1057                if 'beta-q' in varyList:
1058                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1059                if 'sig-0' in varyList:
1060                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1061                if 'sig-1' in varyList:
1062                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1063                if 'sig-q' in varyList:
1064                    dMdv[varyList.index('sig-q')] += dsds2*dervDict['sig']
1065                if 'X' in varyList:
1066                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1067                if 'Y' in varyList:
1068                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
1069                iPeak += 1
1070            except KeyError:        #no more peaks to process
1071                break
1072    return dMdv
1073       
1074def Dict2Values(parmdict, varylist):
1075    '''Use before call to leastsq to setup list of values for the parameters
1076    in parmdict, as selected by key in varylist'''
1077    return [parmdict[key] for key in varylist] 
1078   
1079def Values2Dict(parmdict, varylist, values):
1080    ''' Use after call to leastsq to update the parameter dictionary with
1081    values corresponding to keys in varylist'''
1082    parmdict.update(zip(varylist,values))
1083   
1084def SetBackgroundParms(Background):
1085    'needs a doc string'
1086    if len(Background) == 1:            # fix up old backgrounds
1087        BackGround.append({'nDebye':0,'debyeTerms':[]})
1088    bakType,bakFlag = Background[0][:2]
1089    backVals = Background[0][3:]
1090    backNames = ['Back:'+str(i) for i in range(len(backVals))]
1091    Debye = Background[1]           #also has background peaks stuff
1092    backDict = dict(zip(backNames,backVals))
1093    backVary = []
1094    if bakFlag:
1095        backVary = backNames
1096
1097    backDict['nDebye'] = Debye['nDebye']
1098    debyeDict = {}
1099    debyeList = []
1100    for i in range(Debye['nDebye']):
1101        debyeNames = ['DebyeA:'+str(i),'DebyeR:'+str(i),'DebyeU:'+str(i)]
1102        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1103        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1104    debyeVary = []
1105    for item in debyeList:
1106        if item[1]:
1107            debyeVary.append(item[0])
1108    backDict.update(debyeDict)
1109    backVary += debyeVary
1110
1111    backDict['nPeaks'] = Debye['nPeaks']
1112    peaksDict = {}
1113    peaksList = []
1114    for i in range(Debye['nPeaks']):
1115        peaksNames = ['BkPkpos:'+str(i),'BkPkint:'+str(i),'BkPksig:'+str(i),'BkPkgam:'+str(i)]
1116        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1117        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1118    peaksVary = []
1119    for item in peaksList:
1120        if item[1]:
1121            peaksVary.append(item[0])
1122    backDict.update(peaksDict)
1123    backVary += peaksVary   
1124    return bakType,backDict,backVary
1125   
1126def DoCalibInst(IndexPeaks,Inst):
1127   
1128    def SetInstParms():
1129        dataType = Inst['Type'][0]
1130        insVary = []
1131        insNames = []
1132        insVals = []
1133        for parm in Inst:
1134            insNames.append(parm)
1135            insVals.append(Inst[parm][1])
1136            if parm in ['Lam','difC','difA','difB','Zero',]:
1137                if Inst[parm][2]:
1138                    insVary.append(parm)
1139        instDict = dict(zip(insNames,insVals))
1140        return dataType,instDict,insVary
1141       
1142    def GetInstParms(parmDict,Inst,varyList):
1143        for name in Inst:
1144            Inst[name][1] = parmDict[name]
1145       
1146    def InstPrint(Inst,sigDict):
1147        print 'Instrument Parameters:'
1148        if 'C' in Inst['Type'][0]:
1149            ptfmt = "%12.6f"
1150        else:
1151            ptfmt = "%12.3f"
1152        ptlbls = 'names :'
1153        ptstr =  'values:'
1154        sigstr = 'esds  :'
1155        for parm in Inst:
1156            if parm in  ['Lam','difC','difA','difB','Zero',]:
1157                ptlbls += "%s" % (parm.center(12))
1158                ptstr += ptfmt % (Inst[parm][1])
1159                if parm in sigDict:
1160                    sigstr += ptfmt % (sigDict[parm])
1161                else:
1162                    sigstr += 12*' '
1163        print ptlbls
1164        print ptstr
1165        print sigstr
1166       
1167    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1168        parmDict.update(zip(varyList,values))
1169        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1170
1171    peakPos = []
1172    peakDsp = []
1173    peakWt = []
1174    for peak in IndexPeaks:
1175        if peak[2] and peak[3]:
1176            peakPos.append(peak[0])
1177            peakDsp.append(peak[8])
1178            peakWt.append(1/peak[1])
1179    peakPos = np.array(peakPos)
1180    peakDsp = np.array(peakDsp)
1181    peakWt = np.array(peakWt)
1182    dataType,insDict,insVary = SetInstParms()
1183    parmDict = {}
1184    parmDict.update(insDict)
1185    varyList = insVary
1186    while True:
1187        begin = time.time()
1188        values =  np.array(Dict2Values(parmDict, varyList))
1189        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.0001,
1190            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1191        ncyc = int(result[2]['nfev']/2)
1192        runtime = time.time()-begin   
1193        chisq = np.sum(result[2]['fvec']**2)
1194        Values2Dict(parmDict, varyList, result[0])
1195        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
1196        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(peakPos),' Number of parameters: ',len(varyList)
1197        print 'calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1198        print 'chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)
1199        try:
1200            sig = np.sqrt(np.diag(result[1])*GOF)
1201            if np.any(np.isnan(sig)):
1202                print '*** Least squares aborted - some invalid esds possible ***'
1203            break                   #refinement succeeded - finish up!
1204        except ValueError:          #result[1] is None on singular matrix
1205            print '**** Refinement failed - singular matrix ****'
1206        return
1207       
1208    sigDict = dict(zip(varyList,sig))
1209    GetInstParms(parmDict,Inst,varyList)
1210    InstPrint(Inst,sigDict)
1211           
1212def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
1213    'needs a doc string'
1214       
1215    def GetBackgroundParms(parmList,Background):
1216        iBak = 0
1217        while True:
1218            try:
1219                bakName = 'Back:'+str(iBak)
1220                Background[0][iBak+3] = parmList[bakName]
1221                iBak += 1
1222            except KeyError:
1223                break
1224        iDb = 0
1225        while True:
1226            names = ['DebyeA:','DebyeR:','DebyeU:']
1227            try:
1228                for i,name in enumerate(names):
1229                    val = parmList[name+str(iDb)]
1230                    Background[1]['debyeTerms'][iDb][2*i] = val
1231                iDb += 1
1232            except KeyError:
1233                break
1234        iDb = 0
1235        while True:
1236            names = ['BkPkpos:','BkPkint:','BkPksig:','BkPkgam:']
1237            try:
1238                for i,name in enumerate(names):
1239                    val = parmList[name+str(iDb)]
1240                    Background[1]['peaksList'][iDb][2*i] = val
1241                iDb += 1
1242            except KeyError:
1243                break
1244               
1245    def BackgroundPrint(Background,sigDict):
1246        if Background[0][1]:
1247            print 'Background coefficients for',Background[0][0],'function'
1248            ptfmt = "%12.5f"
1249            ptstr =  'values:'
1250            sigstr = 'esds  :'
1251            for i,back in enumerate(Background[0][3:]):
1252                ptstr += ptfmt % (back)
1253                sigstr += ptfmt % (sigDict['Back:'+str(i)])
1254            print ptstr
1255            print sigstr
1256        else:
1257            print 'Background not refined'
1258        if Background[1]['nDebye']:
1259            parms = ['DebyeA','DebyeR','DebyeU']
1260            print 'Debye diffuse scattering coefficients'
1261            ptfmt = "%12.5f"
1262            names =   'names :'
1263            ptstr =  'values:'
1264            sigstr = 'esds  :'
1265            for item in sigDict:
1266                if 'Debye' in item:
1267                    names += '%12s'%(item)
1268                    sigstr += ptfmt%(sigDict[item])
1269                    parm,id = item.split(':')
1270                    ip = parms.index(parm)
1271                    ptstr += ptfmt%(Background[1]['debyeTerms'][int(id)][2*ip])
1272            print names
1273            print ptstr
1274            print sigstr
1275        if Background[1]['nPeaks']:
1276            parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1277            print 'Peaks in background coefficients'
1278            ptfmt = "%12.5f"
1279            names =   'names :'
1280            ptstr =  'values:'
1281            sigstr = 'esds  :'
1282            for item in sigDict:
1283                if 'BkPk' in item:
1284                    names += '%12s'%(item)
1285                    sigstr += ptfmt%(sigDict[item])
1286                    parm,id = item.split(':')
1287                    ip = parms.index(parm)
1288                    ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
1289            print names
1290            print ptstr
1291            print sigstr
1292                           
1293    def SetInstParms(Inst):
1294        dataType = Inst['Type'][0]
1295        insVary = []
1296        insNames = []
1297        insVals = []
1298        for parm in Inst:
1299            insNames.append(parm)
1300            insVals.append(Inst[parm][1])
1301            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1302                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-q',] and Inst[parm][2]:
1303                    insVary.append(parm)
1304        instDict = dict(zip(insNames,insVals))
1305        instDict['X'] = max(instDict['X'],0.01)
1306        instDict['Y'] = max(instDict['Y'],0.01)
1307        if 'SH/L' in instDict:
1308            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1309        return dataType,instDict,insVary
1310       
1311    def GetInstParms(parmDict,Inst,varyList,Peaks):
1312        for name in Inst:
1313            Inst[name][1] = parmDict[name]
1314        iPeak = 0
1315        while True:
1316            try:
1317                sigName = 'sig'+str(iPeak)
1318                pos = parmDict['pos'+str(iPeak)]
1319                if sigName not in varyList:
1320                    if 'C' in Inst['Type'][0]:
1321                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1322                    else:
1323                        dsp = G2lat.Pos2dsp(Inst,pos)
1324                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1325                gamName = 'gam'+str(iPeak)
1326                if gamName not in varyList:
1327                    if 'C' in Inst['Type'][0]:
1328                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1329                    else:
1330                        dsp = G2lat.Pos2dsp(Inst,pos)
1331                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1332                iPeak += 1
1333            except KeyError:
1334                break
1335       
1336    def InstPrint(Inst,sigDict):
1337        print 'Instrument Parameters:'
1338        ptfmt = "%12.6f"
1339        ptlbls = 'names :'
1340        ptstr =  'values:'
1341        sigstr = 'esds  :'
1342        for parm in Inst:
1343            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1344                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-q',]:
1345                ptlbls += "%s" % (parm.center(12))
1346                ptstr += ptfmt % (Inst[parm][1])
1347                if parm in sigDict:
1348                    sigstr += ptfmt % (sigDict[parm])
1349                else:
1350                    sigstr += 12*' '
1351        print ptlbls
1352        print ptstr
1353        print sigstr
1354
1355    def SetPeaksParms(dataType,Peaks):
1356        peakNames = []
1357        peakVary = []
1358        peakVals = []
1359        if 'C' in dataType:
1360            names = ['pos','int','sig','gam']
1361        else:
1362            names = ['pos','int','alp','bet','sig','gam']
1363        for i,peak in enumerate(Peaks):
1364            for j,name in enumerate(names):
1365                peakVals.append(peak[2*j])
1366                parName = name+str(i)
1367                peakNames.append(parName)
1368                if peak[2*j+1]:
1369                    peakVary.append(parName)
1370        return dict(zip(peakNames,peakVals)),peakVary
1371               
1372    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1373        if 'C' in Inst['Type'][0]:
1374            names = ['pos','int','sig','gam']
1375        else:   #'T'
1376            names = ['pos','int','alp','bet','sig','gam']
1377        for i,peak in enumerate(Peaks):
1378            pos = parmDict['pos'+str(i)]
1379            if 'difC' in Inst:
1380                dsp = pos/Inst['difC'][1]
1381            for j in range(len(names)):
1382                parName = names[j]+str(i)
1383                if parName in varyList:
1384                    peak[2*j] = parmDict[parName]
1385                elif 'alpha' in parName:
1386                    peak[2*j] = parmDict['alpha']/dsp
1387                elif 'beta' in parName:
1388                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1389                elif 'sig' in parName:
1390                    if 'C' in Inst['Type'][0]:
1391                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1392                    else:
1393                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1394                elif 'gam' in parName:
1395                    if 'C' in Inst['Type'][0]:
1396                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1397                    else:
1398                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1399                       
1400    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1401        print 'Peak coefficients:'
1402        if 'C' in dataType:
1403            names = ['pos','int','sig','gam']
1404        else:   #'T'
1405            names = ['pos','int','alp','bet','sig','gam']           
1406        head = 13*' '
1407        for name in names:
1408            head += name.center(10)+'esd'.center(10)
1409        print head
1410        if 'C' in dataType:
1411            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1412        else:
1413            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%10.3f",'bet':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
1414        for i,peak in enumerate(Peaks):
1415            ptstr =  ':'
1416            for j in range(len(names)):
1417                name = names[j]
1418                parName = name+str(i)
1419                ptstr += ptfmt[name] % (parmDict[parName])
1420                if parName in varyList:
1421#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
1422                    ptstr += ptfmt[name] % (sigDict[parName])
1423                else:
1424#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
1425                    ptstr += 10*' '
1426            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1427               
1428    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1429        parmdict.update(zip(varylist,values))
1430        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1431           
1432    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
1433        parmdict.update(zip(varylist,values))
1434        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1435        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1436        if dlg:
1437            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1438            if not GoOn:
1439                return -M           #abort!!
1440        return M
1441       
1442    if controls:
1443        Ftol = controls['min dM/M']
1444        derivType = controls['deriv type']
1445    else:
1446        Ftol = 0.0001
1447        derivType = 'analytic'
1448    if oneCycle:
1449        Ftol = 1.0
1450    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1451    yc *= 0.                            #set calcd ones to zero
1452    yb *= 0.
1453    yd *= 0.
1454    xBeg = np.searchsorted(x,Limits[0])
1455    xFin = np.searchsorted(x,Limits[1])
1456    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1457    dataType,insDict,insVary = SetInstParms(Inst)
1458    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1459    parmDict = {}
1460    parmDict.update(bakDict)
1461    parmDict.update(insDict)
1462    parmDict.update(peakDict)
1463    parmDict['Pdabc'] = []      #dummy Pdabc
1464    parmDict.update(Inst2)      #put in real one if there
1465    varyList = bakVary+insVary+peakVary
1466    while True:
1467        begin = time.time()
1468        values =  np.array(Dict2Values(parmDict, varyList))
1469        if FitPgm == 'LSQ':
1470            try:
1471                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1472                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1473                ncyc = int(result[2]['nfev']/2)
1474            finally:
1475                dlg.Destroy()
1476            runtime = time.time()-begin   
1477            chisq = np.sum(result[2]['fvec']**2)
1478            Values2Dict(parmDict, varyList, result[0])
1479            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1480            GOF = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1481            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1482            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1483            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1484            try:
1485                sig = np.sqrt(np.diag(result[1])*GOF)
1486                if np.any(np.isnan(sig)):
1487                    print '*** Least squares aborted - some invalid esds possible ***'
1488                break                   #refinement succeeded - finish up!
1489            except ValueError:          #result[1] is None on singular matrix
1490                print '**** Refinement failed - singular matrix ****'
1491                Ipvt = result[2]['ipvt']
1492                for i,ipvt in enumerate(Ipvt):
1493                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1494                        print 'Removing parameter: ',varyList[ipvt-1]
1495                        del(varyList[ipvt-1])
1496                        break
1497        elif FitPgm == 'BFGS':
1498            print 'Other program here'
1499            return {}
1500       
1501    sigDict = dict(zip(varyList,sig))
1502    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
1503    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1504    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1505    GetBackgroundParms(parmDict,Background)
1506    BackgroundPrint(Background,sigDict)
1507    GetInstParms(parmDict,Inst,varyList,Peaks)
1508    InstPrint(Inst,sigDict)
1509    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1510    PeaksPrint(dataType,parmDict,sigDict,varyList)
1511    return sigDict
1512
1513def calcIncident(Iparm,xdata):
1514    'needs a doc string'
1515
1516    def IfunAdv(Iparm,xdata):
1517        Itype = Iparm['Itype']
1518        Icoef = Iparm['Icoeff']
1519        DYI = np.ones((12,xdata.shape[0]))
1520        YI = np.ones_like(xdata)*Icoef[0]
1521       
1522        x = xdata/1000.                 #expressions are in ms
1523        if Itype == 'Exponential':
1524            for i in range(1,10,2):
1525                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1526                YI += Icoef[i]*Eterm
1527                DYI[i] *= Eterm
1528                DYI[i+1] *= -Icoef[i]*x**((i+1)/2)           
1529        elif 'Maxwell'in Itype:
1530            Eterm = np.exp(-Icoef[2]/x**2)
1531            DYI[1] = Eterm/x**5
1532            DYI[2] = -Icoef[1]*DYI[1]/x**2
1533            YI += (Icoef[1]*Eterm/x**5)
1534            if 'Exponential' in Itype:
1535                for i in range(3,12,2):
1536                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1537                    YI += Icoef[i]*Eterm
1538                    DYI[i] *= Eterm
1539                    DYI[i+1] *= -Icoef[i]*x**((i+1)/2)
1540            else:   #Chebyschev
1541                T = (2./x)-1.
1542                Ccof = np.ones((12,xdata.shape[0]))
1543                Ccof[1] = T
1544                for i in range(2,12):
1545                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1546                for i in range(1,10):
1547                    YI += Ccof[i]*Icoef[i+2]
1548                    DYI[i+2] =Ccof[i]
1549        return YI,DYI
1550       
1551    Iesd = np.array(Iparm['Iesd'])
1552    Icovar = Iparm['Icovar']
1553    YI,DYI = IfunAdv(Iparm,xdata)
1554    YI = np.where(YI>0,YI,1.)
1555    WYI = np.zeros_like(xdata)
1556    vcov = np.zeros((12,12))
1557    k = 0
1558    for i in range(12):
1559        for j in range(i,12):
1560            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1561            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1562            k += 1
1563    M = np.inner(vcov,DYI.T)
1564    WYI = np.sum(M*DYI,axis=0)
1565    WYI = np.where(WYI>0.,WYI,0.)
1566    return YI,WYI
1567   
1568#testing data
1569NeedTestData = True
1570def TestData():
1571    'needs a doc string'
1572#    global NeedTestData
1573    NeedTestData = False
1574    global bakType
1575    bakType = 'chebyschev'
1576    global xdata
1577    xdata = np.linspace(4.0,40.0,36000)
1578    global parmDict0
1579    parmDict0 = {
1580        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1581        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1582        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1583        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1584        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1585        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1586        }
1587    global parmDict1
1588    parmDict1 = {
1589        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1590        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1591        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1592        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1593        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1594        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1595        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1596        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1597        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1598        }
1599    global parmDict2
1600    parmDict2 = {
1601        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1602        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1603        'Back0':5.,'Back1':-0.02,'Back2':.004,
1604#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1605        }
1606    global varyList
1607    varyList = []
1608
1609def test0():
1610    if NeedTestData: TestData()
1611    msg = 'test '
1612    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1613    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
1614    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1615    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1616    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
1617    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1618   
1619def test1():
1620    if NeedTestData: TestData()
1621    time0 = time.time()
1622    for i in range(100):
1623        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1624    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1625   
1626def test2(name,delt):
1627    if NeedTestData: TestData()
1628    varyList = [name,]
1629    xdata = np.linspace(5.6,5.8,400)
1630    hplot = plotter.add('derivatives test for '+name).gca()
1631    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1632    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1633    parmDict2[name] += delt
1634    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1635    hplot.plot(xdata,(y1-y0)/delt,'r+')
1636   
1637def test3(name,delt):
1638    if NeedTestData: TestData()
1639    names = ['pos','sig','gam','shl']
1640    idx = names.index(name)
1641    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1642    xdata = np.linspace(5.6,5.8,800)
1643    dx = xdata[1]-xdata[0]
1644    hplot = plotter.add('derivatives test for '+name).gca()
1645    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1646    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1647    myDict[name] += delt
1648    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1649    hplot.plot(xdata,(y1-y0)/delt,'r+')
1650
1651if __name__ == '__main__':
1652    import GSASIItestplot as plot
1653    global plotter
1654    plotter = plot.PlotNotebook()
1655#    test0()
1656#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1657    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1658        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1659        test2(name,shft)
1660    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1661        test3(name,shft)
1662    print "OK"
1663    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.