source: trunk/GSASIIpwd.py @ 2357

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

remove import .Pnn file
cleanup CalcRDF

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