source: trunk/GSASIIpwd.py @ 4191

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

put in dummy phase tab & menu for fullrmc
partial implementation of pdf from TOF data; lacks absorption & wavelength dependent scattering length effects.

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