source: trunk/GSASIIpwd.py @ 2651

Last change on this file since 2651 was 2651, checked in by vondreele, 5 years ago

in peakfit add binsperFWHM calc for dach peak fit with warnings if too small

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