source: trunk/GSASIIpwd.py @ 2755

Last change on this file since 2755 was 2755, checked in by vondreele, 5 years ago

extend Export HKLs for powder data to include most all columns in GSASII.py
narrow HorizontalLine? display
add a ReloadSubstances? option in case user changed wavelength in SASD & REFD substances
Add a default unit scatter to substances - only in new projects
add modeling of REFD patterns for x-rays (not yet tested - no data) & CW neutrons (tested)
modify XScattDen to give clearly fo, f' & f" components
add NCScattDen for CW neutrons - uses wave to generate real/imaginary components
fix printing of scattering density units - now correct symbols
remove unused imports from GSASIIsasd.py

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