source: trunk/GSASIIpwd.py @ 4002

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

finish implementation of generating added reflections for Dysnomia. Will currently not show up in reflection list, but output *.pgrid can be viewed by the VESTA program

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