source: trunk/GSASIIpwd.py @ 2493

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

add display of phase weight fraction to DData tab

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