source: trunk/GSASIIpwd.py @ 3000

Last change on this file since 3000 was 3000, checked in by toby, 4 years ago

make the two frame version the trunk as we hit version 3000

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 111.5 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASII powder calculation module*
5==================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2017-08-11 22:34:54 +0000 (Fri, 11 Aug 2017) $
10# $Author: toby $
11# $Revision: 3000 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 3000 2017-08-11 22:34:54Z toby $
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: 3000 $")
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,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 list prevVaryList: Used in sequential refinements to override the
1522      variable list. Defaults as an empty list.
1523    :param bool oneCycle: True if only one cycle of fitting should be performed
1524    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1525      and derivType = controls['deriv type']. If None default values are used.
1526    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1527      Defaults to None, which means no updates are done.
1528    '''
1529    def GetBackgroundParms(parmList,Background):
1530        iBak = 0
1531        while True:
1532            try:
1533                bakName = 'Back;'+str(iBak)
1534                Background[0][iBak+3] = parmList[bakName]
1535                iBak += 1
1536            except KeyError:
1537                break
1538        iDb = 0
1539        while True:
1540            names = ['DebyeA;','DebyeR;','DebyeU;']
1541            try:
1542                for i,name in enumerate(names):
1543                    val = parmList[name+str(iDb)]
1544                    Background[1]['debyeTerms'][iDb][2*i] = val
1545                iDb += 1
1546            except KeyError:
1547                break
1548        iDb = 0
1549        while True:
1550            names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1551            try:
1552                for i,name in enumerate(names):
1553                    val = parmList[name+str(iDb)]
1554                    Background[1]['peaksList'][iDb][2*i] = val
1555                iDb += 1
1556            except KeyError:
1557                break
1558               
1559    def BackgroundPrint(Background,sigDict):
1560        print 'Background coefficients for',Background[0][0],'function'
1561        ptfmt = "%12.5f"
1562        ptstr =  'value: '
1563        sigstr = 'esd  : '
1564        for i,back in enumerate(Background[0][3:]):
1565            ptstr += ptfmt % (back)
1566            if Background[0][1]:
1567                prm = 'Back;'+str(i)
1568                if prm in sigDict:
1569                    sigstr += ptfmt % (sigDict[prm])
1570                else:
1571                    sigstr += " "*12
1572            if len(ptstr) > 75:
1573                print ptstr
1574                if Background[0][1]: print sigstr
1575                ptstr =  'value: '
1576                sigstr = 'esd  : '
1577        if len(ptstr) > 8:
1578            print ptstr
1579            if Background[0][1]: print sigstr
1580
1581        if Background[1]['nDebye']:
1582            parms = ['DebyeA;','DebyeR;','DebyeU;']
1583            print 'Debye diffuse scattering coefficients'
1584            ptfmt = "%12.5f"
1585            print ' term       DebyeA       esd        DebyeR       esd        DebyeU        esd'
1586            for term in range(Background[1]['nDebye']):
1587                line = ' term %d'%(term)
1588                for ip,name in enumerate(parms):
1589                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1590                    if name+str(term) in sigDict:
1591                        line += ptfmt%(sigDict[name+str(term)])
1592                    else:
1593                        line += " "*12
1594                print line
1595        if Background[1]['nPeaks']:
1596            print 'Coefficients for Background Peaks'
1597            ptfmt = "%15.3f"
1598            for j,pl in enumerate(Background[1]['peaksList']):
1599                names =  'peak %3d:'%(j+1)
1600                ptstr =  'values  :'
1601                sigstr = 'esds    :'
1602                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1603                    val = pl[2*i]
1604                    prm = lbl+";"+str(j)
1605                    names += '%15s'%(prm)
1606                    ptstr += ptfmt%(val)
1607                    if prm in sigDict:
1608                        sigstr += ptfmt%(sigDict[prm])
1609                    else:
1610                        sigstr += " "*15
1611                print names
1612                print ptstr
1613                print sigstr
1614                           
1615    def SetInstParms(Inst):
1616        dataType = Inst['Type'][0]
1617        insVary = []
1618        insNames = []
1619        insVals = []
1620        for parm in Inst:
1621            insNames.append(parm)
1622            insVals.append(Inst[parm][1])
1623            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1624                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1625                    insVary.append(parm)
1626        instDict = dict(zip(insNames,insVals))
1627        instDict['X'] = max(instDict['X'],0.01)
1628        instDict['Y'] = max(instDict['Y'],0.01)
1629        if 'SH/L' in instDict:
1630            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1631        return dataType,instDict,insVary
1632       
1633    def GetInstParms(parmDict,Inst,varyList,Peaks):
1634        for name in Inst:
1635            Inst[name][1] = parmDict[name]
1636        iPeak = 0
1637        while True:
1638            try:
1639                sigName = 'sig'+str(iPeak)
1640                pos = parmDict['pos'+str(iPeak)]
1641                if sigName not in varyList:
1642                    if 'C' in Inst['Type'][0]:
1643                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1644                    else:
1645                        dsp = G2lat.Pos2dsp(Inst,pos)
1646                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1647                gamName = 'gam'+str(iPeak)
1648                if gamName not in varyList:
1649                    if 'C' in Inst['Type'][0]:
1650                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1651                    else:
1652                        dsp = G2lat.Pos2dsp(Inst,pos)
1653                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1654                iPeak += 1
1655            except KeyError:
1656                break
1657       
1658    def InstPrint(Inst,sigDict):
1659        print 'Instrument Parameters:'
1660        ptfmt = "%12.6f"
1661        ptlbls = 'names :'
1662        ptstr =  'values:'
1663        sigstr = 'esds  :'
1664        for parm in Inst:
1665            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1666                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1667                ptlbls += "%s" % (parm.center(12))
1668                ptstr += ptfmt % (Inst[parm][1])
1669                if parm in sigDict:
1670                    sigstr += ptfmt % (sigDict[parm])
1671                else:
1672                    sigstr += 12*' '
1673        print ptlbls
1674        print ptstr
1675        print sigstr
1676
1677    def SetPeaksParms(dataType,Peaks):
1678        peakNames = []
1679        peakVary = []
1680        peakVals = []
1681        if 'C' in dataType:
1682            names = ['pos','int','sig','gam']
1683        else:
1684            names = ['pos','int','alp','bet','sig','gam']
1685        for i,peak in enumerate(Peaks):
1686            for j,name in enumerate(names):
1687                peakVals.append(peak[2*j])
1688                parName = name+str(i)
1689                peakNames.append(parName)
1690                if peak[2*j+1]:
1691                    peakVary.append(parName)
1692        return dict(zip(peakNames,peakVals)),peakVary
1693               
1694    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1695        if 'C' in Inst['Type'][0]:
1696            names = ['pos','int','sig','gam']
1697        else:   #'T'
1698            names = ['pos','int','alp','bet','sig','gam']
1699        for i,peak in enumerate(Peaks):
1700            pos = parmDict['pos'+str(i)]
1701            if 'difC' in Inst:
1702                dsp = pos/Inst['difC'][1]
1703            for j in range(len(names)):
1704                parName = names[j]+str(i)
1705                if parName in varyList:
1706                    peak[2*j] = parmDict[parName]
1707                elif 'alpha' in parName:
1708                    peak[2*j] = parmDict['alpha']/dsp
1709                elif 'beta' in parName:
1710                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1711                elif 'sig' in parName:
1712                    if 'C' in Inst['Type'][0]:
1713                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1714                    else:
1715                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1716                elif 'gam' in parName:
1717                    if 'C' in Inst['Type'][0]:
1718                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1719                    else:
1720                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1721                       
1722    def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
1723        print 'Peak coefficients:'
1724        if 'C' in dataType:
1725            names = ['pos','int','sig','gam']
1726        else:   #'T'
1727            names = ['pos','int','alp','bet','sig','gam']           
1728        head = 13*' '
1729        for name in names:
1730            if name in ['alp','bet']:
1731                head += name.center(8)+'esd'.center(8)
1732            else:
1733                head += name.center(10)+'esd'.center(10)
1734        head += 'bins'.center(8)
1735        print head
1736        if 'C' in dataType:
1737            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1738        else:
1739            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1740        for i,peak in enumerate(Peaks):
1741            ptstr =  ':'
1742            for j in range(len(names)):
1743                name = names[j]
1744                parName = name+str(i)
1745                ptstr += ptfmt[name] % (parmDict[parName])
1746                if parName in varyList:
1747                    ptstr += ptfmt[name] % (sigDict[parName])
1748                else:
1749                    if name in ['alp','bet']:
1750                        ptstr += 8*' '
1751                    else:
1752                        ptstr += 10*' '
1753            ptstr += '%9.2f'%(ptsperFW[i])
1754            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1755               
1756    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1757        parmdict.update(zip(varylist,values))
1758        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1759           
1760    def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):       
1761        parmdict.update(zip(varylist,values))
1762        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1763        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1764        if dlg:
1765            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1766            if not GoOn:
1767                return -M           #abort!!
1768        return M
1769       
1770    if controls:
1771        Ftol = controls['min dM/M']
1772    else:
1773        Ftol = 0.0001
1774    if oneCycle:
1775        Ftol = 1.0
1776    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1777    yc *= 0.                            #set calcd ones to zero
1778    yb *= 0.
1779    yd *= 0.
1780    cw = x[1:]-x[:-1]
1781    xBeg = np.searchsorted(x,Limits[0])
1782    xFin = np.searchsorted(x,Limits[1])+1
1783    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1784    dataType,insDict,insVary = SetInstParms(Inst)
1785    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1786    parmDict = {}
1787    parmDict.update(bakDict)
1788    parmDict.update(insDict)
1789    parmDict.update(peakDict)
1790    parmDict['Pdabc'] = []      #dummy Pdabc
1791    parmDict.update(Inst2)      #put in real one if there
1792    if prevVaryList:
1793        varyList = prevVaryList[:]
1794    else:
1795        varyList = bakVary+insVary+peakVary
1796    fullvaryList = varyList[:]
1797    while True:
1798        begin = time.time()
1799        values =  np.array(Dict2Values(parmDict, varyList))
1800        Rvals = {}
1801        badVary = []
1802        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1803               args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1804        ncyc = int(result[2]['nfev']/2)
1805        runtime = time.time()-begin   
1806        chisq = np.sum(result[2]['fvec']**2)
1807        Values2Dict(parmDict, varyList, result[0])
1808        Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1809        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1810        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1811        if ncyc:
1812            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1813        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1814        sig = [0]*len(varyList)
1815        if len(varyList) == 0: break  # if nothing was refined
1816        try:
1817            sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1818            if np.any(np.isnan(sig)):
1819                print '*** Least squares aborted - some invalid esds possible ***'
1820            break                   #refinement succeeded - finish up!
1821        except ValueError:          #result[1] is None on singular matrix
1822            print '**** Refinement failed - singular matrix ****'
1823            Ipvt = result[2]['ipvt']
1824            for i,ipvt in enumerate(Ipvt):
1825                if not np.sum(result[2]['fjac'],axis=1)[i]:
1826                    print 'Removing parameter: ',varyList[ipvt-1]
1827                    badVary.append(varyList[ipvt-1])
1828                    del(varyList[ipvt-1])
1829                    break
1830            else: # nothing removed
1831                break
1832    if dlg: dlg.Destroy()
1833    sigDict = dict(zip(varyList,sig))
1834    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]
1835    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1836    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1837    GetBackgroundParms(parmDict,Background)
1838    if bakVary: BackgroundPrint(Background,sigDict)
1839    GetInstParms(parmDict,Inst,varyList,Peaks)
1840    if insVary: InstPrint(Inst,sigDict)
1841    GetPeaksParms(Inst,parmDict,Peaks,varyList)
1842    binsperFWHM = []
1843    for peak in Peaks:
1844        FWHM = getFWHM(peak[0],Inst)
1845        try:
1846            binsperFWHM.append(FWHM/cw[x.searchsorted(peak[0])])
1847        except IndexError:
1848            binsperFWHM.append(0.)
1849    if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
1850    if len(binsperFWHM):
1851        if min(binsperFWHM) < 1.:
1852            print '*** Warning: calculated peak widths are too narrow to refine profile coefficients ***'
1853            if 'T' in Inst['Type'][0]:
1854                print ' Manually increase sig-0, 1, or 2 in Instrument Parameters'
1855            else:
1856                print ' Manually increase W in Instrument Parameters'
1857        elif min(binsperFWHM) < 4.:
1858            print '*** Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement ***'
1859            print '*** recommended is 6-10; you have %.2f ***'%(min(binsperFWHM))
1860    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
1861   
1862def calcIncident(Iparm,xdata):
1863    'needs a doc string'
1864
1865    def IfunAdv(Iparm,xdata):
1866        Itype = Iparm['Itype']
1867        Icoef = Iparm['Icoeff']
1868        DYI = np.ones((12,xdata.shape[0]))
1869        YI = np.ones_like(xdata)*Icoef[0]
1870       
1871        x = xdata/1000.                 #expressions are in ms
1872        if Itype == 'Exponential':
1873            for i in [1,3,5,7,9]:
1874                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1875                YI += Icoef[i]*Eterm
1876                DYI[i] *= Eterm
1877                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
1878        elif 'Maxwell'in Itype:
1879            Eterm = np.exp(-Icoef[2]/x**2)
1880            DYI[1] = Eterm/x**5
1881            DYI[2] = -Icoef[1]*DYI[1]/x**2
1882            YI += (Icoef[1]*Eterm/x**5)
1883            if 'Exponential' in Itype:
1884                for i in range(3,11,2):
1885                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1886                    YI += Icoef[i]*Eterm
1887                    DYI[i] *= Eterm
1888                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
1889            else:   #Chebyschev
1890                T = (2./x)-1.
1891                Ccof = np.ones((12,xdata.shape[0]))
1892                Ccof[1] = T
1893                for i in range(2,12):
1894                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1895                for i in range(1,10):
1896                    YI += Ccof[i]*Icoef[i+2]
1897                    DYI[i+2] =Ccof[i]
1898        return YI,DYI
1899       
1900    Iesd = np.array(Iparm['Iesd'])
1901    Icovar = Iparm['Icovar']
1902    YI,DYI = IfunAdv(Iparm,xdata)
1903    YI = np.where(YI>0,YI,1.)
1904    WYI = np.zeros_like(xdata)
1905    vcov = np.zeros((12,12))
1906    k = 0
1907    for i in range(12):
1908        for j in range(i,12):
1909            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1910            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1911            k += 1
1912    M = np.inner(vcov,DYI.T)
1913    WYI = np.sum(M*DYI,axis=0)
1914    WYI = np.where(WYI>0.,WYI,0.)
1915    return YI,WYI
1916   
1917################################################################################
1918# Reflectometry calculations
1919################################################################################
1920
1921def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
1922    print 'fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0])
1923   
1924    class RandomDisplacementBounds(object):
1925        """random displacement with bounds"""
1926        def __init__(self, xmin, xmax, stepsize=0.5):
1927            self.xmin = xmin
1928            self.xmax = xmax
1929            self.stepsize = stepsize
1930   
1931        def __call__(self, x):
1932            """take a random step but ensure the new position is within the bounds"""
1933            while True:
1934                # this could be done in a much more clever way, but it will work for example purposes
1935                steps = self.xmax-self.xmin
1936                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
1937                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
1938                    break
1939            return xnew
1940   
1941    def GetModelParms():
1942        parmDict = {}
1943        varyList = []
1944        values = []
1945        bounds = []
1946        parmDict['dQ type'] = data['dQ type']
1947        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
1948        for parm in ['Scale','FltBack']:
1949            parmDict[parm] = data[parm][0]
1950            if data[parm][1]:
1951                varyList.append(parm)
1952                values.append(data[parm][0])
1953                bounds.append(Bounds[parm])
1954        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
1955        parmDict['nLayers'] = len(parmDict['Layer Seq'])
1956        for ilay,layer in enumerate(data['Layers']):
1957            name = layer['Name']
1958            cid = str(ilay)+';'
1959            parmDict[cid+'Name'] = name
1960            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
1961                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
1962                if layer.get(parm,[0.,False])[1]:
1963                    varyList.append(cid+parm)
1964                    value = layer[parm][0]
1965                    values.append(value)
1966                    if value:
1967                        bound = [value*Bfac,value/Bfac]
1968                    else:
1969                        bound = [0.,10.]
1970                    bounds.append(bound)
1971            if name not in ['vacuum','unit scatter']:
1972                parmDict[cid+'rho'] = Substances[name]['Scatt density']
1973                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
1974        return parmDict,varyList,values,bounds
1975   
1976    def SetModelParms():
1977        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
1978        if 'Scale' in varyList:
1979            data['Scale'][0] = parmDict['Scale']
1980            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
1981        print line
1982        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
1983        if 'FltBack' in varyList:
1984            data['FltBack'][0] = parmDict['FltBack']
1985            line += ' esd: %15.3g'%(sigDict['FltBack'])
1986        print line
1987        for ilay,layer in enumerate(data['Layers']):
1988            name = layer['Name']
1989            print ' Parameters for layer: %d %s'%(ilay,name)
1990            cid = str(ilay)+';'
1991            line = ' '
1992            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
1993            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
1994            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
1995                if parm in layer:
1996                    if parm == 'Rough':
1997                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
1998                    else:
1999                        layer[parm][0] = parmDict[cid+parm]
2000                    line += ' %s: %.3f'%(parm,layer[parm][0])
2001                    if cid+parm in varyList:
2002                        line += ' esd: %.3g'%(sigDict[cid+parm])
2003            print line
2004            print line2
2005   
2006    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2007        parmDict.update(zip(varyList,values))
2008        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2009        return M
2010   
2011    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2012        parmDict.update(zip(varyList,values))
2013        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2014        return np.sum(M**2)
2015   
2016    def getREFD(Q,Qsig,parmDict):
2017        Ic = np.ones_like(Q)*parmDict['FltBack']
2018        Scale = parmDict['Scale']
2019        Nlayers = parmDict['nLayers']
2020        Res = parmDict['Res']
2021        depth = np.zeros(Nlayers)
2022        rho = np.zeros(Nlayers)
2023        irho = np.zeros(Nlayers)
2024        sigma = np.zeros(Nlayers)
2025        for ilay,lay in enumerate(parmDict['Layer Seq']):
2026            cid = str(lay)+';'
2027            depth[ilay] = parmDict[cid+'Thick']
2028            sigma[ilay] = parmDict[cid+'Rough']
2029            if parmDict[cid+'Name'] == u'unit scatter':
2030                rho[ilay] = parmDict[cid+'DenMul']
2031                irho[ilay] = parmDict[cid+'iDenMul']
2032            elif 'vacuum' != parmDict[cid+'Name']:
2033                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
2034                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
2035            if cid+'Mag SLD' in parmDict:
2036                rho[ilay] += parmDict[cid+'Mag SLD']
2037        if parmDict['dQ type'] == 'None':
2038            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2039        elif 'const' in parmDict['dQ type']:
2040            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
2041        else:       #dQ/Q in data
2042            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
2043        Ic += AB*Scale
2044        return Ic
2045       
2046    def estimateT0(takestep):
2047        Mmax = -1.e-10
2048        Mmin = 1.e10
2049        for i in range(100):
2050            x0 = takestep(values)
2051            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2052            Mmin = min(M,Mmin)
2053            MMax = max(M,Mmax)
2054        return 1.5*(MMax-Mmin)
2055
2056    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2057    if data.get('2% weight'):
2058        wt = 1./(0.02*Io)**2
2059    Qmin = Limits[1][0]
2060    Qmax = Limits[1][1]
2061    wtFactor = ProfDict['wtFactor']
2062    Bfac = data['Toler']
2063    Ibeg = np.searchsorted(Q,Qmin)
2064    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
2065    Ic[:] = 0
2066    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
2067              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
2068    parmDict,varyList,values,bounds = GetModelParms()
2069    Msg = 'Failed to converge'
2070    if varyList:
2071        if data['Minimizer'] == 'LMLS': 
2072            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
2073                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2074            parmDict.update(zip(varyList,result[0]))
2075            chisq = np.sum(result[2]['fvec']**2)
2076            ncalc = result[2]['nfev']
2077            covM = result[1]
2078            newVals = result[0]
2079        elif data['Minimizer'] == 'Basin Hopping':
2080            xyrng = np.array(bounds).T
2081            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
2082            T0 = estimateT0(take_step)
2083            print ' Estimated temperature: %.3g'%(T0)
2084            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
2085                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
2086                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
2087            chisq = result.fun
2088            ncalc = result.nfev
2089            newVals = result.x
2090            covM = []
2091        elif data['Minimizer'] == 'MC/SA Anneal':
2092            xyrng = np.array(bounds).T
2093            result = G2mth.anneal(sumREFD, values, 
2094                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
2095                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
2096                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
2097                ranRange=0.20,autoRan=False,dlg=None)
2098            newVals = result[0]
2099            parmDict.update(zip(varyList,newVals))
2100            chisq = result[1]
2101            ncalc = result[3]
2102            covM = []
2103            print ' MC/SA final temperature: %.4g'%(result[2])
2104        elif data['Minimizer'] == 'L-BFGS-B':
2105            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
2106                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2107            parmDict.update(zip(varyList,result['x']))
2108            chisq = result.fun
2109            ncalc = result.nfev
2110            newVals = result.x
2111            covM = []
2112    else:   #nothing varied
2113        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2114        chisq = np.sum(M**2)
2115        ncalc = 0
2116        covM = []
2117        sig = []
2118        sigDict = {}
2119        result = []
2120    Rvals = {}
2121    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
2122    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
2123    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
2124    Ib[Ibeg:Ifin] = parmDict['FltBack']
2125    try:
2126        if not len(varyList):
2127            Msg += ' - nothing refined'
2128            raise ValueError
2129        Nans = np.isnan(newVals)
2130        if np.any(Nans):
2131            name = varyList[Nans.nonzero(True)[0]]
2132            Msg += ' Nan result for '+name+'!'
2133            raise ValueError
2134        Negs = np.less_equal(newVals,0.)
2135        if np.any(Negs):
2136            indx = Negs.nonzero()
2137            name = varyList[indx[0][0]]
2138            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
2139                Msg += ' negative coefficient for '+name+'!'
2140                raise ValueError
2141        if len(covM):
2142            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
2143            covMatrix = covM*Rvals['GOF']
2144        else:
2145            sig = np.zeros(len(varyList))
2146            covMatrix = []
2147        sigDict = dict(zip(varyList,sig))
2148        print ' Results of reflectometry data modelling fit:'
2149        print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
2150        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
2151        SetModelParms()
2152        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
2153    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
2154        print Msg
2155        return False,0,0,0,0,0,0,Msg
2156       
2157def makeSLDprofile(data,Substances):
2158   
2159    sq2 = np.sqrt(2.)
2160    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2161    Nlayers = len(laySeq)
2162    laySeq = np.array(laySeq,dtype=int)
2163    interfaces = np.zeros(Nlayers)
2164    rho = np.zeros(Nlayers)
2165    sigma = np.zeros(Nlayers)
2166    name = data['Layers'][0]['Name']
2167    thick = 0.
2168    for ilay,lay in enumerate(laySeq):
2169        layer = data['Layers'][lay]
2170        name = layer['Name']
2171        if 'Thick' in layer:
2172            thick += layer['Thick'][0]
2173            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
2174        if 'Rough' in layer:
2175            sigma[ilay] = max(0.001,layer['Rough'][0])
2176        if name != 'vacuum':
2177            if name == 'unit scatter':
2178                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
2179            else:
2180                rrho = Substances[name]['Scatt density']
2181                irho = Substances[name]['XImag density']
2182                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
2183        if 'Mag SLD' in layer:
2184            rho[ilay] += layer['Mag SLD'][0]
2185    name = data['Layers'][-1]['Name']
2186    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
2187    xr = np.flipud(x)
2188    interfaces[-1] = x[-1]
2189    y = np.ones_like(x)*rho[0]
2190    iBeg = 0
2191    for ilayer in range(Nlayers-1):
2192        delt = rho[ilayer+1]-rho[ilayer]
2193        iPos = np.searchsorted(x,interfaces[ilayer])
2194        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
2195        iBeg = iPos
2196    return x,xr,y   
2197
2198def REFDModelFxn(Profile,Inst,Limits,Substances,data):
2199   
2200    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2201    Qmin = Limits[1][0]
2202    Qmax = Limits[1][1]
2203    iBeg = np.searchsorted(Q,Qmin)
2204    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2205    Ib[:] = data['FltBack'][0]
2206    Ic[:] = 0
2207    Scale = data['Scale'][0]
2208    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2209    Nlayers = len(laySeq)
2210    depth = np.zeros(Nlayers)
2211    rho = np.zeros(Nlayers)
2212    irho = np.zeros(Nlayers)
2213    sigma = np.zeros(Nlayers)
2214    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
2215        layer = data['Layers'][lay]
2216        name = layer['Name']
2217        if 'Thick' in layer:    #skips first & last layers
2218            depth[ilay] = layer['Thick'][0]
2219        if 'Rough' in layer:    #skips first layer
2220            sigma[ilay] = layer['Rough'][0]
2221        if 'unit scatter' == name:
2222            rho[ilay] = layer['DenMul'][0]
2223            irho[ilay] = layer['iDenMul'][0]
2224        else:
2225            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2226            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
2227        if 'Mag SLD' in layer:
2228            rho[ilay] += layer['Mag SLD'][0]
2229    if data['dQ type'] == 'None':
2230        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2231    elif 'const' in data['dQ type']:
2232        res = data['Resolution'][0]/(100.*sateln2)
2233        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
2234    else:       #dQ/Q in data
2235        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
2236    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
2237
2238def abeles(kz, depth, rho, irho=0, sigma=0):
2239    """
2240    Optical matrix form of the reflectivity calculation.
2241    O.S. Heavens, Optical Properties of Thin Solid Films
2242   
2243    Reflectometry as a function of kz for a set of slabs.
2244
2245    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
2246        This is :math:`\\tfrac12 Q_z`.       
2247    :param depth:  float[m] (Ang).
2248        thickness of each layer.  The thickness of the incident medium
2249        and substrate are ignored.
2250    :param rho:  float[n,k] (1e-6/Ang^2)
2251        Real scattering length density for each layer for each kz
2252    :param irho:  float[n,k] (1e-6/Ang^2)
2253        Imaginary scattering length density for each layer for each kz
2254        Note: absorption cross section mu = 2 irho/lambda for neutrons
2255    :param sigma: float[m-1] (Ang)
2256        interfacial roughness.  This is the roughness between a layer
2257        and the previous layer. The sigma array should have m-1 entries.
2258
2259    Slabs are ordered with the surface SLD at index 0 and substrate at
2260    index -1, or reversed if kz < 0.
2261    """
2262    def calc(kz, depth, rho, irho, sigma):
2263        if len(kz) == 0: return kz
2264   
2265        # Complex index of refraction is relative to the incident medium.
2266        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
2267        # in place of kz^2, and ignoring rho_o
2268        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
2269        k = kz
2270   
2271        # According to Heavens, the initial matrix should be [ 1 F; F 1],
2272        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
2273        # multiply versus some coding convenience.
2274        B11 = 1
2275        B22 = 1
2276        B21 = 0
2277        B12 = 0
2278        for i in range(0, len(depth)-1):
2279            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
2280            F = (k - k_next) / (k + k_next)
2281            F *= np.exp(-2*k*k_next*sigma[i]**2)
2282            #print "==== layer",i
2283            #print "kz:", kz
2284            #print "k:", k
2285            #print "k_next:",k_next
2286            #print "F:",F
2287            #print "rho:",rho[:,i+1]
2288            #print "irho:",irho[:,i+1]
2289            #print "d:",depth[i],"sigma:",sigma[i]
2290            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
2291            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
2292            M21 = F*M11
2293            M12 = F*M22
2294            C1 = B11*M11 + B21*M12
2295            C2 = B11*M21 + B21*M22
2296            B11 = C1
2297            B21 = C2
2298            C1 = B12*M11 + B22*M12
2299            C2 = B12*M21 + B22*M22
2300            B12 = C1
2301            B22 = C2
2302            k = k_next
2303   
2304        r = B12/B11
2305        return np.absolute(r)**2
2306
2307    if np.isscalar(kz): kz = np.asarray([kz], 'd')
2308
2309    m = len(depth)
2310
2311    # Make everything into arrays
2312    depth = np.asarray(depth,'d')
2313    rho = np.asarray(rho,'d')
2314    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
2315    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
2316
2317    # Repeat rho,irho columns as needed
2318    if len(rho.shape) == 1:
2319        rho = rho[None,:]
2320        irho = irho[None,:]
2321
2322    return calc(kz, depth, rho, irho, sigma)
2323   
2324def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
2325    y = abeles(kz, depth, rho, irho, sigma)
2326    s = dq/2.
2327    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
2328    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
2329    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
2330    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
2331    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
2332    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
2333    y *= 0.137023
2334    return y
2335       
2336def makeRefdFFT(Limits,Profile):
2337    print 'make fft'
2338    Q,Io = Profile[:2]
2339    Qmin = Limits[1][0]
2340    Qmax = Limits[1][1]
2341    iBeg = np.searchsorted(Q,Qmin)
2342    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2343    Qf = np.linspace(0.,Q[iFin-1],5000)
2344    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
2345    If = QI(Qf)*Qf**4
2346    R = np.linspace(0.,5000.,5000)
2347    F = fft.rfft(If)
2348    return R,F
2349
2350   
2351################################################################################
2352# Stacking fault simulation codes
2353################################################################################
2354
2355def GetStackParms(Layers):
2356   
2357    Parms = []
2358#cell parms
2359    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
2360        Parms.append('cellA')
2361        Parms.append('cellC')
2362    else:
2363        Parms.append('cellA')
2364        Parms.append('cellB')
2365        Parms.append('cellC')
2366        if Layers['Laue'] != 'mmm':
2367            Parms.append('cellG')
2368#Transition parms
2369    for iY in range(len(Layers['Layers'])):
2370        for iX in range(len(Layers['Layers'])):
2371            Parms.append('TransP;%d;%d'%(iY,iX))
2372            Parms.append('TransX;%d;%d'%(iY,iX))
2373            Parms.append('TransY;%d;%d'%(iY,iX))
2374            Parms.append('TransZ;%d;%d'%(iY,iX))
2375    return Parms
2376
2377def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
2378    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
2379   
2380    :param dict Layers: dict with following items
2381
2382      ::
2383
2384       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
2385       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
2386       'Layers':[],'Stacking':[],'Transitions':[]}
2387       
2388    :param str ctrls: controls string to be written on DIFFaX controls.dif file
2389    :param float scale: scale factor
2390    :param dict background: background parameters
2391    :param list limits: min/max 2-theta to be calculated
2392    :param dict inst: instrument parameters dictionary
2393    :param list profile: powder pattern data
2394   
2395    Note that parameters all updated in place   
2396    '''
2397    import atmdata
2398    path = sys.path
2399    for name in path:
2400        if 'bin' in name:
2401            DIFFaX = name+'/DIFFaX.exe'
2402            print ' Execute ',DIFFaX
2403            break
2404    # make form factor file that DIFFaX wants - atom types are GSASII style
2405    sf = open('data.sfc','w')
2406    sf.write('GSASII special form factor file for DIFFaX\n\n')
2407    atTypes = Layers['AtInfo'].keys()
2408    if 'H' not in atTypes:
2409        atTypes.insert(0,'H')
2410    for atType in atTypes:
2411        if atType == 'H': 
2412            blen = -.3741
2413        else:
2414            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
2415        Adat = atmdata.XrayFF[atType]
2416        text = '%4s'%(atType.ljust(4))
2417        for i in range(4):
2418            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
2419        text += '%11.6f%11.6f'%(Adat['fc'],blen)
2420        text += '%3d\n'%(Adat['Z'])
2421        sf.write(text)
2422    sf.close()
2423    #make DIFFaX control.dif file - future use GUI to set some of these flags
2424    cf = open('control.dif','w')
2425    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2426        x0 = profile[0]
2427        iBeg = np.searchsorted(x0,limits[0])
2428        iFin = np.searchsorted(x0,limits[1])+1
2429        if iFin-iBeg > 20000:
2430            iFin = iBeg+20000
2431        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2432        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2433        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
2434    else:
2435        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2436        inst = {'Type':['XSC','XSC',]}
2437    cf.close()
2438    #make DIFFaX data file
2439    df = open('GSASII-DIFFaX.dat','w')
2440    df.write('INSTRUMENTAL\n')
2441    if 'X' in inst['Type'][0]:
2442        df.write('X-RAY\n')
2443    elif 'N' in inst['Type'][0]:
2444        df.write('NEUTRON\n')
2445    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2446        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
2447        U = ateln2*inst['U'][1]/10000.
2448        V = ateln2*inst['V'][1]/10000.
2449        W = ateln2*inst['W'][1]/10000.
2450        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2451        HW = np.sqrt(np.mean(HWHM))
2452    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
2453        if 'Mean' in Layers['selInst']:
2454            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
2455        elif 'Gaussian' in Layers['selInst']:
2456            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
2457        else:
2458            df.write('None\n')
2459    else:
2460        df.write('0.10\nNone\n')
2461    df.write('STRUCTURAL\n')
2462    a,b,c = Layers['Cell'][1:4]
2463    gam = Layers['Cell'][6]
2464    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
2465    laue = Layers['Laue']
2466    if laue == '2/m(ab)':
2467        laue = '2/m(1)'
2468    elif laue == '2/m(c)':
2469        laue = '2/m(2)'
2470    if 'unknown' in Layers['Laue']:
2471        df.write('%s %.3f\n'%(laue,Layers['Toler']))
2472    else:   
2473        df.write('%s\n'%(laue))
2474    df.write('%d\n'%(len(Layers['Layers'])))
2475    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
2476        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
2477    layerNames = []
2478    for layer in Layers['Layers']:
2479        layerNames.append(layer['Name'])
2480    for il,layer in enumerate(Layers['Layers']):
2481        if layer['SameAs']:
2482            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
2483            continue
2484        df.write('LAYER %d\n'%(il+1))
2485        if '-1' in layer['Symm']:
2486            df.write('CENTROSYMMETRIC\n')
2487        else:
2488            df.write('NONE\n')
2489        for ia,atom in enumerate(layer['Atoms']):
2490            [name,atype,x,y,z,frac,Uiso] = atom
2491            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
2492                frac /= 2.
2493            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
2494    df.write('STACKING\n')
2495    df.write('%s\n'%(Layers['Stacking'][0]))
2496    if 'recursive' in Layers['Stacking'][0]:
2497        df.write('%s\n'%Layers['Stacking'][1])
2498    else:
2499        if 'list' in Layers['Stacking'][1]:
2500            Slen = len(Layers['Stacking'][2])
2501            iB = 0
2502            iF = 0
2503            while True:
2504                iF += 68
2505                if iF >= Slen:
2506                    break
2507                iF = min(iF,Slen)
2508                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
2509                iB = iF
2510        else:
2511            df.write('%s\n'%Layers['Stacking'][1])   
2512    df.write('TRANSITIONS\n')
2513    for iY in range(len(Layers['Layers'])):
2514        sumPx = 0.
2515        for iX in range(len(Layers['Layers'])):
2516            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
2517            p = round(p,3)
2518            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
2519            sumPx += p
2520        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
2521            print 'ERROR - Layer probabilities sum to ',sumPx,' DIFFaX will insist it = 1.0'
2522            df.close()
2523            os.remove('data.sfc')
2524            os.remove('control.dif')
2525            os.remove('GSASII-DIFFaX.dat')
2526            return       
2527    df.close()   
2528    time0 = time.time()
2529    try:
2530        subp.call(DIFFaX)
2531    except OSError:
2532        print ' DIFFax.exe is not available for this platform - under development'
2533    print ' DIFFaX time = %.2fs'%(time.time()-time0)
2534    if os.path.exists('GSASII-DIFFaX.spc'):
2535        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
2536        iFin = iBeg+Xpat.shape[1]
2537        bakType,backDict,backVary = SetBackgroundParms(background)
2538        backDict['Lam1'] = G2mth.getWave(inst)
2539        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2540        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
2541        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2542            rv = st.poisson(profile[3][iBeg:iFin])
2543            profile[1][iBeg:iFin] = rv.rvs()
2544            Z = np.ones_like(profile[3][iBeg:iFin])
2545            Z[1::2] *= -1
2546            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2547            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2548        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2549    #cleanup files..
2550        os.remove('GSASII-DIFFaX.spc')
2551    elif os.path.exists('GSASII-DIFFaX.sadp'):
2552        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
2553        Sadp = np.reshape(Sadp,(256,-1))
2554        Layers['Sadp']['Img'] = Sadp
2555        os.remove('GSASII-DIFFaX.sadp')
2556    os.remove('data.sfc')
2557    os.remove('control.dif')
2558    os.remove('GSASII-DIFFaX.dat')
2559   
2560def SetPWDRscan(inst,limits,profile):
2561   
2562    wave = G2mth.getMeanWave(inst)
2563    x0 = profile[0]
2564    iBeg = np.searchsorted(x0,limits[0])
2565    iFin = np.searchsorted(x0,limits[1])
2566    if iFin-iBeg > 20000:
2567        iFin = iBeg+20000
2568    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2569    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
2570    return iFin-iBeg
2571       
2572def SetStackingSF(Layers,debug):
2573# Load scattering factors into DIFFaX arrays
2574    import atmdata
2575    atTypes = Layers['AtInfo'].keys()
2576    aTypes = []
2577    for atype in atTypes:
2578        aTypes.append('%4s'%(atype.ljust(4)))
2579    SFdat = []
2580    for atType in atTypes:
2581        Adat = atmdata.XrayFF[atType]
2582        SF = np.zeros(9)
2583        SF[:8:2] = Adat['fa']
2584        SF[1:8:2] = Adat['fb']
2585        SF[8] = Adat['fc']
2586        SFdat.append(SF)
2587    SFdat = np.array(SFdat)
2588    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
2589   
2590def SetStackingClay(Layers,Type):
2591# Controls
2592    rand.seed()
2593    ranSeed = rand.randint(1,2**16-1)
2594    try:   
2595        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
2596            '6/m','6/mmm'].index(Layers['Laue'])+1
2597    except ValueError:  #for 'unknown'
2598        laueId = -1
2599    if 'SADP' in Type:
2600        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
2601        lmax = int(Layers['Sadp']['Lmax'])
2602    else:
2603        planeId = 0
2604        lmax = 0
2605# Sequences
2606    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
2607    try:
2608        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
2609    except ValueError:
2610        StkParm = -1
2611    if StkParm == 2:    #list
2612        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
2613        Nstk = len(StkSeq)
2614    else:
2615        Nstk = 1
2616        StkSeq = [0,]
2617    if StkParm == -1:
2618        StkParm = int(Layers['Stacking'][1])
2619    Wdth = Layers['Width'][0]
2620    mult = 1
2621    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
2622    LaueSym = Layers['Laue'].ljust(12)
2623    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
2624    return laueId,controls
2625   
2626def SetCellAtoms(Layers):
2627    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
2628# atoms in layers
2629    atTypes = Layers['AtInfo'].keys()
2630    AtomXOU = []
2631    AtomTp = []
2632    LayerSymm = []
2633    LayerNum = []
2634    layerNames = []
2635    Natm = 0
2636    Nuniq = 0
2637    for layer in Layers['Layers']:
2638        layerNames.append(layer['Name'])
2639    for il,layer in enumerate(Layers['Layers']):
2640        if layer['SameAs']:
2641            LayerNum.append(layerNames.index(layer['SameAs'])+1)
2642            continue
2643        else:
2644            LayerNum.append(il+1)
2645            Nuniq += 1
2646        if '-1' in layer['Symm']:
2647            LayerSymm.append(1)
2648        else:
2649            LayerSymm.append(0)
2650        for ia,atom in enumerate(layer['Atoms']):
2651            [name,atype,x,y,z,frac,Uiso] = atom
2652            Natm += 1
2653            AtomTp.append('%4s'%(atype.ljust(4)))
2654            Ta = atTypes.index(atype)+1
2655            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
2656    AtomXOU = np.array(AtomXOU)
2657    Nlayers = len(layerNames)
2658    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
2659    return Nlayers
2660   
2661def SetStackingTrans(Layers,Nlayers):
2662# Transitions
2663    TransX = []
2664    TransP = []
2665    for Ytrans in Layers['Transitions']:
2666        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
2667        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
2668    TransP = np.array(TransP,dtype='float').T
2669    TransX = np.array(TransX,dtype='float')
2670#    GSASIIpath.IPyBreak()
2671    pyx.pygettrans(Nlayers,TransP,TransX)
2672   
2673def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
2674# Scattering factors
2675    SetStackingSF(Layers,debug)
2676# Controls & sequences
2677    laueId,controls = SetStackingClay(Layers,'PWDR')
2678# cell & atoms
2679    Nlayers = SetCellAtoms(Layers)
2680    Volume = Layers['Cell'][7]   
2681# Transitions
2682    SetStackingTrans(Layers,Nlayers)
2683# PWDR scan
2684    Nsteps = SetPWDRscan(inst,limits,profile)
2685# result as Spec
2686    x0 = profile[0]
2687    profile[3] = np.zeros(len(profile[0]))
2688    profile[4] = np.zeros(len(profile[0]))
2689    profile[5] = np.zeros(len(profile[0]))
2690    iBeg = np.searchsorted(x0,limits[0])
2691    iFin = np.searchsorted(x0,limits[1])+1
2692    if iFin-iBeg > 20000:
2693        iFin = iBeg+20000
2694    Nspec = 20001       
2695    spec = np.zeros(Nspec,dtype='double')   
2696    time0 = time.time()
2697    pyx.pygetspc(controls,Nspec,spec)
2698    print ' GETSPC time = %.2fs'%(time.time()-time0)
2699    time0 = time.time()
2700    U = ateln2*inst['U'][1]/10000.
2701    V = ateln2*inst['V'][1]/10000.
2702    W = ateln2*inst['W'][1]/10000.
2703    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2704    HW = np.sqrt(np.mean(HWHM))
2705    BrdSpec = np.zeros(Nsteps)
2706    if 'Mean' in Layers['selInst']:
2707        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
2708    elif 'Gaussian' in Layers['selInst']:
2709        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
2710    else:
2711        BrdSpec = spec[:Nsteps]
2712    BrdSpec /= Volume
2713    iFin = iBeg+Nsteps
2714    bakType,backDict,backVary = SetBackgroundParms(background)
2715    backDict['Lam1'] = G2mth.getWave(inst)
2716    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2717    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
2718    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2719        try:
2720            rv = st.poisson(profile[3][iBeg:iFin])
2721            profile[1][iBeg:iFin] = rv.rvs()
2722        except ValueError:
2723            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
2724        Z = np.ones_like(profile[3][iBeg:iFin])
2725        Z[1::2] *= -1
2726        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2727        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2728    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2729    print ' Broadening time = %.2fs'%(time.time()-time0)
2730   
2731def CalcStackingSADP(Layers,debug):
2732   
2733# Scattering factors
2734    SetStackingSF(Layers,debug)
2735# Controls & sequences
2736    laueId,controls = SetStackingClay(Layers,'SADP')
2737# cell & atoms
2738    Nlayers = SetCellAtoms(Layers)   
2739# Transitions
2740    SetStackingTrans(Layers,Nlayers)
2741# result as Sadp
2742    Nspec = 20001       
2743    spec = np.zeros(Nspec,dtype='double')   
2744    time0 = time.time()
2745    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
2746    Sapd = np.zeros((256,256))
2747    iB = 0
2748    for i in range(hkLim):
2749        iF = iB+Nblk
2750        p1 = 127+int(i*Incr)
2751        p2 = 128-int(i*Incr)
2752        if Nblk == 128:
2753            if i:
2754                Sapd[128:,p1] = spec[iB:iF]
2755                Sapd[:128,p1] = spec[iF:iB:-1]
2756            Sapd[128:,p2] = spec[iB:iF]
2757            Sapd[:128,p2] = spec[iF:iB:-1]
2758        else:
2759            if i:
2760                Sapd[:,p1] = spec[iB:iF]
2761            Sapd[:,p2] = spec[iB:iF]
2762        iB += Nblk
2763    Layers['Sadp']['Img'] = Sapd
2764    print ' GETSAD time = %.2fs'%(time.time()-time0)
2765#    GSASIIpath.IPyBreak()
2766   
2767#testing data
2768NeedTestData = True
2769def TestData():
2770    'needs a doc string'
2771#    global NeedTestData
2772    global bakType
2773    bakType = 'chebyschev'
2774    global xdata
2775    xdata = np.linspace(4.0,40.0,36000)
2776    global parmDict0
2777    parmDict0 = {
2778        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
2779        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
2780        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
2781        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
2782        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
2783        'Back0':5.384,'Back1':-0.015,'Back2':.004,
2784        }
2785    global parmDict1
2786    parmDict1 = {
2787        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
2788        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
2789        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
2790        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
2791        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
2792        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
2793        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
2794        'Back0':36.897,'Back1':-0.508,'Back2':.006,
2795        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2796        }
2797    global parmDict2
2798    parmDict2 = {
2799        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
2800        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
2801        'Back0':5.,'Back1':-0.02,'Back2':.004,
2802#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2803        }
2804    global varyList
2805    varyList = []
2806
2807def test0():
2808    if NeedTestData: TestData()
2809    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
2810    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
2811    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
2812    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
2813    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
2814    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
2815   
2816def test1():
2817    if NeedTestData: TestData()
2818    time0 = time.time()
2819    for i in range(100):
2820        getPeakProfile(parmDict1,xdata,varyList,bakType)
2821    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
2822   
2823def test2(name,delt):
2824    if NeedTestData: TestData()
2825    varyList = [name,]
2826    xdata = np.linspace(5.6,5.8,400)
2827    hplot = plotter.add('derivatives test for '+name).gca()
2828    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
2829    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2830    parmDict2[name] += delt
2831    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2832    hplot.plot(xdata,(y1-y0)/delt,'r+')
2833   
2834def test3(name,delt):
2835    if NeedTestData: TestData()
2836    names = ['pos','sig','gam','shl']
2837    idx = names.index(name)
2838    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
2839    xdata = np.linspace(5.6,5.8,800)
2840    dx = xdata[1]-xdata[0]
2841    hplot = plotter.add('derivatives test for '+name).gca()
2842    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
2843    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2844    myDict[name] += delt
2845    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2846    hplot.plot(xdata,(y1-y0)/delt,'r+')
2847
2848if __name__ == '__main__':
2849    import GSASIItestplot as plot
2850    global plotter
2851    plotter = plot.PlotNotebook()
2852#    test0()
2853#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
2854    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
2855        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
2856        test2(name,shft)
2857    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
2858        test3(name,shft)
2859    print "OK"
2860    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.