source: trunk/GSASIIpwd.py @ 3712

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

add new configuration item: PDF_Rmax maximum radius for G(r) calculations; rarely changed by user
Change PDF GUI Rmax to be for the max r for G(r) plot - no longer r max for calculation
This speeds up optimization
Fix problems with PDF peak table - principally the rowlabels needed to be str not int - crashed wx 4.0 & gave bizarre effects in wx 3.

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