source: trunk/GSASIIpwd.py @ 4080

Last change on this file since 4080 was 4080, checked in by vondreele, 2 years ago

fix to ellipsoid size math & GUI revisions

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 125.9 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASII powder calculation module*
5==================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2019-08-07 16:35:30 +0000 (Wed, 07 Aug 2019) $
10# $Author: vondreele $
11# $Revision: 4080 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 4080 2019-08-07 16:35:30Z 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: 4080 $")
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):
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    for h,k,l,d in HKL:
1151        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1152        if ext and 'MagSpGrp' in SGData:
1153            ext = G2spc.checkMagextc([h,k,l],SGData)
1154        if not ext:
1155            if Inst == None:
1156                HKLs.append([h,k,l,d,0,-1])
1157            else:
1158                HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1159    return np.array(HKLs)
1160
1161def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1162    'needs a doc string'
1163    HKLs = []
1164    vec = np.array(Vec)
1165    vstar = np.sqrt(G2lat.calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1166    dvec = 1./(maxH*vstar+1./dmin)
1167    HKL = G2lat.GenHLaue(dvec,SGData,A)       
1168    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1169    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1170    ifMag = False
1171    if 'MagSpGrp' in SGData:
1172        ifMag = True
1173    for h,k,l,d in HKL:
1174        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1175        if not ext and d >= dmin:
1176            HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1177        for dH in SSdH:
1178            if dH:
1179                DH = SSdH[dH]
1180                H = [h+DH[0],k+DH[1],l+DH[2]]
1181                d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1182                if d >= dmin:
1183                    HKLM = np.array([h,k,l,dH])
1184                    if G2spc.checkSSextc(HKLM,SSGData) or ifMag:
1185                        HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])   
1186    return G2lat.sortHKLd(HKLs,True,True,True)
1187
1188def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
1189    'Computes the profile for a powder pattern'
1190   
1191    yb = getBackground('',parmDict,bakType,dataType,xdata)[0]
1192    yc = np.zeros_like(yb)
1193    cw = np.diff(xdata)
1194    cw = np.append(cw,cw[-1])
1195    if 'C' in dataType:
1196        shl = max(parmDict['SH/L'],0.002)
1197        Ka2 = False
1198        if 'Lam1' in parmDict.keys():
1199            Ka2 = True
1200            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1201            kRatio = parmDict['I(L2)/I(L1)']
1202        iPeak = 0
1203        while True:
1204            try:
1205                pos = parmDict['pos'+str(iPeak)]
1206                tth = (pos-parmDict['Zero'])
1207                intens = parmDict['int'+str(iPeak)]
1208                sigName = 'sig'+str(iPeak)
1209                if sigName in varyList:
1210                    sig = parmDict[sigName]
1211                else:
1212                    sig = G2mth.getCWsig(parmDict,tth)
1213                sig = max(sig,0.001)          #avoid neg sigma^2
1214                gamName = 'gam'+str(iPeak)
1215                if gamName in varyList:
1216                    gam = parmDict[gamName]
1217                else:
1218                    gam = G2mth.getCWgam(parmDict,tth)
1219                gam = max(gam,0.001)             #avoid neg gamma
1220                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1221                iBeg = np.searchsorted(xdata,pos-fmin)
1222                iFin = np.searchsorted(xdata,pos+fmin)
1223                if not iBeg+iFin:       #peak below low limit
1224                    iPeak += 1
1225                    continue
1226                elif not iBeg-iFin:     #peak above high limit
1227                    return yb+yc
1228                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1229                if Ka2:
1230                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1231                    iBeg = np.searchsorted(xdata,pos2-fmin)
1232                    iFin = np.searchsorted(xdata,pos2+fmin)
1233                    if iBeg-iFin:
1234                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1235                iPeak += 1
1236            except KeyError:        #no more peaks to process
1237                return yb+yc
1238    else:
1239        Pdabc = parmDict['Pdabc']
1240        difC = parmDict['difC']
1241        iPeak = 0
1242        while True:
1243            try:
1244                pos = parmDict['pos'+str(iPeak)]               
1245                tof = pos-parmDict['Zero']
1246                dsp = tof/difC
1247                intens = parmDict['int'+str(iPeak)]
1248                alpName = 'alp'+str(iPeak)
1249                if alpName in varyList:
1250                    alp = parmDict[alpName]
1251                else:
1252                    if len(Pdabc):
1253                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1254                    else:
1255                        alp = G2mth.getTOFalpha(parmDict,dsp)
1256                alp = max(0.1,alp)
1257                betName = 'bet'+str(iPeak)
1258                if betName in varyList:
1259                    bet = parmDict[betName]
1260                else:
1261                    if len(Pdabc):
1262                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1263                    else:
1264                        bet = G2mth.getTOFbeta(parmDict,dsp)
1265                bet = max(0.0001,bet)
1266                sigName = 'sig'+str(iPeak)
1267                if sigName in varyList:
1268                    sig = parmDict[sigName]
1269                else:
1270                    sig = G2mth.getTOFsig(parmDict,dsp)
1271                gamName = 'gam'+str(iPeak)
1272                if gamName in varyList:
1273                    gam = parmDict[gamName]
1274                else:
1275                    gam = G2mth.getTOFgamma(parmDict,dsp)
1276                gam = max(gam,0.001)             #avoid neg gamma
1277                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1278                iBeg = np.searchsorted(xdata,pos-fmin)
1279                iFin = np.searchsorted(xdata,pos+fmax)
1280                lenX = len(xdata)
1281                if not iBeg:
1282                    iFin = np.searchsorted(xdata,pos+fmax)
1283                elif iBeg == lenX:
1284                    iFin = iBeg
1285                else:
1286                    iFin = np.searchsorted(xdata,pos+fmax)
1287                if not iBeg+iFin:       #peak below low limit
1288                    iPeak += 1
1289                    continue
1290                elif not iBeg-iFin:     #peak above high limit
1291                    return yb+yc
1292                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1293                iPeak += 1
1294            except KeyError:        #no more peaks to process
1295                return yb+yc
1296           
1297def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
1298    'needs a doc string'
1299# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1300    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1301    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,dataType,xdata)
1302    if 'Back;0' in varyList:            #background derivs are in front if present
1303        dMdv[0:len(dMdb)] = dMdb
1304    names = ['DebyeA','DebyeR','DebyeU']
1305    for name in varyList:
1306        if 'Debye' in name:
1307            parm,Id = name.split(';')
1308            ip = names.index(parm)
1309            dMdv[varyList.index(name)] = dMddb[3*int(Id)+ip]
1310    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1311    for name in varyList:
1312        if 'BkPk' in name:
1313            parm,Id = name.split(';')
1314            ip = names.index(parm)
1315            dMdv[varyList.index(name)] = dMdpk[4*int(Id)+ip]
1316    cw = np.diff(xdata)
1317    cw = np.append(cw,cw[-1])
1318    if 'C' in dataType:
1319        shl = max(parmDict['SH/L'],0.002)
1320        Ka2 = False
1321        if 'Lam1' in parmDict.keys():
1322            Ka2 = True
1323            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1324            kRatio = parmDict['I(L2)/I(L1)']
1325        iPeak = 0
1326        while True:
1327            try:
1328                pos = parmDict['pos'+str(iPeak)]
1329                tth = (pos-parmDict['Zero'])
1330                intens = parmDict['int'+str(iPeak)]
1331                sigName = 'sig'+str(iPeak)
1332                if sigName in varyList:
1333                    sig = parmDict[sigName]
1334                    dsdU = dsdV = dsdW = 0
1335                else:
1336                    sig = G2mth.getCWsig(parmDict,tth)
1337                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1338                sig = max(sig,0.001)          #avoid neg sigma
1339                gamName = 'gam'+str(iPeak)
1340                if gamName in varyList:
1341                    gam = parmDict[gamName]
1342                    dgdX = dgdY = dgdZ = 0
1343                else:
1344                    gam = G2mth.getCWgam(parmDict,tth)
1345                    dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1346                gam = max(gam,0.001)             #avoid neg gamma
1347                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1348                iBeg = np.searchsorted(xdata,pos-fmin)
1349                iFin = np.searchsorted(xdata,pos+fmin)
1350                if not iBeg+iFin:       #peak below low limit
1351                    iPeak += 1
1352                    continue
1353                elif not iBeg-iFin:     #peak above high limit
1354                    break
1355                dMdpk = np.zeros(shape=(6,len(xdata)))
1356                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1357                for i in range(1,5):
1358                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
1359                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
1360                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1361                if Ka2:
1362                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1363                    iBeg = np.searchsorted(xdata,pos2-fmin)
1364                    iFin = np.searchsorted(xdata,pos2+fmin)
1365                    if iBeg-iFin:
1366                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1367                        for i in range(1,5):
1368                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
1369                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
1370                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
1371                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1372                for parmName in ['pos','int','sig','gam']:
1373                    try:
1374                        idx = varyList.index(parmName+str(iPeak))
1375                        dMdv[idx] = dervDict[parmName]
1376                    except ValueError:
1377                        pass
1378                if 'U' in varyList:
1379                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1380                if 'V' in varyList:
1381                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1382                if 'W' in varyList:
1383                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1384                if 'X' in varyList:
1385                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1386                if 'Y' in varyList:
1387                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1388                if 'Z' in varyList:
1389                    dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1390                if 'SH/L' in varyList:
1391                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
1392                if 'I(L2)/I(L1)' in varyList:
1393                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1394                iPeak += 1
1395            except KeyError:        #no more peaks to process
1396                break
1397    else:
1398        Pdabc = parmDict['Pdabc']
1399        difC = parmDict['difC']
1400        iPeak = 0
1401        while True:
1402            try:
1403                pos = parmDict['pos'+str(iPeak)]               
1404                tof = pos-parmDict['Zero']
1405                dsp = tof/difC
1406                intens = parmDict['int'+str(iPeak)]
1407                alpName = 'alp'+str(iPeak)
1408                if alpName in varyList:
1409                    alp = parmDict[alpName]
1410                else:
1411                    if len(Pdabc):
1412                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1413                        dada0 = 0
1414                    else:
1415                        alp = G2mth.getTOFalpha(parmDict,dsp)
1416                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1417                betName = 'bet'+str(iPeak)
1418                if betName in varyList:
1419                    bet = parmDict[betName]
1420                else:
1421                    if len(Pdabc):
1422                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1423                        dbdb0 = dbdb1 = dbdb2 = 0
1424                    else:
1425                        bet = G2mth.getTOFbeta(parmDict,dsp)
1426                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1427                sigName = 'sig'+str(iPeak)
1428                if sigName in varyList:
1429                    sig = parmDict[sigName]
1430                    dsds0 = dsds1 = dsds2 = dsds3 = 0
1431                else:
1432                    sig = G2mth.getTOFsig(parmDict,dsp)
1433                    dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1434                gamName = 'gam'+str(iPeak)
1435                if gamName in varyList:
1436                    gam = parmDict[gamName]
1437                    dsdX = dsdY = dsdZ = 0
1438                else:
1439                    gam = G2mth.getTOFgamma(parmDict,dsp)
1440                    dsdX,dsdY,dsdZ = G2mth.getTOFgammaDeriv(dsp)
1441                gam = max(gam,0.001)             #avoid neg gamma
1442                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1443                iBeg = np.searchsorted(xdata,pos-fmin)
1444                lenX = len(xdata)
1445                if not iBeg:
1446                    iFin = np.searchsorted(xdata,pos+fmax)
1447                elif iBeg == lenX:
1448                    iFin = iBeg
1449                else:
1450                    iFin = np.searchsorted(xdata,pos+fmax)
1451                if not iBeg+iFin:       #peak below low limit
1452                    iPeak += 1
1453                    continue
1454                elif not iBeg-iFin:     #peak above high limit
1455                    break
1456                dMdpk = np.zeros(shape=(7,len(xdata)))
1457                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1458                for i in range(1,6):
1459                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1460                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1461                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1462                for parmName in ['pos','int','alp','bet','sig','gam']:
1463                    try:
1464                        idx = varyList.index(parmName+str(iPeak))
1465                        dMdv[idx] = dervDict[parmName]
1466                    except ValueError:
1467                        pass
1468                if 'alpha' in varyList:
1469                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1470                if 'beta-0' in varyList:
1471                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1472                if 'beta-1' in varyList:
1473                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1474                if 'beta-q' in varyList:
1475                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1476                if 'sig-0' in varyList:
1477                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1478                if 'sig-1' in varyList:
1479                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1480                if 'sig-2' in varyList:
1481                    dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1482                if 'sig-q' in varyList:
1483                    dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1484                if 'X' in varyList:
1485                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1486                if 'Y' in varyList:
1487                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']
1488                if 'Z' in varyList:
1489                    dMdv[varyList.index('Z')] += dsdZ*dervDict['gam']
1490                iPeak += 1
1491            except KeyError:        #no more peaks to process
1492                break
1493    return dMdv
1494       
1495def Dict2Values(parmdict, varylist):
1496    '''Use before call to leastsq to setup list of values for the parameters
1497    in parmdict, as selected by key in varylist'''
1498    return [parmdict[key] for key in varylist] 
1499   
1500def Values2Dict(parmdict, varylist, values):
1501    ''' Use after call to leastsq to update the parameter dictionary with
1502    values corresponding to keys in varylist'''
1503    parmdict.update(zip(varylist,values))
1504   
1505def SetBackgroundParms(Background):
1506    'Loads background parameters into dicts/lists to create varylist & parmdict'
1507    if len(Background) == 1:            # fix up old backgrounds
1508        Background.append({'nDebye':0,'debyeTerms':[]})
1509    bakType,bakFlag = Background[0][:2]
1510    backVals = Background[0][3:]
1511    backNames = ['Back;'+str(i) for i in range(len(backVals))]
1512    Debye = Background[1]           #also has background peaks stuff
1513    backDict = dict(zip(backNames,backVals))
1514    backVary = []
1515    if bakFlag:
1516        backVary = backNames
1517
1518    backDict['nDebye'] = Debye['nDebye']
1519    debyeDict = {}
1520    debyeList = []
1521    for i in range(Debye['nDebye']):
1522        debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1523        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1524        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1525    debyeVary = []
1526    for item in debyeList:
1527        if item[1]:
1528            debyeVary.append(item[0])
1529    backDict.update(debyeDict)
1530    backVary += debyeVary
1531
1532    backDict['nPeaks'] = Debye['nPeaks']
1533    peaksDict = {}
1534    peaksList = []
1535    for i in range(Debye['nPeaks']):
1536        peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1537        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1538        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1539    peaksVary = []
1540    for item in peaksList:
1541        if item[1]:
1542            peaksVary.append(item[0])
1543    backDict.update(peaksDict)
1544    backVary += peaksVary
1545    return bakType,backDict,backVary
1546   
1547def DoCalibInst(IndexPeaks,Inst):
1548   
1549    def SetInstParms():
1550        dataType = Inst['Type'][0]
1551        insVary = []
1552        insNames = []
1553        insVals = []
1554        for parm in Inst:
1555            insNames.append(parm)
1556            insVals.append(Inst[parm][1])
1557            if parm in ['Lam','difC','difA','difB','Zero',]:
1558                if Inst[parm][2]:
1559                    insVary.append(parm)
1560        instDict = dict(zip(insNames,insVals))
1561        return dataType,instDict,insVary
1562       
1563    def GetInstParms(parmDict,Inst,varyList):
1564        for name in Inst:
1565            Inst[name][1] = parmDict[name]
1566       
1567    def InstPrint(Inst,sigDict):
1568        print ('Instrument Parameters:')
1569        if 'C' in Inst['Type'][0]:
1570            ptfmt = "%12.6f"
1571        else:
1572            ptfmt = "%12.3f"
1573        ptlbls = 'names :'
1574        ptstr =  'values:'
1575        sigstr = 'esds  :'
1576        for parm in Inst:
1577            if parm in  ['Lam','difC','difA','difB','Zero',]:
1578                ptlbls += "%s" % (parm.center(12))
1579                ptstr += ptfmt % (Inst[parm][1])
1580                if parm in sigDict:
1581                    sigstr += ptfmt % (sigDict[parm])
1582                else:
1583                    sigstr += 12*' '
1584        print (ptlbls)
1585        print (ptstr)
1586        print (sigstr)
1587       
1588    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1589        parmDict.update(zip(varyList,values))
1590        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1591
1592    peakPos = []
1593    peakDsp = []
1594    peakWt = []
1595    for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1596        if peak[2] and peak[3] and sig > 0.:
1597            peakPos.append(peak[0])
1598            peakDsp.append(peak[-1])    #d-calc
1599#            peakWt.append(peak[-1]**2/sig**2)   #weight by d**2
1600            peakWt.append(1./(sig*peak[-1]))   #
1601    peakPos = np.array(peakPos)
1602    peakDsp = np.array(peakDsp)
1603    peakWt = np.array(peakWt)
1604    dataType,insDict,insVary = SetInstParms()
1605    parmDict = {}
1606    parmDict.update(insDict)
1607    varyList = insVary
1608    if not len(varyList):
1609        G2fil.G2Print ('**** ERROR - nothing to refine! ****')
1610        return False
1611    while True:
1612        begin = time.time()
1613        values =  np.array(Dict2Values(parmDict, varyList))
1614        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1615            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1616        ncyc = int(result[2]['nfev']/2)
1617        runtime = time.time()-begin   
1618        chisq = np.sum(result[2]['fvec']**2)
1619        Values2Dict(parmDict, varyList, result[0])
1620        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
1621        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(peakPos),len(varyList)))
1622        G2fil.G2Print ('calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
1623        G2fil.G2Print ('chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF))
1624        try:
1625            sig = np.sqrt(np.diag(result[1])*GOF)
1626            if np.any(np.isnan(sig)):
1627                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***')
1628            break                   #refinement succeeded - finish up!
1629        except ValueError:          #result[1] is None on singular matrix
1630            G2fil.G2Print ('**** Refinement failed - singular matrix ****')
1631       
1632    sigDict = dict(zip(varyList,sig))
1633    GetInstParms(parmDict,Inst,varyList)
1634    InstPrint(Inst,sigDict)
1635    return True
1636           
1637def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
1638    '''Called to perform a peak fit, refining the selected items in the peak
1639    table as well as selected items in the background.
1640
1641    :param str FitPgm: type of fit to perform. At present this is ignored.
1642    :param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1643      four values followed by a refine flag where the values are: position, intensity,
1644      sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1645      tree entry, dict item "peaks"
1646    :param list Background: describes the background. List with two items.
1647      Item 0 specifies a background model and coefficients. Item 1 is a dict.
1648      From the Histogram/Background tree entry.
1649    :param list Limits: min and max x-value to use
1650    :param dict Inst: Instrument parameters
1651    :param dict Inst2: more Instrument parameters
1652    :param numpy.array data: a 5xn array. data[0] is the x-values,
1653      data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1654      calc, background and difference intensities, respectively.
1655    :param array fixback: fixed background values
1656    :param list prevVaryList: Used in sequential refinements to override the
1657      variable list. Defaults as an empty list.
1658    :param bool oneCycle: True if only one cycle of fitting should be performed
1659    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1660      and derivType = controls['deriv type']. If None default values are used.
1661    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1662      Defaults to None, which means no updates are done.
1663    '''
1664    def GetBackgroundParms(parmList,Background):
1665        iBak = 0
1666        while True:
1667            try:
1668                bakName = 'Back;'+str(iBak)
1669                Background[0][iBak+3] = parmList[bakName]
1670                iBak += 1
1671            except KeyError:
1672                break
1673        iDb = 0
1674        while True:
1675            names = ['DebyeA;','DebyeR;','DebyeU;']
1676            try:
1677                for i,name in enumerate(names):
1678                    val = parmList[name+str(iDb)]
1679                    Background[1]['debyeTerms'][iDb][2*i] = val
1680                iDb += 1
1681            except KeyError:
1682                break
1683        iDb = 0
1684        while True:
1685            names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1686            try:
1687                for i,name in enumerate(names):
1688                    val = parmList[name+str(iDb)]
1689                    Background[1]['peaksList'][iDb][2*i] = val
1690                iDb += 1
1691            except KeyError:
1692                break
1693               
1694    def BackgroundPrint(Background,sigDict):
1695        print ('Background coefficients for '+Background[0][0]+' function')
1696        ptfmt = "%12.5f"
1697        ptstr =  'value: '
1698        sigstr = 'esd  : '
1699        for i,back in enumerate(Background[0][3:]):
1700            ptstr += ptfmt % (back)
1701            if Background[0][1]:
1702                prm = 'Back;'+str(i)
1703                if prm in sigDict:
1704                    sigstr += ptfmt % (sigDict[prm])
1705                else:
1706                    sigstr += " "*12
1707            if len(ptstr) > 75:
1708                print (ptstr)
1709                if Background[0][1]: print (sigstr)
1710                ptstr =  'value: '
1711                sigstr = 'esd  : '
1712        if len(ptstr) > 8:
1713            print (ptstr)
1714            if Background[0][1]: print (sigstr)
1715
1716        if Background[1]['nDebye']:
1717            parms = ['DebyeA;','DebyeR;','DebyeU;']
1718            print ('Debye diffuse scattering coefficients')
1719            ptfmt = "%12.5f"
1720            print (' term       DebyeA       esd        DebyeR       esd        DebyeU        esd')
1721            for term in range(Background[1]['nDebye']):
1722                line = ' term %d'%(term)
1723                for ip,name in enumerate(parms):
1724                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1725                    if name+str(term) in sigDict:
1726                        line += ptfmt%(sigDict[name+str(term)])
1727                    else:
1728                        line += " "*12
1729                print (line)
1730        if Background[1]['nPeaks']:
1731            print ('Coefficients for Background Peaks')
1732            ptfmt = "%15.3f"
1733            for j,pl in enumerate(Background[1]['peaksList']):
1734                names =  'peak %3d:'%(j+1)
1735                ptstr =  'values  :'
1736                sigstr = 'esds    :'
1737                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1738                    val = pl[2*i]
1739                    prm = lbl+";"+str(j)
1740                    names += '%15s'%(prm)
1741                    ptstr += ptfmt%(val)
1742                    if prm in sigDict:
1743                        sigstr += ptfmt%(sigDict[prm])
1744                    else:
1745                        sigstr += " "*15
1746                print (names)
1747                print (ptstr)
1748                print (sigstr)
1749                           
1750    def SetInstParms(Inst):
1751        dataType = Inst['Type'][0]
1752        insVary = []
1753        insNames = []
1754        insVals = []
1755        for parm in Inst:
1756            insNames.append(parm)
1757            insVals.append(Inst[parm][1])
1758            if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
1759                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1760                    insVary.append(parm)
1761        instDict = dict(zip(insNames,insVals))
1762#        instDict['X'] = max(instDict['X'],0.01)
1763#        instDict['Y'] = max(instDict['Y'],0.01)
1764        if 'SH/L' in instDict:
1765            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1766        return dataType,instDict,insVary
1767       
1768    def GetInstParms(parmDict,Inst,varyList,Peaks):
1769        for name in Inst:
1770            Inst[name][1] = parmDict[name]
1771        iPeak = 0
1772        while True:
1773            try:
1774                sigName = 'sig'+str(iPeak)
1775                pos = parmDict['pos'+str(iPeak)]
1776                if sigName not in varyList:
1777                    if 'C' in Inst['Type'][0]:
1778                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1779                    else:
1780                        dsp = G2lat.Pos2dsp(Inst,pos)
1781                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1782                gamName = 'gam'+str(iPeak)
1783                if gamName not in varyList:
1784                    if 'C' in Inst['Type'][0]:
1785                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1786                    else:
1787                        dsp = G2lat.Pos2dsp(Inst,pos)
1788                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1789                iPeak += 1
1790            except KeyError:
1791                break
1792       
1793    def InstPrint(Inst,sigDict):
1794        print ('Instrument Parameters:')
1795        ptfmt = "%12.6f"
1796        ptlbls = 'names :'
1797        ptstr =  'values:'
1798        sigstr = 'esds  :'
1799        for parm in Inst:
1800            if parm in  ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
1801                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1802                ptlbls += "%s" % (parm.center(12))
1803                ptstr += ptfmt % (Inst[parm][1])
1804                if parm in sigDict:
1805                    sigstr += ptfmt % (sigDict[parm])
1806                else:
1807                    sigstr += 12*' '
1808        print (ptlbls)
1809        print (ptstr)
1810        print (sigstr)
1811
1812    def SetPeaksParms(dataType,Peaks):
1813        peakNames = []
1814        peakVary = []
1815        peakVals = []
1816        if 'C' in dataType:
1817            names = ['pos','int','sig','gam']
1818        else:
1819            names = ['pos','int','alp','bet','sig','gam']
1820        for i,peak in enumerate(Peaks):
1821            for j,name in enumerate(names):
1822                peakVals.append(peak[2*j])
1823                parName = name+str(i)
1824                peakNames.append(parName)
1825                if peak[2*j+1]:
1826                    peakVary.append(parName)
1827        return dict(zip(peakNames,peakVals)),peakVary
1828               
1829    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1830        if 'C' in Inst['Type'][0]:
1831            names = ['pos','int','sig','gam']
1832        else:   #'T'
1833            names = ['pos','int','alp','bet','sig','gam']
1834        for i,peak in enumerate(Peaks):
1835            pos = parmDict['pos'+str(i)]
1836            if 'difC' in Inst:
1837                dsp = pos/Inst['difC'][1]
1838            for j in range(len(names)):
1839                parName = names[j]+str(i)
1840                if parName in varyList:
1841                    peak[2*j] = parmDict[parName]
1842                elif 'alpha' in parName:
1843                    peak[2*j] = parmDict['alpha']/dsp
1844                elif 'beta' in parName:
1845                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1846                elif 'sig' in parName:
1847                    if 'C' in Inst['Type'][0]:
1848                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1849                    else:
1850                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1851                elif 'gam' in parName:
1852                    if 'C' in Inst['Type'][0]:
1853                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1854                    else:
1855                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1856                       
1857    def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
1858        print ('Peak coefficients:')
1859        if 'C' in dataType:
1860            names = ['pos','int','sig','gam']
1861        else:   #'T'
1862            names = ['pos','int','alp','bet','sig','gam']           
1863        head = 13*' '
1864        for name in names:
1865            if name in ['alp','bet']:
1866                head += name.center(8)+'esd'.center(8)
1867            else:
1868                head += name.center(10)+'esd'.center(10)
1869        head += 'bins'.center(8)
1870        print (head)
1871        if 'C' in dataType:
1872            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1873        else:
1874            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1875        for i,peak in enumerate(Peaks):
1876            ptstr =  ':'
1877            for j in range(len(names)):
1878                name = names[j]
1879                parName = name+str(i)
1880                ptstr += ptfmt[name] % (parmDict[parName])
1881                if parName in varyList:
1882                    ptstr += ptfmt[name] % (sigDict[parName])
1883                else:
1884                    if name in ['alp','bet']:
1885                        ptstr += 8*' '
1886                    else:
1887                        ptstr += 10*' '
1888            ptstr += '%9.2f'%(ptsperFW[i])
1889            print ('%s'%(('Peak'+str(i+1)).center(8)),ptstr)
1890               
1891    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1892        parmdict.update(zip(varylist,values))
1893        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1894           
1895    def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):       
1896        parmdict.update(zip(varylist,values))
1897        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1898        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1899        if dlg:
1900            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1901            if not GoOn:
1902                return -M           #abort!!
1903        return M
1904       
1905    if controls:
1906        Ftol = controls['min dM/M']
1907    else:
1908        Ftol = 0.0001
1909    if oneCycle:
1910        Ftol = 1.0
1911    x,y,w,yc,yb,yd = data   #these are numpy arrays - remove masks!
1912    if fixback is None:
1913        fixback = np.zeros_like(y)
1914    yc *= 0.                            #set calcd ones to zero
1915    yb *= 0.
1916    yd *= 0.
1917    cw = x[1:]-x[:-1]
1918    xBeg = np.searchsorted(x,Limits[0])
1919    xFin = np.searchsorted(x,Limits[1])+1
1920    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1921    dataType,insDict,insVary = SetInstParms(Inst)
1922    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1923    parmDict = {}
1924    parmDict.update(bakDict)
1925    parmDict.update(insDict)
1926    parmDict.update(peakDict)
1927    parmDict['Pdabc'] = []      #dummy Pdabc
1928    parmDict.update(Inst2)      #put in real one if there
1929    if prevVaryList:
1930        varyList = prevVaryList[:]
1931    else:
1932        varyList = bakVary+insVary+peakVary
1933    fullvaryList = varyList[:]
1934    while True:
1935        begin = time.time()
1936        values =  np.array(Dict2Values(parmDict, varyList))
1937        Rvals = {}
1938        badVary = []
1939        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1940               args=(x[xBeg:xFin],(y+fixback)[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1941        ncyc = int(result[2]['nfev']/2)
1942        runtime = time.time()-begin   
1943        chisq = np.sum(result[2]['fvec']**2)
1944        Values2Dict(parmDict, varyList, result[0])
1945        Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*(y+fixback)[xBeg:xFin]**2))*100.      #to %
1946        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1947        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList)))
1948        if ncyc:
1949            G2fil.G2Print ('fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
1950        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
1951        sig = [0]*len(varyList)
1952        if len(varyList) == 0: break  # if nothing was refined
1953        try:
1954            sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1955            if np.any(np.isnan(sig)):
1956                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***')
1957            break                   #refinement succeeded - finish up!
1958        except ValueError:          #result[1] is None on singular matrix
1959            G2fil.G2Print ('**** Refinement failed - singular matrix ****')
1960            Ipvt = result[2]['ipvt']
1961            for i,ipvt in enumerate(Ipvt):
1962                if not np.sum(result[2]['fjac'],axis=1)[i]:
1963                    G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
1964                    badVary.append(varyList[ipvt-1])
1965                    del(varyList[ipvt-1])
1966                    break
1967            else: # nothing removed
1968                break
1969    if dlg: dlg.Destroy()
1970    sigDict = dict(zip(varyList,sig))
1971    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]-fixback[xBeg:xFin]
1972    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)-fixback[xBeg:xFin]
1973    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1974    GetBackgroundParms(parmDict,Background)
1975    if bakVary: BackgroundPrint(Background,sigDict)
1976    GetInstParms(parmDict,Inst,varyList,Peaks)
1977    if insVary: InstPrint(Inst,sigDict)
1978    GetPeaksParms(Inst,parmDict,Peaks,varyList)
1979    binsperFWHM = []
1980    for peak in Peaks:
1981        FWHM = getFWHM(peak[0],Inst)
1982        try:
1983            binsperFWHM.append(FWHM/cw[x.searchsorted(peak[0])])
1984        except IndexError:
1985            binsperFWHM.append(0.)
1986    if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
1987    if len(binsperFWHM):
1988        if min(binsperFWHM) < 1.:
1989            G2fil.G2Print ('*** Warning: calculated peak widths are too narrow to refine profile coefficients ***')
1990            if 'T' in Inst['Type'][0]:
1991                G2fil.G2Print (' Manually increase sig-0, 1, or 2 in Instrument Parameters')
1992            else:
1993                G2fil.G2Print (' Manually increase W in Instrument Parameters')
1994        elif min(binsperFWHM) < 4.:
1995            G2fil.G2Print ('*** Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement ***')
1996            G2fil.G2Print ('*** recommended is 6-10; you have %.2f ***'%(min(binsperFWHM)))
1997    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
1998   
1999def calcIncident(Iparm,xdata):
2000    'needs a doc string'
2001
2002    def IfunAdv(Iparm,xdata):
2003        Itype = Iparm['Itype']
2004        Icoef = Iparm['Icoeff']
2005        DYI = np.ones((12,xdata.shape[0]))
2006        YI = np.ones_like(xdata)*Icoef[0]
2007       
2008        x = xdata/1000.                 #expressions are in ms
2009        if Itype == 'Exponential':
2010            for i in [1,3,5,7,9]:
2011                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
2012                YI += Icoef[i]*Eterm
2013                DYI[i] *= Eterm
2014                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
2015        elif 'Maxwell'in Itype:
2016            Eterm = np.exp(-Icoef[2]/x**2)
2017            DYI[1] = Eterm/x**5
2018            DYI[2] = -Icoef[1]*DYI[1]/x**2
2019            YI += (Icoef[1]*Eterm/x**5)
2020            if 'Exponential' in Itype:
2021                for i in range(3,11,2):
2022                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
2023                    YI += Icoef[i]*Eterm
2024                    DYI[i] *= Eterm
2025                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
2026            else:   #Chebyschev
2027                T = (2./x)-1.
2028                Ccof = np.ones((12,xdata.shape[0]))
2029                Ccof[1] = T
2030                for i in range(2,12):
2031                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
2032                for i in range(1,10):
2033                    YI += Ccof[i]*Icoef[i+2]
2034                    DYI[i+2] =Ccof[i]
2035        return YI,DYI
2036       
2037    Iesd = np.array(Iparm['Iesd'])
2038    Icovar = Iparm['Icovar']
2039    YI,DYI = IfunAdv(Iparm,xdata)
2040    YI = np.where(YI>0,YI,1.)
2041    WYI = np.zeros_like(xdata)
2042    vcov = np.zeros((12,12))
2043    k = 0
2044    for i in range(12):
2045        for j in range(i,12):
2046            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
2047            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
2048            k += 1
2049    M = np.inner(vcov,DYI.T)
2050    WYI = np.sum(M*DYI,axis=0)
2051    WYI = np.where(WYI>0.,WYI,0.)
2052    return YI,WYI
2053   
2054################################################################################
2055# Reflectometry calculations
2056################################################################################
2057
2058def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
2059    G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
2060   
2061    class RandomDisplacementBounds(object):
2062        """random displacement with bounds"""
2063        def __init__(self, xmin, xmax, stepsize=0.5):
2064            self.xmin = xmin
2065            self.xmax = xmax
2066            self.stepsize = stepsize
2067   
2068        def __call__(self, x):
2069            """take a random step but ensure the new position is within the bounds"""
2070            while True:
2071                # this could be done in a much more clever way, but it will work for example purposes
2072                steps = self.xmax-self.xmin
2073                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
2074                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
2075                    break
2076            return xnew
2077   
2078    def GetModelParms():
2079        parmDict = {}
2080        varyList = []
2081        values = []
2082        bounds = []
2083        parmDict['dQ type'] = data['dQ type']
2084        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
2085        for parm in ['Scale','FltBack']:
2086            parmDict[parm] = data[parm][0]
2087            if data[parm][1]:
2088                varyList.append(parm)
2089                values.append(data[parm][0])
2090                bounds.append(Bounds[parm])
2091        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
2092        parmDict['nLayers'] = len(parmDict['Layer Seq'])
2093        for ilay,layer in enumerate(data['Layers']):
2094            name = layer['Name']
2095            cid = str(ilay)+';'
2096            parmDict[cid+'Name'] = name
2097            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2098                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
2099                if layer.get(parm,[0.,False])[1]:
2100                    varyList.append(cid+parm)
2101                    value = layer[parm][0]
2102                    values.append(value)
2103                    if value:
2104                        bound = [value*Bfac,value/Bfac]
2105                    else:
2106                        bound = [0.,10.]
2107                    bounds.append(bound)
2108            if name not in ['vacuum','unit scatter']:
2109                parmDict[cid+'rho'] = Substances[name]['Scatt density']
2110                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
2111        return parmDict,varyList,values,bounds
2112   
2113    def SetModelParms():
2114        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
2115        if 'Scale' in varyList:
2116            data['Scale'][0] = parmDict['Scale']
2117            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
2118        G2fil.G2Print (line)
2119        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
2120        if 'FltBack' in varyList:
2121            data['FltBack'][0] = parmDict['FltBack']
2122            line += ' esd: %15.3g'%(sigDict['FltBack'])
2123        G2fil.G2Print (line)
2124        for ilay,layer in enumerate(data['Layers']):
2125            name = layer['Name']
2126            G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
2127            cid = str(ilay)+';'
2128            line = ' '
2129            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
2130            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
2131            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2132                if parm in layer:
2133                    if parm == 'Rough':
2134                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
2135                    else:
2136                        layer[parm][0] = parmDict[cid+parm]
2137                    line += ' %s: %.3f'%(parm,layer[parm][0])
2138                    if cid+parm in varyList:
2139                        line += ' esd: %.3g'%(sigDict[cid+parm])
2140            G2fil.G2Print (line)
2141            G2fil.G2Print (line2)
2142   
2143    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2144        parmDict.update(zip(varyList,values))
2145        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2146        return M
2147   
2148    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2149        parmDict.update(zip(varyList,values))
2150        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2151        return np.sum(M**2)
2152   
2153    def getREFD(Q,Qsig,parmDict):
2154        Ic = np.ones_like(Q)*parmDict['FltBack']
2155        Scale = parmDict['Scale']
2156        Nlayers = parmDict['nLayers']
2157        Res = parmDict['Res']
2158        depth = np.zeros(Nlayers)
2159        rho = np.zeros(Nlayers)
2160        irho = np.zeros(Nlayers)
2161        sigma = np.zeros(Nlayers)
2162        for ilay,lay in enumerate(parmDict['Layer Seq']):
2163            cid = str(lay)+';'
2164            depth[ilay] = parmDict[cid+'Thick']
2165            sigma[ilay] = parmDict[cid+'Rough']
2166            if parmDict[cid+'Name'] == u'unit scatter':
2167                rho[ilay] = parmDict[cid+'DenMul']
2168                irho[ilay] = parmDict[cid+'iDenMul']
2169            elif 'vacuum' != parmDict[cid+'Name']:
2170                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
2171                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
2172            if cid+'Mag SLD' in parmDict:
2173                rho[ilay] += parmDict[cid+'Mag SLD']
2174        if parmDict['dQ type'] == 'None':
2175            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2176        elif 'const' in parmDict['dQ type']:
2177            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
2178        else:       #dQ/Q in data
2179            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
2180        Ic += AB*Scale
2181        return Ic
2182       
2183    def estimateT0(takestep):
2184        Mmax = -1.e-10
2185        Mmin = 1.e10
2186        for i in range(100):
2187            x0 = takestep(values)
2188            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2189            Mmin = min(M,Mmin)
2190            MMax = max(M,Mmax)
2191        return 1.5*(MMax-Mmin)
2192
2193    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2194    if data.get('2% weight'):
2195        wt = 1./(0.02*Io)**2
2196    Qmin = Limits[1][0]
2197    Qmax = Limits[1][1]
2198    wtFactor = ProfDict['wtFactor']
2199    Bfac = data['Toler']
2200    Ibeg = np.searchsorted(Q,Qmin)
2201    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
2202    Ic[:] = 0
2203    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
2204              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
2205    parmDict,varyList,values,bounds = GetModelParms()
2206    Msg = 'Failed to converge'
2207    if varyList:
2208        if data['Minimizer'] == 'LMLS': 
2209            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
2210                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2211            parmDict.update(zip(varyList,result[0]))
2212            chisq = np.sum(result[2]['fvec']**2)
2213            ncalc = result[2]['nfev']
2214            covM = result[1]
2215            newVals = result[0]
2216        elif data['Minimizer'] == 'Basin Hopping':
2217            xyrng = np.array(bounds).T
2218            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
2219            T0 = estimateT0(take_step)
2220            G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
2221            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
2222                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
2223                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
2224            chisq = result.fun
2225            ncalc = result.nfev
2226            newVals = result.x
2227            covM = []
2228        elif data['Minimizer'] == 'MC/SA Anneal':
2229            xyrng = np.array(bounds).T
2230            result = G2mth.anneal(sumREFD, values, 
2231                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
2232                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
2233                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
2234                ranRange=0.20,autoRan=False,dlg=None)
2235            newVals = result[0]
2236            parmDict.update(zip(varyList,newVals))
2237            chisq = result[1]
2238            ncalc = result[3]
2239            covM = []
2240            G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
2241        elif data['Minimizer'] == 'L-BFGS-B':
2242            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
2243                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2244            parmDict.update(zip(varyList,result['x']))
2245            chisq = result.fun
2246            ncalc = result.nfev
2247            newVals = result.x
2248            covM = []
2249    else:   #nothing varied
2250        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2251        chisq = np.sum(M**2)
2252        ncalc = 0
2253        covM = []
2254        sig = []
2255        sigDict = {}
2256        result = []
2257    Rvals = {}
2258    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
2259    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
2260    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
2261    Ib[Ibeg:Ifin] = parmDict['FltBack']
2262    try:
2263        if not len(varyList):
2264            Msg += ' - nothing refined'
2265            raise ValueError
2266        Nans = np.isnan(newVals)
2267        if np.any(Nans):
2268            name = varyList[Nans.nonzero(True)[0]]
2269            Msg += ' Nan result for '+name+'!'
2270            raise ValueError
2271        Negs = np.less_equal(newVals,0.)
2272        if np.any(Negs):
2273            indx = Negs.nonzero()
2274            name = varyList[indx[0][0]]
2275            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
2276                Msg += ' negative coefficient for '+name+'!'
2277                raise ValueError
2278        if len(covM):
2279            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
2280            covMatrix = covM*Rvals['GOF']
2281        else:
2282            sig = np.zeros(len(varyList))
2283            covMatrix = []
2284        sigDict = dict(zip(varyList,sig))
2285        G2fil.G2Print (' Results of reflectometry data modelling fit:')
2286        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
2287        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2288        SetModelParms()
2289        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
2290    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
2291        G2fil.G2Print (Msg)
2292        return False,0,0,0,0,0,0,Msg
2293       
2294def makeSLDprofile(data,Substances):
2295   
2296    sq2 = np.sqrt(2.)
2297    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2298    Nlayers = len(laySeq)
2299    laySeq = np.array(laySeq,dtype=int)
2300    interfaces = np.zeros(Nlayers)
2301    rho = np.zeros(Nlayers)
2302    sigma = np.zeros(Nlayers)
2303    name = data['Layers'][0]['Name']
2304    thick = 0.
2305    for ilay,lay in enumerate(laySeq):
2306        layer = data['Layers'][lay]
2307        name = layer['Name']
2308        if 'Thick' in layer:
2309            thick += layer['Thick'][0]
2310            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
2311        if 'Rough' in layer:
2312            sigma[ilay] = max(0.001,layer['Rough'][0])
2313        if name != 'vacuum':
2314            if name == 'unit scatter':
2315                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
2316            else:
2317                rrho = Substances[name]['Scatt density']
2318                irho = Substances[name]['XImag density']
2319                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
2320        if 'Mag SLD' in layer:
2321            rho[ilay] += layer['Mag SLD'][0]
2322    name = data['Layers'][-1]['Name']
2323    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
2324    xr = np.flipud(x)
2325    interfaces[-1] = x[-1]
2326    y = np.ones_like(x)*rho[0]
2327    iBeg = 0
2328    for ilayer in range(Nlayers-1):
2329        delt = rho[ilayer+1]-rho[ilayer]
2330        iPos = np.searchsorted(x,interfaces[ilayer])
2331        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
2332        iBeg = iPos
2333    return x,xr,y   
2334
2335def REFDModelFxn(Profile,Inst,Limits,Substances,data):
2336   
2337    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2338    Qmin = Limits[1][0]
2339    Qmax = Limits[1][1]
2340    iBeg = np.searchsorted(Q,Qmin)
2341    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2342    Ib[:] = data['FltBack'][0]
2343    Ic[:] = 0
2344    Scale = data['Scale'][0]
2345    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2346    Nlayers = len(laySeq)
2347    depth = np.zeros(Nlayers)
2348    rho = np.zeros(Nlayers)
2349    irho = np.zeros(Nlayers)
2350    sigma = np.zeros(Nlayers)
2351    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
2352        layer = data['Layers'][lay]
2353        name = layer['Name']
2354        if 'Thick' in layer:    #skips first & last layers
2355            depth[ilay] = layer['Thick'][0]
2356        if 'Rough' in layer:    #skips first layer
2357            sigma[ilay] = layer['Rough'][0]
2358        if 'unit scatter' == name:
2359            rho[ilay] = layer['DenMul'][0]
2360            irho[ilay] = layer['iDenMul'][0]
2361        else:
2362            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2363            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
2364        if 'Mag SLD' in layer:
2365            rho[ilay] += layer['Mag SLD'][0]
2366    if data['dQ type'] == 'None':
2367        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2368    elif 'const' in data['dQ type']:
2369        res = data['Resolution'][0]/(100.*sateln2)
2370        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
2371    else:       #dQ/Q in data
2372        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
2373    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
2374
2375def abeles(kz, depth, rho, irho=0, sigma=0):
2376    """
2377    Optical matrix form of the reflectivity calculation.
2378    O.S. Heavens, Optical Properties of Thin Solid Films
2379   
2380    Reflectometry as a function of kz for a set of slabs.
2381
2382    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
2383        This is :math:`\\tfrac12 Q_z`.       
2384    :param depth:  float[m] (Ang).
2385        thickness of each layer.  The thickness of the incident medium
2386        and substrate are ignored.
2387    :param rho:  float[n,k] (1e-6/Ang^2)
2388        Real scattering length density for each layer for each kz
2389    :param irho:  float[n,k] (1e-6/Ang^2)
2390        Imaginary scattering length density for each layer for each kz
2391        Note: absorption cross section mu = 2 irho/lambda for neutrons
2392    :param sigma: float[m-1] (Ang)
2393        interfacial roughness.  This is the roughness between a layer
2394        and the previous layer. The sigma array should have m-1 entries.
2395
2396    Slabs are ordered with the surface SLD at index 0 and substrate at
2397    index -1, or reversed if kz < 0.
2398    """
2399    def calc(kz, depth, rho, irho, sigma):
2400        if len(kz) == 0: return kz
2401   
2402        # Complex index of refraction is relative to the incident medium.
2403        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
2404        # in place of kz^2, and ignoring rho_o
2405        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
2406        k = kz
2407   
2408        # According to Heavens, the initial matrix should be [ 1 F; F 1],
2409        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
2410        # multiply versus some coding convenience.
2411        B11 = 1
2412        B22 = 1
2413        B21 = 0
2414        B12 = 0
2415        for i in range(0, len(depth)-1):
2416            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
2417            F = (k - k_next) / (k + k_next)
2418            F *= np.exp(-2*k*k_next*sigma[i]**2)
2419            #print "==== layer",i
2420            #print "kz:", kz
2421            #print "k:", k
2422            #print "k_next:",k_next
2423            #print "F:",F
2424            #print "rho:",rho[:,i+1]
2425            #print "irho:",irho[:,i+1]
2426            #print "d:",depth[i],"sigma:",sigma[i]
2427            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
2428            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
2429            M21 = F*M11
2430            M12 = F*M22
2431            C1 = B11*M11 + B21*M12
2432            C2 = B11*M21 + B21*M22
2433            B11 = C1
2434            B21 = C2
2435            C1 = B12*M11 + B22*M12
2436            C2 = B12*M21 + B22*M22
2437            B12 = C1
2438            B22 = C2
2439            k = k_next
2440   
2441        r = B12/B11
2442        return np.absolute(r)**2
2443
2444    if np.isscalar(kz): kz = np.asarray([kz], 'd')
2445
2446    m = len(depth)
2447
2448    # Make everything into arrays
2449    depth = np.asarray(depth,'d')
2450    rho = np.asarray(rho,'d')
2451    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
2452    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
2453
2454    # Repeat rho,irho columns as needed
2455    if len(rho.shape) == 1:
2456        rho = rho[None,:]
2457        irho = irho[None,:]
2458
2459    return calc(kz, depth, rho, irho, sigma)
2460   
2461def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
2462    y = abeles(kz, depth, rho, irho, sigma)
2463    s = dq/2.
2464    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
2465    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
2466    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
2467    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
2468    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
2469    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
2470    y *= 0.137023
2471    return y
2472       
2473def makeRefdFFT(Limits,Profile):
2474    G2fil.G2Print ('make fft')
2475    Q,Io = Profile[:2]
2476    Qmin = Limits[1][0]
2477    Qmax = Limits[1][1]
2478    iBeg = np.searchsorted(Q,Qmin)
2479    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2480    Qf = np.linspace(0.,Q[iFin-1],5000)
2481    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
2482    If = QI(Qf)*Qf**4
2483    R = np.linspace(0.,5000.,5000)
2484    F = fft.rfft(If)
2485    return R,F
2486
2487   
2488################################################################################
2489# Stacking fault simulation codes
2490################################################################################
2491
2492def GetStackParms(Layers):
2493   
2494    Parms = []
2495#cell parms
2496    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
2497        Parms.append('cellA')
2498        Parms.append('cellC')
2499    else:
2500        Parms.append('cellA')
2501        Parms.append('cellB')
2502        Parms.append('cellC')
2503        if Layers['Laue'] != 'mmm':
2504            Parms.append('cellG')
2505#Transition parms
2506    for iY in range(len(Layers['Layers'])):
2507        for iX in range(len(Layers['Layers'])):
2508            Parms.append('TransP;%d;%d'%(iY,iX))
2509            Parms.append('TransX;%d;%d'%(iY,iX))
2510            Parms.append('TransY;%d;%d'%(iY,iX))
2511            Parms.append('TransZ;%d;%d'%(iY,iX))
2512    return Parms
2513
2514def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
2515    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
2516   
2517    :param dict Layers: dict with following items
2518
2519      ::
2520
2521       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
2522       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
2523       'Layers':[],'Stacking':[],'Transitions':[]}
2524       
2525    :param str ctrls: controls string to be written on DIFFaX controls.dif file
2526    :param float scale: scale factor
2527    :param dict background: background parameters
2528    :param list limits: min/max 2-theta to be calculated
2529    :param dict inst: instrument parameters dictionary
2530    :param list profile: powder pattern data
2531   
2532    Note that parameters all updated in place   
2533    '''
2534    import atmdata
2535    path = sys.path
2536    for name in path:
2537        if 'bin' in name:
2538            DIFFaX = name+'/DIFFaX.exe'
2539            G2fil.G2Print (' Execute '+DIFFaX)
2540            break
2541    # make form factor file that DIFFaX wants - atom types are GSASII style
2542    sf = open('data.sfc','w')
2543    sf.write('GSASII special form factor file for DIFFaX\n\n')
2544    atTypes = list(Layers['AtInfo'].keys())
2545    if 'H' not in atTypes:
2546        atTypes.insert(0,'H')
2547    for atType in atTypes:
2548        if atType == 'H': 
2549            blen = -.3741
2550        else:
2551            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
2552        Adat = atmdata.XrayFF[atType]
2553        text = '%4s'%(atType.ljust(4))
2554        for i in range(4):
2555            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
2556        text += '%11.6f%11.6f'%(Adat['fc'],blen)
2557        text += '%3d\n'%(Adat['Z'])
2558        sf.write(text)
2559    sf.close()
2560    #make DIFFaX control.dif file - future use GUI to set some of these flags
2561    cf = open('control.dif','w')
2562    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2563        x0 = profile[0]
2564        iBeg = np.searchsorted(x0,limits[0])
2565        iFin = np.searchsorted(x0,limits[1])+1
2566        if iFin-iBeg > 20000:
2567            iFin = iBeg+20000
2568        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2569        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2570        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
2571    else:
2572        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2573        inst = {'Type':['XSC','XSC',]}
2574    cf.close()
2575    #make DIFFaX data file
2576    df = open('GSASII-DIFFaX.dat','w')
2577    df.write('INSTRUMENTAL\n')
2578    if 'X' in inst['Type'][0]:
2579        df.write('X-RAY\n')
2580    elif 'N' in inst['Type'][0]:
2581        df.write('NEUTRON\n')
2582    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2583        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
2584        U = ateln2*inst['U'][1]/10000.
2585        V = ateln2*inst['V'][1]/10000.
2586        W = ateln2*inst['W'][1]/10000.
2587        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2588        HW = np.sqrt(np.mean(HWHM))
2589    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
2590        if 'Mean' in Layers['selInst']:
2591            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
2592        elif 'Gaussian' in Layers['selInst']:
2593            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
2594        else:
2595            df.write('None\n')
2596    else:
2597        df.write('0.10\nNone\n')
2598    df.write('STRUCTURAL\n')
2599    a,b,c = Layers['Cell'][1:4]
2600    gam = Layers['Cell'][6]
2601    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
2602    laue = Layers['Laue']
2603    if laue == '2/m(ab)':
2604        laue = '2/m(1)'
2605    elif laue == '2/m(c)':
2606        laue = '2/m(2)'
2607    if 'unknown' in Layers['Laue']:
2608        df.write('%s %.3f\n'%(laue,Layers['Toler']))
2609    else:   
2610        df.write('%s\n'%(laue))
2611    df.write('%d\n'%(len(Layers['Layers'])))
2612    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
2613        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
2614    layerNames = []
2615    for layer in Layers['Layers']:
2616        layerNames.append(layer['Name'])
2617    for il,layer in enumerate(Layers['Layers']):
2618        if layer['SameAs']:
2619            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
2620            continue
2621        df.write('LAYER %d\n'%(il+1))
2622        if '-1' in layer['Symm']:
2623            df.write('CENTROSYMMETRIC\n')
2624        else:
2625            df.write('NONE\n')
2626        for ia,atom in enumerate(layer['Atoms']):
2627            [name,atype,x,y,z,frac,Uiso] = atom
2628            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
2629                frac /= 2.
2630            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
2631    df.write('STACKING\n')
2632    df.write('%s\n'%(Layers['Stacking'][0]))
2633    if 'recursive' in Layers['Stacking'][0]:
2634        df.write('%s\n'%Layers['Stacking'][1])
2635    else:
2636        if 'list' in Layers['Stacking'][1]:
2637            Slen = len(Layers['Stacking'][2])
2638            iB = 0
2639            iF = 0
2640            while True:
2641                iF += 68
2642                if iF >= Slen:
2643                    break
2644                iF = min(iF,Slen)
2645                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
2646                iB = iF
2647        else:
2648            df.write('%s\n'%Layers['Stacking'][1])   
2649    df.write('TRANSITIONS\n')
2650    for iY in range(len(Layers['Layers'])):
2651        sumPx = 0.
2652        for iX in range(len(Layers['Layers'])):
2653            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
2654            p = round(p,3)
2655            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
2656            sumPx += p
2657        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
2658            G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
2659            df.close()
2660            os.remove('data.sfc')
2661            os.remove('control.dif')
2662            os.remove('GSASII-DIFFaX.dat')
2663            return       
2664    df.close()   
2665    time0 = time.time()
2666    try:
2667        subp.call(DIFFaX)
2668    except OSError:
2669        G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
2670    G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
2671    if os.path.exists('GSASII-DIFFaX.spc'):
2672        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
2673        iFin = iBeg+Xpat.shape[1]
2674        bakType,backDict,backVary = SetBackgroundParms(background)
2675        backDict['Lam1'] = G2mth.getWave(inst)
2676        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2677        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
2678        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2679            rv = st.poisson(profile[3][iBeg:iFin])
2680            profile[1][iBeg:iFin] = rv.rvs()
2681            Z = np.ones_like(profile[3][iBeg:iFin])
2682            Z[1::2] *= -1
2683            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2684            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2685        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2686    #cleanup files..
2687        os.remove('GSASII-DIFFaX.spc')
2688    elif os.path.exists('GSASII-DIFFaX.sadp'):
2689        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
2690        Sadp = np.reshape(Sadp,(256,-1))
2691        Layers['Sadp']['Img'] = Sadp
2692        os.remove('GSASII-DIFFaX.sadp')
2693    os.remove('data.sfc')
2694    os.remove('control.dif')
2695    os.remove('GSASII-DIFFaX.dat')
2696   
2697def SetPWDRscan(inst,limits,profile):
2698   
2699    wave = G2mth.getMeanWave(inst)
2700    x0 = profile[0]
2701    iBeg = np.searchsorted(x0,limits[0])
2702    iFin = np.searchsorted(x0,limits[1])
2703    if iFin-iBeg > 20000:
2704        iFin = iBeg+20000
2705    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2706    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
2707    return iFin-iBeg
2708       
2709def SetStackingSF(Layers,debug):
2710# Load scattering factors into DIFFaX arrays
2711    import atmdata
2712    atTypes = Layers['AtInfo'].keys()
2713    aTypes = []
2714    for atype in atTypes:
2715        aTypes.append('%4s'%(atype.ljust(4)))
2716    SFdat = []
2717    for atType in atTypes:
2718        Adat = atmdata.XrayFF[atType]
2719        SF = np.zeros(9)
2720        SF[:8:2] = Adat['fa']
2721        SF[1:8:2] = Adat['fb']
2722        SF[8] = Adat['fc']
2723        SFdat.append(SF)
2724    SFdat = np.array(SFdat)
2725    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
2726   
2727def SetStackingClay(Layers,Type):
2728# Controls
2729    rand.seed()
2730    ranSeed = rand.randint(1,2**16-1)
2731    try:   
2732        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
2733            '6/m','6/mmm'].index(Layers['Laue'])+1
2734    except ValueError:  #for 'unknown'
2735        laueId = -1
2736    if 'SADP' in Type:
2737        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
2738        lmax = int(Layers['Sadp']['Lmax'])
2739    else:
2740        planeId = 0
2741        lmax = 0
2742# Sequences
2743    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
2744    try:
2745        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
2746    except ValueError:
2747        StkParm = -1
2748    if StkParm == 2:    #list
2749        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
2750        Nstk = len(StkSeq)
2751    else:
2752        Nstk = 1
2753        StkSeq = [0,]
2754    if StkParm == -1:
2755        StkParm = int(Layers['Stacking'][1])
2756    Wdth = Layers['Width'][0]
2757    mult = 1
2758    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
2759    LaueSym = Layers['Laue'].ljust(12)
2760    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
2761    return laueId,controls
2762   
2763def SetCellAtoms(Layers):
2764    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
2765# atoms in layers
2766    atTypes = list(Layers['AtInfo'].keys())
2767    AtomXOU = []
2768    AtomTp = []
2769    LayerSymm = []
2770    LayerNum = []
2771    layerNames = []
2772    Natm = 0
2773    Nuniq = 0
2774    for layer in Layers['Layers']:
2775        layerNames.append(layer['Name'])
2776    for il,layer in enumerate(Layers['Layers']):
2777        if layer['SameAs']:
2778            LayerNum.append(layerNames.index(layer['SameAs'])+1)
2779            continue
2780        else:
2781            LayerNum.append(il+1)
2782            Nuniq += 1
2783        if '-1' in layer['Symm']:
2784            LayerSymm.append(1)
2785        else:
2786            LayerSymm.append(0)
2787        for ia,atom in enumerate(layer['Atoms']):
2788            [name,atype,x,y,z,frac,Uiso] = atom
2789            Natm += 1
2790            AtomTp.append('%4s'%(atype.ljust(4)))
2791            Ta = atTypes.index(atype)+1
2792            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
2793    AtomXOU = np.array(AtomXOU)
2794    Nlayers = len(layerNames)
2795    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
2796    return Nlayers
2797   
2798def SetStackingTrans(Layers,Nlayers):
2799# Transitions
2800    TransX = []
2801    TransP = []
2802    for Ytrans in Layers['Transitions']:
2803        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
2804        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
2805    TransP = np.array(TransP,dtype='float').T
2806    TransX = np.array(TransX,dtype='float')
2807#    GSASIIpath.IPyBreak()
2808    pyx.pygettrans(Nlayers,TransP,TransX)
2809   
2810def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
2811# Scattering factors
2812    SetStackingSF(Layers,debug)
2813# Controls & sequences
2814    laueId,controls = SetStackingClay(Layers,'PWDR')
2815# cell & atoms
2816    Nlayers = SetCellAtoms(Layers)
2817    Volume = Layers['Cell'][7]   
2818# Transitions
2819    SetStackingTrans(Layers,Nlayers)
2820# PWDR scan
2821    Nsteps = SetPWDRscan(inst,limits,profile)
2822# result as Spec
2823    x0 = profile[0]
2824    profile[3] = np.zeros(len(profile[0]))
2825    profile[4] = np.zeros(len(profile[0]))
2826    profile[5] = np.zeros(len(profile[0]))
2827    iBeg = np.searchsorted(x0,limits[0])
2828    iFin = np.searchsorted(x0,limits[1])+1
2829    if iFin-iBeg > 20000:
2830        iFin = iBeg+20000
2831    Nspec = 20001       
2832    spec = np.zeros(Nspec,dtype='double')   
2833    time0 = time.time()
2834    pyx.pygetspc(controls,Nspec,spec)
2835    G2fil.G2Print (' GETSPC time = %.2fs'%(time.time()-time0))
2836    time0 = time.time()
2837    U = ateln2*inst['U'][1]/10000.
2838    V = ateln2*inst['V'][1]/10000.
2839    W = ateln2*inst['W'][1]/10000.
2840    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2841    HW = np.sqrt(np.mean(HWHM))
2842    BrdSpec = np.zeros(Nsteps)
2843    if 'Mean' in Layers['selInst']:
2844        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
2845    elif 'Gaussian' in Layers['selInst']:
2846        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
2847    else:
2848        BrdSpec = spec[:Nsteps]
2849    BrdSpec /= Volume
2850    iFin = iBeg+Nsteps
2851    bakType,backDict,backVary = SetBackgroundParms(background)
2852    backDict['Lam1'] = G2mth.getWave(inst)
2853    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2854    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
2855    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2856        try:
2857            rv = st.poisson(profile[3][iBeg:iFin])
2858            profile[1][iBeg:iFin] = rv.rvs()
2859        except ValueError:
2860            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
2861        Z = np.ones_like(profile[3][iBeg:iFin])
2862        Z[1::2] *= -1
2863        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2864        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2865    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2866    G2fil.G2Print (' Broadening time = %.2fs'%(time.time()-time0))
2867   
2868def CalcStackingSADP(Layers,debug):
2869   
2870# Scattering factors
2871    SetStackingSF(Layers,debug)
2872# Controls & sequences
2873    laueId,controls = SetStackingClay(Layers,'SADP')
2874# cell & atoms
2875    Nlayers = SetCellAtoms(Layers)   
2876# Transitions
2877    SetStackingTrans(Layers,Nlayers)
2878# result as Sadp
2879    Nspec = 20001       
2880    spec = np.zeros(Nspec,dtype='double')   
2881    time0 = time.time()
2882    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
2883    Sapd = np.zeros((256,256))
2884    iB = 0
2885    for i in range(hkLim):
2886        iF = iB+Nblk
2887        p1 = 127+int(i*Incr)
2888        p2 = 128-int(i*Incr)
2889        if Nblk == 128:
2890            if i:
2891                Sapd[128:,p1] = spec[iB:iF]
2892                Sapd[:128,p1] = spec[iF:iB:-1]
2893            Sapd[128:,p2] = spec[iB:iF]
2894            Sapd[:128,p2] = spec[iF:iB:-1]
2895        else:
2896            if i:
2897                Sapd[:,p1] = spec[iB:iF]
2898            Sapd[:,p2] = spec[iB:iF]
2899        iB += Nblk
2900    Layers['Sadp']['Img'] = Sapd
2901    G2fil.G2Print (' GETSAD time = %.2fs'%(time.time()-time0))
2902   
2903###############################################################################
2904#### Maximum Entropy Method - Dysnomia
2905###############################################################################
2906   
2907def makePRFfile(data,MEMtype):
2908    ''' makes Dysnomia .prf control file from Dysnomia GUI controls
2909   
2910    :param dict data: GSAS-II phase data
2911    :param int MEMtype: 1 for neutron data with negative scattering lengths
2912                        0 otherwise
2913    :returns str: name of Dysnomia control file
2914    '''
2915
2916    generalData = data['General']
2917    pName = generalData['Name'].replace(' ','_')
2918    DysData = data['Dysnomia']
2919    prfName = pName+'.prf'
2920    prf = open(prfName,'w')
2921    prf.write('$PREFERENCES\n')
2922    prf.write(pName+'.mem\n') #or .fos?
2923    prf.write(pName+'.out\n')
2924    prf.write(pName+'.pgrid\n')
2925    prf.write(pName+'.fba\n')
2926    prf.write(pName+'_eps.raw\n')
2927    prf.write('%d\n'%MEMtype)
2928    if DysData['DenStart'] == 'uniform':
2929        prf.write('0\n')
2930    else:
2931        prf.write('1\n')
2932    if DysData['Optimize'] == 'ZSPA':
2933        prf.write('0\n')
2934    else:
2935        prf.write('1\n')
2936    prf.write('1\n')
2937    if DysData['Lagrange'][0] == 'user':
2938        prf.write('0\n')
2939    else:
2940        prf.write('1\n')
2941    prf.write('%.4f %d\n'%(DysData['Lagrange'][1],DysData['wt pwr']))
2942    prf.write('%.3f\n'%DysData['Lagrange'][2])
2943    prf.write('%.2f\n'%DysData['E_factor'])
2944    prf.write('1\n')
2945    prf.write('0\n')
2946    prf.write('%d\n'%DysData['Ncyc'])
2947    prf.write('1\n')
2948    prf.write('1 0 0 0 0 0 0 0\n')
2949    if DysData['prior'] == 'uniform':
2950        prf.write('0\n')
2951    else:
2952        prf.write('1\n')
2953    prf.close()
2954    return prfName
2955
2956def makeMEMfile(data,reflData,MEMtype,DYSNOMIA):
2957    ''' make Dysnomia .mem file of reflection data, etc.
2958
2959    :param dict data: GSAS-II phase data
2960    :param list reflData: GSAS-II reflection data
2961    :param int MEMtype: 1 for neutron data with negative scattering lengths
2962                        0 otherwise
2963    :param str DYSNOMIA: path to dysnomia.exe
2964    '''
2965   
2966    DysData = data['Dysnomia']
2967    generalData = data['General']
2968    cell = generalData['Cell'][1:7]
2969    A = G2lat.cell2A(cell)
2970    SGData = generalData['SGData']
2971    pName = generalData['Name'].replace(' ','_')
2972    memName = pName+'.mem'
2973    Map = generalData['Map']
2974    Type = Map['Type']
2975    UseList = Map['RefList']
2976    mem = open(memName,'w')
2977    mem.write('%s\n'%(generalData['Name']+' from '+UseList[0]))
2978    a,b,c,alp,bet,gam = cell
2979    mem.write('%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n'%(a,b,c,alp,bet,gam))
2980    mem.write('      0.0000000      0.0000000     -1    0    0    0     P\n')   #dummy PO stuff
2981    SGSym = generalData['SGData']['SpGrp']
2982    try:
2983        SGId = G2spc.spgbyNum.index(SGSym)
2984    except ValueError:
2985        return False
2986    org = 1
2987    if SGSym in G2spc.spg2origins:
2988        org = 2
2989    mapsize = Map['rho'].shape
2990    sumZ = 0.
2991    sumpos = 0.
2992    sumneg = 0.
2993    mem.write('%5d%5d%5d%5d%5d\n'%(SGId,org,mapsize[0],mapsize[1],mapsize[2]))
2994    for atm in generalData['NoAtoms']:
2995        Nat = generalData['NoAtoms'][atm]
2996        AtInfo = G2elem.GetAtomInfo(atm)
2997        sumZ += Nat*AtInfo['Z']
2998        isotope = generalData['Isotope'][atm]
2999        blen = generalData['Isotopes'][atm][isotope]['SL'][0]
3000        if blen < 0.:
3001            sumneg += blen*Nat
3002        else:
3003            sumpos += blen*Nat
3004    if 'X' in Type:
3005        mem.write('%10.2f  0.001\n'%sumZ)
3006    elif 'N' in Type and MEMtype:
3007        mem.write('%10.3f%10.3f 0.001\n'%(sumpos,sumneg))
3008    else:
3009        mem.write('%10.3f 0.001\n'%sumpos)
3010       
3011    dmin = DysData['MEMdmin']
3012    refSet = G2lat.GenHLaue(dmin,SGData,A)      #list of h,k,l,d
3013    refDict = {'%d %d %d'%(ref[0],ref[1],ref[2]):ref for ref in refSet}
3014       
3015    refs = []
3016    prevpos = 0.
3017    for ref in reflData:
3018        if ref[3] < 0:
3019            continue
3020        if 'T' in Type:
3021            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,x,x,x,prfo = ref[:16]
3022            FWHM = getgamFW(gam,sig)
3023        else:
3024            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,prfo = ref[:13]
3025            g = gam/100.    #centideg -> deg
3026            s = np.sqrt(max(sig,0.0001))/100.   #var -> sig in deg
3027            FWHM = getgamFW(g,s)
3028        delt = pos-prevpos
3029        refs.append([h,k,l,mult,pos,FWHM,Fobs,phase,delt])
3030        prevpos = pos
3031           
3032    ovlp = DysData['overlap']
3033    refs1 = []
3034    refs2 = []
3035    nref2 = 0
3036    iref = 0
3037    Nref = len(refs)
3038    start = False
3039    while iref < Nref-1:
3040        if refs[iref+1][-1] < ovlp*refs[iref][5]:
3041            if refs[iref][-1] > ovlp*refs[iref][5]:
3042                refs2.append([])
3043                start = True
3044            if nref2 == len(refs2):
3045                refs2.append([])
3046            refs2[nref2].append(refs[iref])
3047        else:
3048            if start:
3049                refs2[nref2].append(refs[iref])
3050                start = False
3051                nref2 += 1
3052            else:
3053                refs1.append(refs[iref])
3054        iref += 1
3055    if start:
3056        refs2[nref2].append(refs[iref])
3057    else:
3058        refs1.append(refs[iref])
3059   
3060    mem.write('%5d\n'%len(refs1))
3061    for ref in refs1:
3062        h,k,l = ref[:3]
3063        hkl = '%d %d %d'%(h,k,l)
3064        if hkl in refDict:
3065            del refDict[hkl]
3066        Fobs = np.sqrt(ref[6])
3067        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)))
3068    while True:
3069        if not len(refs2[-1]):
3070            del refs2[-1]
3071        else:
3072            break
3073    mem.write('%5d\n'%len(refs2))
3074    for iref2,ref2 in enumerate(refs2):
3075        mem.write('#%5d\n'%iref2)
3076        mem.write('%5d\n'%len(ref2))
3077        Gsum = 0.
3078        Msum = 0
3079        for ref in ref2:
3080            Gsum += ref[6]*ref[3]
3081            Msum += ref[3]
3082        G = np.sqrt(Gsum/Msum)
3083        h,k,l = ref2[0][:3]
3084        hkl = '%d %d %d'%(h,k,l)
3085        if hkl in refDict:
3086            del refDict[hkl]
3087        mem.write('%5d%5d%5d%10.3f%10.3f%5d\n'%(h,k,l,G,max(0.01*G,0.1),ref2[0][3]))
3088        for ref in ref2[1:]:
3089            h,k,l,m = ref[:4]
3090            mem.write('%5d%5d%5d%5d\n'%(h,k,l,m))
3091            hkl = '%d %d %d'%(h,k,l)
3092            if hkl in refDict:
3093                del refDict[hkl]
3094    if len(refDict):
3095        mem.write('%d\n'%len(refDict))
3096        for hkl in list(refDict.keys()):
3097            h,k,l = refDict[hkl][:3]
3098            mem.write('%5d%5d%5d\n'%(h,k,l))
3099    else:
3100        mem.write('0\n')
3101    mem.close()
3102    return True
3103
3104def MEMupdateReflData(prfName,data,reflData):
3105    ''' Update reflection data with new Fosq, phase result from Dysnomia
3106
3107    :param str prfName: phase.mem file name
3108    :param list reflData: GSAS-II reflection data
3109    '''
3110   
3111    generalData = data['General']
3112    cell = generalData['Cell'][1:7]
3113    A = G2lat.cell2A(cell)
3114    reflDict = {}
3115    newRefs = []
3116    for iref,ref in enumerate(reflData):
3117        if ref[3] > 0:
3118            newRefs.append(ref)
3119            reflDict[hash('%5d%5d%5d'%(ref[0],ref[1],ref[2]))] = iref
3120    fbaName = os.path.splitext(prfName)[0]+'.fba'
3121    try:
3122        fba = open(fbaName,'r')
3123    except FileNotFoundError:
3124        return False
3125    fba.readline()
3126    Nref = int(fba.readline()[:-1])
3127    fbalines = fba.readlines()
3128    for line in fbalines[:Nref]:
3129        info = line.split()
3130        h = int(info[0])
3131        k = int(info[1])
3132        l = int(info[2])
3133        FoR = float(info[3])
3134        FoI = float(info[4])
3135        Fosq = FoR**2+FoI**2
3136        phase = npatan2d(FoI,FoR)
3137        try:
3138            refId = reflDict[hash('%5d%5d%5d'%(h,k,l))]
3139        except KeyError:    #added reflections at end skipped
3140            d = float(1/np.sqrt(G2lat.calc_rDsq([h,k,l],A)))
3141            newRefs.append([h,k,l,-1,d,0.,0.01,1.0,Fosq,Fosq,phase,1.0,1.0,1.0,1.0])
3142            continue
3143        newRefs[refId][8] = Fosq
3144        newRefs[refId][10] = phase
3145    newRefs = np.array(newRefs)
3146    return True,newRefs
3147   
3148#testing data
3149NeedTestData = True
3150def TestData():
3151    'needs a doc string'
3152#    global NeedTestData
3153    global bakType
3154    bakType = 'chebyschev'
3155    global xdata
3156    xdata = np.linspace(4.0,40.0,36000)
3157    global parmDict0
3158    parmDict0 = {
3159        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
3160        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
3161        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
3162        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
3163        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'Z':0.0,'SH/L':0.002,
3164        'Back0':5.384,'Back1':-0.015,'Back2':.004,
3165        }
3166    global parmDict1
3167    parmDict1 = {
3168        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
3169        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
3170        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
3171        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
3172        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
3173        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
3174        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'Z':0.0,'SH/L':0.002,
3175        'Back0':36.897,'Back1':-0.508,'Back2':.006,
3176        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
3177        }
3178    global parmDict2
3179    parmDict2 = {
3180        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
3181        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'Z':0.0,'SH/L':0.02,
3182        'Back0':5.,'Back1':-0.02,'Back2':.004,
3183#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
3184        }
3185    global varyList
3186    varyList = []
3187
3188def test0():
3189    if NeedTestData: TestData()
3190    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
3191    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
3192    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
3193    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
3194    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
3195    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
3196   
3197def test1():
3198    if NeedTestData: TestData()
3199    time0 = time.time()
3200    for i in range(100):
3201        getPeakProfile(parmDict1,xdata,varyList,bakType)
3202    G2fil.G2Print ('100+6*Ka1-2 peaks=1200 peaks %.2f'%time.time()-time0)
3203   
3204def test2(name,delt):
3205    if NeedTestData: TestData()
3206    varyList = [name,]
3207    xdata = np.linspace(5.6,5.8,400)
3208    hplot = plotter.add('derivatives test for '+name).gca()
3209    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
3210    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
3211    parmDict2[name] += delt
3212    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
3213    hplot.plot(xdata,(y1-y0)/delt,'r+')
3214   
3215def test3(name,delt):
3216    if NeedTestData: TestData()
3217    names = ['pos','sig','gam','shl']
3218    idx = names.index(name)
3219    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
3220    xdata = np.linspace(5.6,5.8,800)
3221    dx = xdata[1]-xdata[0]
3222    hplot = plotter.add('derivatives test for '+name).gca()
3223    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
3224    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
3225    myDict[name] += delt
3226    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
3227    hplot.plot(xdata,(y1-y0)/delt,'r+')
3228
3229if __name__ == '__main__':
3230    import GSASIItestplot as plot
3231    global plotter
3232    plotter = plot.PlotNotebook()
3233#    test0()
3234#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)']:
3235    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
3236        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['Z',0.0001],['SH/L',0.00005]]:
3237        test2(name,shft)
3238    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
3239        test3(name,shft)
3240    G2fil.G2Print ("OK")
3241    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.