source: trunk/GSASIIpwd.py @ 2522

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

Do Ivo's changes & some more work on mag derivatives.

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