source: trunk/GSASIIpwd.py @ 2561

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

add Rmax to PDF Controls; default = 100.; will not be exact
PDF always generates 5000 points for 0-Rmax; independent of chosen Qmax or image binning.
Trap attempt to calculate PDF without chemical formula - now get error message.

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