source: trunk/GSASIIpwd.py @ 2688

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

sequential PDF peak fitting now implemented

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