source: trunk/GSASIIpwd.py @ 3136

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

make GSAS-II python 3.6 compliant & preserve python 2.7 use;changes:
do from future import division, print_function for all GSAS-II py sources
all menu items revised to be py 2.7/3.6 compliant
all wx.OPEN --> wx.FD_OPEN in file dialogs
all integer divides (typically for image pixel math) made explicit with ; ambiguous ones made floats as appropriate
all print "stuff" --> print (stuff)
all print >> pFile,'stuff' --> pFile.writeCIFtemplate('stuff')
all read file opens made explicit 'r' or 'rb'
all cPickle imports made for py2.7 or 3.6 as cPickle or _pickle; test for '2' platform.version_tuple[0] for py 2.7
define cPickleload to select load(fp) or load(fp,encoding='latin-1') for loading gpx files; provides cross compatibility between py 2.7/3.6 gpx files
make dict.keys() as explicit list(dict.keys()) as needed (NB: possible source of remaining py3.6 bugs)
make zip(a,b) as explicit list(zip(a,b)) as needed (NB: possible source of remaining py3.6 bugs)
select unichr/chr according test for '2' platform.version_tuple[0] for py 2.7 (G2pwdGUI * G2plot) for special characters
select wg.EVT_GRID_CELL_CHANGE (classic) or wg.EVT_GRID_CELL_CHANGED (phoenix) in grid Bind
maxint --> maxsize; used in random number stuff
raise Exception,"stuff" --> raise Exception("stuff")
wx 'classic' sizer.DeleteWindows?() or 'phoenix' sizer.Clear(True)
wx 'classic' SetToolTipString?(text) or 'phoenix' SetToolTip?(wx.ToolTip?(text)); define SetToolTipString?(self,text) to handle the choice in plots
status.SetFields? --> status.SetStatusText?
'classic' AddSimpleTool? or 'phoenix' self.AddTool? for plot toolbar; Bind different as well
define GetItemPydata? as it doesn't exist in wx 'phoenix'
allow python versions 2.7 & 3.6 to run GSAS-II
Bind override commented out - no logging capability (NB: remove all logging code?)
all import ContentsValidator? open filename & test if valid then close; filepointer removed from Reader
binary importers (mostly images) test for 'byte' type & convert as needed to satisfy py 3.6 str/byte rules

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