source: trunk/GSASIIpwd.py @ 2838

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

add a bit of printing to LS & show no. of penalty fxns

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