source: trunk/GSASIIpwd.py @ 2761

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

fix parameter bugs for reflectometry
add more substances

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