source: trunk/GSASIIpwd.py @ 1459

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

add instrument parameters (flight path & detector 2-theta) needed for TOF
rework reflection list for TOF
change default diff curve & reflection marker offsets
change weighting to instrument constants calibration to be 1/esd2 from peak fit positions - works a lot better
1st shot at TOF Rietveld refinement with derivatives - need to be checked for correctness (some are wrong)

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