source: trunk/GSASIIpwd.py @ 2660

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

fix issues with new pdf creation & fully trap missing elements problem
progress on pdf peak display, selection & limit changes
remove Add from background, dark & container PDF stuff

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