source: trunk/GSASIIpwd.py @ 2546

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

cleanup of all GSASII main routines - remove unused variables, dead code, etc.

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