source: trunk/GSASIIpwd.py @ 3435

Last change on this file since 3435 was 3435, checked in by vondreele, 3 years ago

make new routine GetSpGrpfromUser? & use it inGeneral & Transform fo space group input
modify getHKLpeak to check for magnetic space group extinctions - passes hkl that is allowed by either.
modify G2pwdGUI to retain SGData in Unit Cell data (in ssopts dict)

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