source: trunk/GSASIIpwd.py @ 4142

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

fix panalytical reflectometry mporter - now converts 2theta to q
some fixes to reflectometry stuff

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