source: trunk/GSASIIpwd.py @ 4195

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

add a find utility to GSASIIfiles for finding a file within a directory tree
faster version of ExpandCell? & define SwapItems? in GSASIIlattice
changes to RMCProfile interface

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 133.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-03 09:00:37 +0000 (Tue, 03 Dec 2019) $
10# $Author: vondreele $
11# $Revision: 4195 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 4195 2019-12-03 09:00:37Z 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: 4195 $")
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 = parmDict[pfx+'BkPkint;'+str(iD)]
932            pkS = parmDict[pfx+'BkPksig;'+str(iD)]
933            pkG = parmDict[pfx+'BkPkgam;'+str(iD)]
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 = parmDict[hfx+'BkPkint;'+str(iD)]
1065            pkS = parmDict[hfx+'BkPksig;'+str(iD)]
1066            pkG = parmDict[hfx+'BkPkgam;'+str(iD)]
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    Atoms = newPhase['Atoms']
2153    Cell = newPhase['General']['Cell'][1:7]
2154    fname = Name+'.rmc6f'
2155    fl = open(fname,'w')
2156    fl.write('(Version 6f format configuration file)\n')
2157    for item in Meta:
2158        fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2159    fl.write('%-35s%3d%3d%3d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2160    fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2161            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2162    A,B = G2lat. cell2AB(Cell)
2163    fl.write('Lattice vectors (Ang):\n')
2164    for i in [0,1,2]:
2165        fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2166    fl.write('Atoms (fractional coordinates):\n')
2167    Natm = np.core.defchararray.count(np.array(Atcodes),'+')
2168    Natm = np.count_nonzero(Natm-1)
2169    nat = 0
2170    for atm in Atseq:
2171        for iat,atom in enumerate(Atoms):
2172            if atom[1] == atm:
2173                nat += 1
2174                atcode = Atcodes[iat].split(':')
2175                cell = [0,0,0]
2176                if '+' in atcode[1]:
2177                    cell = eval(atcode[1].split('+')[1])
2178                fl.write('%6d%4s  [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(       
2179                        nat,atom[1],atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2180    fl.close()
2181    return fname
2182
2183def MakeBragg(G2frame,Name,Phase,PWId):
2184    PWDdata = G2frame.GetPWDRdatafromTree(PWId)
2185    generalData = Phase['General']
2186    Vol = generalData['Cell'][7]
2187    Data = PWDdata['Data']
2188    Inst = PWDdata['Instrument Parameters'][0]
2189    Bank = int(Inst['Bank'][1])
2190    Sample = PWDdata['Sample Parameters']
2191    Scale = Sample['Scale'][0]
2192    Limits = PWDdata['Limits'][1]
2193    Ibeg = np.searchsorted(Data[0],Limits[0])
2194    Ifin = np.searchsorted(Data[0],Limits[1])+1
2195    fname = Name+'.bragg'
2196    fl = open(fname,'w')
2197    fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-1,Bank,Scale,Vol))
2198    if 'T' in Inst['Type'][0]:
2199        fl.write('%12s%12s\n'%('   TOF,ms','  I(obs)'))
2200        for i in range(Ibeg,Ifin-1):
2201            fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2202    else:
2203        fl.write('%12s%12s\n'%('   2-theta, deg','  I(obs)'))
2204        for i in range(Ibeg,Ifin-1):
2205            fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))       
2206    fl.close()
2207    return fname
2208   
2209################################################################################
2210#### Reflectometry calculations
2211################################################################################
2212
2213def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
2214    G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
2215   
2216    class RandomDisplacementBounds(object):
2217        """random displacement with bounds"""
2218        def __init__(self, xmin, xmax, stepsize=0.5):
2219            self.xmin = xmin
2220            self.xmax = xmax
2221            self.stepsize = stepsize
2222   
2223        def __call__(self, x):
2224            """take a random step but ensure the new position is within the bounds"""
2225            while True:
2226                # this could be done in a much more clever way, but it will work for example purposes
2227                steps = self.xmax-self.xmin
2228                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
2229                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
2230                    break
2231            return xnew
2232   
2233    def GetModelParms():
2234        parmDict = {}
2235        varyList = []
2236        values = []
2237        bounds = []
2238        parmDict['dQ type'] = data['dQ type']
2239        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
2240        for parm in ['Scale','FltBack']:
2241            parmDict[parm] = data[parm][0]
2242            if data[parm][1]:
2243                varyList.append(parm)
2244                values.append(data[parm][0])
2245                bounds.append(Bounds[parm])
2246        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
2247        parmDict['nLayers'] = len(parmDict['Layer Seq'])
2248        for ilay,layer in enumerate(data['Layers']):
2249            name = layer['Name']
2250            cid = str(ilay)+';'
2251            parmDict[cid+'Name'] = name
2252            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2253                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
2254                if layer.get(parm,[0.,False])[1]:
2255                    varyList.append(cid+parm)
2256                    value = layer[parm][0]
2257                    values.append(value)
2258                    if value:
2259                        bound = [value*Bfac,value/Bfac]
2260                    else:
2261                        bound = [0.,10.]
2262                    bounds.append(bound)
2263            if name not in ['vacuum','unit scatter']:
2264                parmDict[cid+'rho'] = Substances[name]['Scatt density']
2265                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
2266        return parmDict,varyList,values,bounds
2267   
2268    def SetModelParms():
2269        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
2270        if 'Scale' in varyList:
2271            data['Scale'][0] = parmDict['Scale']
2272            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
2273        G2fil.G2Print (line)
2274        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
2275        if 'FltBack' in varyList:
2276            data['FltBack'][0] = parmDict['FltBack']
2277            line += ' esd: %15.3g'%(sigDict['FltBack'])
2278        G2fil.G2Print (line)
2279        for ilay,layer in enumerate(data['Layers']):
2280            name = layer['Name']
2281            G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
2282            cid = str(ilay)+';'
2283            line = ' '
2284            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
2285            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
2286            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2287                if parm in layer:
2288                    if parm == 'Rough':
2289                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
2290                    else:
2291                        layer[parm][0] = parmDict[cid+parm]
2292                    line += ' %s: %.3f'%(parm,layer[parm][0])
2293                    if cid+parm in varyList:
2294                        line += ' esd: %.3g'%(sigDict[cid+parm])
2295            G2fil.G2Print (line)
2296            G2fil.G2Print (line2)
2297   
2298    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2299        parmDict.update(zip(varyList,values))
2300        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2301        return M
2302   
2303    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2304        parmDict.update(zip(varyList,values))
2305        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2306        return np.sum(M**2)
2307   
2308    def getREFD(Q,Qsig,parmDict):
2309        Ic = np.ones_like(Q)*parmDict['FltBack']
2310        Scale = parmDict['Scale']
2311        Nlayers = parmDict['nLayers']
2312        Res = parmDict['Res']
2313        depth = np.zeros(Nlayers)
2314        rho = np.zeros(Nlayers)
2315        irho = np.zeros(Nlayers)
2316        sigma = np.zeros(Nlayers)
2317        for ilay,lay in enumerate(parmDict['Layer Seq']):
2318            cid = str(lay)+';'
2319            depth[ilay] = parmDict[cid+'Thick']
2320            sigma[ilay] = parmDict[cid+'Rough']
2321            if parmDict[cid+'Name'] == u'unit scatter':
2322                rho[ilay] = parmDict[cid+'DenMul']
2323                irho[ilay] = parmDict[cid+'iDenMul']
2324            elif 'vacuum' != parmDict[cid+'Name']:
2325                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
2326                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
2327            if cid+'Mag SLD' in parmDict:
2328                rho[ilay] += parmDict[cid+'Mag SLD']
2329        if parmDict['dQ type'] == 'None':
2330            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2331        elif 'const' in parmDict['dQ type']:
2332            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
2333        else:       #dQ/Q in data
2334            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
2335        Ic += AB*Scale
2336        return Ic
2337       
2338    def estimateT0(takestep):
2339        Mmax = -1.e-10
2340        Mmin = 1.e10
2341        for i in range(100):
2342            x0 = takestep(values)
2343            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2344            Mmin = min(M,Mmin)
2345            MMax = max(M,Mmax)
2346        return 1.5*(MMax-Mmin)
2347
2348    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2349    if data.get('2% weight'):
2350        wt = 1./(0.02*Io)**2
2351    Qmin = Limits[1][0]
2352    Qmax = Limits[1][1]
2353    wtFactor = ProfDict['wtFactor']
2354    Bfac = data['Toler']
2355    Ibeg = np.searchsorted(Q,Qmin)
2356    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
2357    Ic[:] = 0
2358    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
2359              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
2360    parmDict,varyList,values,bounds = GetModelParms()
2361    Msg = 'Failed to converge'
2362    if varyList:
2363        if data['Minimizer'] == 'LMLS': 
2364            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
2365                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2366            parmDict.update(zip(varyList,result[0]))
2367            chisq = np.sum(result[2]['fvec']**2)
2368            ncalc = result[2]['nfev']
2369            covM = result[1]
2370            newVals = result[0]
2371        elif data['Minimizer'] == 'Basin Hopping':
2372            xyrng = np.array(bounds).T
2373            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
2374            T0 = estimateT0(take_step)
2375            G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
2376            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
2377                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
2378                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
2379            chisq = result.fun
2380            ncalc = result.nfev
2381            newVals = result.x
2382            covM = []
2383        elif data['Minimizer'] == 'MC/SA Anneal':
2384            xyrng = np.array(bounds).T
2385            result = G2mth.anneal(sumREFD, values, 
2386                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
2387                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
2388                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
2389                ranRange=0.20,autoRan=False,dlg=None)
2390            newVals = result[0]
2391            parmDict.update(zip(varyList,newVals))
2392            chisq = result[1]
2393            ncalc = result[3]
2394            covM = []
2395            G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
2396        elif data['Minimizer'] == 'L-BFGS-B':
2397            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
2398                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2399            parmDict.update(zip(varyList,result['x']))
2400            chisq = result.fun
2401            ncalc = result.nfev
2402            newVals = result.x
2403            covM = []
2404    else:   #nothing varied
2405        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2406        chisq = np.sum(M**2)
2407        ncalc = 0
2408        covM = []
2409        sig = []
2410        sigDict = {}
2411        result = []
2412    Rvals = {}
2413    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
2414    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
2415    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
2416    Ib[Ibeg:Ifin] = parmDict['FltBack']
2417    try:
2418        if not len(varyList):
2419            Msg += ' - nothing refined'
2420            raise ValueError
2421        Nans = np.isnan(newVals)
2422        if np.any(Nans):
2423            name = varyList[Nans.nonzero(True)[0]]
2424            Msg += ' Nan result for '+name+'!'
2425            raise ValueError
2426        Negs = np.less_equal(newVals,0.)
2427        if np.any(Negs):
2428            indx = Negs.nonzero()
2429            name = varyList[indx[0][0]]
2430            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
2431                Msg += ' negative coefficient for '+name+'!'
2432                raise ValueError
2433        if len(covM):
2434            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
2435            covMatrix = covM*Rvals['GOF']
2436        else:
2437            sig = np.zeros(len(varyList))
2438            covMatrix = []
2439        sigDict = dict(zip(varyList,sig))
2440        G2fil.G2Print (' Results of reflectometry data modelling fit:')
2441        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
2442        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2443        SetModelParms()
2444        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
2445    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
2446        G2fil.G2Print (Msg)
2447        return False,0,0,0,0,0,0,Msg
2448       
2449def makeSLDprofile(data,Substances):
2450   
2451    sq2 = np.sqrt(2.)
2452    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2453    Nlayers = len(laySeq)
2454    laySeq = np.array(laySeq,dtype=int)
2455    interfaces = np.zeros(Nlayers)
2456    rho = np.zeros(Nlayers)
2457    sigma = np.zeros(Nlayers)
2458    name = data['Layers'][0]['Name']
2459    thick = 0.
2460    for ilay,lay in enumerate(laySeq):
2461        layer = data['Layers'][lay]
2462        name = layer['Name']
2463        if 'Thick' in layer:
2464            thick += layer['Thick'][0]
2465            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
2466        if 'Rough' in layer:
2467            sigma[ilay] = max(0.001,layer['Rough'][0])
2468        if name != 'vacuum':
2469            if name == 'unit scatter':
2470                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
2471            else:
2472                rrho = Substances[name]['Scatt density']
2473                irho = Substances[name]['XImag density']
2474                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
2475        if 'Mag SLD' in layer:
2476            rho[ilay] += layer['Mag SLD'][0]
2477    name = data['Layers'][-1]['Name']
2478    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
2479    xr = np.flipud(x)
2480    interfaces[-1] = x[-1]
2481    y = np.ones_like(x)*rho[0]
2482    iBeg = 0
2483    for ilayer in range(Nlayers-1):
2484        delt = rho[ilayer+1]-rho[ilayer]
2485        iPos = np.searchsorted(x,interfaces[ilayer])
2486        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
2487        iBeg = iPos
2488    return x,xr,y   
2489
2490def REFDModelFxn(Profile,Inst,Limits,Substances,data):
2491   
2492    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2493    Qmin = Limits[1][0]
2494    Qmax = Limits[1][1]
2495    iBeg = np.searchsorted(Q,Qmin)
2496    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2497    Ib[:] = data['FltBack'][0]
2498    Ic[:] = 0
2499    Scale = data['Scale'][0]
2500    if data['Layer Seq'] == []:
2501        return
2502    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2503    Nlayers = len(laySeq)
2504    depth = np.zeros(Nlayers)
2505    rho = np.zeros(Nlayers)
2506    irho = np.zeros(Nlayers)
2507    sigma = np.zeros(Nlayers)
2508    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
2509        layer = data['Layers'][lay]
2510        name = layer['Name']
2511        if 'Thick' in layer:    #skips first & last layers
2512            depth[ilay] = layer['Thick'][0]
2513        if 'Rough' in layer:    #skips first layer
2514            sigma[ilay] = layer['Rough'][0]
2515        if 'unit scatter' == name:
2516            rho[ilay] = layer['DenMul'][0]
2517            irho[ilay] = layer['iDenMul'][0]
2518        else:
2519            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2520            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
2521        if 'Mag SLD' in layer:
2522            rho[ilay] += layer['Mag SLD'][0]
2523    if data['dQ type'] == 'None':
2524        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2525    elif 'const' in data['dQ type']:
2526        res = data['Resolution'][0]/(100.*sateln2)
2527        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
2528    else:       #dQ/Q in data
2529        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
2530    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
2531
2532def abeles(kz, depth, rho, irho=0, sigma=0):
2533    """
2534    Optical matrix form of the reflectivity calculation.
2535    O.S. Heavens, Optical Properties of Thin Solid Films
2536   
2537    Reflectometry as a function of kz for a set of slabs.
2538
2539    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
2540        This is :math:`\\tfrac12 Q_z`.       
2541    :param depth:  float[m] (Ang).
2542        thickness of each layer.  The thickness of the incident medium
2543        and substrate are ignored.
2544    :param rho:  float[n,k] (1e-6/Ang^2)
2545        Real scattering length density for each layer for each kz
2546    :param irho:  float[n,k] (1e-6/Ang^2)
2547        Imaginary scattering length density for each layer for each kz
2548        Note: absorption cross section mu = 2 irho/lambda for neutrons
2549    :param sigma: float[m-1] (Ang)
2550        interfacial roughness.  This is the roughness between a layer
2551        and the previous layer. The sigma array should have m-1 entries.
2552
2553    Slabs are ordered with the surface SLD at index 0 and substrate at
2554    index -1, or reversed if kz < 0.
2555    """
2556    def calc(kz, depth, rho, irho, sigma):
2557        if len(kz) == 0: return kz
2558   
2559        # Complex index of refraction is relative to the incident medium.
2560        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
2561        # in place of kz^2, and ignoring rho_o
2562        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
2563        k = kz
2564   
2565        # According to Heavens, the initial matrix should be [ 1 F; F 1],
2566        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
2567        # multiply versus some coding convenience.
2568        B11 = 1
2569        B22 = 1
2570        B21 = 0
2571        B12 = 0
2572        for i in range(0, len(depth)-1):
2573            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
2574            F = (k - k_next) / (k + k_next)
2575            F *= np.exp(-2*k*k_next*sigma[i]**2)
2576            #print "==== layer",i
2577            #print "kz:", kz
2578            #print "k:", k
2579            #print "k_next:",k_next
2580            #print "F:",F
2581            #print "rho:",rho[:,i+1]
2582            #print "irho:",irho[:,i+1]
2583            #print "d:",depth[i],"sigma:",sigma[i]
2584            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
2585            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
2586            M21 = F*M11
2587            M12 = F*M22
2588            C1 = B11*M11 + B21*M12
2589            C2 = B11*M21 + B21*M22
2590            B11 = C1
2591            B21 = C2
2592            C1 = B12*M11 + B22*M12
2593            C2 = B12*M21 + B22*M22
2594            B12 = C1
2595            B22 = C2
2596            k = k_next
2597   
2598        r = B12/B11
2599        return np.absolute(r)**2
2600
2601    if np.isscalar(kz): kz = np.asarray([kz], 'd')
2602
2603    m = len(depth)
2604
2605    # Make everything into arrays
2606    depth = np.asarray(depth,'d')
2607    rho = np.asarray(rho,'d')
2608    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
2609    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
2610
2611    # Repeat rho,irho columns as needed
2612    if len(rho.shape) == 1:
2613        rho = rho[None,:]
2614        irho = irho[None,:]
2615
2616    return calc(kz, depth, rho, irho, sigma)
2617   
2618def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
2619    y = abeles(kz, depth, rho, irho, sigma)
2620    s = dq/2.
2621    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
2622    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
2623    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
2624    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
2625    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
2626    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
2627    y *= 0.137023
2628    return y
2629       
2630def makeRefdFFT(Limits,Profile):
2631    G2fil.G2Print ('make fft')
2632    Q,Io = Profile[:2]
2633    Qmin = Limits[1][0]
2634    Qmax = Limits[1][1]
2635    iBeg = np.searchsorted(Q,Qmin)
2636    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2637    Qf = np.linspace(0.,Q[iFin-1],5000)
2638    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
2639    If = QI(Qf)*Qf**4
2640    R = np.linspace(0.,5000.,5000)
2641    F = fft.rfft(If)
2642    return R,F
2643
2644   
2645################################################################################
2646#### Stacking fault simulation codes
2647################################################################################
2648
2649def GetStackParms(Layers):
2650   
2651    Parms = []
2652#cell parms
2653    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
2654        Parms.append('cellA')
2655        Parms.append('cellC')
2656    else:
2657        Parms.append('cellA')
2658        Parms.append('cellB')
2659        Parms.append('cellC')
2660        if Layers['Laue'] != 'mmm':
2661            Parms.append('cellG')
2662#Transition parms
2663    for iY in range(len(Layers['Layers'])):
2664        for iX in range(len(Layers['Layers'])):
2665            Parms.append('TransP;%d;%d'%(iY,iX))
2666            Parms.append('TransX;%d;%d'%(iY,iX))
2667            Parms.append('TransY;%d;%d'%(iY,iX))
2668            Parms.append('TransZ;%d;%d'%(iY,iX))
2669    return Parms
2670
2671def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
2672    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
2673   
2674    :param dict Layers: dict with following items
2675
2676      ::
2677
2678       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
2679       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
2680       'Layers':[],'Stacking':[],'Transitions':[]}
2681       
2682    :param str ctrls: controls string to be written on DIFFaX controls.dif file
2683    :param float scale: scale factor
2684    :param dict background: background parameters
2685    :param list limits: min/max 2-theta to be calculated
2686    :param dict inst: instrument parameters dictionary
2687    :param list profile: powder pattern data
2688   
2689    Note that parameters all updated in place   
2690    '''
2691    import atmdata
2692    path = sys.path
2693    for name in path:
2694        if 'bin' in name:
2695            DIFFaX = name+'/DIFFaX.exe'
2696            G2fil.G2Print (' Execute '+DIFFaX)
2697            break
2698    # make form factor file that DIFFaX wants - atom types are GSASII style
2699    sf = open('data.sfc','w')
2700    sf.write('GSASII special form factor file for DIFFaX\n\n')
2701    atTypes = list(Layers['AtInfo'].keys())
2702    if 'H' not in atTypes:
2703        atTypes.insert(0,'H')
2704    for atType in atTypes:
2705        if atType == 'H': 
2706            blen = -.3741
2707        else:
2708            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
2709        Adat = atmdata.XrayFF[atType]
2710        text = '%4s'%(atType.ljust(4))
2711        for i in range(4):
2712            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
2713        text += '%11.6f%11.6f'%(Adat['fc'],blen)
2714        text += '%3d\n'%(Adat['Z'])
2715        sf.write(text)
2716    sf.close()
2717    #make DIFFaX control.dif file - future use GUI to set some of these flags
2718    cf = open('control.dif','w')
2719    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2720        x0 = profile[0]
2721        iBeg = np.searchsorted(x0,limits[0])
2722        iFin = np.searchsorted(x0,limits[1])+1
2723        if iFin-iBeg > 20000:
2724            iFin = iBeg+20000
2725        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2726        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2727        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
2728    else:
2729        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2730        inst = {'Type':['XSC','XSC',]}
2731    cf.close()
2732    #make DIFFaX data file
2733    df = open('GSASII-DIFFaX.dat','w')
2734    df.write('INSTRUMENTAL\n')
2735    if 'X' in inst['Type'][0]:
2736        df.write('X-RAY\n')
2737    elif 'N' in inst['Type'][0]:
2738        df.write('NEUTRON\n')
2739    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2740        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
2741        U = ateln2*inst['U'][1]/10000.
2742        V = ateln2*inst['V'][1]/10000.
2743        W = ateln2*inst['W'][1]/10000.
2744        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2745        HW = np.sqrt(np.mean(HWHM))
2746    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
2747        if 'Mean' in Layers['selInst']:
2748            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
2749        elif 'Gaussian' in Layers['selInst']:
2750            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
2751        else:
2752            df.write('None\n')
2753    else:
2754        df.write('0.10\nNone\n')
2755    df.write('STRUCTURAL\n')
2756    a,b,c = Layers['Cell'][1:4]
2757    gam = Layers['Cell'][6]
2758    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
2759    laue = Layers['Laue']
2760    if laue == '2/m(ab)':
2761        laue = '2/m(1)'
2762    elif laue == '2/m(c)':
2763        laue = '2/m(2)'
2764    if 'unknown' in Layers['Laue']:
2765        df.write('%s %.3f\n'%(laue,Layers['Toler']))
2766    else:   
2767        df.write('%s\n'%(laue))
2768    df.write('%d\n'%(len(Layers['Layers'])))
2769    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
2770        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
2771    layerNames = []
2772    for layer in Layers['Layers']:
2773        layerNames.append(layer['Name'])
2774    for il,layer in enumerate(Layers['Layers']):
2775        if layer['SameAs']:
2776            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
2777            continue
2778        df.write('LAYER %d\n'%(il+1))
2779        if '-1' in layer['Symm']:
2780            df.write('CENTROSYMMETRIC\n')
2781        else:
2782            df.write('NONE\n')
2783        for ia,atom in enumerate(layer['Atoms']):
2784            [name,atype,x,y,z,frac,Uiso] = atom
2785            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
2786                frac /= 2.
2787            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
2788    df.write('STACKING\n')
2789    df.write('%s\n'%(Layers['Stacking'][0]))
2790    if 'recursive' in Layers['Stacking'][0]:
2791        df.write('%s\n'%Layers['Stacking'][1])
2792    else:
2793        if 'list' in Layers['Stacking'][1]:
2794            Slen = len(Layers['Stacking'][2])
2795            iB = 0
2796            iF = 0
2797            while True:
2798                iF += 68
2799                if iF >= Slen:
2800                    break
2801                iF = min(iF,Slen)
2802                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
2803                iB = iF
2804        else:
2805            df.write('%s\n'%Layers['Stacking'][1])   
2806    df.write('TRANSITIONS\n')
2807    for iY in range(len(Layers['Layers'])):
2808        sumPx = 0.
2809        for iX in range(len(Layers['Layers'])):
2810            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
2811            p = round(p,3)
2812            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
2813            sumPx += p
2814        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
2815            G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
2816            df.close()
2817            os.remove('data.sfc')
2818            os.remove('control.dif')
2819            os.remove('GSASII-DIFFaX.dat')
2820            return       
2821    df.close()   
2822    time0 = time.time()
2823    try:
2824        subp.call(DIFFaX)
2825    except OSError:
2826        G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
2827    G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
2828    if os.path.exists('GSASII-DIFFaX.spc'):
2829        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
2830        iFin = iBeg+Xpat.shape[1]
2831        bakType,backDict,backVary = SetBackgroundParms(background)
2832        backDict['Lam1'] = G2mth.getWave(inst)
2833        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2834        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
2835        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2836            rv = st.poisson(profile[3][iBeg:iFin])
2837            profile[1][iBeg:iFin] = rv.rvs()
2838            Z = np.ones_like(profile[3][iBeg:iFin])
2839            Z[1::2] *= -1
2840            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2841            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2842        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2843    #cleanup files..
2844        os.remove('GSASII-DIFFaX.spc')
2845    elif os.path.exists('GSASII-DIFFaX.sadp'):
2846        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
2847        Sadp = np.reshape(Sadp,(256,-1))
2848        Layers['Sadp']['Img'] = Sadp
2849        os.remove('GSASII-DIFFaX.sadp')
2850    os.remove('data.sfc')
2851    os.remove('control.dif')
2852    os.remove('GSASII-DIFFaX.dat')
2853   
2854def SetPWDRscan(inst,limits,profile):
2855   
2856    wave = G2mth.getMeanWave(inst)
2857    x0 = profile[0]
2858    iBeg = np.searchsorted(x0,limits[0])
2859    iFin = np.searchsorted(x0,limits[1])
2860    if iFin-iBeg > 20000:
2861        iFin = iBeg+20000
2862    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2863    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
2864    return iFin-iBeg
2865       
2866def SetStackingSF(Layers,debug):
2867# Load scattering factors into DIFFaX arrays
2868    import atmdata
2869    atTypes = Layers['AtInfo'].keys()
2870    aTypes = []
2871    for atype in atTypes:
2872        aTypes.append('%4s'%(atype.ljust(4)))
2873    SFdat = []
2874    for atType in atTypes:
2875        Adat = atmdata.XrayFF[atType]
2876        SF = np.zeros(9)
2877        SF[:8:2] = Adat['fa']
2878        SF[1:8:2] = Adat['fb']
2879        SF[8] = Adat['fc']
2880        SFdat.append(SF)
2881    SFdat = np.array(SFdat)
2882    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
2883   
2884def SetStackingClay(Layers,Type):
2885# Controls
2886    rand.seed()
2887    ranSeed = rand.randint(1,2**16-1)
2888    try:   
2889        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
2890            '6/m','6/mmm'].index(Layers['Laue'])+1
2891    except ValueError:  #for 'unknown'
2892        laueId = -1
2893    if 'SADP' in Type:
2894        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
2895        lmax = int(Layers['Sadp']['Lmax'])
2896    else:
2897        planeId = 0
2898        lmax = 0
2899# Sequences
2900    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
2901    try:
2902        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
2903    except ValueError:
2904        StkParm = -1
2905    if StkParm == 2:    #list
2906        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
2907        Nstk = len(StkSeq)
2908    else:
2909        Nstk = 1
2910        StkSeq = [0,]
2911    if StkParm == -1:
2912        StkParm = int(Layers['Stacking'][1])
2913    Wdth = Layers['Width'][0]
2914    mult = 1
2915    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
2916    LaueSym = Layers['Laue'].ljust(12)
2917    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
2918    return laueId,controls
2919   
2920def SetCellAtoms(Layers):
2921    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
2922# atoms in layers
2923    atTypes = list(Layers['AtInfo'].keys())
2924    AtomXOU = []
2925    AtomTp = []
2926    LayerSymm = []
2927    LayerNum = []
2928    layerNames = []
2929    Natm = 0
2930    Nuniq = 0
2931    for layer in Layers['Layers']:
2932        layerNames.append(layer['Name'])
2933    for il,layer in enumerate(Layers['Layers']):
2934        if layer['SameAs']:
2935            LayerNum.append(layerNames.index(layer['SameAs'])+1)
2936            continue
2937        else:
2938            LayerNum.append(il+1)
2939            Nuniq += 1
2940        if '-1' in layer['Symm']:
2941            LayerSymm.append(1)
2942        else:
2943            LayerSymm.append(0)
2944        for ia,atom in enumerate(layer['Atoms']):
2945            [name,atype,x,y,z,frac,Uiso] = atom
2946            Natm += 1
2947            AtomTp.append('%4s'%(atype.ljust(4)))
2948            Ta = atTypes.index(atype)+1
2949            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
2950    AtomXOU = np.array(AtomXOU)
2951    Nlayers = len(layerNames)
2952    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
2953    return Nlayers
2954   
2955def SetStackingTrans(Layers,Nlayers):
2956# Transitions
2957    TransX = []
2958    TransP = []
2959    for Ytrans in Layers['Transitions']:
2960        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
2961        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
2962    TransP = np.array(TransP,dtype='float').T
2963    TransX = np.array(TransX,dtype='float')
2964#    GSASIIpath.IPyBreak()
2965    pyx.pygettrans(Nlayers,TransP,TransX)
2966   
2967def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
2968# Scattering factors
2969    SetStackingSF(Layers,debug)
2970# Controls & sequences
2971    laueId,controls = SetStackingClay(Layers,'PWDR')
2972# cell & atoms
2973    Nlayers = SetCellAtoms(Layers)
2974    Volume = Layers['Cell'][7]   
2975# Transitions
2976    SetStackingTrans(Layers,Nlayers)
2977# PWDR scan
2978    Nsteps = SetPWDRscan(inst,limits,profile)
2979# result as Spec
2980    x0 = profile[0]
2981    profile[3] = np.zeros(len(profile[0]))
2982    profile[4] = np.zeros(len(profile[0]))
2983    profile[5] = np.zeros(len(profile[0]))
2984    iBeg = np.searchsorted(x0,limits[0])
2985    iFin = np.searchsorted(x0,limits[1])+1
2986    if iFin-iBeg > 20000:
2987        iFin = iBeg+20000
2988    Nspec = 20001       
2989    spec = np.zeros(Nspec,dtype='double')   
2990    time0 = time.time()
2991    pyx.pygetspc(controls,Nspec,spec)
2992    G2fil.G2Print (' GETSPC time = %.2fs'%(time.time()-time0))
2993    time0 = time.time()
2994    U = ateln2*inst['U'][1]/10000.
2995    V = ateln2*inst['V'][1]/10000.
2996    W = ateln2*inst['W'][1]/10000.
2997    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2998    HW = np.sqrt(np.mean(HWHM))
2999    BrdSpec = np.zeros(Nsteps)
3000    if 'Mean' in Layers['selInst']:
3001        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
3002    elif 'Gaussian' in Layers['selInst']:
3003        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
3004    else:
3005        BrdSpec = spec[:Nsteps]
3006    BrdSpec /= Volume
3007    iFin = iBeg+Nsteps
3008    bakType,backDict,backVary = SetBackgroundParms(background)
3009    backDict['Lam1'] = G2mth.getWave(inst)
3010    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
3011    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
3012    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
3013        try:
3014            rv = st.poisson(profile[3][iBeg:iFin])
3015            profile[1][iBeg:iFin] = rv.rvs()
3016        except ValueError:
3017            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
3018        Z = np.ones_like(profile[3][iBeg:iFin])
3019        Z[1::2] *= -1
3020        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3021        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3022    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3023    G2fil.G2Print (' Broadening time = %.2fs'%(time.time()-time0))
3024   
3025def CalcStackingSADP(Layers,debug):
3026   
3027# Scattering factors
3028    SetStackingSF(Layers,debug)
3029# Controls & sequences
3030    laueId,controls = SetStackingClay(Layers,'SADP')
3031# cell & atoms
3032    Nlayers = SetCellAtoms(Layers)   
3033# Transitions
3034    SetStackingTrans(Layers,Nlayers)
3035# result as Sadp
3036    Nspec = 20001       
3037    spec = np.zeros(Nspec,dtype='double')   
3038    time0 = time.time()
3039    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
3040    Sapd = np.zeros((256,256))
3041    iB = 0
3042    for i in range(hkLim):
3043        iF = iB+Nblk
3044        p1 = 127+int(i*Incr)
3045        p2 = 128-int(i*Incr)
3046        if Nblk == 128:
3047            if i:
3048                Sapd[128:,p1] = spec[iB:iF]
3049                Sapd[:128,p1] = spec[iF:iB:-1]
3050            Sapd[128:,p2] = spec[iB:iF]
3051            Sapd[:128,p2] = spec[iF:iB:-1]
3052        else:
3053            if i:
3054                Sapd[:,p1] = spec[iB:iF]
3055            Sapd[:,p2] = spec[iB:iF]
3056        iB += Nblk
3057    Layers['Sadp']['Img'] = Sapd
3058    G2fil.G2Print (' GETSAD time = %.2fs'%(time.time()-time0))
3059   
3060###############################################################################
3061#### Maximum Entropy Method - Dysnomia
3062###############################################################################
3063   
3064def makePRFfile(data,MEMtype):
3065    ''' makes Dysnomia .prf control file from Dysnomia GUI controls
3066   
3067    :param dict data: GSAS-II phase data
3068    :param int MEMtype: 1 for neutron data with negative scattering lengths
3069                        0 otherwise
3070    :returns str: name of Dysnomia control file
3071    '''
3072
3073    generalData = data['General']
3074    pName = generalData['Name'].replace(' ','_')
3075    DysData = data['Dysnomia']
3076    prfName = pName+'.prf'
3077    prf = open(prfName,'w')
3078    prf.write('$PREFERENCES\n')
3079    prf.write(pName+'.mem\n') #or .fos?
3080    prf.write(pName+'.out\n')
3081    prf.write(pName+'.pgrid\n')
3082    prf.write(pName+'.fba\n')
3083    prf.write(pName+'_eps.raw\n')
3084    prf.write('%d\n'%MEMtype)
3085    if DysData['DenStart'] == 'uniform':
3086        prf.write('0\n')
3087    else:
3088        prf.write('1\n')
3089    if DysData['Optimize'] == 'ZSPA':
3090        prf.write('0\n')
3091    else:
3092        prf.write('1\n')
3093    prf.write('1\n')
3094    if DysData['Lagrange'][0] == 'user':
3095        prf.write('0\n')
3096    else:
3097        prf.write('1\n')
3098    prf.write('%.4f %d\n'%(DysData['Lagrange'][1],DysData['wt pwr']))
3099    prf.write('%.3f\n'%DysData['Lagrange'][2])
3100    prf.write('%.2f\n'%DysData['E_factor'])
3101    prf.write('1\n')
3102    prf.write('0\n')
3103    prf.write('%d\n'%DysData['Ncyc'])
3104    prf.write('1\n')
3105    prf.write('1 0 0 0 0 0 0 0\n')
3106    if DysData['prior'] == 'uniform':
3107        prf.write('0\n')
3108    else:
3109        prf.write('1\n')
3110    prf.close()
3111    return prfName
3112
3113def makeMEMfile(data,reflData,MEMtype,DYSNOMIA):
3114    ''' make Dysnomia .mem file of reflection data, etc.
3115
3116    :param dict data: GSAS-II phase data
3117    :param list reflData: GSAS-II reflection data
3118    :param int MEMtype: 1 for neutron data with negative scattering lengths
3119                        0 otherwise
3120    :param str DYSNOMIA: path to dysnomia.exe
3121    '''
3122   
3123    DysData = data['Dysnomia']
3124    generalData = data['General']
3125    cell = generalData['Cell'][1:7]
3126    A = G2lat.cell2A(cell)
3127    SGData = generalData['SGData']
3128    pName = generalData['Name'].replace(' ','_')
3129    memName = pName+'.mem'
3130    Map = generalData['Map']
3131    Type = Map['Type']
3132    UseList = Map['RefList']
3133    mem = open(memName,'w')
3134    mem.write('%s\n'%(generalData['Name']+' from '+UseList[0]))
3135    a,b,c,alp,bet,gam = cell
3136    mem.write('%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n'%(a,b,c,alp,bet,gam))
3137    mem.write('      0.0000000      0.0000000     -1    0    0    0     P\n')   #dummy PO stuff
3138    SGSym = generalData['SGData']['SpGrp']
3139    try:
3140        SGId = G2spc.spgbyNum.index(SGSym)
3141    except ValueError:
3142        return False
3143    org = 1
3144    if SGSym in G2spc.spg2origins:
3145        org = 2
3146    mapsize = Map['rho'].shape
3147    sumZ = 0.
3148    sumpos = 0.
3149    sumneg = 0.
3150    mem.write('%5d%5d%5d%5d%5d\n'%(SGId,org,mapsize[0],mapsize[1],mapsize[2]))
3151    for atm in generalData['NoAtoms']:
3152        Nat = generalData['NoAtoms'][atm]
3153        AtInfo = G2elem.GetAtomInfo(atm)
3154        sumZ += Nat*AtInfo['Z']
3155        isotope = generalData['Isotope'][atm]
3156        blen = generalData['Isotopes'][atm][isotope]['SL'][0]
3157        if blen < 0.:
3158            sumneg += blen*Nat
3159        else:
3160            sumpos += blen*Nat
3161    if 'X' in Type:
3162        mem.write('%10.2f  0.001\n'%sumZ)
3163    elif 'N' in Type and MEMtype:
3164        mem.write('%10.3f%10.3f 0.001\n'%(sumpos,sumneg))
3165    else:
3166        mem.write('%10.3f 0.001\n'%sumpos)
3167       
3168    dmin = DysData['MEMdmin']
3169    TOFlam = 2.0*dmin*npsind(80.0)
3170    refSet = G2lat.GenHLaue(dmin,SGData,A)      #list of h,k,l,d
3171    refDict = {'%d %d %d'%(ref[0],ref[1],ref[2]):ref for ref in refSet}
3172       
3173    refs = []
3174    prevpos = 0.
3175    for ref in reflData:
3176        if ref[3] < 0:
3177            continue
3178        if 'T' in Type:
3179            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,x,x,x,prfo = ref[:16]
3180            s = np.sqrt(max(sig,0.0001))   #var -> sig in deg
3181            FWHM = getgamFW(gam,s)
3182            if dsp < dmin:
3183                continue
3184            theta = npasind(TOFlam/(2.*dsp))
3185            FWHM *= nptand(theta)/pos
3186            pos = 2.*theta
3187        else:
3188            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,prfo = ref[:13]
3189            g = gam/100.    #centideg -> deg
3190            s = np.sqrt(max(sig,0.0001))/100.   #var -> sig in deg
3191            FWHM = getgamFW(g,s)
3192        delt = pos-prevpos
3193        refs.append([h,k,l,mult,pos,FWHM,Fobs,phase,delt])
3194        prevpos = pos
3195           
3196    ovlp = DysData['overlap']
3197    refs1 = []
3198    refs2 = []
3199    nref2 = 0
3200    iref = 0
3201    Nref = len(refs)
3202    start = False
3203    while iref < Nref-1:
3204        if refs[iref+1][-1] < ovlp*refs[iref][5]:
3205            if refs[iref][-1] > ovlp*refs[iref][5]:
3206                refs2.append([])
3207                start = True
3208            if nref2 == len(refs2):
3209                refs2.append([])
3210            refs2[nref2].append(refs[iref])
3211        else:
3212            if start:
3213                refs2[nref2].append(refs[iref])
3214                start = False
3215                nref2 += 1
3216            else:
3217                refs1.append(refs[iref])
3218        iref += 1
3219    if start:
3220        refs2[nref2].append(refs[iref])
3221    else:
3222        refs1.append(refs[iref])
3223   
3224    mem.write('%5d\n'%len(refs1))
3225    for ref in refs1:
3226        h,k,l = ref[:3]
3227        hkl = '%d %d %d'%(h,k,l)
3228        if hkl in refDict:
3229            del refDict[hkl]
3230        Fobs = np.sqrt(ref[6])
3231        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)))
3232    while True and nref2:
3233        if not len(refs2[-1]):
3234            del refs2[-1]
3235        else:
3236            break
3237    mem.write('%5d\n'%len(refs2))
3238    for iref2,ref2 in enumerate(refs2):
3239        mem.write('#%5d\n'%iref2)
3240        mem.write('%5d\n'%len(ref2))
3241        Gsum = 0.
3242        Msum = 0
3243        for ref in ref2:
3244            Gsum += ref[6]*ref[3]
3245            Msum += ref[3]
3246        G = np.sqrt(Gsum/Msum)
3247        h,k,l = ref2[0][:3]
3248        hkl = '%d %d %d'%(h,k,l)
3249        if hkl in refDict:
3250            del refDict[hkl]
3251        mem.write('%5d%5d%5d%10.3f%10.3f%5d\n'%(h,k,l,G,max(0.01*G,0.1),ref2[0][3]))
3252        for ref in ref2[1:]:
3253            h,k,l,m = ref[:4]
3254            mem.write('%5d%5d%5d%5d\n'%(h,k,l,m))
3255            hkl = '%d %d %d'%(h,k,l)
3256            if hkl in refDict:
3257                del refDict[hkl]
3258    if len(refDict):
3259        mem.write('%d\n'%len(refDict))
3260        for hkl in list(refDict.keys()):
3261            h,k,l = refDict[hkl][:3]
3262            mem.write('%5d%5d%5d\n'%(h,k,l))
3263    else:
3264        mem.write('0\n')
3265    mem.close()
3266    return True
3267
3268def MEMupdateReflData(prfName,data,reflData):
3269    ''' Update reflection data with new Fosq, phase result from Dysnomia
3270
3271    :param str prfName: phase.mem file name
3272    :param list reflData: GSAS-II reflection data
3273    '''
3274   
3275    generalData = data['General']
3276    Map = generalData['Map']
3277    Type = Map['Type']
3278    cell = generalData['Cell'][1:7]
3279    A = G2lat.cell2A(cell)
3280    reflDict = {}
3281    newRefs = []
3282    for iref,ref in enumerate(reflData):
3283        if ref[3] > 0:
3284            newRefs.append(ref)
3285            reflDict[hash('%5d%5d%5d'%(ref[0],ref[1],ref[2]))] = iref
3286    fbaName = os.path.splitext(prfName)[0]+'.fba'
3287    try:
3288        fba = open(fbaName,'r')
3289    except FileNotFoundError:
3290        return False
3291    fba.readline()
3292    Nref = int(fba.readline()[:-1])
3293    fbalines = fba.readlines()
3294    for line in fbalines[:Nref]:
3295        info = line.split()
3296        h = int(info[0])
3297        k = int(info[1])
3298        l = int(info[2])
3299        FoR = float(info[3])
3300        FoI = float(info[4])
3301        Fosq = FoR**2+FoI**2
3302        phase = npatan2d(FoI,FoR)
3303        try:
3304            refId = reflDict[hash('%5d%5d%5d'%(h,k,l))]
3305        except KeyError:    #added reflections at end skipped
3306            d = float(1/np.sqrt(G2lat.calc_rDsq([h,k,l],A)))
3307            if 'T' in Type:
3308                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])
3309            else:
3310                newRefs.append([h,k,l,-1,d,0.,0.01,1.0,Fosq,Fosq,phase,1.0,1.0,1.0,1.0])
3311            continue
3312        newRefs[refId][8] = Fosq
3313        newRefs[refId][10] = phase
3314    newRefs = np.array(newRefs)
3315    return True,newRefs
3316   
3317#### testing data
3318NeedTestData = True
3319def TestData():
3320    'needs a doc string'
3321#    global NeedTestData
3322    global bakType
3323    bakType = 'chebyschev'
3324    global xdata
3325    xdata = np.linspace(4.0,40.0,36000)
3326    global parmDict0
3327    parmDict0 = {
3328        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
3329        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
3330        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
3331        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
3332        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'Z':0.0,'SH/L':0.002,
3333        'Back0':5.384,'Back1':-0.015,'Back2':.004,
3334        }
3335    global parmDict1
3336    parmDict1 = {
3337        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
3338        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
3339        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
3340        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
3341        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
3342        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
3343        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'Z':0.0,'SH/L':0.002,
3344        'Back0':36.897,'Back1':-0.508,'Back2':.006,
3345        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
3346        }
3347    global parmDict2
3348    parmDict2 = {
3349        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
3350        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'Z':0.0,'SH/L':0.02,
3351        'Back0':5.,'Back1':-0.02,'Back2':.004,
3352#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
3353        }
3354    global varyList
3355    varyList = []
3356
3357def test0():
3358    if NeedTestData: TestData()
3359    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
3360    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
3361    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
3362    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
3363    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
3364    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
3365   
3366def test1():
3367    if NeedTestData: TestData()
3368    time0 = time.time()
3369    for i in range(100):
3370        getPeakProfile(parmDict1,xdata,varyList,bakType)
3371    G2fil.G2Print ('100+6*Ka1-2 peaks=1200 peaks %.2f'%time.time()-time0)
3372   
3373def test2(name,delt):
3374    if NeedTestData: TestData()
3375    varyList = [name,]
3376    xdata = np.linspace(5.6,5.8,400)
3377    hplot = plotter.add('derivatives test for '+name).gca()
3378    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
3379    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
3380    parmDict2[name] += delt
3381    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
3382    hplot.plot(xdata,(y1-y0)/delt,'r+')
3383   
3384def test3(name,delt):
3385    if NeedTestData: TestData()
3386    names = ['pos','sig','gam','shl']
3387    idx = names.index(name)
3388    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
3389    xdata = np.linspace(5.6,5.8,800)
3390    dx = xdata[1]-xdata[0]
3391    hplot = plotter.add('derivatives test for '+name).gca()
3392    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
3393    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
3394    myDict[name] += delt
3395    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
3396    hplot.plot(xdata,(y1-y0)/delt,'r+')
3397
3398if __name__ == '__main__':
3399    import GSASIItestplot as plot
3400    global plotter
3401    plotter = plot.PlotNotebook()
3402#    test0()
3403#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)']:
3404    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
3405        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['Z',0.0001],['SH/L',0.00005]]:
3406        test2(name,shft)
3407    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
3408        test3(name,shft)
3409    G2fil.G2Print ("OK")
3410    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.