source: trunk/GSASIIpwd.py @ 2356

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

revision to RDF calculation/plot

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