source: trunk/GSASIIpwd.py @ 2206

Last change on this file since 2206 was 2206, checked in by vondreele, 6 years ago

find errors in stacking fortran - replace both binwin directories
principal problem - transition probability matrix transposed in G2pwd
& HW needed to be sqrt(HW)
setup a debug mode for stacking fault stuff
G2plot - comment if page.Context lines - caused problems

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