source: trunk/GSASIIpwd.py @ 2412

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

fixes to PDF stuff

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