source: trunk/GSASIIpwd.py @ 3444

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

fix problem of masked xdata values as limits to linspace for background functions & derivatives; use the data instead of the masked data

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