source: trunk/GSASIIpwd.py @ 4021

Last change on this file since 4021 was 4021, checked in by toby, 3 years ago

implement filter for screen messages; start to replace print() with G2Print(); scripting changes for PDF

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