source: trunk/GSASIIpwd.py @ 3561

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

add Print operators option to SGMessageBox for showing operators on console
fix weighted diff plot for when background file is used for SASD & REFL data
Add obs-back & calc-back options to RFD plotting (Background)
fix bug in Sample parms when a substance was deleted (SASD)
fix bugs in LoadUntCells? & ImportUnitCells? - bad defaults
new routines in G2spc: TextOs? - makes formatted list of all symmetry operators & Text2MT: convert sym. op text to matrix form; works for either order ('x+1/2' or '1/2+x' styles)

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