source: trunk/GSASIIpwd.py @ 2600

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

remove comment about mag structures under construction
set powder peak widths to be 50*FWHM for TOF & 75/50 for CW (bigger on low side for axial divergence)
Tune up density map color rendering - make "red" really orange (brighter) & both orange/green more intense.

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