source: trunk/GSASIIpwd.py @ 3164

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

fix neutron resonant scattering factor calcs. - neglected constant term
put tyr/except around import pypowder in G2pwd & G2math - to allow testing

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