source: trunk/GSASIIpwd.py @ 4199

Last change on this file since 4199 was 4199, checked in by vondreele, 23 months ago

enhance MakeRMC6f for RMCProfile

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