source: trunk/GSASIIpwd.py @ 3157

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

add Lorentzian term 'Z' to CW profile instrument parameters - constant term
change sig-q function to d*sig-q from sig-q/d2 for TOF profile instrument parameters
add filter to SumDialog? for PWDR & IMG data

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