source: trunk/GSASIIpwd.py @ 2681

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

fix missing 'Oblique' problem in G2imgGUI
fix inner 2-theta limit issue in OnTransferAngles?
PDF Peaks can now clear all peaks
PDF peak picking now gives pos,mag & sig=0.085 as default
calculated PDF peak fit curve now plotted on G(R) from PDF Peaks
fitting of PDF peaks implemented - works fine; needs text output

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