source: trunk/GSASIIpwd.py @ 3906

Last change on this file since 3906 was 3906, checked in by toby, 4 years ago

implement fixed backgound file in Riteveld Fit; show Rietveld & peak fits with bkg file added to calc bkg, not subtracted

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