source: trunk/GSASIIpwd.py @ 4200

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

Add visualization of RMCProfile results, i.e. G(R), S(Q), partials, etc.
MOre enhancements for RMCProfile stuff

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