source: trunk/GSASIIpwd.py @ 3774

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

fix super indexing problem in transposeHKLF
fix reflection generation for incommensurate mag case in G2lattice & G2pwd
clean up non Fourier modulation calcs & remove analytic derivative stuff (now numeric)
fix uij derivative bug
work on incommensurate magnetic sturcture factors - not working yet
clean up of testDeriv - better choices for delt & force reflection regeneration before each derivative test

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