source: trunk/GSASIIpwd.py @ 4193

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

fix comment lines so easier to see in spyder outline
use a dummy powder import class for rd in the Load from defaults for instrument parameters

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