source: trunk/GSASIIpwd.py @ 3231

Last change on this file since 3231 was 3231, checked in by vondreele, 4 years ago

fix bug in IndexPeaks? display
fix import pypowder issue in G2math
fix ssSpaceGroup call in G2obj
add a dmax to Pawley controls; eliminate reflections with d>dmax (useful for TOF)
change reflist to be sorted as found in tree; used for Fourier, Charge flip & MCSA controls
add a buttonHandler to PlotXYZ to retrieve cursor position on generic 2D plots
use buttonHandler to get 2D modulation vector from plot & try it for indexing
fix a TOF modulation bug in getPowderProfile

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