source: trunk/GSASIIpwd.py @ 2777

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

implement import of getPDFx G(R) files NB: these have no matching PWDR entries
replace all scipy.fft with numpy.fft
add a plot SLD button for reflectometry
for PDF Peaks - Atom elements from periodic table - Bond No. still not working so doesn't really matter

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 109.5 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASII powder calculation module*
5==================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2017-04-12 20:12:45 +0000 (Wed, 12 Apr 2017) $
10# $Author: vondreele $
11# $Revision: 2777 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 2777 2017-04-12 20:12:45Z 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
26import numpy.fft as fft
27import scipy.interpolate as si
28import scipy.stats as st
29import scipy.optimize as so
30import scipy.special as sp
31
32import GSASIIpath
33GSASIIpath.SetVersionNumber("$Revision: 2777 $")
34import GSASIIlattice as G2lat
35import GSASIIspc as G2spc
36import GSASIIElem as G2elem
37import GSASIImath as G2mth
38import pypowder as pyd
39try:
40    import pydiffax as pyx
41except ImportError:
42    print 'pydiffax is not available for this platform - under develpment'
43
44   
45# trig functions in degrees
46tand = lambda x: math.tan(x*math.pi/180.)
47atand = lambda x: 180.*math.atan(x)/math.pi
48atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
49cosd = lambda x: math.cos(x*math.pi/180.)
50acosd = lambda x: 180.*math.acos(x)/math.pi
51rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
52#numpy versions
53npsind = lambda x: np.sin(x*np.pi/180.)
54npasind = lambda x: 180.*np.arcsin(x)/math.pi
55npcosd = lambda x: np.cos(x*math.pi/180.)
56npacosd = lambda x: 180.*np.arccos(x)/math.pi
57nptand = lambda x: np.tan(x*math.pi/180.)
58npatand = lambda x: 180.*np.arctan(x)/np.pi
59npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
60npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave    #=d*
61npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)    #=2pi*d*
62ateln2 = 8.0*math.log(2.0)
63nxs = np.newaxis
64
65################################################################################
66#### Powder utilities
67################################################################################
68
69def PhaseWtSum(G2frame,histo):
70    '''
71    Calculate sum of phase mass*shase fraction for PWDR data (exclude magnetic phases)
72   
73    :param G2frame: GSASII main frame structure
74    :param str histo: histogram name
75    :returns sum(scale*mass) for phases in histo
76    '''
77    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
78    wtSum = 0.0
79    for phase in Phases:
80        if Phases[phase]['General']['Type'] != 'magnetic':
81            if histo in Phases[phase]['Histograms']:
82                mass = Phases[phase]['General']['Mass']
83                phFr = Phases[phase]['Histograms'][histo]['Scale'][0]
84                wtSum += mass*phFr
85    return wtSum
86   
87################################################################################
88#GSASII pdf calculation routines
89################################################################################
90       
91def Transmission(Geometry,Abs,Diam):
92    '''
93    Calculate sample transmission
94
95    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
96    :param float Abs: absorption coeff in cm-1
97    :param float Diam: sample thickness/diameter in mm
98    '''
99    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
100        MuR = Abs*Diam/20.0
101        if MuR <= 3.0:
102            T0 = 16/(3.*math.pi)
103            T1 = -0.045780
104            T2 = -0.02489
105            T3 = 0.003045
106            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
107            if T < -20.:
108                return 2.06e-9
109            else:
110                return math.exp(T)
111        else:
112            T1 = 1.433902
113            T2 = 0.013869+0.337894
114            T3 = 1.933433+1.163198
115            T4 = 0.044365-0.04259
116            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
117            return T/100.
118    elif 'plate' in Geometry:
119        MuR = Abs*Diam/10.
120        return math.exp(-MuR)
121    elif 'Bragg' in Geometry:
122        return 0.0
123       
124def SurfaceRough(SRA,SRB,Tth):
125    ''' Suortti (J. Appl. Cryst, 5,325-331, 1972) surface roughness correction
126    :param float SRA: Suortti surface roughness parameter
127    :param float SRB: Suortti surface roughness parameter
128    :param float Tth: 2-theta(deg) - can be numpy array
129   
130    '''
131    sth = npsind(Tth/2.)
132    T1 = np.exp(-SRB/sth)
133    T2 = SRA+(1.-SRA)*np.exp(-SRB)
134    return (SRA+(1.-SRA)*T1)/T2
135   
136def SurfaceRoughDerv(SRA,SRB,Tth):
137    ''' Suortti surface roughness correction derivatives
138    :param float SRA: Suortti surface roughness parameter (dimensionless)
139    :param float SRB: Suortti surface roughness parameter (dimensionless)
140    :param float Tth: 2-theta(deg) - can be numpy array
141    :return list: [dydSRA,dydSRB] derivatives to be used for intensity derivative
142    '''
143    sth = npsind(Tth/2.)
144    T1 = np.exp(-SRB/sth)
145    T2 = SRA+(1.-SRA)*np.exp(-SRB)
146    Trans = (SRA+(1.-SRA)*T1)/T2
147    dydSRA = ((1.-T1)*T2-(1.-np.exp(-SRB))*Trans)/T2**2
148    dydSRB = ((SRA-1.)*T1*T2/sth-Trans*(SRA-T2))/T2**2
149    return [dydSRA,dydSRB]
150
151def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
152    '''Calculate sample absorption
153    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
154    :param float MuR: absorption coeff * sample thickness/2 or radius
155    :param Tth: 2-theta scattering angle - can be numpy array
156    :param float Phi: flat plate tilt angle - future
157    :param float Psi: flat plate tilt axis - future
158    '''
159   
160    def muRunder3(MuR,Sth2):
161        T0 = 16.0/(3.*np.pi)
162        T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
163            0.109561*np.sqrt(Sth2)-26.04556
164        T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
165            1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
166        T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
167        Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
168        return np.exp(Trns)
169   
170    def muRover3(MuR,Sth2):
171        T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
172            10.02088*Sth2**3-3.36778*Sth2**4
173        T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
174            (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
175        T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
176            0.13576*np.sqrt(Sth2)+1.163198
177        T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
178        Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
179        return Trns/100.
180       
181    Sth2 = npsind(Tth/2.0)**2
182    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
183        if 'array' in str(type(MuR)):
184            MuRSTh2 = np.concatenate((MuR,Sth2))
185            AbsCr = np.where(MuRSTh2[0]<=3.0,muRunder3(MuRSTh2[0],MuRSTh2[1]),muRover3(MuRSTh2[0],MuRSTh2[1]))
186            return AbsCr
187        else:
188            if MuR <= 3.0:
189                return muRunder3(MuR,Sth2)
190            else:
191                return muRover3(MuR,Sth2)
192    elif 'Bragg' in Geometry:
193        return 1.0
194    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
195        # and only defined for 2theta < 90
196        MuT = 2.*MuR
197        T1 = np.exp(-MuT)
198        T2 = np.exp(-MuT/npcosd(Tth))
199        Tb = MuT-MuT/npcosd(Tth)
200        return (T2-T1)/Tb
201    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
202        MuT = 2.*MuR
203        cth = npcosd(Tth/2.0)
204        return np.exp(-MuT/cth)/cth
205       
206def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
207    'needs a doc string'
208    dA = 0.001
209    AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
210    if MuR:
211        AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
212        return (AbsP-AbsM)/(2.0*dA)
213    else:
214        return (AbsP-1.)/dA
215       
216def Polarization(Pola,Tth,Azm=0.0):
217    """   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
218
219    :param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
220    :param Azm: azimuthal angle e.g. 0.0 in plane of polarization
221    :param Tth: 2-theta scattering angle - can be numpy array
222      which (if either) of these is "right"?
223    :return: (pola, dpdPola)
224      * pola = ((1-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \
225        (1-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
226      * dpdPola: derivative needed for least squares
227
228    """
229    pola = ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
230        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
231    dpdPola = -npsind(Tth)**2*(npsind(Azm)**2-npcosd(Azm)**2)
232    return pola,dpdPola
233   
234def Oblique(ObCoeff,Tth):
235    'currently assumes detector is normal to beam'
236    if ObCoeff:
237        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
238    else:
239        return 1.0
240               
241def Ruland(RulCoff,wave,Q,Compton):
242    'needs a doc string'
243    C = 2.9978e8
244    D = 1.5e-3
245    hmc = 0.024262734687
246    sinth2 = (Q*wave/(4.0*np.pi))**2
247    dlam = (wave**2)*Compton*Q/C
248    dlam_c = 2.0*hmc*sinth2-D*wave**2
249    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
250   
251def LorchWeight(Q):
252    'needs a doc string'
253    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
254           
255def GetAsfMean(ElList,Sthl2):
256    '''Calculate various scattering factor terms for PDF calcs
257
258    :param dict ElList: element dictionary contains scattering factor coefficients, etc.
259    :param np.array Sthl2: numpy array of sin theta/lambda squared values
260    :returns: mean(f^2), mean(f)^2, mean(compton)
261    '''
262    sumNoAtoms = 0.0
263    FF = np.zeros_like(Sthl2)
264    FF2 = np.zeros_like(Sthl2)
265    CF = np.zeros_like(Sthl2)
266    for El in ElList:
267        sumNoAtoms += ElList[El]['FormulaNo']
268    for El in ElList:
269        el = ElList[El]
270        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
271        cf = G2elem.ComptonFac(el,Sthl2)
272        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
273        FF2 += ff2*el['FormulaNo']/sumNoAtoms
274        CF += cf*el['FormulaNo']/sumNoAtoms
275    return FF2,FF**2,CF
276   
277def GetNumDensity(ElList,Vol):
278    'needs a doc string'
279    sumNoAtoms = 0.0
280    for El in ElList:
281        sumNoAtoms += ElList[El]['FormulaNo']
282    return sumNoAtoms/Vol
283           
284def CalcPDF(data,inst,limits,xydata):
285    '''Computes I(Q), S(Q) & G(r) from Sample, Bkg, etc. diffraction patterns loaded into
286    dict xydata; results are placed in xydata.
287    Calculation parameters are found in dicts data and inst and list limits.
288    The return value is at present an empty list.
289    '''
290    auxPlot = []
291    Ibeg = np.searchsorted(xydata['Sample'][1][0],limits[0])
292    Ifin = np.searchsorted(xydata['Sample'][1][0],limits[1])+1
293    #subtract backgrounds - if any & use PWDR limits
294#    GSASIIpath.IPyBreak()
295    IofQ = copy.deepcopy(xydata['Sample'])
296    IofQ[1] = np.array(IofQ[1])[:,Ibeg:Ifin]
297    if data['Sample Bkg.']['Name']:
298        IofQ[1][1] += xydata['Sample Bkg.'][1][1][Ibeg:Ifin]*data['Sample Bkg.']['Mult']
299    if data['Container']['Name']:
300        xycontainer = xydata['Container'][1][1]*data['Container']['Mult']
301        if data['Container Bkg.']['Name']:
302            xycontainer += xydata['Container Bkg.'][1][1][Ibeg:Ifin]*data['Container Bkg.']['Mult']
303        IofQ[1][1] += xycontainer[Ibeg:Ifin]
304    data['IofQmin'] = IofQ[1][1][-1]
305    IofQ[1][1] -= data.get('Flat Bkg',0.)
306    #get element data & absorption coeff.
307    ElList = data['ElList']
308    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
309    #Apply angle dependent corrections
310    Tth = IofQ[1][0]
311    MuR = Abs*data['Diam']/20.0
312    IofQ[1][1] /= Absorb(data['Geometry'],MuR,Tth)
313    if 'X' in inst['Type'][0]:
314        IofQ[1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
315    if data['DetType'] == 'Image plate':
316        IofQ[1][1] *= Oblique(data['ObliqCoeff'],Tth)
317    XY = IofQ[1]   
318    #convert to Q
319#    nQpoints = len(XY[0])     #points for Q interpolation
320    nQpoints = 5000
321    if 'C' in inst['Type'][0]:
322        wave = G2mth.getWave(inst)
323        minQ = npT2q(Tth[0],wave)
324        maxQ = npT2q(Tth[-1],wave)   
325        Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
326        dq = Qpoints[1]-Qpoints[0]
327        XY[0] = npT2q(XY[0],wave)
328    elif 'T' in inst['Type'][0]:
329        difC = inst['difC'][1]
330        minQ = 2.*np.pi*difC/Tth[-1]
331        maxQ = 2.*np.pi*difC/Tth[0]
332        Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
333        dq = Qpoints[1]-Qpoints[0]
334        XY[0] = 2.*np.pi*difC/XY[0]
335    Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0])    #interpolate I(Q)
336    Qdata -= np.min(Qdata)*data['BackRatio']
337   
338    qLimits = data['QScaleLim']
339    minQ = np.searchsorted(Qpoints,qLimits[0])
340    maxQ = np.searchsorted(Qpoints,qLimits[1])+1
341    newdata = []
342    if len(IofQ) < 3:
343        xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],'']
344    else:
345        xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],IofQ[2]]
346    for item in xydata['IofQ'][1]:
347        newdata.append(item[:maxQ])
348    xydata['IofQ'][1] = newdata
349   
350    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
351    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
352    Q = xydata['SofQ'][1][0]
353#    auxPlot.append([Q,np.copy(CF),'CF-unCorr'])
354    ruland = Ruland(data['Ruland'],wave,Q,CF)
355#    auxPlot.append([Q,ruland,'Ruland'])     
356    CF *= ruland
357#    auxPlot.append([Q,CF,'CF-Corr'])
358    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
359    xydata['SofQ'][1][1] *= scale
360    xydata['SofQ'][1][1] -= CF
361    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
362    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
363    xydata['SofQ'][1][1] *= scale
364    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
365    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
366    if data['Lorch']:
367        xydata['FofQ'][1][1] *= LorchWeight(Q)   
368    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
369    nR = len(xydata['GofR'][1][1])
370    mul = int(round(2.*np.pi*nR/(data.get('Rmax',100.)*qLimits[1])))
371    xydata['GofR'][1][0] = 2.*np.pi*np.linspace(0,nR,nR,endpoint=True)/(mul*qLimits[1])
372    xydata['GofR'][1][1] = -dq*np.imag(fft.fft(xydata['FofQ'][1][1],mul*nR)[:nR])
373    if data.get('noRing',True):
374        xydata['GofR'][1][1] = np.where(xydata['GofR'][1][0]<0.5,0.,xydata['GofR'][1][1])
375    return auxPlot
376   
377def PDFPeakFit(peaks,data):
378    rs2pi = 1./np.sqrt(2*np.pi)
379   
380    def MakeParms(peaks):
381        varyList = []
382        parmDict = {'slope':peaks['Background'][1][1]}
383        if peaks['Background'][2]:
384            varyList.append('slope')
385        for i,peak in enumerate(peaks['Peaks']):
386            parmDict['PDFpos;'+str(i)] = peak[0]
387            parmDict['PDFmag;'+str(i)] = peak[1]
388            parmDict['PDFsig;'+str(i)] = peak[2]
389            if 'P' in peak[3]:
390                varyList.append('PDFpos;'+str(i))
391            if 'M' in peak[3]:
392                varyList.append('PDFmag;'+str(i))
393            if 'S' in peak[3]:
394                varyList.append('PDFsig;'+str(i))
395        return parmDict,varyList
396       
397    def SetParms(peaks,parmDict,varyList):
398        if 'slope' in varyList:
399            peaks['Background'][1][1] = parmDict['slope']
400        for i,peak in enumerate(peaks['Peaks']):
401            if 'PDFpos;'+str(i) in varyList:
402                peak[0] = parmDict['PDFpos;'+str(i)]
403            if 'PDFmag;'+str(i) in varyList:
404                peak[1] = parmDict['PDFmag;'+str(i)]
405            if 'PDFsig;'+str(i) in varyList:
406                peak[2] = parmDict['PDFsig;'+str(i)]
407       
408   
409    def CalcPDFpeaks(parmdict,Xdata):
410        Z = parmDict['slope']*Xdata
411        ipeak = 0
412        while True:
413            try:
414                pos = parmdict['PDFpos;'+str(ipeak)]
415                mag = parmdict['PDFmag;'+str(ipeak)]
416                wid = parmdict['PDFsig;'+str(ipeak)]
417                wid2 = 2.*wid**2
418                Z += mag*rs2pi*np.exp(-(Xdata-pos)**2/wid2)/wid
419                ipeak += 1
420            except KeyError:        #no more peaks to process
421                return Z
422               
423    def errPDFProfile(values,xdata,ydata,parmdict,varylist):       
424        parmdict.update(zip(varylist,values))
425        M = CalcPDFpeaks(parmdict,xdata)-ydata
426        return M
427           
428    newpeaks = copy.copy(peaks)
429    iBeg = np.searchsorted(data[1][0],newpeaks['Limits'][0])
430    iFin = np.searchsorted(data[1][0],newpeaks['Limits'][1])+1
431    X = data[1][0][iBeg:iFin]
432    Y = data[1][1][iBeg:iFin]
433    parmDict,varyList = MakeParms(peaks)
434    if not len(varyList):
435        print ' Nothing varied'
436        return newpeaks,None,None,None,None,None
437   
438    Rvals = {}
439    values =  np.array(Dict2Values(parmDict, varyList))
440    result = so.leastsq(errPDFProfile,values,full_output=True,ftol=0.0001,
441           args=(X,Y,parmDict,varyList))
442    chisq = np.sum(result[2]['fvec']**2)
443    Values2Dict(parmDict, varyList, result[0])
444    SetParms(peaks,parmDict,varyList)
445    Rvals['Rwp'] = np.sqrt(chisq/np.sum(Y**2))*100.      #to %
446    chisq = np.sum(result[2]['fvec']**2)/(len(X)-len(values))   #reduced chi^2 = M/(Nobs-Nvar)
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.signal as signal
454    auxPlot = []
455    if 'C' in inst['Type'][0]:
456        Tth = pwddata[0]
457        wave = G2mth.getWave(inst)
458        minQ = npT2q(Tth[0],wave)
459        maxQ = npT2q(Tth[-1],wave)
460        powQ = npT2q(Tth,wave) 
461    elif 'T' in inst['Type'][0]:
462        TOF = pwddata[0]
463        difC = inst['difC'][1]
464        minQ = 2.*np.pi*difC/TOF[-1]
465        maxQ = 2.*np.pi*difC/TOF[0]
466        powQ = 2.*np.pi*difC/TOF
467    piDQ = np.pi/(maxQ-minQ)
468    Qpoints = np.linspace(minQ,maxQ,len(pwddata[0]),endpoint=True)
469    if RDFcontrols['UseObsCalc']:
470        Qdata = si.griddata(powQ,pwddata[1]-pwddata[3],Qpoints,method=RDFcontrols['Smooth'],fill_value=0.)
471    else:
472        Qdata = si.griddata(powQ,pwddata[1]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
473    Qdata *= np.sin((Qpoints-minQ)*piDQ)/piDQ
474    Qdata *= 0.5*np.sqrt(Qpoints)       #Qbin normalization
475#    GSASIIpath.IPyBreak()
476    dq = Qpoints[1]-Qpoints[0]
477    nR = len(Qdata)
478    R = 0.5*np.pi*np.linspace(0,nR,nR)/(4.*maxQ)
479    iFin = np.searchsorted(R,RDFcontrols['maxR'])+1
480    bBut,aBut = signal.butter(4,0.01)
481    Qsmooth = signal.filtfilt(bBut,aBut,Qdata)
482#    auxPlot.append([Qpoints,Qdata,'interpolate:'+RDFcontrols['Smooth']])
483#    auxPlot.append([Qpoints,Qsmooth,'interpolate:'+RDFcontrols['Smooth']])
484    DofR = dq*np.imag(fft.fft(Qsmooth,16*nR)[:nR])
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.fft(z)
694        Df = fft.ifft(Z.prod(axis=0)).real
695    else:
696        z = np.column_stack([Norm,Cauchy]).T
697        Z = fft.fft(z)
698        Df = fft.fftshift(fft.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])+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# Reflectometry calculations
1913################################################################################
1914
1915def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
1916    print 'fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0])
1917   
1918    class RandomDisplacementBounds(object):
1919        """random displacement with bounds"""
1920        def __init__(self, xmin, xmax, stepsize=0.5):
1921            self.xmin = xmin
1922            self.xmax = xmax
1923            self.stepsize = stepsize
1924   
1925        def __call__(self, x):
1926            """take a random step but ensure the new position is within the bounds"""
1927            while True:
1928                # this could be done in a much more clever way, but it will work for example purposes
1929                steps = self.xmax-self.xmin
1930                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
1931                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
1932                    break
1933            return xnew
1934   
1935    def GetModelParms():
1936        parmDict = {}
1937        varyList = []
1938        values = []
1939        bounds = []
1940        parmDict['Res'] = data['Resolution'][0]/(100.*np.sqrt(ateln2))     #% FWHM-->decimal sig
1941        for parm in ['Scale','FltBack']:
1942            parmDict[parm] = data[parm][0]
1943            if data[parm][1]:
1944                varyList.append(parm)
1945                values.append(data[parm][0])
1946                bounds.append(Bounds[parm])
1947        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
1948        parmDict['nLayers'] = len(parmDict['Layer Seq'])
1949        for ilay,layer in enumerate(data['Layers']):
1950            name = layer['Name']
1951            cid = str(ilay)+';'
1952            parmDict[cid+'Name'] = name
1953            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
1954                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
1955                if layer.get(parm,[0.,False])[1]:
1956                    varyList.append(cid+parm)
1957                    value = layer[parm][0]
1958                    values.append(value)
1959                    if value:
1960                        bound = [value*Bfac,value/Bfac]
1961                    else:
1962                        bound = [0.,10.]
1963                    bounds.append(bound)
1964            if name not in ['vacuum','unit scatter']:
1965                parmDict[cid+'rho'] = Substances[name]['Scatt density']
1966                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
1967        return parmDict,varyList,values,bounds
1968   
1969    def SetModelParms():
1970        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
1971        if 'Scale' in varyList:
1972            data['Scale'][0] = parmDict['Scale']
1973            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
1974        print line
1975        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
1976        if 'FltBack' in varyList:
1977            data['FltBack'][0] = parmDict['FltBack']
1978            line += ' esd: %15.3g'%(sigDict['FltBack'])
1979        print line
1980        for ilay,layer in enumerate(data['Layers']):
1981            name = layer['Name']
1982            print ' Parameters for layer: %d %s'%(ilay,name)
1983            cid = str(ilay)+';'
1984            line = ' '
1985            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
1986            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)**parmDict[cid+'DenMul'])
1987            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
1988                if parm in layer:
1989                    layer[parm][0] = parmDict[cid+parm]
1990                    line += ' %s: %.3f'%(parm,layer[parm][0])
1991                    if cid+parm in varyList:
1992                        line += ' esd: %.3g'%(sigDict[cid+parm])
1993            print line
1994            print line2
1995   
1996    def calcREFD(values,Q,Io,wt,parmDict,varyList):
1997        parmDict.update(zip(varyList,values))
1998        M = np.sqrt(wt)*(getREFD(Q,parmDict)-Io)
1999        return M
2000   
2001    def sumREFD(values,Q,Io,wt,parmDict,varyList):
2002        parmDict.update(zip(varyList,values))
2003        M = np.sqrt(wt)*(getREFD(Q,parmDict)-Io)
2004        return np.sum(M**2)
2005   
2006    def getREFD(Q,parmDict):
2007        Ic = np.ones_like(Q)*parmDict['FltBack']
2008        Scale = parmDict['Scale']
2009        Nlayers = parmDict['nLayers']
2010        Res = parmDict['Res']
2011        depth = np.zeros(Nlayers)
2012        rho = np.zeros(Nlayers)
2013        irho = np.zeros(Nlayers)
2014        sigma = np.zeros(Nlayers)
2015        for ilay,lay in enumerate(parmDict['Layer Seq']):
2016            cid = str(lay)+';'
2017            depth[ilay] = parmDict[cid+'Thick']
2018            sigma[ilay] = parmDict[cid+'Rough']
2019            if parmDict[cid+'Name'] == u'unit scatter':
2020                rho[ilay] = parmDict[cid+'DenMul']
2021                irho[ilay] = parmDict[cid+'iDenMul']
2022            elif 'vacuum' != parmDict[cid+'Name']:
2023                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
2024                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
2025            if cid+'Mag SLD' in parmDict:
2026                rho[ilay] += parmDict[cid+'Mag SLD']
2027        A,B = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2028        Ic += (A**2+B**2)*Scale
2029        #TODO: smear Ic by instrument resolution Qsig
2030        return Ic
2031       
2032        def Smear(f,w,z,dq):
2033            y = f(w,z)
2034            s = dq/ateln2
2035            y += 0.1354*(f(w,z+2*s)+f(w,z-2*s))
2036            y += 0.24935*(f(w,z-1.666667*s)+f(w,z+1.666667*s)) 
2037            y += 0.4111*(f(w,z-1.333333*s)+f(w,z+1.333333*s)) 
2038            y += 0.60653*(f(w,z-s) +f(w,z+s))
2039            y += 0.80074*(f(w,z-0.6666667*s)+f(w,z+0.6666667*s))
2040            y += 0.94596*(f(w,z-0.3333333*s)+f(w,z+0.3333333*s))
2041            y *= 0.137023
2042            return y
2043       
2044    def estimateT0(takestep):
2045        Mmax = -1.e-10
2046        Mmin = 1.e10
2047        for i in range(100):
2048            x0 = takestep(values)
2049            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)
2050            Mmin = min(M,Mmin)
2051            MMax = max(M,Mmax)
2052        return 1.5*(MMax-Mmin)
2053
2054    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2055    if data.get('2% weight'):
2056        wt = 1./(0.02*Io)**2
2057    Qmin = Limits[1][0]
2058    Qmax = Limits[1][1]
2059    wtFactor = ProfDict['wtFactor']
2060    Bfac = data['Toler']
2061    Ibeg = np.searchsorted(Q,Qmin)
2062    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
2063    Ic[:] = 0
2064    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
2065              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
2066    parmDict,varyList,values,bounds = GetModelParms()
2067    Msg = 'Failed to converge'
2068    if varyList:
2069        if data['Minimizer'] == 'LMLS': 
2070            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
2071                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList))
2072            parmDict.update(zip(varyList,result[0]))
2073            chisq = np.sum(result[2]['fvec']**2)
2074            ncalc = result[2]['nfev']
2075            covM = result[1]
2076            newVals = result[0]
2077        elif data['Minimizer'] == 'Basin Hopping':
2078            xyrng = np.array(bounds).T
2079            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
2080            T0 = estimateT0(take_step)
2081            print ' Estimated temperature: %.3g'%(T0)
2082            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
2083                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
2084                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)})
2085            chisq = result.fun
2086            ncalc = result.nfev
2087            newVals = result.x
2088            covM = []
2089        elif data['Minimizer'] == 'MC/SA Anneal':
2090            xyrng = np.array(bounds).T
2091            result = G2mth.anneal(sumREFD, values, args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList),
2092                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
2093                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
2094                ranRange=0.20,autoRan=False,dlg=None)
2095            newVals = result[0]
2096            parmDict.update(zip(varyList,newVals))
2097            chisq = result[1]
2098            ncalc = result[3]
2099            covM = []
2100            print ' MC/SA final temperature: %.4g'%(result[2])
2101        elif data['Minimizer'] == 'L-BFGS-B':
2102            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
2103                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList))
2104            parmDict.update(zip(varyList,result['x']))
2105            chisq = result.fun
2106            ncalc = result.nfev
2107            newVals = result.x
2108            covM = []
2109    else:   #nothing varied
2110        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)
2111        chisq = np.sum(M**2)
2112        ncalc = 0
2113        covM = []
2114        sig = []
2115        sigDict = {}
2116        result = []
2117    Rvals = {}
2118    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
2119    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
2120    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],parmDict)
2121    Ib[Ibeg:Ifin] = parmDict['FltBack']
2122    try:
2123        if not len(varyList):
2124            Msg += ' - nothing refined'
2125            raise ValueError
2126        Nans = np.isnan(newVals)
2127        if np.any(Nans):
2128            name = varyList[Nans.nonzero(True)[0]]
2129            Msg += ' Nan result for '+name+'!'
2130            raise ValueError
2131        Negs = np.less_equal(newVals,0.)
2132        if np.any(Negs):
2133            indx = Negs.nonzero()
2134            name = varyList[indx[0][0]]
2135            if name != 'FltBack' and 'Mag SLD' not in name:
2136                Msg += ' negative coefficient for '+name+'!'
2137                raise ValueError
2138        if len(covM):
2139            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
2140            covMatrix = covM*Rvals['GOF']
2141        else:
2142            sig = np.zeros(len(varyList))
2143            covMatrix = []
2144        sigDict = dict(zip(varyList,sig))
2145        print ' Results of reflectometry data modelling fit:'
2146        print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
2147        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
2148        SetModelParms()
2149        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
2150    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
2151        print Msg
2152        return False,0,0,0,0,0,0,Msg
2153       
2154def makeSLDprofile(data,Substances):
2155   
2156    sq2 = np.sqrt(2.)
2157    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2158    Nlayers = len(laySeq)
2159    laySeq = np.array(laySeq,dtype=int)
2160    interfaces = np.zeros(Nlayers)
2161    rho = np.zeros(Nlayers)
2162    sigma = np.zeros(Nlayers)
2163    name = data['Layers'][0]['Name']
2164    thick = 0.
2165    for ilay,lay in enumerate(laySeq):
2166        layer = data['Layers'][lay]
2167        name = layer['Name']
2168        if 'Thick' in layer:
2169            thick += layer['Thick'][0]
2170            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
2171        if 'Rough' in layer:
2172            sigma[ilay] = max(0.001,layer['Rough'][0])
2173        if name != 'vacuum':
2174            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2175        if 'Mag SLD' in layer:
2176            rho[ilay] += layer['Mag SLD'][0]
2177    name = data['Layers'][-1]['Name']
2178    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
2179    xr = np.flipud(x)
2180    interfaces[-1] = x[-1]
2181    y = np.ones_like(x)*rho[0]
2182    iBeg = 0
2183    for ilayer in range(Nlayers-1):
2184        delt = rho[ilayer+1]-rho[ilayer]
2185        iPos = np.searchsorted(x,interfaces[ilayer])
2186        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
2187        iBeg = iPos
2188    return x,xr,y   
2189
2190def REFDModelFxn(Profile,Inst,Limits,Substances,data):
2191   
2192    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2193    Qmin = Limits[1][0]
2194    Qmax = Limits[1][1]
2195    iBeg = np.searchsorted(Q,Qmin)
2196    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2197    Ib[:] = data['FltBack'][0]
2198    Ic[:] = 0
2199    Scale = data['Scale'][0]
2200    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2201    Nlayers = len(laySeq)
2202    depth = np.zeros(Nlayers)
2203    rho = np.zeros(Nlayers)
2204    irho = np.zeros(Nlayers)
2205    sigma = np.zeros(Nlayers)
2206    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
2207        layer = data['Layers'][lay]
2208        name = layer['Name']
2209        if 'Thick' in layer:    #skips first & last layers
2210            depth[ilay] = layer['Thick'][0]
2211        if 'Rough' in layer:    #skips first layer
2212            sigma[ilay] = layer['Rough'][0]
2213        if 'unit scatter' == name:
2214            rho[ilay] = layer['DenMul'][0]
2215            irho[ilay] = layer['iDenMul'][0]
2216        else:
2217            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2218            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
2219        if 'Mag SLD' in layer:
2220            rho[ilay] += layer['Mag SLD'][0]
2221    A,B = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2222    Ic[iBeg:iFin] = (A**2+B**2)*Scale+Ib[iBeg:iFin]
2223
2224def abeles(kz, depth, rho, irho=0, sigma=0):
2225    """
2226    Optical matrix form of the reflectivity calculation.
2227    O.S. Heavens, Optical Properties of Thin Solid Films
2228   
2229    Reflectometry as a function of kz for a set of slabs.
2230
2231    :Parameters:
2232
2233    *kz* : float[n] | |1/Ang|
2234        Scattering vector $2\pi\sin(\theta)/\lambda$. This is $\tfrac12 Q_z$.
2235    *depth* :  float[m] | |Ang|
2236        thickness of each layer.  The thickness of the incident medium
2237        and substrate are ignored.
2238    *rho*, *irho* :  float[n,k] | |1e-6/Ang^2|
2239        real and imaginary scattering length density for each layer for each kz
2240        Note: absorption cross section mu = 2 irho/lambda for neutrons
2241    *sigma* : float[m-1] | |Ang|
2242        interfacial roughness.  This is the roughness between a layer
2243        and the previous layer. The sigma array should have m-1 entries.
2244
2245    Slabs are ordered with the surface SLD at index 0 and substrate at
2246    index -1, or reversed if kz < 0.
2247    """
2248    def calc(kz, depth, rho, irho, sigma):
2249        if len(kz) == 0: return kz
2250   
2251        # Complex index of refraction is relative to the incident medium.
2252        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
2253        # in place of kz^2, and ignoring rho_o
2254        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
2255        k = kz
2256   
2257        # According to Heavens, the initial matrix should be [ 1 F; F 1],
2258        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
2259        # multiply versus some coding convenience.
2260        B11 = 1
2261        B22 = 1
2262        B21 = 0
2263        B12 = 0
2264        for i in range(0, len(depth)-1):
2265            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
2266            F = (k - k_next) / (k + k_next)
2267            F *= np.exp(-2*k*k_next*sigma[i]**2)
2268            #print "==== layer",i
2269            #print "kz:", kz
2270            #print "k:", k
2271            #print "k_next:",k_next
2272            #print "F:",F
2273            #print "rho:",rho[:,i+1]
2274            #print "irho:",irho[:,i+1]
2275            #print "d:",depth[i],"sigma:",sigma[i]
2276            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
2277            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
2278            M21 = F*M11
2279            M12 = F*M22
2280            C1 = B11*M11 + B21*M12
2281            C2 = B11*M21 + B21*M22
2282            B11 = C1
2283            B21 = C2
2284            C1 = B12*M11 + B22*M12
2285            C2 = B12*M21 + B22*M22
2286            B12 = C1
2287            B22 = C2
2288            k = k_next
2289   
2290        r = B12/B11
2291        return np.real(r),np.imag(r)
2292
2293    if np.isscalar(kz): kz = np.asarray([kz], 'd')
2294
2295    m = len(depth)
2296
2297    # Make everything into arrays
2298    depth = np.asarray(depth,'d')
2299    rho = np.asarray(rho,'d')
2300    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
2301    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
2302
2303    # Repeat rho,irho columns as needed
2304    if len(rho.shape) == 1:
2305        rho = rho[None,:]
2306        irho = irho[None,:]
2307
2308    return calc(kz, depth, rho, irho, sigma)
2309   
2310def makeRefdFFT(Limits,Profile):
2311    print 'make fft'
2312    Q,Io = Profile[:2]
2313    Qmin = Limits[1][0]
2314    Qmax = Limits[1][1]
2315    iBeg = np.searchsorted(Q,Qmin)
2316    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2317    Qf = np.linspace(0.,Q[iFin-1],5000)
2318    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
2319    If = QI(Qf)*Qf**4
2320    R = np.linspace(0.,5000.,5000)
2321    F = fft.rfft(If)
2322    return R,F
2323
2324   
2325################################################################################
2326# Stacking fault simulation codes
2327################################################################################
2328
2329def GetStackParms(Layers):
2330   
2331    Parms = []
2332#cell parms
2333    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
2334        Parms.append('cellA')
2335        Parms.append('cellC')
2336    else:
2337        Parms.append('cellA')
2338        Parms.append('cellB')
2339        Parms.append('cellC')
2340        if Layers['Laue'] != 'mmm':
2341            Parms.append('cellG')
2342#Transition parms
2343    for iY in range(len(Layers['Layers'])):
2344        for iX in range(len(Layers['Layers'])):
2345            Parms.append('TransP;%d;%d'%(iY,iX))
2346            Parms.append('TransX;%d;%d'%(iY,iX))
2347            Parms.append('TransY;%d;%d'%(iY,iX))
2348            Parms.append('TransZ;%d;%d'%(iY,iX))
2349    return Parms
2350
2351def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
2352    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
2353   
2354    param: Layers dict: 'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
2355                        'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
2356                        'Layers':[],'Stacking':[],'Transitions':[]}
2357    param: ctrls string: controls string to be written on DIFFaX controls.dif file
2358    param: scale float: scale factor
2359    param: background dict: background parameters
2360    param: limits list: min/max 2-theta to be calculated
2361    param: inst dict: instrument parameters dictionary
2362    param: profile list: powder pattern data
2363   
2364    all updated in place   
2365    '''
2366    import atmdata
2367    path = sys.path
2368    for name in path:
2369        if 'bin' in name:
2370            DIFFaX = name+'/DIFFaX.exe'
2371            print ' Execute ',DIFFaX
2372            break
2373    # make form factor file that DIFFaX wants - atom types are GSASII style
2374    sf = open('data.sfc','w')
2375    sf.write('GSASII special form factor file for DIFFaX\n\n')
2376    atTypes = Layers['AtInfo'].keys()
2377    if 'H' not in atTypes:
2378        atTypes.insert(0,'H')
2379    for atType in atTypes:
2380        if atType == 'H': 
2381            blen = -.3741
2382        else:
2383            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
2384        Adat = atmdata.XrayFF[atType]
2385        text = '%4s'%(atType.ljust(4))
2386        for i in range(4):
2387            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
2388        text += '%11.6f%11.6f'%(Adat['fc'],blen)
2389        text += '%3d\n'%(Adat['Z'])
2390        sf.write(text)
2391    sf.close()
2392    #make DIFFaX control.dif file - future use GUI to set some of these flags
2393    cf = open('control.dif','w')
2394    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2395        x0 = profile[0]
2396        iBeg = np.searchsorted(x0,limits[0])
2397        iFin = np.searchsorted(x0,limits[1])+1
2398        if iFin-iBeg > 20000:
2399            iFin = iBeg+20000
2400        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2401        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2402        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
2403    else:
2404        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2405        inst = {'Type':['XSC','XSC',]}
2406    cf.close()
2407    #make DIFFaX data file
2408    df = open('GSASII-DIFFaX.dat','w')
2409    df.write('INSTRUMENTAL\n')
2410    if 'X' in inst['Type'][0]:
2411        df.write('X-RAY\n')
2412    elif 'N' in inst['Type'][0]:
2413        df.write('NEUTRON\n')
2414    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2415        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
2416        U = ateln2*inst['U'][1]/10000.
2417        V = ateln2*inst['V'][1]/10000.
2418        W = ateln2*inst['W'][1]/10000.
2419        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2420        HW = np.sqrt(np.mean(HWHM))
2421    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
2422        if 'Mean' in Layers['selInst']:
2423            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
2424        elif 'Gaussian' in Layers['selInst']:
2425            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
2426        else:
2427            df.write('None\n')
2428    else:
2429        df.write('0.10\nNone\n')
2430    df.write('STRUCTURAL\n')
2431    a,b,c = Layers['Cell'][1:4]
2432    gam = Layers['Cell'][6]
2433    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
2434    laue = Layers['Laue']
2435    if laue == '2/m(ab)':
2436        laue = '2/m(1)'
2437    elif laue == '2/m(c)':
2438        laue = '2/m(2)'
2439    if 'unknown' in Layers['Laue']:
2440        df.write('%s %.3f\n'%(laue,Layers['Toler']))
2441    else:   
2442        df.write('%s\n'%(laue))
2443    df.write('%d\n'%(len(Layers['Layers'])))
2444    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
2445        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
2446    layerNames = []
2447    for layer in Layers['Layers']:
2448        layerNames.append(layer['Name'])
2449    for il,layer in enumerate(Layers['Layers']):
2450        if layer['SameAs']:
2451            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
2452            continue
2453        df.write('LAYER %d\n'%(il+1))
2454        if '-1' in layer['Symm']:
2455            df.write('CENTROSYMMETRIC\n')
2456        else:
2457            df.write('NONE\n')
2458        for ia,atom in enumerate(layer['Atoms']):
2459            [name,atype,x,y,z,frac,Uiso] = atom
2460            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
2461                frac /= 2.
2462            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
2463    df.write('STACKING\n')
2464    df.write('%s\n'%(Layers['Stacking'][0]))
2465    if 'recursive' in Layers['Stacking'][0]:
2466        df.write('%s\n'%Layers['Stacking'][1])
2467    else:
2468        if 'list' in Layers['Stacking'][1]:
2469            Slen = len(Layers['Stacking'][2])
2470            iB = 0
2471            iF = 0
2472            while True:
2473                iF += 68
2474                if iF >= Slen:
2475                    break
2476                iF = min(iF,Slen)
2477                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
2478                iB = iF
2479        else:
2480            df.write('%s\n'%Layers['Stacking'][1])   
2481    df.write('TRANSITIONS\n')
2482    for iY in range(len(Layers['Layers'])):
2483        sumPx = 0.
2484        for iX in range(len(Layers['Layers'])):
2485            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
2486            p = round(p,3)
2487            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
2488            sumPx += p
2489        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
2490            print 'ERROR - Layer probabilities sum to ',sumPx,' DIFFaX will insist it = 1.0'
2491            df.close()
2492            os.remove('data.sfc')
2493            os.remove('control.dif')
2494            os.remove('GSASII-DIFFaX.dat')
2495            return       
2496    df.close()   
2497    time0 = time.time()
2498    try:
2499        subp.call(DIFFaX)
2500    except OSError:
2501        print ' DIFFax.exe is not available for this platform - under development'
2502    print ' DIFFaX time = %.2fs'%(time.time()-time0)
2503    if os.path.exists('GSASII-DIFFaX.spc'):
2504        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
2505        iFin = iBeg+Xpat.shape[1]
2506        bakType,backDict,backVary = SetBackgroundParms(background)
2507        backDict['Lam1'] = G2mth.getWave(inst)
2508        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2509        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
2510        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2511            rv = st.poisson(profile[3][iBeg:iFin])
2512            profile[1][iBeg:iFin] = rv.rvs()
2513            Z = np.ones_like(profile[3][iBeg:iFin])
2514            Z[1::2] *= -1
2515            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2516            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2517        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2518    #cleanup files..
2519        os.remove('GSASII-DIFFaX.spc')
2520    elif os.path.exists('GSASII-DIFFaX.sadp'):
2521        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
2522        Sadp = np.reshape(Sadp,(256,-1))
2523        Layers['Sadp']['Img'] = Sadp
2524        os.remove('GSASII-DIFFaX.sadp')
2525    os.remove('data.sfc')
2526    os.remove('control.dif')
2527    os.remove('GSASII-DIFFaX.dat')
2528   
2529def SetPWDRscan(inst,limits,profile):
2530   
2531    wave = G2mth.getMeanWave(inst)
2532    x0 = profile[0]
2533    iBeg = np.searchsorted(x0,limits[0])
2534    iFin = np.searchsorted(x0,limits[1])+1
2535    if iFin-iBeg > 20000:
2536        iFin = iBeg+20000
2537    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2538    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
2539    return iFin-iBeg
2540       
2541def SetStackingSF(Layers,debug):
2542# Load scattering factors into DIFFaX arrays
2543    import atmdata
2544    atTypes = Layers['AtInfo'].keys()
2545    aTypes = []
2546    for atype in atTypes:
2547        aTypes.append('%4s'%(atype.ljust(4)))
2548    SFdat = []
2549    for atType in atTypes:
2550        Adat = atmdata.XrayFF[atType]
2551        SF = np.zeros(9)
2552        SF[:8:2] = Adat['fa']
2553        SF[1:8:2] = Adat['fb']
2554        SF[8] = Adat['fc']
2555        SFdat.append(SF)
2556    SFdat = np.array(SFdat)
2557    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
2558   
2559def SetStackingClay(Layers,Type):
2560# Controls
2561    rand.seed()
2562    ranSeed = rand.randint(1,2**16-1)
2563    try:   
2564        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
2565            '6/m','6/mmm'].index(Layers['Laue'])+1
2566    except ValueError:  #for 'unknown'
2567        laueId = -1
2568    if 'SADP' in Type:
2569        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
2570        lmax = int(Layers['Sadp']['Lmax'])
2571    else:
2572        planeId = 0
2573        lmax = 0
2574# Sequences
2575    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
2576    try:
2577        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
2578    except ValueError:
2579        StkParm = -1
2580    if StkParm == 2:    #list
2581        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
2582        Nstk = len(StkSeq)
2583    else:
2584        Nstk = 1
2585        StkSeq = [0,]
2586    if StkParm == -1:
2587        StkParm = int(Layers['Stacking'][1])
2588    Wdth = Layers['Width'][0]
2589    mult = 1
2590    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
2591    LaueSym = Layers['Laue'].ljust(12)
2592    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
2593    return laueId,controls
2594   
2595def SetCellAtoms(Layers):
2596    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
2597# atoms in layers
2598    atTypes = Layers['AtInfo'].keys()
2599    AtomXOU = []
2600    AtomTp = []
2601    LayerSymm = []
2602    LayerNum = []
2603    layerNames = []
2604    Natm = 0
2605    Nuniq = 0
2606    for layer in Layers['Layers']:
2607        layerNames.append(layer['Name'])
2608    for il,layer in enumerate(Layers['Layers']):
2609        if layer['SameAs']:
2610            LayerNum.append(layerNames.index(layer['SameAs'])+1)
2611            continue
2612        else:
2613            LayerNum.append(il+1)
2614            Nuniq += 1
2615        if '-1' in layer['Symm']:
2616            LayerSymm.append(1)
2617        else:
2618            LayerSymm.append(0)
2619        for ia,atom in enumerate(layer['Atoms']):
2620            [name,atype,x,y,z,frac,Uiso] = atom
2621            Natm += 1
2622            AtomTp.append('%4s'%(atype.ljust(4)))
2623            Ta = atTypes.index(atype)+1
2624            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
2625    AtomXOU = np.array(AtomXOU)
2626    Nlayers = len(layerNames)
2627    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
2628    return Nlayers
2629   
2630def SetStackingTrans(Layers,Nlayers):
2631# Transitions
2632    TransX = []
2633    TransP = []
2634    for Ytrans in Layers['Transitions']:
2635        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
2636        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
2637    TransP = np.array(TransP,dtype='float').T
2638    TransX = np.array(TransX,dtype='float')
2639#    GSASIIpath.IPyBreak()
2640    pyx.pygettrans(Nlayers,TransP,TransX)
2641   
2642def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
2643# Scattering factors
2644    SetStackingSF(Layers,debug)
2645# Controls & sequences
2646    laueId,controls = SetStackingClay(Layers,'PWDR')
2647# cell & atoms
2648    Nlayers = SetCellAtoms(Layers)
2649    Volume = Layers['Cell'][7]   
2650# Transitions
2651    SetStackingTrans(Layers,Nlayers)
2652# PWDR scan
2653    Nsteps = SetPWDRscan(inst,limits,profile)
2654# result as Spec
2655    x0 = profile[0]
2656    profile[3] = np.zeros(len(profile[0]))
2657    profile[4] = np.zeros(len(profile[0]))
2658    profile[5] = np.zeros(len(profile[0]))
2659    iBeg = np.searchsorted(x0,limits[0])
2660    iFin = np.searchsorted(x0,limits[1])+1
2661    if iFin-iBeg > 20000:
2662        iFin = iBeg+20000
2663    Nspec = 20001       
2664    spec = np.zeros(Nspec,dtype='double')   
2665    time0 = time.time()
2666    pyx.pygetspc(controls,Nspec,spec)
2667    print ' GETSPC time = %.2fs'%(time.time()-time0)
2668    time0 = time.time()
2669    U = ateln2*inst['U'][1]/10000.
2670    V = ateln2*inst['V'][1]/10000.
2671    W = ateln2*inst['W'][1]/10000.
2672    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2673    HW = np.sqrt(np.mean(HWHM))
2674    BrdSpec = np.zeros(Nsteps)
2675    if 'Mean' in Layers['selInst']:
2676        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
2677    elif 'Gaussian' in Layers['selInst']:
2678        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
2679    else:
2680        BrdSpec = spec[:Nsteps]
2681    BrdSpec /= Volume
2682    iFin = iBeg+Nsteps
2683    bakType,backDict,backVary = SetBackgroundParms(background)
2684    backDict['Lam1'] = G2mth.getWave(inst)
2685    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2686    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
2687    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2688        rv = st.poisson(profile[3][iBeg:iFin])
2689        profile[1][iBeg:iFin] = rv.rvs()
2690        Z = np.ones_like(profile[3][iBeg:iFin])
2691        Z[1::2] *= -1
2692        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2693        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2694    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2695    print ' Broadening time = %.2fs'%(time.time()-time0)
2696   
2697def CalcStackingSADP(Layers,debug):
2698   
2699# Scattering factors
2700    SetStackingSF(Layers,debug)
2701# Controls & sequences
2702    laueId,controls = SetStackingClay(Layers,'SADP')
2703# cell & atoms
2704    Nlayers = SetCellAtoms(Layers)   
2705# Transitions
2706    SetStackingTrans(Layers,Nlayers)
2707# result as Sadp
2708    Nspec = 20001       
2709    spec = np.zeros(Nspec,dtype='double')   
2710    time0 = time.time()
2711    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
2712    Sapd = np.zeros((256,256))
2713    iB = 0
2714    for i in range(hkLim):
2715        iF = iB+Nblk
2716        p1 = 127+int(i*Incr)
2717        p2 = 128-int(i*Incr)
2718        if Nblk == 128:
2719            if i:
2720                Sapd[128:,p1] = spec[iB:iF]
2721                Sapd[:128,p1] = spec[iF:iB:-1]
2722            Sapd[128:,p2] = spec[iB:iF]
2723            Sapd[:128,p2] = spec[iF:iB:-1]
2724        else:
2725            if i:
2726                Sapd[:,p1] = spec[iB:iF]
2727            Sapd[:,p2] = spec[iB:iF]
2728        iB += Nblk
2729    Layers['Sadp']['Img'] = Sapd
2730    print ' GETSAD time = %.2fs'%(time.time()-time0)
2731#    GSASIIpath.IPyBreak()
2732   
2733#testing data
2734NeedTestData = True
2735def TestData():
2736    'needs a doc string'
2737#    global NeedTestData
2738    global bakType
2739    bakType = 'chebyschev'
2740    global xdata
2741    xdata = np.linspace(4.0,40.0,36000)
2742    global parmDict0
2743    parmDict0 = {
2744        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
2745        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
2746        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
2747        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
2748        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
2749        'Back0':5.384,'Back1':-0.015,'Back2':.004,
2750        }
2751    global parmDict1
2752    parmDict1 = {
2753        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
2754        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
2755        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
2756        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
2757        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
2758        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
2759        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
2760        'Back0':36.897,'Back1':-0.508,'Back2':.006,
2761        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2762        }
2763    global parmDict2
2764    parmDict2 = {
2765        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
2766        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
2767        'Back0':5.,'Back1':-0.02,'Back2':.004,
2768#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2769        }
2770    global varyList
2771    varyList = []
2772
2773def test0():
2774    if NeedTestData: TestData()
2775    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
2776    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
2777    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
2778    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
2779    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
2780    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
2781   
2782def test1():
2783    if NeedTestData: TestData()
2784    time0 = time.time()
2785    for i in range(100):
2786        getPeakProfile(parmDict1,xdata,varyList,bakType)
2787    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
2788   
2789def test2(name,delt):
2790    if NeedTestData: TestData()
2791    varyList = [name,]
2792    xdata = np.linspace(5.6,5.8,400)
2793    hplot = plotter.add('derivatives test for '+name).gca()
2794    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
2795    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2796    parmDict2[name] += delt
2797    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2798    hplot.plot(xdata,(y1-y0)/delt,'r+')
2799   
2800def test3(name,delt):
2801    if NeedTestData: TestData()
2802    names = ['pos','sig','gam','shl']
2803    idx = names.index(name)
2804    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
2805    xdata = np.linspace(5.6,5.8,800)
2806    dx = xdata[1]-xdata[0]
2807    hplot = plotter.add('derivatives test for '+name).gca()
2808    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
2809    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2810    myDict[name] += delt
2811    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2812    hplot.plot(xdata,(y1-y0)/delt,'r+')
2813
2814if __name__ == '__main__':
2815    import GSASIItestplot as plot
2816    global plotter
2817    plotter = plot.PlotNotebook()
2818#    test0()
2819#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
2820    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
2821        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
2822        test2(name,shft)
2823    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
2824        test3(name,shft)
2825    print "OK"
2826    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.