source: trunk/GSASIIpwd.py @ 4102

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

add transfer of flat background to images corrected for 1/dist2
correct error in recalibrate; eliminate duplicate reflections in e.g. silicon from genHKLpeak

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