source: trunk/GSASIIpwd.py @ 3991

Last change on this file since 3991 was 3991, checked in by vondreele, 2 years ago

process Dysnomia reflection table result saving new Fosq & phase in ReflData?. Then compute new fourier map.

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