source: trunk/GSASIIpwd.py @ 3070

Last change on this file since 3070 was 3070, checked in by vondreele, 4 years ago

implement use of background pattern for PWDR data.
Used for plot & peak fitting; not Rietveld refinement
fix background parm copy issue - needed deepcopy

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