source: trunk/GSASIIpwd.py @ 3778

Last change on this file since 3778 was 3778, checked in by vondreele, 3 years ago

correct calculation of FWHM for TOF data & fix the alp values in the exported peaks table
a fix for incommensurate mag str fctr calcs

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 113.4 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-01-10 20:10:58 +0000 (Thu, 10 Jan 2019) $
10# $Author: vondreele $
11# $Revision: 3778 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 3778 2019-01-10 20:10:58Z vondreele $
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: 3778 $")
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):
719    'needs a doc string'
720    if 'T' in dataType:
721        q = 2.*np.pi*parmDict[pfx+'difC']/xdata
722    elif 'C' in dataType:
723        wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
724        q = npT2q(xdata,wave)
725    yb = np.zeros_like(xdata)
726    nBak = 0
727    cw = np.diff(xdata)
728    cw = np.append(cw,cw[-1])
729    sumBk = [0.,0.,0]
730    while True:
731        key = pfx+'Back;'+str(nBak)
732        if key in parmDict:
733            nBak += 1
734        else:
735            break
736#empirical functions
737    if bakType in ['chebyschev','cosine']:
738        dt = xdata[-1]-xdata[0]   
739        for iBak in range(nBak):
740            key = pfx+'Back;'+str(iBak)
741            if bakType == 'chebyschev':
742                ybi = parmDict[key]*(-1.+2.*(xdata-xdata[0])/dt)**iBak
743            elif bakType == 'cosine':
744                ybi = parmDict[key]*npcosd(180.*xdata*iBak/xdata[-1])
745            yb += ybi
746        sumBk[0] = np.sum(yb)
747    elif bakType in ['Q^2 power series','Q^-2 power series']:
748        QT = 1.
749        yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
750        for iBak in range(nBak-1):
751            key = pfx+'Back;'+str(iBak+1)
752            if '-2' in bakType:
753                QT *= (iBak+1)*q**-2
754            else:
755                QT *= q**2/(iBak+1)
756            yb += QT*parmDict[key]
757        sumBk[0] = np.sum(yb)
758    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
759        if nBak == 1:
760            yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
761        elif nBak == 2:
762            dX = xdata[-1]-xdata[0]
763            T2 = (xdata-xdata[0])/dX
764            T1 = 1.0-T2
765            yb = parmDict[pfx+'Back;0']*T1+parmDict[pfx+'Back;1']*T2
766        else:
767            xnomask = ma.getdata(xdata)
768            xmin,xmax = xnomask[0],xnomask[-1]
769            if bakType == 'lin interpolate':
770                bakPos = np.linspace(xmin,xmax,nBak,True)
771            elif bakType == 'inv interpolate':
772                bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
773            elif bakType == 'log interpolate':
774                bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
775            bakPos[0] = xmin
776            bakPos[-1] = xmax
777            bakVals = np.zeros(nBak)
778            for i in range(nBak):
779                bakVals[i] = parmDict[pfx+'Back;'+str(i)]
780            bakInt = si.interp1d(bakPos,bakVals,'linear')
781            yb = bakInt(ma.getdata(xdata))
782        sumBk[0] = np.sum(yb)
783#Debye function       
784    if pfx+'difC' in parmDict:
785        ff = 1.
786    else:       
787        try:
788            wave = parmDict[pfx+'Lam']
789        except KeyError:
790            wave = parmDict[pfx+'Lam1']
791        SQ = (q/(4.*np.pi))**2
792        FF = G2elem.GetFormFactorCoeff('Si')[0]
793        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
794    iD = 0       
795    while True:
796        try:
797            dbA = parmDict[pfx+'DebyeA;'+str(iD)]
798            dbR = parmDict[pfx+'DebyeR;'+str(iD)]
799            dbU = parmDict[pfx+'DebyeU;'+str(iD)]
800            ybi = ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
801            yb += ybi
802            sumBk[1] += np.sum(ybi)
803            iD += 1       
804        except KeyError:
805            break
806#peaks
807    iD = 0
808    while True:
809        try:
810            pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
811            pkI = parmDict[pfx+'BkPkint;'+str(iD)]
812            pkS = parmDict[pfx+'BkPksig;'+str(iD)]
813            pkG = parmDict[pfx+'BkPkgam;'+str(iD)]
814            if 'C' in dataType:
815                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
816            else: #'T'OF
817                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
818            iBeg = np.searchsorted(xdata,pkP-fmin)
819            iFin = np.searchsorted(xdata,pkP+fmax)
820            lenX = len(xdata)
821            if not iBeg:
822                iFin = np.searchsorted(xdata,pkP+fmax)
823            elif iBeg == lenX:
824                iFin = iBeg
825            else:
826                iFin = np.searchsorted(xdata,pkP+fmax)
827            if 'C' in dataType:
828                ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])
829                yb[iBeg:iFin] += ybi
830            else:   #'T'OF
831                ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
832                yb[iBeg:iFin] += ybi
833            sumBk[2] += np.sum(ybi)
834            iD += 1       
835        except KeyError:
836            break
837        except ValueError:
838            print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****')
839            break       
840    return yb,sumBk
841   
842def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata):
843    'needs a doc string'
844    if 'T' in dataType:
845        q = 2.*np.pi*parmDict[hfx+'difC']/xdata
846    elif 'C' in dataType:
847        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
848        q = 2.*np.pi*npsind(xdata/2.)/wave
849    nBak = 0
850    while True:
851        key = hfx+'Back;'+str(nBak)
852        if key in parmDict:
853            nBak += 1
854        else:
855            break
856    dydb = np.zeros(shape=(nBak,len(xdata)))
857    dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
858    dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
859    cw = np.diff(xdata)
860    cw = np.append(cw,cw[-1])
861
862    if bakType in ['chebyschev','cosine']:
863        dt = xdata[-1]-xdata[0]   
864        for iBak in range(nBak):   
865            if bakType == 'chebyschev':
866                dydb[iBak] = (-1.+2.*(xdata-xdata[0])/dt)**iBak
867            elif bakType == 'cosine':
868                dydb[iBak] = npcosd(180.*xdata*iBak/xdata[-1])
869    elif bakType in ['Q^2 power series','Q^-2 power series']:
870        QT = 1.
871        dydb[0] = np.ones_like(xdata)
872        for iBak in range(nBak-1):
873            if '-2' in bakType:
874                QT *= (iBak+1)*q**-2
875            else:
876                QT *= q**2/(iBak+1)
877            dydb[iBak+1] = QT
878    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
879        if nBak == 1:
880            dydb[0] = np.ones_like(xdata)
881        elif nBak == 2:
882            dX = xdata[-1]-xdata[0]
883            T2 = (xdata-xdata[0])/dX
884            T1 = 1.0-T2
885            dydb = [T1,T2]
886        else:
887            xnomask = ma.getdata(xdata)
888            xmin,xmax = xnomask[0],xnomask[-1]
889            if bakType == 'lin interpolate':
890                bakPos = np.linspace(xmin,xmax,nBak,True)
891            elif bakType == 'inv interpolate':
892                bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
893            elif bakType == 'log interpolate':
894                bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
895            bakPos[0] = xmin
896            bakPos[-1] = xmax
897            for i,pos in enumerate(bakPos):
898                if i == 0:
899                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
900                elif i == len(bakPos)-1:
901                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
902                else:
903                    dydb[i] = np.where(xdata>bakPos[i],
904                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
905                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
906    if hfx+'difC' in parmDict:
907        ff = 1.
908    else:
909        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
910        q = npT2q(xdata,wave)
911        SQ = (q/(4*np.pi))**2
912        FF = G2elem.GetFormFactorCoeff('Si')[0]
913        ff = np.array(G2elem.ScatFac(FF,SQ)[0])*np.pi**2    #needs pi^2~10. for cw data (why?)
914    iD = 0       
915    while True:
916        try:
917            if hfx+'difC' in parmDict:
918                q = 2*np.pi*parmDict[hfx+'difC']/xdata
919            dbA = parmDict[hfx+'DebyeA;'+str(iD)]
920            dbR = parmDict[hfx+'DebyeR;'+str(iD)]
921            dbU = parmDict[hfx+'DebyeU;'+str(iD)]
922            sqr = np.sin(q*dbR)/(q*dbR)
923            cqr = np.cos(q*dbR)
924            temp = np.exp(-dbU*q**2)
925            dyddb[3*iD] = ff*sqr*temp
926            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(dbR)
927            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2
928            iD += 1
929        except KeyError:
930            break
931    iD = 0
932    while True:
933        try:
934            pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
935            pkI = parmDict[hfx+'BkPkint;'+str(iD)]
936            pkS = parmDict[hfx+'BkPksig;'+str(iD)]
937            pkG = parmDict[hfx+'BkPkgam;'+str(iD)]
938            if 'C' in dataType:
939                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
940            else: #'T'OF
941                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
942            iBeg = np.searchsorted(xdata,pkP-fmin)
943            iFin = np.searchsorted(xdata,pkP+fmax)
944            lenX = len(xdata)
945            if not iBeg:
946                iFin = np.searchsorted(xdata,pkP+fmax)
947            elif iBeg == lenX:
948                iFin = iBeg
949            else:
950                iFin = np.searchsorted(xdata,pkP+fmax)
951            if 'C' in dataType:
952                Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
953                dydpk[4*iD][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdp
954                dydpk[4*iD+1][iBeg:iFin] += 100.*cw[iBeg:iFin]*Df
955                dydpk[4*iD+2][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFds
956                dydpk[4*iD+3][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdg
957            else:   #'T'OF
958                Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
959                dydpk[4*iD][iBeg:iFin] += pkI*dFdp
960                dydpk[4*iD+1][iBeg:iFin] += Df
961                dydpk[4*iD+2][iBeg:iFin] += pkI*dFds
962                dydpk[4*iD+3][iBeg:iFin] += pkI*dFdg
963            iD += 1       
964        except KeyError:
965            break
966        except ValueError:
967            print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****')
968            break       
969    return dydb,dyddb,dydpk
970
971#use old fortran routine
972def getFCJVoigt3(pos,sig,gam,shl,xdata):
973    '''Compute the Finger-Cox-Jepcoat modified Pseudo-Voigt function for a
974    CW powder peak in external Fortran routine
975    '''
976    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
977#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
978    Df /= np.sum(Df)
979    return Df
980
981def getdFCJVoigt3(pos,sig,gam,shl,xdata):
982    '''Compute analytic derivatives the Finger-Cox-Jepcoat modified Pseudo-Voigt
983    function for a CW powder peak
984    '''
985    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
986#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
987    return Df,dFdp,dFds,dFdg,dFdsh
988
989def getPsVoigt(pos,sig,gam,xdata):
990    'needs a doc string'
991   
992    Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
993    Df /= np.sum(Df)
994    return Df
995
996def getdPsVoigt(pos,sig,gam,xdata):
997    'needs a doc string'
998   
999    Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
1000    return Df,dFdp,dFds,dFdg
1001
1002def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
1003    'needs a doc string'
1004    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1005    Df /= np.sum(Df)
1006    return Df 
1007   
1008def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
1009    'needs a doc string'
1010    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1011    return Df,dFdp,dFda,dFdb,dFds,dFdg   
1012
1013def ellipseSize(H,Sij,GB):
1014    'needs a doc string'
1015    HX = np.inner(H.T,GB)
1016    lenHX = np.sqrt(np.sum(HX**2))
1017    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
1018    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
1019    lenR = np.sqrt(np.sum(R**2))
1020    return lenR
1021
1022def ellipseSizeDerv(H,Sij,GB):
1023    'needs a doc string'
1024    lenR = ellipseSize(H,Sij,GB)
1025    delt = 0.001
1026    dRdS = np.zeros(6)
1027    for i in range(6):
1028        Sij[i] -= delt
1029        lenM = ellipseSize(H,Sij,GB)
1030        Sij[i] += 2.*delt
1031        lenP = ellipseSize(H,Sij,GB)
1032        Sij[i] -= delt
1033        dRdS[i] = (lenP-lenM)/(2.*delt)
1034    return lenR,dRdS
1035
1036def getHKLpeak(dmin,SGData,A,Inst=None):
1037    '''
1038    Generates allowed by symmetry reflections with d >= dmin
1039    NB: GenHKLf & checkMagextc return True for extinct reflections
1040
1041    :param dmin:  minimum d-spacing
1042    :param SGData: space group data obtained from SpcGroup
1043    :param A: lattice parameter terms A1-A6
1044    :param Inst: instrument parameter info
1045    :returns: HKLs: np.array hkl, etc for allowed reflections
1046
1047    '''
1048    HKL = G2lat.GenHLaue(dmin,SGData,A)       
1049    HKLs = []
1050    for h,k,l,d in HKL:
1051        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1052        if ext and 'MagSpGrp' in SGData:
1053            ext = G2spc.checkMagextc([h,k,l],SGData)
1054        if not ext:
1055            if Inst == None:
1056                HKLs.append([h,k,l,d,0,-1])
1057            else:
1058                HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1059    return np.array(HKLs)
1060
1061def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1062    'needs a doc string'
1063    HKLs = []
1064    vec = np.array(Vec)
1065    vstar = np.sqrt(G2lat.calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1066    dvec = 1./(maxH*vstar+1./dmin)
1067    HKL = G2lat.GenHLaue(dvec,SGData,A)       
1068    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1069    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1070    ifMag = False
1071    if 'MagSpGrp' in SGData:
1072        ifMag = True
1073    for h,k,l,d in HKL:
1074        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1075        if not ext and d >= dmin:
1076            HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1077        for dH in SSdH:
1078            if dH:
1079                DH = SSdH[dH]
1080                H = [h+DH[0],k+DH[1],l+DH[2]]
1081                d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1082                if d >= dmin:
1083                    HKLM = np.array([h,k,l,dH])
1084                    if G2spc.checkSSextc(HKLM,SSGData) or ifMag:
1085                        HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])   
1086    return G2lat.sortHKLd(HKLs,True,True,True)
1087
1088def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
1089    'needs a doc string'
1090   
1091    yb = getBackground('',parmDict,bakType,dataType,xdata)[0]
1092    yc = np.zeros_like(yb)
1093    cw = np.diff(xdata)
1094    cw = np.append(cw,cw[-1])
1095    if 'C' in dataType:
1096        shl = max(parmDict['SH/L'],0.002)
1097        Ka2 = False
1098        if 'Lam1' in parmDict.keys():
1099            Ka2 = True
1100            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1101            kRatio = parmDict['I(L2)/I(L1)']
1102        iPeak = 0
1103        while True:
1104            try:
1105                pos = parmDict['pos'+str(iPeak)]
1106                tth = (pos-parmDict['Zero'])
1107                intens = parmDict['int'+str(iPeak)]
1108                sigName = 'sig'+str(iPeak)
1109                if sigName in varyList:
1110                    sig = parmDict[sigName]
1111                else:
1112                    sig = G2mth.getCWsig(parmDict,tth)
1113                sig = max(sig,0.001)          #avoid neg sigma^2
1114                gamName = 'gam'+str(iPeak)
1115                if gamName in varyList:
1116                    gam = parmDict[gamName]
1117                else:
1118                    gam = G2mth.getCWgam(parmDict,tth)
1119                gam = max(gam,0.001)             #avoid neg gamma
1120                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1121                iBeg = np.searchsorted(xdata,pos-fmin)
1122                iFin = np.searchsorted(xdata,pos+fmin)
1123                if not iBeg+iFin:       #peak below low limit
1124                    iPeak += 1
1125                    continue
1126                elif not iBeg-iFin:     #peak above high limit
1127                    return yb+yc
1128                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1129                if Ka2:
1130                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1131                    iBeg = np.searchsorted(xdata,pos2-fmin)
1132                    iFin = np.searchsorted(xdata,pos2+fmin)
1133                    if iBeg-iFin:
1134                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1135                iPeak += 1
1136            except KeyError:        #no more peaks to process
1137                return yb+yc
1138    else:
1139        Pdabc = parmDict['Pdabc']
1140        difC = parmDict['difC']
1141        iPeak = 0
1142        while True:
1143            try:
1144                pos = parmDict['pos'+str(iPeak)]               
1145                tof = pos-parmDict['Zero']
1146                dsp = tof/difC
1147                intens = parmDict['int'+str(iPeak)]
1148                alpName = 'alp'+str(iPeak)
1149                if alpName in varyList:
1150                    alp = parmDict[alpName]
1151                else:
1152                    if len(Pdabc):
1153                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1154                    else:
1155                        alp = G2mth.getTOFalpha(parmDict,dsp)
1156                alp = max(0.1,alp)
1157                betName = 'bet'+str(iPeak)
1158                if betName in varyList:
1159                    bet = parmDict[betName]
1160                else:
1161                    if len(Pdabc):
1162                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1163                    else:
1164                        bet = G2mth.getTOFbeta(parmDict,dsp)
1165                bet = max(0.0001,bet)
1166                sigName = 'sig'+str(iPeak)
1167                if sigName in varyList:
1168                    sig = parmDict[sigName]
1169                else:
1170                    sig = G2mth.getTOFsig(parmDict,dsp)
1171                gamName = 'gam'+str(iPeak)
1172                if gamName in varyList:
1173                    gam = parmDict[gamName]
1174                else:
1175                    gam = G2mth.getTOFgamma(parmDict,dsp)
1176                gam = max(gam,0.001)             #avoid neg gamma
1177                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1178                iBeg = np.searchsorted(xdata,pos-fmin)
1179                iFin = np.searchsorted(xdata,pos+fmax)
1180                lenX = len(xdata)
1181                if not iBeg:
1182                    iFin = np.searchsorted(xdata,pos+fmax)
1183                elif iBeg == lenX:
1184                    iFin = iBeg
1185                else:
1186                    iFin = np.searchsorted(xdata,pos+fmax)
1187                if not iBeg+iFin:       #peak below low limit
1188                    iPeak += 1
1189                    continue
1190                elif not iBeg-iFin:     #peak above high limit
1191                    return yb+yc
1192                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1193                iPeak += 1
1194            except KeyError:        #no more peaks to process
1195                return yb+yc
1196           
1197def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
1198    'needs a doc string'
1199# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1200    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1201    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,dataType,xdata)
1202    if 'Back;0' in varyList:            #background derivs are in front if present
1203        dMdv[0:len(dMdb)] = dMdb
1204    names = ['DebyeA','DebyeR','DebyeU']
1205    for name in varyList:
1206        if 'Debye' in name:
1207            parm,id = name.split(';')
1208            ip = names.index(parm)
1209            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
1210    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1211    for name in varyList:
1212        if 'BkPk' in name:
1213            parm,id = name.split(';')
1214            ip = names.index(parm)
1215            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
1216    cw = np.diff(xdata)
1217    cw = np.append(cw,cw[-1])
1218    if 'C' in dataType:
1219        shl = max(parmDict['SH/L'],0.002)
1220        Ka2 = False
1221        if 'Lam1' in parmDict.keys():
1222            Ka2 = True
1223            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1224            kRatio = parmDict['I(L2)/I(L1)']
1225        iPeak = 0
1226        while True:
1227            try:
1228                pos = parmDict['pos'+str(iPeak)]
1229                tth = (pos-parmDict['Zero'])
1230                intens = parmDict['int'+str(iPeak)]
1231                sigName = 'sig'+str(iPeak)
1232                if sigName in varyList:
1233                    sig = parmDict[sigName]
1234                    dsdU = dsdV = dsdW = 0
1235                else:
1236                    sig = G2mth.getCWsig(parmDict,tth)
1237                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1238                sig = max(sig,0.001)          #avoid neg sigma
1239                gamName = 'gam'+str(iPeak)
1240                if gamName in varyList:
1241                    gam = parmDict[gamName]
1242                    dgdX = dgdY = dgdZ = 0
1243                else:
1244                    gam = G2mth.getCWgam(parmDict,tth)
1245                    dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1246                gam = max(gam,0.001)             #avoid neg gamma
1247                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1248                iBeg = np.searchsorted(xdata,pos-fmin)
1249                iFin = np.searchsorted(xdata,pos+fmin)
1250                if not iBeg+iFin:       #peak below low limit
1251                    iPeak += 1
1252                    continue
1253                elif not iBeg-iFin:     #peak above high limit
1254                    break
1255                dMdpk = np.zeros(shape=(6,len(xdata)))
1256                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1257                for i in range(1,5):
1258                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
1259                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
1260                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1261                if Ka2:
1262                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1263                    iBeg = np.searchsorted(xdata,pos2-fmin)
1264                    iFin = np.searchsorted(xdata,pos2+fmin)
1265                    if iBeg-iFin:
1266                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1267                        for i in range(1,5):
1268                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
1269                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
1270                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
1271                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1272                for parmName in ['pos','int','sig','gam']:
1273                    try:
1274                        idx = varyList.index(parmName+str(iPeak))
1275                        dMdv[idx] = dervDict[parmName]
1276                    except ValueError:
1277                        pass
1278                if 'U' in varyList:
1279                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1280                if 'V' in varyList:
1281                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1282                if 'W' in varyList:
1283                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1284                if 'X' in varyList:
1285                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1286                if 'Y' in varyList:
1287                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1288                if 'Z' in varyList:
1289                    dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1290                if 'SH/L' in varyList:
1291                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
1292                if 'I(L2)/I(L1)' in varyList:
1293                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1294                iPeak += 1
1295            except KeyError:        #no more peaks to process
1296                break
1297    else:
1298        Pdabc = parmDict['Pdabc']
1299        difC = parmDict['difC']
1300        iPeak = 0
1301        while True:
1302            try:
1303                pos = parmDict['pos'+str(iPeak)]               
1304                tof = pos-parmDict['Zero']
1305                dsp = tof/difC
1306                intens = parmDict['int'+str(iPeak)]
1307                alpName = 'alp'+str(iPeak)
1308                if alpName in varyList:
1309                    alp = parmDict[alpName]
1310                else:
1311                    if len(Pdabc):
1312                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1313                        dada0 = 0
1314                    else:
1315                        alp = G2mth.getTOFalpha(parmDict,dsp)
1316                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1317                betName = 'bet'+str(iPeak)
1318                if betName in varyList:
1319                    bet = parmDict[betName]
1320                else:
1321                    if len(Pdabc):
1322                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1323                        dbdb0 = dbdb1 = dbdb2 = 0
1324                    else:
1325                        bet = G2mth.getTOFbeta(parmDict,dsp)
1326                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1327                sigName = 'sig'+str(iPeak)
1328                if sigName in varyList:
1329                    sig = parmDict[sigName]
1330                    dsds0 = dsds1 = dsds2 = dsds3 = 0
1331                else:
1332                    sig = G2mth.getTOFsig(parmDict,dsp)
1333                    dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1334                gamName = 'gam'+str(iPeak)
1335                if gamName in varyList:
1336                    gam = parmDict[gamName]
1337                    dsdX = dsdY = dsdZ = 0
1338                else:
1339                    gam = G2mth.getTOFgamma(parmDict,dsp)
1340                    dsdX,dsdY,dsdZ = G2mth.getTOFgammaDeriv(dsp)
1341                gam = max(gam,0.001)             #avoid neg gamma
1342                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1343                iBeg = np.searchsorted(xdata,pos-fmin)
1344                lenX = len(xdata)
1345                if not iBeg:
1346                    iFin = np.searchsorted(xdata,pos+fmax)
1347                elif iBeg == lenX:
1348                    iFin = iBeg
1349                else:
1350                    iFin = np.searchsorted(xdata,pos+fmax)
1351                if not iBeg+iFin:       #peak below low limit
1352                    iPeak += 1
1353                    continue
1354                elif not iBeg-iFin:     #peak above high limit
1355                    break
1356                dMdpk = np.zeros(shape=(7,len(xdata)))
1357                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1358                for i in range(1,6):
1359                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1360                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1361                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1362                for parmName in ['pos','int','alp','bet','sig','gam']:
1363                    try:
1364                        idx = varyList.index(parmName+str(iPeak))
1365                        dMdv[idx] = dervDict[parmName]
1366                    except ValueError:
1367                        pass
1368                if 'alpha' in varyList:
1369                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1370                if 'beta-0' in varyList:
1371                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1372                if 'beta-1' in varyList:
1373                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1374                if 'beta-q' in varyList:
1375                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1376                if 'sig-0' in varyList:
1377                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1378                if 'sig-1' in varyList:
1379                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1380                if 'sig-2' in varyList:
1381                    dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1382                if 'sig-q' in varyList:
1383                    dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1384                if 'X' in varyList:
1385                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1386                if 'Y' in varyList:
1387                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']
1388                if 'Z' in varyList:
1389                    dMdv[varyList.index('Z')] += dsdZ*dervDict['gam']
1390                iPeak += 1
1391            except KeyError:        #no more peaks to process
1392                break
1393    return dMdv
1394       
1395def Dict2Values(parmdict, varylist):
1396    '''Use before call to leastsq to setup list of values for the parameters
1397    in parmdict, as selected by key in varylist'''
1398    return [parmdict[key] for key in varylist] 
1399   
1400def Values2Dict(parmdict, varylist, values):
1401    ''' Use after call to leastsq to update the parameter dictionary with
1402    values corresponding to keys in varylist'''
1403    parmdict.update(zip(varylist,values))
1404   
1405def SetBackgroundParms(Background):
1406    'needs a doc string'
1407    if len(Background) == 1:            # fix up old backgrounds
1408        Background.append({'nDebye':0,'debyeTerms':[]})
1409    bakType,bakFlag = Background[0][:2]
1410    backVals = Background[0][3:]
1411    backNames = ['Back;'+str(i) for i in range(len(backVals))]
1412    Debye = Background[1]           #also has background peaks stuff
1413    backDict = dict(zip(backNames,backVals))
1414    backVary = []
1415    if bakFlag:
1416        backVary = backNames
1417
1418    backDict['nDebye'] = Debye['nDebye']
1419    debyeDict = {}
1420    debyeList = []
1421    for i in range(Debye['nDebye']):
1422        debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1423        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1424        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1425    debyeVary = []
1426    for item in debyeList:
1427        if item[1]:
1428            debyeVary.append(item[0])
1429    backDict.update(debyeDict)
1430    backVary += debyeVary
1431
1432    backDict['nPeaks'] = Debye['nPeaks']
1433    peaksDict = {}
1434    peaksList = []
1435    for i in range(Debye['nPeaks']):
1436        peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1437        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1438        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1439    peaksVary = []
1440    for item in peaksList:
1441        if item[1]:
1442            peaksVary.append(item[0])
1443    backDict.update(peaksDict)
1444    backVary += peaksVary   
1445    return bakType,backDict,backVary
1446   
1447def DoCalibInst(IndexPeaks,Inst):
1448   
1449    def SetInstParms():
1450        dataType = Inst['Type'][0]
1451        insVary = []
1452        insNames = []
1453        insVals = []
1454        for parm in Inst:
1455            insNames.append(parm)
1456            insVals.append(Inst[parm][1])
1457            if parm in ['Lam','difC','difA','difB','Zero',]:
1458                if Inst[parm][2]:
1459                    insVary.append(parm)
1460        instDict = dict(zip(insNames,insVals))
1461        return dataType,instDict,insVary
1462       
1463    def GetInstParms(parmDict,Inst,varyList):
1464        for name in Inst:
1465            Inst[name][1] = parmDict[name]
1466       
1467    def InstPrint(Inst,sigDict):
1468        print ('Instrument Parameters:')
1469        if 'C' in Inst['Type'][0]:
1470            ptfmt = "%12.6f"
1471        else:
1472            ptfmt = "%12.3f"
1473        ptlbls = 'names :'
1474        ptstr =  'values:'
1475        sigstr = 'esds  :'
1476        for parm in Inst:
1477            if parm in  ['Lam','difC','difA','difB','Zero',]:
1478                ptlbls += "%s" % (parm.center(12))
1479                ptstr += ptfmt % (Inst[parm][1])
1480                if parm in sigDict:
1481                    sigstr += ptfmt % (sigDict[parm])
1482                else:
1483                    sigstr += 12*' '
1484        print (ptlbls)
1485        print (ptstr)
1486        print (sigstr)
1487       
1488    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1489        parmDict.update(zip(varyList,values))
1490        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1491
1492    peakPos = []
1493    peakDsp = []
1494    peakWt = []
1495    for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1496        if peak[2] and peak[3] and sig > 0.:
1497            peakPos.append(peak[0])
1498            peakDsp.append(peak[-1])    #d-calc
1499#            peakWt.append(peak[-1]**2/sig**2)   #weight by d**2
1500            peakWt.append(1./(sig*peak[-1]))   #
1501    peakPos = np.array(peakPos)
1502    peakDsp = np.array(peakDsp)
1503    peakWt = np.array(peakWt)
1504    dataType,insDict,insVary = SetInstParms()
1505    parmDict = {}
1506    parmDict.update(insDict)
1507    varyList = insVary
1508    if not len(varyList):
1509        print ('**** ERROR - nothing to refine! ****')
1510        return False
1511    while True:
1512        begin = time.time()
1513        values =  np.array(Dict2Values(parmDict, varyList))
1514        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1515            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1516        ncyc = int(result[2]['nfev']/2)
1517        runtime = time.time()-begin   
1518        chisq = np.sum(result[2]['fvec']**2)
1519        Values2Dict(parmDict, varyList, result[0])
1520        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
1521        print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(peakPos),len(varyList)))
1522        print ('calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
1523        print ('chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF))
1524        try:
1525            sig = np.sqrt(np.diag(result[1])*GOF)
1526            if np.any(np.isnan(sig)):
1527                print ('*** Least squares aborted - some invalid esds possible ***')
1528            break                   #refinement succeeded - finish up!
1529        except ValueError:          #result[1] is None on singular matrix
1530            print ('**** Refinement failed - singular matrix ****')
1531       
1532    sigDict = dict(zip(varyList,sig))
1533    GetInstParms(parmDict,Inst,varyList)
1534    InstPrint(Inst,sigDict)
1535    return True
1536           
1537def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
1538    '''Called to perform a peak fit, refining the selected items in the peak
1539    table as well as selected items in the background.
1540
1541    :param str FitPgm: type of fit to perform. At present this is ignored.
1542    :param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1543      four values followed by a refine flag where the values are: position, intensity,
1544      sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1545      tree entry, dict item "peaks"
1546    :param list Background: describes the background. List with two items.
1547      Item 0 specifies a background model and coefficients. Item 1 is a dict.
1548      From the Histogram/Background tree entry.
1549    :param list Limits: min and max x-value to use
1550    :param dict Inst: Instrument parameters
1551    :param dict Inst2: more Instrument parameters
1552    :param numpy.array data: a 5xn array. data[0] is the x-values,
1553      data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1554      calc, background and difference intensities, respectively.
1555    :param array fixback: fixed background values
1556    :param list prevVaryList: Used in sequential refinements to override the
1557      variable list. Defaults as an empty list.
1558    :param bool oneCycle: True if only one cycle of fitting should be performed
1559    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1560      and derivType = controls['deriv type']. If None default values are used.
1561    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1562      Defaults to None, which means no updates are done.
1563    '''
1564    def GetBackgroundParms(parmList,Background):
1565        iBak = 0
1566        while True:
1567            try:
1568                bakName = 'Back;'+str(iBak)
1569                Background[0][iBak+3] = parmList[bakName]
1570                iBak += 1
1571            except KeyError:
1572                break
1573        iDb = 0
1574        while True:
1575            names = ['DebyeA;','DebyeR;','DebyeU;']
1576            try:
1577                for i,name in enumerate(names):
1578                    val = parmList[name+str(iDb)]
1579                    Background[1]['debyeTerms'][iDb][2*i] = val
1580                iDb += 1
1581            except KeyError:
1582                break
1583        iDb = 0
1584        while True:
1585            names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1586            try:
1587                for i,name in enumerate(names):
1588                    val = parmList[name+str(iDb)]
1589                    Background[1]['peaksList'][iDb][2*i] = val
1590                iDb += 1
1591            except KeyError:
1592                break
1593               
1594    def BackgroundPrint(Background,sigDict):
1595        print ('Background coefficients for '+Background[0][0]+' function')
1596        ptfmt = "%12.5f"
1597        ptstr =  'value: '
1598        sigstr = 'esd  : '
1599        for i,back in enumerate(Background[0][3:]):
1600            ptstr += ptfmt % (back)
1601            if Background[0][1]:
1602                prm = 'Back;'+str(i)
1603                if prm in sigDict:
1604                    sigstr += ptfmt % (sigDict[prm])
1605                else:
1606                    sigstr += " "*12
1607            if len(ptstr) > 75:
1608                print (ptstr)
1609                if Background[0][1]: print (sigstr)
1610                ptstr =  'value: '
1611                sigstr = 'esd  : '
1612        if len(ptstr) > 8:
1613            print (ptstr)
1614            if Background[0][1]: print (sigstr)
1615
1616        if Background[1]['nDebye']:
1617            parms = ['DebyeA;','DebyeR;','DebyeU;']
1618            print ('Debye diffuse scattering coefficients')
1619            ptfmt = "%12.5f"
1620            print (' term       DebyeA       esd        DebyeR       esd        DebyeU        esd')
1621            for term in range(Background[1]['nDebye']):
1622                line = ' term %d'%(term)
1623                for ip,name in enumerate(parms):
1624                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1625                    if name+str(term) in sigDict:
1626                        line += ptfmt%(sigDict[name+str(term)])
1627                    else:
1628                        line += " "*12
1629                print (line)
1630        if Background[1]['nPeaks']:
1631            print ('Coefficients for Background Peaks')
1632            ptfmt = "%15.3f"
1633            for j,pl in enumerate(Background[1]['peaksList']):
1634                names =  'peak %3d:'%(j+1)
1635                ptstr =  'values  :'
1636                sigstr = 'esds    :'
1637                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1638                    val = pl[2*i]
1639                    prm = lbl+";"+str(j)
1640                    names += '%15s'%(prm)
1641                    ptstr += ptfmt%(val)
1642                    if prm in sigDict:
1643                        sigstr += ptfmt%(sigDict[prm])
1644                    else:
1645                        sigstr += " "*15
1646                print (names)
1647                print (ptstr)
1648                print (sigstr)
1649                           
1650    def SetInstParms(Inst):
1651        dataType = Inst['Type'][0]
1652        insVary = []
1653        insNames = []
1654        insVals = []
1655        for parm in Inst:
1656            insNames.append(parm)
1657            insVals.append(Inst[parm][1])
1658            if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
1659                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1660                    insVary.append(parm)
1661        instDict = dict(zip(insNames,insVals))
1662#        instDict['X'] = max(instDict['X'],0.01)
1663#        instDict['Y'] = max(instDict['Y'],0.01)
1664        if 'SH/L' in instDict:
1665            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1666        return dataType,instDict,insVary
1667       
1668    def GetInstParms(parmDict,Inst,varyList,Peaks):
1669        for name in Inst:
1670            Inst[name][1] = parmDict[name]
1671        iPeak = 0
1672        while True:
1673            try:
1674                sigName = 'sig'+str(iPeak)
1675                pos = parmDict['pos'+str(iPeak)]
1676                if sigName not in varyList:
1677                    if 'C' in Inst['Type'][0]:
1678                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1679                    else:
1680                        dsp = G2lat.Pos2dsp(Inst,pos)
1681                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1682                gamName = 'gam'+str(iPeak)
1683                if gamName not in varyList:
1684                    if 'C' in Inst['Type'][0]:
1685                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1686                    else:
1687                        dsp = G2lat.Pos2dsp(Inst,pos)
1688                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1689                iPeak += 1
1690            except KeyError:
1691                break
1692       
1693    def InstPrint(Inst,sigDict):
1694        print ('Instrument Parameters:')
1695        ptfmt = "%12.6f"
1696        ptlbls = 'names :'
1697        ptstr =  'values:'
1698        sigstr = 'esds  :'
1699        for parm in Inst:
1700            if parm in  ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
1701                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1702                ptlbls += "%s" % (parm.center(12))
1703                ptstr += ptfmt % (Inst[parm][1])
1704                if parm in sigDict:
1705                    sigstr += ptfmt % (sigDict[parm])
1706                else:
1707                    sigstr += 12*' '
1708        print (ptlbls)
1709        print (ptstr)
1710        print (sigstr)
1711
1712    def SetPeaksParms(dataType,Peaks):
1713        peakNames = []
1714        peakVary = []
1715        peakVals = []
1716        if 'C' in dataType:
1717            names = ['pos','int','sig','gam']
1718        else:
1719            names = ['pos','int','alp','bet','sig','gam']
1720        for i,peak in enumerate(Peaks):
1721            for j,name in enumerate(names):
1722                peakVals.append(peak[2*j])
1723                parName = name+str(i)
1724                peakNames.append(parName)
1725                if peak[2*j+1]:
1726                    peakVary.append(parName)
1727        return dict(zip(peakNames,peakVals)),peakVary
1728               
1729    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1730        if 'C' in Inst['Type'][0]:
1731            names = ['pos','int','sig','gam']
1732        else:   #'T'
1733            names = ['pos','int','alp','bet','sig','gam']
1734        for i,peak in enumerate(Peaks):
1735            pos = parmDict['pos'+str(i)]
1736            if 'difC' in Inst:
1737                dsp = pos/Inst['difC'][1]
1738            for j in range(len(names)):
1739                parName = names[j]+str(i)
1740                if parName in varyList:
1741                    peak[2*j] = parmDict[parName]
1742                elif 'alpha' in parName:
1743                    peak[2*j] = parmDict['alpha']/dsp
1744                elif 'beta' in parName:
1745                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1746                elif 'sig' in parName:
1747                    if 'C' in Inst['Type'][0]:
1748                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1749                    else:
1750                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1751                elif 'gam' in parName:
1752                    if 'C' in Inst['Type'][0]:
1753                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1754                    else:
1755                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1756                       
1757    def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
1758        print ('Peak coefficients:')
1759        if 'C' in dataType:
1760            names = ['pos','int','sig','gam']
1761        else:   #'T'
1762            names = ['pos','int','alp','bet','sig','gam']           
1763        head = 13*' '
1764        for name in names:
1765            if name in ['alp','bet']:
1766                head += name.center(8)+'esd'.center(8)
1767            else:
1768                head += name.center(10)+'esd'.center(10)
1769        head += 'bins'.center(8)
1770        print (head)
1771        if 'C' in dataType:
1772            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1773        else:
1774            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1775        for i,peak in enumerate(Peaks):
1776            ptstr =  ':'
1777            for j in range(len(names)):
1778                name = names[j]
1779                parName = name+str(i)
1780                ptstr += ptfmt[name] % (parmDict[parName])
1781                if parName in varyList:
1782                    ptstr += ptfmt[name] % (sigDict[parName])
1783                else:
1784                    if name in ['alp','bet']:
1785                        ptstr += 8*' '
1786                    else:
1787                        ptstr += 10*' '
1788            ptstr += '%9.2f'%(ptsperFW[i])
1789            print ('%s'%(('Peak'+str(i+1)).center(8)),ptstr)
1790               
1791    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1792        parmdict.update(zip(varylist,values))
1793        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1794           
1795    def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):       
1796        parmdict.update(zip(varylist,values))
1797        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1798        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1799        if dlg:
1800            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1801            if not GoOn:
1802                return -M           #abort!!
1803        return M
1804       
1805    if controls:
1806        Ftol = controls['min dM/M']
1807    else:
1808        Ftol = 0.0001
1809    if oneCycle:
1810        Ftol = 1.0
1811    x,y,w,yc,yb,yd = data   #these are numpy arrays - remove masks!
1812    yc *= 0.                            #set calcd ones to zero
1813    yb *= 0.
1814    yd *= 0.
1815    cw = x[1:]-x[:-1]
1816    xBeg = np.searchsorted(x,Limits[0])
1817    xFin = np.searchsorted(x,Limits[1])+1
1818    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1819    dataType,insDict,insVary = SetInstParms(Inst)
1820    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1821    parmDict = {}
1822    parmDict.update(bakDict)
1823    parmDict.update(insDict)
1824    parmDict.update(peakDict)
1825    parmDict['Pdabc'] = []      #dummy Pdabc
1826    parmDict.update(Inst2)      #put in real one if there
1827    if prevVaryList:
1828        varyList = prevVaryList[:]
1829    else:
1830        varyList = bakVary+insVary+peakVary
1831    fullvaryList = varyList[:]
1832    while True:
1833        begin = time.time()
1834        values =  np.array(Dict2Values(parmDict, varyList))
1835        Rvals = {}
1836        badVary = []
1837        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1838               args=(x[xBeg:xFin],(y+fixback)[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1839        ncyc = int(result[2]['nfev']/2)
1840        runtime = time.time()-begin   
1841        chisq = np.sum(result[2]['fvec']**2)
1842        Values2Dict(parmDict, varyList, result[0])
1843        Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*(y+fixback)[xBeg:xFin]**2))*100.      #to %
1844        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1845        print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList)))
1846        if ncyc:
1847            print ('fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
1848        print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
1849        sig = [0]*len(varyList)
1850        if len(varyList) == 0: break  # if nothing was refined
1851        try:
1852            sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1853            if np.any(np.isnan(sig)):
1854                print ('*** Least squares aborted - some invalid esds possible ***')
1855            break                   #refinement succeeded - finish up!
1856        except ValueError:          #result[1] is None on singular matrix
1857            print ('**** Refinement failed - singular matrix ****')
1858            Ipvt = result[2]['ipvt']
1859            for i,ipvt in enumerate(Ipvt):
1860                if not np.sum(result[2]['fjac'],axis=1)[i]:
1861                    print ('Removing parameter: '+varyList[ipvt-1])
1862                    badVary.append(varyList[ipvt-1])
1863                    del(varyList[ipvt-1])
1864                    break
1865            else: # nothing removed
1866                break
1867    if dlg: dlg.Destroy()
1868    sigDict = dict(zip(varyList,sig))
1869    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]
1870    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1871    yd[xBeg:xFin] = (y+fixback)[xBeg:xFin]-yc[xBeg:xFin]
1872    GetBackgroundParms(parmDict,Background)
1873    if bakVary: BackgroundPrint(Background,sigDict)
1874    GetInstParms(parmDict,Inst,varyList,Peaks)
1875    if insVary: InstPrint(Inst,sigDict)
1876    GetPeaksParms(Inst,parmDict,Peaks,varyList)
1877    binsperFWHM = []
1878    for peak in Peaks:
1879        FWHM = getFWHM(peak[0],Inst)
1880        try:
1881            binsperFWHM.append(FWHM/cw[x.searchsorted(peak[0])])
1882        except IndexError:
1883            binsperFWHM.append(0.)
1884    if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
1885    if len(binsperFWHM):
1886        if min(binsperFWHM) < 1.:
1887            print ('*** Warning: calculated peak widths are too narrow to refine profile coefficients ***')
1888            if 'T' in Inst['Type'][0]:
1889                print (' Manually increase sig-0, 1, or 2 in Instrument Parameters')
1890            else:
1891                print (' Manually increase W in Instrument Parameters')
1892        elif min(binsperFWHM) < 4.:
1893            print ('*** Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement ***')
1894            print ('*** recommended is 6-10; you have %.2f ***'%(min(binsperFWHM)))
1895    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
1896   
1897def calcIncident(Iparm,xdata):
1898    'needs a doc string'
1899
1900    def IfunAdv(Iparm,xdata):
1901        Itype = Iparm['Itype']
1902        Icoef = Iparm['Icoeff']
1903        DYI = np.ones((12,xdata.shape[0]))
1904        YI = np.ones_like(xdata)*Icoef[0]
1905       
1906        x = xdata/1000.                 #expressions are in ms
1907        if Itype == 'Exponential':
1908            for i in [1,3,5,7,9]:
1909                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1910                YI += Icoef[i]*Eterm
1911                DYI[i] *= Eterm
1912                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
1913        elif 'Maxwell'in Itype:
1914            Eterm = np.exp(-Icoef[2]/x**2)
1915            DYI[1] = Eterm/x**5
1916            DYI[2] = -Icoef[1]*DYI[1]/x**2
1917            YI += (Icoef[1]*Eterm/x**5)
1918            if 'Exponential' in Itype:
1919                for i in range(3,11,2):
1920                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1921                    YI += Icoef[i]*Eterm
1922                    DYI[i] *= Eterm
1923                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
1924            else:   #Chebyschev
1925                T = (2./x)-1.
1926                Ccof = np.ones((12,xdata.shape[0]))
1927                Ccof[1] = T
1928                for i in range(2,12):
1929                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1930                for i in range(1,10):
1931                    YI += Ccof[i]*Icoef[i+2]
1932                    DYI[i+2] =Ccof[i]
1933        return YI,DYI
1934       
1935    Iesd = np.array(Iparm['Iesd'])
1936    Icovar = Iparm['Icovar']
1937    YI,DYI = IfunAdv(Iparm,xdata)
1938    YI = np.where(YI>0,YI,1.)
1939    WYI = np.zeros_like(xdata)
1940    vcov = np.zeros((12,12))
1941    k = 0
1942    for i in range(12):
1943        for j in range(i,12):
1944            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1945            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1946            k += 1
1947    M = np.inner(vcov,DYI.T)
1948    WYI = np.sum(M*DYI,axis=0)
1949    WYI = np.where(WYI>0.,WYI,0.)
1950    return YI,WYI
1951   
1952################################################################################
1953# Reflectometry calculations
1954################################################################################
1955
1956def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
1957    print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
1958   
1959    class RandomDisplacementBounds(object):
1960        """random displacement with bounds"""
1961        def __init__(self, xmin, xmax, stepsize=0.5):
1962            self.xmin = xmin
1963            self.xmax = xmax
1964            self.stepsize = stepsize
1965   
1966        def __call__(self, x):
1967            """take a random step but ensure the new position is within the bounds"""
1968            while True:
1969                # this could be done in a much more clever way, but it will work for example purposes
1970                steps = self.xmax-self.xmin
1971                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
1972                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
1973                    break
1974            return xnew
1975   
1976    def GetModelParms():
1977        parmDict = {}
1978        varyList = []
1979        values = []
1980        bounds = []
1981        parmDict['dQ type'] = data['dQ type']
1982        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
1983        for parm in ['Scale','FltBack']:
1984            parmDict[parm] = data[parm][0]
1985            if data[parm][1]:
1986                varyList.append(parm)
1987                values.append(data[parm][0])
1988                bounds.append(Bounds[parm])
1989        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
1990        parmDict['nLayers'] = len(parmDict['Layer Seq'])
1991        for ilay,layer in enumerate(data['Layers']):
1992            name = layer['Name']
1993            cid = str(ilay)+';'
1994            parmDict[cid+'Name'] = name
1995            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
1996                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
1997                if layer.get(parm,[0.,False])[1]:
1998                    varyList.append(cid+parm)
1999                    value = layer[parm][0]
2000                    values.append(value)
2001                    if value:
2002                        bound = [value*Bfac,value/Bfac]
2003                    else:
2004                        bound = [0.,10.]
2005                    bounds.append(bound)
2006            if name not in ['vacuum','unit scatter']:
2007                parmDict[cid+'rho'] = Substances[name]['Scatt density']
2008                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
2009        return parmDict,varyList,values,bounds
2010   
2011    def SetModelParms():
2012        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
2013        if 'Scale' in varyList:
2014            data['Scale'][0] = parmDict['Scale']
2015            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
2016        print (line)
2017        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
2018        if 'FltBack' in varyList:
2019            data['FltBack'][0] = parmDict['FltBack']
2020            line += ' esd: %15.3g'%(sigDict['FltBack'])
2021        print (line)
2022        for ilay,layer in enumerate(data['Layers']):
2023            name = layer['Name']
2024            print (' Parameters for layer: %d %s'%(ilay,name))
2025            cid = str(ilay)+';'
2026            line = ' '
2027            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
2028            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
2029            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2030                if parm in layer:
2031                    if parm == 'Rough':
2032                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
2033                    else:
2034                        layer[parm][0] = parmDict[cid+parm]
2035                    line += ' %s: %.3f'%(parm,layer[parm][0])
2036                    if cid+parm in varyList:
2037                        line += ' esd: %.3g'%(sigDict[cid+parm])
2038            print (line)
2039            print (line2)
2040   
2041    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2042        parmDict.update(zip(varyList,values))
2043        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2044        return M
2045   
2046    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2047        parmDict.update(zip(varyList,values))
2048        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2049        return np.sum(M**2)
2050   
2051    def getREFD(Q,Qsig,parmDict):
2052        Ic = np.ones_like(Q)*parmDict['FltBack']
2053        Scale = parmDict['Scale']
2054        Nlayers = parmDict['nLayers']
2055        Res = parmDict['Res']
2056        depth = np.zeros(Nlayers)
2057        rho = np.zeros(Nlayers)
2058        irho = np.zeros(Nlayers)
2059        sigma = np.zeros(Nlayers)
2060        for ilay,lay in enumerate(parmDict['Layer Seq']):
2061            cid = str(lay)+';'
2062            depth[ilay] = parmDict[cid+'Thick']
2063            sigma[ilay] = parmDict[cid+'Rough']
2064            if parmDict[cid+'Name'] == u'unit scatter':
2065                rho[ilay] = parmDict[cid+'DenMul']
2066                irho[ilay] = parmDict[cid+'iDenMul']
2067            elif 'vacuum' != parmDict[cid+'Name']:
2068                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
2069                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
2070            if cid+'Mag SLD' in parmDict:
2071                rho[ilay] += parmDict[cid+'Mag SLD']
2072        if parmDict['dQ type'] == 'None':
2073            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2074        elif 'const' in parmDict['dQ type']:
2075            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
2076        else:       #dQ/Q in data
2077            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
2078        Ic += AB*Scale
2079        return Ic
2080       
2081    def estimateT0(takestep):
2082        Mmax = -1.e-10
2083        Mmin = 1.e10
2084        for i in range(100):
2085            x0 = takestep(values)
2086            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2087            Mmin = min(M,Mmin)
2088            MMax = max(M,Mmax)
2089        return 1.5*(MMax-Mmin)
2090
2091    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2092    if data.get('2% weight'):
2093        wt = 1./(0.02*Io)**2
2094    Qmin = Limits[1][0]
2095    Qmax = Limits[1][1]
2096    wtFactor = ProfDict['wtFactor']
2097    Bfac = data['Toler']
2098    Ibeg = np.searchsorted(Q,Qmin)
2099    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
2100    Ic[:] = 0
2101    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
2102              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
2103    parmDict,varyList,values,bounds = GetModelParms()
2104    Msg = 'Failed to converge'
2105    if varyList:
2106        if data['Minimizer'] == 'LMLS': 
2107            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
2108                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2109            parmDict.update(zip(varyList,result[0]))
2110            chisq = np.sum(result[2]['fvec']**2)
2111            ncalc = result[2]['nfev']
2112            covM = result[1]
2113            newVals = result[0]
2114        elif data['Minimizer'] == 'Basin Hopping':
2115            xyrng = np.array(bounds).T
2116            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
2117            T0 = estimateT0(take_step)
2118            print (' Estimated temperature: %.3g'%(T0))
2119            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
2120                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
2121                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
2122            chisq = result.fun
2123            ncalc = result.nfev
2124            newVals = result.x
2125            covM = []
2126        elif data['Minimizer'] == 'MC/SA Anneal':
2127            xyrng = np.array(bounds).T
2128            result = G2mth.anneal(sumREFD, values, 
2129                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
2130                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
2131                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
2132                ranRange=0.20,autoRan=False,dlg=None)
2133            newVals = result[0]
2134            parmDict.update(zip(varyList,newVals))
2135            chisq = result[1]
2136            ncalc = result[3]
2137            covM = []
2138            print (' MC/SA final temperature: %.4g'%(result[2]))
2139        elif data['Minimizer'] == 'L-BFGS-B':
2140            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
2141                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2142            parmDict.update(zip(varyList,result['x']))
2143            chisq = result.fun
2144            ncalc = result.nfev
2145            newVals = result.x
2146            covM = []
2147    else:   #nothing varied
2148        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2149        chisq = np.sum(M**2)
2150        ncalc = 0
2151        covM = []
2152        sig = []
2153        sigDict = {}
2154        result = []
2155    Rvals = {}
2156    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
2157    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
2158    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
2159    Ib[Ibeg:Ifin] = parmDict['FltBack']
2160    try:
2161        if not len(varyList):
2162            Msg += ' - nothing refined'
2163            raise ValueError
2164        Nans = np.isnan(newVals)
2165        if np.any(Nans):
2166            name = varyList[Nans.nonzero(True)[0]]
2167            Msg += ' Nan result for '+name+'!'
2168            raise ValueError
2169        Negs = np.less_equal(newVals,0.)
2170        if np.any(Negs):
2171            indx = Negs.nonzero()
2172            name = varyList[indx[0][0]]
2173            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
2174                Msg += ' negative coefficient for '+name+'!'
2175                raise ValueError
2176        if len(covM):
2177            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
2178            covMatrix = covM*Rvals['GOF']
2179        else:
2180            sig = np.zeros(len(varyList))
2181            covMatrix = []
2182        sigDict = dict(zip(varyList,sig))
2183        print (' Results of reflectometry data modelling fit:')
2184        print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
2185        print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2186        SetModelParms()
2187        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
2188    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
2189        print (Msg)
2190        return False,0,0,0,0,0,0,Msg
2191       
2192def makeSLDprofile(data,Substances):
2193   
2194    sq2 = np.sqrt(2.)
2195    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2196    Nlayers = len(laySeq)
2197    laySeq = np.array(laySeq,dtype=int)
2198    interfaces = np.zeros(Nlayers)
2199    rho = np.zeros(Nlayers)
2200    sigma = np.zeros(Nlayers)
2201    name = data['Layers'][0]['Name']
2202    thick = 0.
2203    for ilay,lay in enumerate(laySeq):
2204        layer = data['Layers'][lay]
2205        name = layer['Name']
2206        if 'Thick' in layer:
2207            thick += layer['Thick'][0]
2208            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
2209        if 'Rough' in layer:
2210            sigma[ilay] = max(0.001,layer['Rough'][0])
2211        if name != 'vacuum':
2212            if name == 'unit scatter':
2213                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
2214            else:
2215                rrho = Substances[name]['Scatt density']
2216                irho = Substances[name]['XImag density']
2217                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
2218        if 'Mag SLD' in layer:
2219            rho[ilay] += layer['Mag SLD'][0]
2220    name = data['Layers'][-1]['Name']
2221    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
2222    xr = np.flipud(x)
2223    interfaces[-1] = x[-1]
2224    y = np.ones_like(x)*rho[0]
2225    iBeg = 0
2226    for ilayer in range(Nlayers-1):
2227        delt = rho[ilayer+1]-rho[ilayer]
2228        iPos = np.searchsorted(x,interfaces[ilayer])
2229        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
2230        iBeg = iPos
2231    return x,xr,y   
2232
2233def REFDModelFxn(Profile,Inst,Limits,Substances,data):
2234   
2235    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2236    Qmin = Limits[1][0]
2237    Qmax = Limits[1][1]
2238    iBeg = np.searchsorted(Q,Qmin)
2239    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2240    Ib[:] = data['FltBack'][0]
2241    Ic[:] = 0
2242    Scale = data['Scale'][0]
2243    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2244    Nlayers = len(laySeq)
2245    depth = np.zeros(Nlayers)
2246    rho = np.zeros(Nlayers)
2247    irho = np.zeros(Nlayers)
2248    sigma = np.zeros(Nlayers)
2249    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
2250        layer = data['Layers'][lay]
2251        name = layer['Name']
2252        if 'Thick' in layer:    #skips first & last layers
2253            depth[ilay] = layer['Thick'][0]
2254        if 'Rough' in layer:    #skips first layer
2255            sigma[ilay] = layer['Rough'][0]
2256        if 'unit scatter' == name:
2257            rho[ilay] = layer['DenMul'][0]
2258            irho[ilay] = layer['iDenMul'][0]
2259        else:
2260            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2261            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
2262        if 'Mag SLD' in layer:
2263            rho[ilay] += layer['Mag SLD'][0]
2264    if data['dQ type'] == 'None':
2265        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2266    elif 'const' in data['dQ type']:
2267        res = data['Resolution'][0]/(100.*sateln2)
2268        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
2269    else:       #dQ/Q in data
2270        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
2271    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
2272
2273def abeles(kz, depth, rho, irho=0, sigma=0):
2274    """
2275    Optical matrix form of the reflectivity calculation.
2276    O.S. Heavens, Optical Properties of Thin Solid Films
2277   
2278    Reflectometry as a function of kz for a set of slabs.
2279
2280    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
2281        This is :math:`\\tfrac12 Q_z`.       
2282    :param depth:  float[m] (Ang).
2283        thickness of each layer.  The thickness of the incident medium
2284        and substrate are ignored.
2285    :param rho:  float[n,k] (1e-6/Ang^2)
2286        Real scattering length density for each layer for each kz
2287    :param irho:  float[n,k] (1e-6/Ang^2)
2288        Imaginary scattering length density for each layer for each kz
2289        Note: absorption cross section mu = 2 irho/lambda for neutrons
2290    :param sigma: float[m-1] (Ang)
2291        interfacial roughness.  This is the roughness between a layer
2292        and the previous layer. The sigma array should have m-1 entries.
2293
2294    Slabs are ordered with the surface SLD at index 0 and substrate at
2295    index -1, or reversed if kz < 0.
2296    """
2297    def calc(kz, depth, rho, irho, sigma):
2298        if len(kz) == 0: return kz
2299   
2300        # Complex index of refraction is relative to the incident medium.
2301        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
2302        # in place of kz^2, and ignoring rho_o
2303        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
2304        k = kz
2305   
2306        # According to Heavens, the initial matrix should be [ 1 F; F 1],
2307        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
2308        # multiply versus some coding convenience.
2309        B11 = 1
2310        B22 = 1
2311        B21 = 0
2312        B12 = 0
2313        for i in range(0, len(depth)-1):
2314            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
2315            F = (k - k_next) / (k + k_next)
2316            F *= np.exp(-2*k*k_next*sigma[i]**2)
2317            #print "==== layer",i
2318            #print "kz:", kz
2319            #print "k:", k
2320            #print "k_next:",k_next
2321            #print "F:",F
2322            #print "rho:",rho[:,i+1]
2323            #print "irho:",irho[:,i+1]
2324            #print "d:",depth[i],"sigma:",sigma[i]
2325            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
2326            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
2327            M21 = F*M11
2328            M12 = F*M22
2329            C1 = B11*M11 + B21*M12
2330            C2 = B11*M21 + B21*M22
2331            B11 = C1
2332            B21 = C2
2333            C1 = B12*M11 + B22*M12
2334            C2 = B12*M21 + B22*M22
2335            B12 = C1
2336            B22 = C2
2337            k = k_next
2338   
2339        r = B12/B11
2340        return np.absolute(r)**2
2341
2342    if np.isscalar(kz): kz = np.asarray([kz], 'd')
2343
2344    m = len(depth)
2345
2346    # Make everything into arrays
2347    depth = np.asarray(depth,'d')
2348    rho = np.asarray(rho,'d')
2349    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
2350    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
2351
2352    # Repeat rho,irho columns as needed
2353    if len(rho.shape) == 1:
2354        rho = rho[None,:]
2355        irho = irho[None,:]
2356
2357    return calc(kz, depth, rho, irho, sigma)
2358   
2359def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
2360    y = abeles(kz, depth, rho, irho, sigma)
2361    s = dq/2.
2362    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
2363    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
2364    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
2365    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
2366    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
2367    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
2368    y *= 0.137023
2369    return y
2370       
2371def makeRefdFFT(Limits,Profile):
2372    print ('make fft')
2373    Q,Io = Profile[:2]
2374    Qmin = Limits[1][0]
2375    Qmax = Limits[1][1]
2376    iBeg = np.searchsorted(Q,Qmin)
2377    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2378    Qf = np.linspace(0.,Q[iFin-1],5000)
2379    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
2380    If = QI(Qf)*Qf**4
2381    R = np.linspace(0.,5000.,5000)
2382    F = fft.rfft(If)
2383    return R,F
2384
2385   
2386################################################################################
2387# Stacking fault simulation codes
2388################################################################################
2389
2390def GetStackParms(Layers):
2391   
2392    Parms = []
2393#cell parms
2394    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
2395        Parms.append('cellA')
2396        Parms.append('cellC')
2397    else:
2398        Parms.append('cellA')
2399        Parms.append('cellB')
2400        Parms.append('cellC')
2401        if Layers['Laue'] != 'mmm':
2402            Parms.append('cellG')
2403#Transition parms
2404    for iY in range(len(Layers['Layers'])):
2405        for iX in range(len(Layers['Layers'])):
2406            Parms.append('TransP;%d;%d'%(iY,iX))
2407            Parms.append('TransX;%d;%d'%(iY,iX))
2408            Parms.append('TransY;%d;%d'%(iY,iX))
2409            Parms.append('TransZ;%d;%d'%(iY,iX))
2410    return Parms
2411
2412def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
2413    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
2414   
2415    :param dict Layers: dict with following items
2416
2417      ::
2418
2419       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
2420       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
2421       'Layers':[],'Stacking':[],'Transitions':[]}
2422       
2423    :param str ctrls: controls string to be written on DIFFaX controls.dif file
2424    :param float scale: scale factor
2425    :param dict background: background parameters
2426    :param list limits: min/max 2-theta to be calculated
2427    :param dict inst: instrument parameters dictionary
2428    :param list profile: powder pattern data
2429   
2430    Note that parameters all updated in place   
2431    '''
2432    import atmdata
2433    path = sys.path
2434    for name in path:
2435        if 'bin' in name:
2436            DIFFaX = name+'/DIFFaX.exe'
2437            print (' Execute '+DIFFaX)
2438            break
2439    # make form factor file that DIFFaX wants - atom types are GSASII style
2440    sf = open('data.sfc','w')
2441    sf.write('GSASII special form factor file for DIFFaX\n\n')
2442    atTypes = list(Layers['AtInfo'].keys())
2443    if 'H' not in atTypes:
2444        atTypes.insert(0,'H')
2445    for atType in atTypes:
2446        if atType == 'H': 
2447            blen = -.3741
2448        else:
2449            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
2450        Adat = atmdata.XrayFF[atType]
2451        text = '%4s'%(atType.ljust(4))
2452        for i in range(4):
2453            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
2454        text += '%11.6f%11.6f'%(Adat['fc'],blen)
2455        text += '%3d\n'%(Adat['Z'])
2456        sf.write(text)
2457    sf.close()
2458    #make DIFFaX control.dif file - future use GUI to set some of these flags
2459    cf = open('control.dif','w')
2460    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2461        x0 = profile[0]
2462        iBeg = np.searchsorted(x0,limits[0])
2463        iFin = np.searchsorted(x0,limits[1])+1
2464        if iFin-iBeg > 20000:
2465            iFin = iBeg+20000
2466        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2467        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2468        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
2469    else:
2470        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2471        inst = {'Type':['XSC','XSC',]}
2472    cf.close()
2473    #make DIFFaX data file
2474    df = open('GSASII-DIFFaX.dat','w')
2475    df.write('INSTRUMENTAL\n')
2476    if 'X' in inst['Type'][0]:
2477        df.write('X-RAY\n')
2478    elif 'N' in inst['Type'][0]:
2479        df.write('NEUTRON\n')
2480    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2481        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
2482        U = ateln2*inst['U'][1]/10000.
2483        V = ateln2*inst['V'][1]/10000.
2484        W = ateln2*inst['W'][1]/10000.
2485        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2486        HW = np.sqrt(np.mean(HWHM))
2487    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
2488        if 'Mean' in Layers['selInst']:
2489            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
2490        elif 'Gaussian' in Layers['selInst']:
2491            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
2492        else:
2493            df.write('None\n')
2494    else:
2495        df.write('0.10\nNone\n')
2496    df.write('STRUCTURAL\n')
2497    a,b,c = Layers['Cell'][1:4]
2498    gam = Layers['Cell'][6]
2499    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
2500    laue = Layers['Laue']
2501    if laue == '2/m(ab)':
2502        laue = '2/m(1)'
2503    elif laue == '2/m(c)':
2504        laue = '2/m(2)'
2505    if 'unknown' in Layers['Laue']:
2506        df.write('%s %.3f\n'%(laue,Layers['Toler']))
2507    else:   
2508        df.write('%s\n'%(laue))
2509    df.write('%d\n'%(len(Layers['Layers'])))
2510    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
2511        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
2512    layerNames = []
2513    for layer in Layers['Layers']:
2514        layerNames.append(layer['Name'])
2515    for il,layer in enumerate(Layers['Layers']):
2516        if layer['SameAs']:
2517            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
2518            continue
2519        df.write('LAYER %d\n'%(il+1))
2520        if '-1' in layer['Symm']:
2521            df.write('CENTROSYMMETRIC\n')
2522        else:
2523            df.write('NONE\n')
2524        for ia,atom in enumerate(layer['Atoms']):
2525            [name,atype,x,y,z,frac,Uiso] = atom
2526            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
2527                frac /= 2.
2528            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
2529    df.write('STACKING\n')
2530    df.write('%s\n'%(Layers['Stacking'][0]))
2531    if 'recursive' in Layers['Stacking'][0]:
2532        df.write('%s\n'%Layers['Stacking'][1])
2533    else:
2534        if 'list' in Layers['Stacking'][1]:
2535            Slen = len(Layers['Stacking'][2])
2536            iB = 0
2537            iF = 0
2538            while True:
2539                iF += 68
2540                if iF >= Slen:
2541                    break
2542                iF = min(iF,Slen)
2543                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
2544                iB = iF
2545        else:
2546            df.write('%s\n'%Layers['Stacking'][1])   
2547    df.write('TRANSITIONS\n')
2548    for iY in range(len(Layers['Layers'])):
2549        sumPx = 0.
2550        for iX in range(len(Layers['Layers'])):
2551            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
2552            p = round(p,3)
2553            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
2554            sumPx += p
2555        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
2556            print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
2557            df.close()
2558            os.remove('data.sfc')
2559            os.remove('control.dif')
2560            os.remove('GSASII-DIFFaX.dat')
2561            return       
2562    df.close()   
2563    time0 = time.time()
2564    try:
2565        subp.call(DIFFaX)
2566    except OSError:
2567        print (' DIFFax.exe is not available for this platform - under development')
2568    print (' DIFFaX time = %.2fs'%(time.time()-time0))
2569    if os.path.exists('GSASII-DIFFaX.spc'):
2570        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
2571        iFin = iBeg+Xpat.shape[1]
2572        bakType,backDict,backVary = SetBackgroundParms(background)
2573        backDict['Lam1'] = G2mth.getWave(inst)
2574        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2575        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
2576        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2577            rv = st.poisson(profile[3][iBeg:iFin])
2578            profile[1][iBeg:iFin] = rv.rvs()
2579            Z = np.ones_like(profile[3][iBeg:iFin])
2580            Z[1::2] *= -1
2581            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2582            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2583        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2584    #cleanup files..
2585        os.remove('GSASII-DIFFaX.spc')
2586    elif os.path.exists('GSASII-DIFFaX.sadp'):
2587        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
2588        Sadp = np.reshape(Sadp,(256,-1))
2589        Layers['Sadp']['Img'] = Sadp
2590        os.remove('GSASII-DIFFaX.sadp')
2591    os.remove('data.sfc')
2592    os.remove('control.dif')
2593    os.remove('GSASII-DIFFaX.dat')
2594   
2595def SetPWDRscan(inst,limits,profile):
2596   
2597    wave = G2mth.getMeanWave(inst)
2598    x0 = profile[0]
2599    iBeg = np.searchsorted(x0,limits[0])
2600    iFin = np.searchsorted(x0,limits[1])
2601    if iFin-iBeg > 20000:
2602        iFin = iBeg+20000
2603    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2604    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
2605    return iFin-iBeg
2606       
2607def SetStackingSF(Layers,debug):
2608# Load scattering factors into DIFFaX arrays
2609    import atmdata
2610    atTypes = Layers['AtInfo'].keys()
2611    aTypes = []
2612    for atype in atTypes:
2613        aTypes.append('%4s'%(atype.ljust(4)))
2614    SFdat = []
2615    for atType in atTypes:
2616        Adat = atmdata.XrayFF[atType]
2617        SF = np.zeros(9)
2618        SF[:8:2] = Adat['fa']
2619        SF[1:8:2] = Adat['fb']
2620        SF[8] = Adat['fc']
2621        SFdat.append(SF)
2622    SFdat = np.array(SFdat)
2623    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
2624   
2625def SetStackingClay(Layers,Type):
2626# Controls
2627    rand.seed()
2628    ranSeed = rand.randint(1,2**16-1)
2629    try:   
2630        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
2631            '6/m','6/mmm'].index(Layers['Laue'])+1
2632    except ValueError:  #for 'unknown'
2633        laueId = -1
2634    if 'SADP' in Type:
2635        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
2636        lmax = int(Layers['Sadp']['Lmax'])
2637    else:
2638        planeId = 0
2639        lmax = 0
2640# Sequences
2641    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
2642    try:
2643        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
2644    except ValueError:
2645        StkParm = -1
2646    if StkParm == 2:    #list
2647        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
2648        Nstk = len(StkSeq)
2649    else:
2650        Nstk = 1
2651        StkSeq = [0,]
2652    if StkParm == -1:
2653        StkParm = int(Layers['Stacking'][1])
2654    Wdth = Layers['Width'][0]
2655    mult = 1
2656    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
2657    LaueSym = Layers['Laue'].ljust(12)
2658    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
2659    return laueId,controls
2660   
2661def SetCellAtoms(Layers):
2662    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
2663# atoms in layers
2664    atTypes = list(Layers['AtInfo'].keys())
2665    AtomXOU = []
2666    AtomTp = []
2667    LayerSymm = []
2668    LayerNum = []
2669    layerNames = []
2670    Natm = 0
2671    Nuniq = 0
2672    for layer in Layers['Layers']:
2673        layerNames.append(layer['Name'])
2674    for il,layer in enumerate(Layers['Layers']):
2675        if layer['SameAs']:
2676            LayerNum.append(layerNames.index(layer['SameAs'])+1)
2677            continue
2678        else:
2679            LayerNum.append(il+1)
2680            Nuniq += 1
2681        if '-1' in layer['Symm']:
2682            LayerSymm.append(1)
2683        else:
2684            LayerSymm.append(0)
2685        for ia,atom in enumerate(layer['Atoms']):
2686            [name,atype,x,y,z,frac,Uiso] = atom
2687            Natm += 1
2688            AtomTp.append('%4s'%(atype.ljust(4)))
2689            Ta = atTypes.index(atype)+1
2690            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
2691    AtomXOU = np.array(AtomXOU)
2692    Nlayers = len(layerNames)
2693    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
2694    return Nlayers
2695   
2696def SetStackingTrans(Layers,Nlayers):
2697# Transitions
2698    TransX = []
2699    TransP = []
2700    for Ytrans in Layers['Transitions']:
2701        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
2702        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
2703    TransP = np.array(TransP,dtype='float').T
2704    TransX = np.array(TransX,dtype='float')
2705#    GSASIIpath.IPyBreak()
2706    pyx.pygettrans(Nlayers,TransP,TransX)
2707   
2708def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
2709# Scattering factors
2710    SetStackingSF(Layers,debug)
2711# Controls & sequences
2712    laueId,controls = SetStackingClay(Layers,'PWDR')
2713# cell & atoms
2714    Nlayers = SetCellAtoms(Layers)
2715    Volume = Layers['Cell'][7]   
2716# Transitions
2717    SetStackingTrans(Layers,Nlayers)
2718# PWDR scan
2719    Nsteps = SetPWDRscan(inst,limits,profile)
2720# result as Spec
2721    x0 = profile[0]
2722    profile[3] = np.zeros(len(profile[0]))
2723    profile[4] = np.zeros(len(profile[0]))
2724    profile[5] = np.zeros(len(profile[0]))
2725    iBeg = np.searchsorted(x0,limits[0])
2726    iFin = np.searchsorted(x0,limits[1])+1
2727    if iFin-iBeg > 20000:
2728        iFin = iBeg+20000
2729    Nspec = 20001       
2730    spec = np.zeros(Nspec,dtype='double')   
2731    time0 = time.time()
2732    pyx.pygetspc(controls,Nspec,spec)
2733    print (' GETSPC time = %.2fs'%(time.time()-time0))
2734    time0 = time.time()
2735    U = ateln2*inst['U'][1]/10000.
2736    V = ateln2*inst['V'][1]/10000.
2737    W = ateln2*inst['W'][1]/10000.
2738    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2739    HW = np.sqrt(np.mean(HWHM))
2740    BrdSpec = np.zeros(Nsteps)
2741    if 'Mean' in Layers['selInst']:
2742        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
2743    elif 'Gaussian' in Layers['selInst']:
2744        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
2745    else:
2746        BrdSpec = spec[:Nsteps]
2747    BrdSpec /= Volume
2748    iFin = iBeg+Nsteps
2749    bakType,backDict,backVary = SetBackgroundParms(background)
2750    backDict['Lam1'] = G2mth.getWave(inst)
2751    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2752    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
2753    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2754        try:
2755            rv = st.poisson(profile[3][iBeg:iFin])
2756            profile[1][iBeg:iFin] = rv.rvs()
2757        except ValueError:
2758            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
2759        Z = np.ones_like(profile[3][iBeg:iFin])
2760        Z[1::2] *= -1
2761        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2762        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2763    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2764    print (' Broadening time = %.2fs'%(time.time()-time0))
2765   
2766def CalcStackingSADP(Layers,debug):
2767   
2768# Scattering factors
2769    SetStackingSF(Layers,debug)
2770# Controls & sequences
2771    laueId,controls = SetStackingClay(Layers,'SADP')
2772# cell & atoms
2773    Nlayers = SetCellAtoms(Layers)   
2774# Transitions
2775    SetStackingTrans(Layers,Nlayers)
2776# result as Sadp
2777    Nspec = 20001       
2778    spec = np.zeros(Nspec,dtype='double')   
2779    time0 = time.time()
2780    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
2781    Sapd = np.zeros((256,256))
2782    iB = 0
2783    for i in range(hkLim):
2784        iF = iB+Nblk
2785        p1 = 127+int(i*Incr)
2786        p2 = 128-int(i*Incr)
2787        if Nblk == 128:
2788            if i:
2789                Sapd[128:,p1] = spec[iB:iF]
2790                Sapd[:128,p1] = spec[iF:iB:-1]
2791            Sapd[128:,p2] = spec[iB:iF]
2792            Sapd[:128,p2] = spec[iF:iB:-1]
2793        else:
2794            if i:
2795                Sapd[:,p1] = spec[iB:iF]
2796            Sapd[:,p2] = spec[iB:iF]
2797        iB += Nblk
2798    Layers['Sadp']['Img'] = Sapd
2799    print (' GETSAD time = %.2fs'%(time.time()-time0))
2800#    GSASIIpath.IPyBreak()
2801   
2802#testing data
2803NeedTestData = True
2804def TestData():
2805    'needs a doc string'
2806#    global NeedTestData
2807    global bakType
2808    bakType = 'chebyschev'
2809    global xdata
2810    xdata = np.linspace(4.0,40.0,36000)
2811    global parmDict0
2812    parmDict0 = {
2813        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
2814        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
2815        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
2816        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
2817        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'Z':0.0,'SH/L':0.002,
2818        'Back0':5.384,'Back1':-0.015,'Back2':.004,
2819        }
2820    global parmDict1
2821    parmDict1 = {
2822        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
2823        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
2824        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
2825        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
2826        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
2827        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
2828        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'Z':0.0,'SH/L':0.002,
2829        'Back0':36.897,'Back1':-0.508,'Back2':.006,
2830        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2831        }
2832    global parmDict2
2833    parmDict2 = {
2834        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
2835        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'Z':0.0,'SH/L':0.02,
2836        'Back0':5.,'Back1':-0.02,'Back2':.004,
2837#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2838        }
2839    global varyList
2840    varyList = []
2841
2842def test0():
2843    if NeedTestData: TestData()
2844    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
2845    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
2846    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
2847    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
2848    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
2849    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
2850   
2851def test1():
2852    if NeedTestData: TestData()
2853    time0 = time.time()
2854    for i in range(100):
2855        getPeakProfile(parmDict1,xdata,varyList,bakType)
2856    print ('100+6*Ka1-2 peaks=1200 peaks %.2f'%time.time()-time0)
2857   
2858def test2(name,delt):
2859    if NeedTestData: TestData()
2860    varyList = [name,]
2861    xdata = np.linspace(5.6,5.8,400)
2862    hplot = plotter.add('derivatives test for '+name).gca()
2863    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
2864    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2865    parmDict2[name] += delt
2866    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2867    hplot.plot(xdata,(y1-y0)/delt,'r+')
2868   
2869def test3(name,delt):
2870    if NeedTestData: TestData()
2871    names = ['pos','sig','gam','shl']
2872    idx = names.index(name)
2873    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
2874    xdata = np.linspace(5.6,5.8,800)
2875    dx = xdata[1]-xdata[0]
2876    hplot = plotter.add('derivatives test for '+name).gca()
2877    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
2878    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2879    myDict[name] += delt
2880    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2881    hplot.plot(xdata,(y1-y0)/delt,'r+')
2882
2883if __name__ == '__main__':
2884    import GSASIItestplot as plot
2885    global plotter
2886    plotter = plot.PlotNotebook()
2887#    test0()
2888#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)']:
2889    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
2890        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['Z',0.0001],['SH/L',0.00005]]:
2891        test2(name,shft)
2892    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
2893        test3(name,shft)
2894    print ("OK")
2895    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.