source: trunk/GSASIIpwd.py @ 4192

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

fix OnFileReopen?; broke when some file were deleted beforehand
Some preliminaries for rmc stuff
allow use of default instrument parameters upon Load command - same idea; cancel brings up default list choice

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