source: trunk/GSASIIpwd.py @ 2808

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

small fixes to reflectometry calcs.
fix Hessian restraint bug - wrong sign for Vec
add theta, refl, (sig) import option for reflectometry

  • 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-04-25 19:01:51 +0000 (Tue, 25 Apr 2017) $
10# $Author: vondreele $
11# $Revision: 2808 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 2808 2017-04-25 19:01:51Z 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: 2808 $")
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   
957    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
958#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
959    Df /= np.sum(Df)
960    return Df
961
962def getdFCJVoigt3(pos,sig,gam,shl,xdata):
963    'needs a doc string'
964   
965    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
966#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
967    return Df,dFdp,dFds,dFdg,dFdsh
968
969def getPsVoigt(pos,sig,gam,xdata):
970    'needs a doc string'
971   
972    Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
973    Df /= np.sum(Df)
974    return Df
975
976def getdPsVoigt(pos,sig,gam,xdata):
977    'needs a doc string'
978   
979    Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
980    return Df,dFdp,dFds,dFdg
981
982def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
983    'needs a doc string'
984    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
985    Df /= np.sum(Df)
986    return Df 
987   
988def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
989    'needs a doc string'
990    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
991    return Df,dFdp,dFda,dFdb,dFds,dFdg   
992
993def ellipseSize(H,Sij,GB):
994    'needs a doc string'
995    HX = np.inner(H.T,GB)
996    lenHX = np.sqrt(np.sum(HX**2))
997    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
998    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
999    lenR = np.sqrt(np.sum(R**2))
1000    return lenR
1001
1002def ellipseSizeDerv(H,Sij,GB):
1003    'needs a doc string'
1004    lenR = ellipseSize(H,Sij,GB)
1005    delt = 0.001
1006    dRdS = np.zeros(6)
1007    for i in range(6):
1008        Sij[i] -= delt
1009        lenM = ellipseSize(H,Sij,GB)
1010        Sij[i] += 2.*delt
1011        lenP = ellipseSize(H,Sij,GB)
1012        Sij[i] -= delt
1013        dRdS[i] = (lenP-lenM)/(2.*delt)
1014    return lenR,dRdS
1015
1016def getHKLpeak(dmin,SGData,A,Inst=None):
1017    'needs a doc string'
1018    HKL = G2lat.GenHLaue(dmin,SGData,A)       
1019    HKLs = []
1020    for h,k,l,d in HKL:
1021        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1022        if not ext:
1023            if Inst == None:
1024                HKLs.append([h,k,l,d,0,-1])
1025            else:
1026                HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1027    return HKLs
1028
1029def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1030    'needs a doc string'
1031    HKLs = []
1032    vec = np.array(Vec)
1033    vstar = np.sqrt(G2lat.calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1034    dvec = 1./(maxH*vstar+1./dmin)
1035    HKL = G2lat.GenHLaue(dvec,SGData,A)       
1036    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1037    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1038    for h,k,l,d in HKL:
1039        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1040        if not ext and d >= dmin:
1041            HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1042        for dH in SSdH:
1043            if dH:
1044                DH = SSdH[dH]
1045                H = [h+DH[0],k+DH[1],l+DH[2]]
1046                d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1047                if d >= dmin:
1048                    HKLM = np.array([h,k,l,dH])
1049                    if G2spc.checkSSextc(HKLM,SSGData):
1050                        HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])   
1051#    GSASIIpath.IPyBreak()
1052    return G2lat.sortHKLd(HKLs,True,True,True)
1053
1054def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
1055    'needs a doc string'
1056   
1057    yb = getBackground('',parmDict,bakType,dataType,xdata)[0]
1058    yc = np.zeros_like(yb)
1059    cw = np.diff(xdata)
1060    cw = np.append(cw,cw[-1])
1061    if 'C' in dataType:
1062        shl = max(parmDict['SH/L'],0.002)
1063        Ka2 = False
1064        if 'Lam1' in parmDict.keys():
1065            Ka2 = True
1066            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1067            kRatio = parmDict['I(L2)/I(L1)']
1068        iPeak = 0
1069        while True:
1070            try:
1071                pos = parmDict['pos'+str(iPeak)]
1072                tth = (pos-parmDict['Zero'])
1073                intens = parmDict['int'+str(iPeak)]
1074                sigName = 'sig'+str(iPeak)
1075                if sigName in varyList:
1076                    sig = parmDict[sigName]
1077                else:
1078                    sig = G2mth.getCWsig(parmDict,tth)
1079                sig = max(sig,0.001)          #avoid neg sigma^2
1080                gamName = 'gam'+str(iPeak)
1081                if gamName in varyList:
1082                    gam = parmDict[gamName]
1083                else:
1084                    gam = G2mth.getCWgam(parmDict,tth)
1085                gam = max(gam,0.001)             #avoid neg gamma
1086                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1087                iBeg = np.searchsorted(xdata,pos-fmin)
1088                iFin = np.searchsorted(xdata,pos+fmin)
1089                if not iBeg+iFin:       #peak below low limit
1090                    iPeak += 1
1091                    continue
1092                elif not iBeg-iFin:     #peak above high limit
1093                    return yb+yc
1094                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1095                if Ka2:
1096                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1097                    iBeg = np.searchsorted(xdata,pos2-fmin)
1098                    iFin = np.searchsorted(xdata,pos2+fmin)
1099                    if iBeg-iFin:
1100                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1101                iPeak += 1
1102            except KeyError:        #no more peaks to process
1103                return yb+yc
1104    else:
1105        Pdabc = parmDict['Pdabc']
1106        difC = parmDict['difC']
1107        iPeak = 0
1108        while True:
1109            try:
1110                pos = parmDict['pos'+str(iPeak)]               
1111                tof = pos-parmDict['Zero']
1112                dsp = tof/difC
1113                intens = parmDict['int'+str(iPeak)]
1114                alpName = 'alp'+str(iPeak)
1115                if alpName in varyList:
1116                    alp = parmDict[alpName]
1117                else:
1118                    if len(Pdabc):
1119                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1120                    else:
1121                        alp = G2mth.getTOFalpha(parmDict,dsp)
1122                alp = max(0.0001,alp)
1123                betName = 'bet'+str(iPeak)
1124                if betName in varyList:
1125                    bet = parmDict[betName]
1126                else:
1127                    if len(Pdabc):
1128                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1129                    else:
1130                        bet = G2mth.getTOFbeta(parmDict,dsp)
1131                bet = max(0.0001,bet)
1132                sigName = 'sig'+str(iPeak)
1133                if sigName in varyList:
1134                    sig = parmDict[sigName]
1135                else:
1136                    sig = G2mth.getTOFsig(parmDict,dsp)
1137                gamName = 'gam'+str(iPeak)
1138                if gamName in varyList:
1139                    gam = parmDict[gamName]
1140                else:
1141                    gam = G2mth.getTOFgamma(parmDict,dsp)
1142                gam = max(gam,0.001)             #avoid neg gamma
1143                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1144                iBeg = np.searchsorted(xdata,pos-fmin)
1145                iFin = np.searchsorted(xdata,pos+fmax)
1146                lenX = len(xdata)
1147                if not iBeg:
1148                    iFin = np.searchsorted(xdata,pos+fmax)
1149                elif iBeg == lenX:
1150                    iFin = iBeg
1151                else:
1152                    iFin = np.searchsorted(xdata,pos+fmax)
1153                if not iBeg+iFin:       #peak below low limit
1154                    iPeak += 1
1155                    continue
1156                elif not iBeg-iFin:     #peak above high limit
1157                    return yb+yc
1158                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1159                iPeak += 1
1160            except KeyError:        #no more peaks to process
1161                return yb+yc
1162           
1163def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
1164    'needs a doc string'
1165# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1166    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1167    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,dataType,xdata)
1168    if 'Back;0' in varyList:            #background derivs are in front if present
1169        dMdv[0:len(dMdb)] = dMdb
1170    names = ['DebyeA','DebyeR','DebyeU']
1171    for name in varyList:
1172        if 'Debye' in name:
1173            parm,id = name.split(';')
1174            ip = names.index(parm)
1175            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
1176    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1177    for name in varyList:
1178        if 'BkPk' in name:
1179            parm,id = name.split(';')
1180            ip = names.index(parm)
1181            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
1182    cw = np.diff(xdata)
1183    cw = np.append(cw,cw[-1])
1184    if 'C' in dataType:
1185        shl = max(parmDict['SH/L'],0.002)
1186        Ka2 = False
1187        if 'Lam1' in parmDict.keys():
1188            Ka2 = True
1189            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1190            kRatio = parmDict['I(L2)/I(L1)']
1191        iPeak = 0
1192        while True:
1193            try:
1194                pos = parmDict['pos'+str(iPeak)]
1195                tth = (pos-parmDict['Zero'])
1196                intens = parmDict['int'+str(iPeak)]
1197                sigName = 'sig'+str(iPeak)
1198                if sigName in varyList:
1199                    sig = parmDict[sigName]
1200                    dsdU = dsdV = dsdW = 0
1201                else:
1202                    sig = G2mth.getCWsig(parmDict,tth)
1203                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1204                sig = max(sig,0.001)          #avoid neg sigma
1205                gamName = 'gam'+str(iPeak)
1206                if gamName in varyList:
1207                    gam = parmDict[gamName]
1208                    dgdX = dgdY = 0
1209                else:
1210                    gam = G2mth.getCWgam(parmDict,tth)
1211                    dgdX,dgdY = G2mth.getCWgamDeriv(tth)
1212                gam = max(gam,0.001)             #avoid neg gamma
1213                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1214                iBeg = np.searchsorted(xdata,pos-fmin)
1215                iFin = np.searchsorted(xdata,pos+fmin)
1216                if not iBeg+iFin:       #peak below low limit
1217                    iPeak += 1
1218                    continue
1219                elif not iBeg-iFin:     #peak above high limit
1220                    break
1221                dMdpk = np.zeros(shape=(6,len(xdata)))
1222                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1223                for i in range(1,5):
1224                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
1225                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
1226                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1227                if Ka2:
1228                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1229                    iBeg = np.searchsorted(xdata,pos2-fmin)
1230                    iFin = np.searchsorted(xdata,pos2+fmin)
1231                    if iBeg-iFin:
1232                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1233                        for i in range(1,5):
1234                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
1235                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
1236                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
1237                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1238                for parmName in ['pos','int','sig','gam']:
1239                    try:
1240                        idx = varyList.index(parmName+str(iPeak))
1241                        dMdv[idx] = dervDict[parmName]
1242                    except ValueError:
1243                        pass
1244                if 'U' in varyList:
1245                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1246                if 'V' in varyList:
1247                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1248                if 'W' in varyList:
1249                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1250                if 'X' in varyList:
1251                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1252                if 'Y' in varyList:
1253                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1254                if 'SH/L' in varyList:
1255                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
1256                if 'I(L2)/I(L1)' in varyList:
1257                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1258                iPeak += 1
1259            except KeyError:        #no more peaks to process
1260                break
1261    else:
1262        Pdabc = parmDict['Pdabc']
1263        difC = parmDict['difC']
1264        iPeak = 0
1265        while True:
1266            try:
1267                pos = parmDict['pos'+str(iPeak)]               
1268                tof = pos-parmDict['Zero']
1269                dsp = tof/difC
1270                intens = parmDict['int'+str(iPeak)]
1271                alpName = 'alp'+str(iPeak)
1272                if alpName in varyList:
1273                    alp = parmDict[alpName]
1274                else:
1275                    if len(Pdabc):
1276                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1277                        dada0 = 0
1278                    else:
1279                        alp = G2mth.getTOFalpha(parmDict,dsp)
1280                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1281                betName = 'bet'+str(iPeak)
1282                if betName in varyList:
1283                    bet = parmDict[betName]
1284                else:
1285                    if len(Pdabc):
1286                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1287                        dbdb0 = dbdb1 = dbdb2 = 0
1288                    else:
1289                        bet = G2mth.getTOFbeta(parmDict,dsp)
1290                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1291                sigName = 'sig'+str(iPeak)
1292                if sigName in varyList:
1293                    sig = parmDict[sigName]
1294                    dsds0 = dsds1 = dsds2 = dsds3 = 0
1295                else:
1296                    sig = G2mth.getTOFsig(parmDict,dsp)
1297                    dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1298                gamName = 'gam'+str(iPeak)
1299                if gamName in varyList:
1300                    gam = parmDict[gamName]
1301                    dsdX = dsdY = 0
1302                else:
1303                    gam = G2mth.getTOFgamma(parmDict,dsp)
1304                    dsdX,dsdY = G2mth.getTOFgammaDeriv(dsp)
1305                gam = max(gam,0.001)             #avoid neg gamma
1306                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1307                iBeg = np.searchsorted(xdata,pos-fmin)
1308                lenX = len(xdata)
1309                if not iBeg:
1310                    iFin = np.searchsorted(xdata,pos+fmax)
1311                elif iBeg == lenX:
1312                    iFin = iBeg
1313                else:
1314                    iFin = np.searchsorted(xdata,pos+fmax)
1315                if not iBeg+iFin:       #peak below low limit
1316                    iPeak += 1
1317                    continue
1318                elif not iBeg-iFin:     #peak above high limit
1319                    break
1320                dMdpk = np.zeros(shape=(7,len(xdata)))
1321                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1322                for i in range(1,6):
1323                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1324                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1325                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1326                for parmName in ['pos','int','alp','bet','sig','gam']:
1327                    try:
1328                        idx = varyList.index(parmName+str(iPeak))
1329                        dMdv[idx] = dervDict[parmName]
1330                    except ValueError:
1331                        pass
1332                if 'alpha' in varyList:
1333                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1334                if 'beta-0' in varyList:
1335                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1336                if 'beta-1' in varyList:
1337                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1338                if 'beta-q' in varyList:
1339                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1340                if 'sig-0' in varyList:
1341                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1342                if 'sig-1' in varyList:
1343                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1344                if 'sig-2' in varyList:
1345                    dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1346                if 'sig-q' in varyList:
1347                    dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1348                if 'X' in varyList:
1349                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1350                if 'Y' in varyList:
1351                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
1352                iPeak += 1
1353            except KeyError:        #no more peaks to process
1354                break
1355    return dMdv
1356       
1357def Dict2Values(parmdict, varylist):
1358    '''Use before call to leastsq to setup list of values for the parameters
1359    in parmdict, as selected by key in varylist'''
1360    return [parmdict[key] for key in varylist] 
1361   
1362def Values2Dict(parmdict, varylist, values):
1363    ''' Use after call to leastsq to update the parameter dictionary with
1364    values corresponding to keys in varylist'''
1365    parmdict.update(zip(varylist,values))
1366   
1367def SetBackgroundParms(Background):
1368    'needs a doc string'
1369    if len(Background) == 1:            # fix up old backgrounds
1370        Background.append({'nDebye':0,'debyeTerms':[]})
1371    bakType,bakFlag = Background[0][:2]
1372    backVals = Background[0][3:]
1373    backNames = ['Back;'+str(i) for i in range(len(backVals))]
1374    Debye = Background[1]           #also has background peaks stuff
1375    backDict = dict(zip(backNames,backVals))
1376    backVary = []
1377    if bakFlag:
1378        backVary = backNames
1379
1380    backDict['nDebye'] = Debye['nDebye']
1381    debyeDict = {}
1382    debyeList = []
1383    for i in range(Debye['nDebye']):
1384        debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1385        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1386        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1387    debyeVary = []
1388    for item in debyeList:
1389        if item[1]:
1390            debyeVary.append(item[0])
1391    backDict.update(debyeDict)
1392    backVary += debyeVary
1393
1394    backDict['nPeaks'] = Debye['nPeaks']
1395    peaksDict = {}
1396    peaksList = []
1397    for i in range(Debye['nPeaks']):
1398        peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1399        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1400        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1401    peaksVary = []
1402    for item in peaksList:
1403        if item[1]:
1404            peaksVary.append(item[0])
1405    backDict.update(peaksDict)
1406    backVary += peaksVary   
1407    return bakType,backDict,backVary
1408   
1409def DoCalibInst(IndexPeaks,Inst):
1410   
1411    def SetInstParms():
1412        dataType = Inst['Type'][0]
1413        insVary = []
1414        insNames = []
1415        insVals = []
1416        for parm in Inst:
1417            insNames.append(parm)
1418            insVals.append(Inst[parm][1])
1419            if parm in ['Lam','difC','difA','difB','Zero',]:
1420                if Inst[parm][2]:
1421                    insVary.append(parm)
1422        instDict = dict(zip(insNames,insVals))
1423        return dataType,instDict,insVary
1424       
1425    def GetInstParms(parmDict,Inst,varyList):
1426        for name in Inst:
1427            Inst[name][1] = parmDict[name]
1428       
1429    def InstPrint(Inst,sigDict):
1430        print 'Instrument Parameters:'
1431        if 'C' in Inst['Type'][0]:
1432            ptfmt = "%12.6f"
1433        else:
1434            ptfmt = "%12.3f"
1435        ptlbls = 'names :'
1436        ptstr =  'values:'
1437        sigstr = 'esds  :'
1438        for parm in Inst:
1439            if parm in  ['Lam','difC','difA','difB','Zero',]:
1440                ptlbls += "%s" % (parm.center(12))
1441                ptstr += ptfmt % (Inst[parm][1])
1442                if parm in sigDict:
1443                    sigstr += ptfmt % (sigDict[parm])
1444                else:
1445                    sigstr += 12*' '
1446        print ptlbls
1447        print ptstr
1448        print sigstr
1449       
1450    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1451        parmDict.update(zip(varyList,values))
1452        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1453
1454    peakPos = []
1455    peakDsp = []
1456    peakWt = []
1457    for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1458        if peak[2] and peak[3] and sig > 0.:
1459            peakPos.append(peak[0])
1460            peakDsp.append(peak[-1])    #d-calc
1461            peakWt.append(peak[-1]**2/sig**2)   #weight by d**2
1462    peakPos = np.array(peakPos)
1463    peakDsp = np.array(peakDsp)
1464    peakWt = np.array(peakWt)
1465    dataType,insDict,insVary = SetInstParms()
1466    parmDict = {}
1467    parmDict.update(insDict)
1468    varyList = insVary
1469    if not len(varyList):
1470        print '**** ERROR - nothing to refine! ****'
1471        return False
1472    while True:
1473        begin = time.time()
1474        values =  np.array(Dict2Values(parmDict, varyList))
1475        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1476            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1477        ncyc = int(result[2]['nfev']/2)
1478        runtime = time.time()-begin   
1479        chisq = np.sum(result[2]['fvec']**2)
1480        Values2Dict(parmDict, varyList, result[0])
1481        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
1482        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(peakPos),' Number of parameters: ',len(varyList)
1483        print 'calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1484        print 'chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)
1485        try:
1486            sig = np.sqrt(np.diag(result[1])*GOF)
1487            if np.any(np.isnan(sig)):
1488                print '*** Least squares aborted - some invalid esds possible ***'
1489            break                   #refinement succeeded - finish up!
1490        except ValueError:          #result[1] is None on singular matrix
1491            print '**** Refinement failed - singular matrix ****'
1492       
1493    sigDict = dict(zip(varyList,sig))
1494    GetInstParms(parmDict,Inst,varyList)
1495    InstPrint(Inst,sigDict)
1496    return True
1497           
1498def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
1499    '''Called to perform a peak fit, refining the selected items in the peak
1500    table as well as selected items in the background.
1501
1502    :param str FitPgm: type of fit to perform. At present this is ignored.
1503    :param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1504      four values followed by a refine flag where the values are: position, intensity,
1505      sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1506      tree entry, dict item "peaks"
1507    :param list Background: describes the background. List with two items.
1508      Item 0 specifies a background model and coefficients. Item 1 is a dict.
1509      From the Histogram/Background tree entry.
1510    :param list Limits: min and max x-value to use
1511    :param dict Inst: Instrument parameters
1512    :param dict Inst2: more Instrument parameters
1513    :param numpy.array data: a 5xn array. data[0] is the x-values,
1514      data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1515      calc, background and difference intensities, respectively.
1516    :param list prevVaryList: Used in sequential refinements to override the
1517      variable list. Defaults as an empty list.
1518    :param bool oneCycle: True if only one cycle of fitting should be performed
1519    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1520      and derivType = controls['deriv type']. If None default values are used.
1521    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1522      Defaults to None, which means no updates are done.
1523    '''
1524    def GetBackgroundParms(parmList,Background):
1525        iBak = 0
1526        while True:
1527            try:
1528                bakName = 'Back;'+str(iBak)
1529                Background[0][iBak+3] = parmList[bakName]
1530                iBak += 1
1531            except KeyError:
1532                break
1533        iDb = 0
1534        while True:
1535            names = ['DebyeA;','DebyeR;','DebyeU;']
1536            try:
1537                for i,name in enumerate(names):
1538                    val = parmList[name+str(iDb)]
1539                    Background[1]['debyeTerms'][iDb][2*i] = val
1540                iDb += 1
1541            except KeyError:
1542                break
1543        iDb = 0
1544        while True:
1545            names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1546            try:
1547                for i,name in enumerate(names):
1548                    val = parmList[name+str(iDb)]
1549                    Background[1]['peaksList'][iDb][2*i] = val
1550                iDb += 1
1551            except KeyError:
1552                break
1553               
1554    def BackgroundPrint(Background,sigDict):
1555        print 'Background coefficients for',Background[0][0],'function'
1556        ptfmt = "%12.5f"
1557        ptstr =  'value: '
1558        sigstr = 'esd  : '
1559        for i,back in enumerate(Background[0][3:]):
1560            ptstr += ptfmt % (back)
1561            if Background[0][1]:
1562                prm = 'Back;'+str(i)
1563                if prm in sigDict:
1564                    sigstr += ptfmt % (sigDict[prm])
1565                else:
1566                    sigstr += " "*12
1567            if len(ptstr) > 75:
1568                print ptstr
1569                if Background[0][1]: print sigstr
1570                ptstr =  'value: '
1571                sigstr = 'esd  : '
1572        if len(ptstr) > 8:
1573            print ptstr
1574            if Background[0][1]: print sigstr
1575
1576        if Background[1]['nDebye']:
1577            parms = ['DebyeA;','DebyeR;','DebyeU;']
1578            print 'Debye diffuse scattering coefficients'
1579            ptfmt = "%12.5f"
1580            print ' term       DebyeA       esd        DebyeR       esd        DebyeU        esd'
1581            for term in range(Background[1]['nDebye']):
1582                line = ' term %d'%(term)
1583                for ip,name in enumerate(parms):
1584                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1585                    if name+str(term) in sigDict:
1586                        line += ptfmt%(sigDict[name+str(term)])
1587                    else:
1588                        line += " "*12
1589                print line
1590        if Background[1]['nPeaks']:
1591            print 'Coefficients for Background Peaks'
1592            ptfmt = "%15.3f"
1593            for j,pl in enumerate(Background[1]['peaksList']):
1594                names =  'peak %3d:'%(j+1)
1595                ptstr =  'values  :'
1596                sigstr = 'esds    :'
1597                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1598                    val = pl[2*i]
1599                    prm = lbl+";"+str(j)
1600                    names += '%15s'%(prm)
1601                    ptstr += ptfmt%(val)
1602                    if prm in sigDict:
1603                        sigstr += ptfmt%(sigDict[prm])
1604                    else:
1605                        sigstr += " "*15
1606                print names
1607                print ptstr
1608                print sigstr
1609                           
1610    def SetInstParms(Inst):
1611        dataType = Inst['Type'][0]
1612        insVary = []
1613        insNames = []
1614        insVals = []
1615        for parm in Inst:
1616            insNames.append(parm)
1617            insVals.append(Inst[parm][1])
1618            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1619                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1620                    insVary.append(parm)
1621        instDict = dict(zip(insNames,insVals))
1622        instDict['X'] = max(instDict['X'],0.01)
1623        instDict['Y'] = max(instDict['Y'],0.01)
1624        if 'SH/L' in instDict:
1625            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1626        return dataType,instDict,insVary
1627       
1628    def GetInstParms(parmDict,Inst,varyList,Peaks):
1629        for name in Inst:
1630            Inst[name][1] = parmDict[name]
1631        iPeak = 0
1632        while True:
1633            try:
1634                sigName = 'sig'+str(iPeak)
1635                pos = parmDict['pos'+str(iPeak)]
1636                if sigName not in varyList:
1637                    if 'C' in Inst['Type'][0]:
1638                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1639                    else:
1640                        dsp = G2lat.Pos2dsp(Inst,pos)
1641                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1642                gamName = 'gam'+str(iPeak)
1643                if gamName not in varyList:
1644                    if 'C' in Inst['Type'][0]:
1645                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1646                    else:
1647                        dsp = G2lat.Pos2dsp(Inst,pos)
1648                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1649                iPeak += 1
1650            except KeyError:
1651                break
1652       
1653    def InstPrint(Inst,sigDict):
1654        print 'Instrument Parameters:'
1655        ptfmt = "%12.6f"
1656        ptlbls = 'names :'
1657        ptstr =  'values:'
1658        sigstr = 'esds  :'
1659        for parm in Inst:
1660            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1661                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1662                ptlbls += "%s" % (parm.center(12))
1663                ptstr += ptfmt % (Inst[parm][1])
1664                if parm in sigDict:
1665                    sigstr += ptfmt % (sigDict[parm])
1666                else:
1667                    sigstr += 12*' '
1668        print ptlbls
1669        print ptstr
1670        print sigstr
1671
1672    def SetPeaksParms(dataType,Peaks):
1673        peakNames = []
1674        peakVary = []
1675        peakVals = []
1676        if 'C' in dataType:
1677            names = ['pos','int','sig','gam']
1678        else:
1679            names = ['pos','int','alp','bet','sig','gam']
1680        for i,peak in enumerate(Peaks):
1681            for j,name in enumerate(names):
1682                peakVals.append(peak[2*j])
1683                parName = name+str(i)
1684                peakNames.append(parName)
1685                if peak[2*j+1]:
1686                    peakVary.append(parName)
1687        return dict(zip(peakNames,peakVals)),peakVary
1688               
1689    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1690        if 'C' in Inst['Type'][0]:
1691            names = ['pos','int','sig','gam']
1692        else:   #'T'
1693            names = ['pos','int','alp','bet','sig','gam']
1694        for i,peak in enumerate(Peaks):
1695            pos = parmDict['pos'+str(i)]
1696            if 'difC' in Inst:
1697                dsp = pos/Inst['difC'][1]
1698            for j in range(len(names)):
1699                parName = names[j]+str(i)
1700                if parName in varyList:
1701                    peak[2*j] = parmDict[parName]
1702                elif 'alpha' in parName:
1703                    peak[2*j] = parmDict['alpha']/dsp
1704                elif 'beta' in parName:
1705                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1706                elif 'sig' in parName:
1707                    if 'C' in Inst['Type'][0]:
1708                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1709                    else:
1710                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1711                elif 'gam' in parName:
1712                    if 'C' in Inst['Type'][0]:
1713                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1714                    else:
1715                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1716                       
1717    def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
1718        print 'Peak coefficients:'
1719        if 'C' in dataType:
1720            names = ['pos','int','sig','gam']
1721        else:   #'T'
1722            names = ['pos','int','alp','bet','sig','gam']           
1723        head = 13*' '
1724        for name in names:
1725            if name in ['alp','bet']:
1726                head += name.center(8)+'esd'.center(8)
1727            else:
1728                head += name.center(10)+'esd'.center(10)
1729        head += 'bins'.center(8)
1730        print head
1731        if 'C' in dataType:
1732            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1733        else:
1734            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1735        for i,peak in enumerate(Peaks):
1736            ptstr =  ':'
1737            for j in range(len(names)):
1738                name = names[j]
1739                parName = name+str(i)
1740                ptstr += ptfmt[name] % (parmDict[parName])
1741                if parName in varyList:
1742                    ptstr += ptfmt[name] % (sigDict[parName])
1743                else:
1744                    if name in ['alp','bet']:
1745                        ptstr += 8*' '
1746                    else:
1747                        ptstr += 10*' '
1748            ptstr += '%9.2f'%(ptsperFW[i])
1749            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1750               
1751    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1752        parmdict.update(zip(varylist,values))
1753        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1754           
1755    def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):       
1756        parmdict.update(zip(varylist,values))
1757        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1758        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1759        if dlg:
1760            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1761            if not GoOn:
1762                return -M           #abort!!
1763        return M
1764       
1765    if controls:
1766        Ftol = controls['min dM/M']
1767    else:
1768        Ftol = 0.0001
1769    if oneCycle:
1770        Ftol = 1.0
1771    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1772    yc *= 0.                            #set calcd ones to zero
1773    yb *= 0.
1774    yd *= 0.
1775    cw = x[1:]-x[:-1]
1776    xBeg = np.searchsorted(x,Limits[0])
1777    xFin = np.searchsorted(x,Limits[1])+1
1778    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1779    dataType,insDict,insVary = SetInstParms(Inst)
1780    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1781    parmDict = {}
1782    parmDict.update(bakDict)
1783    parmDict.update(insDict)
1784    parmDict.update(peakDict)
1785    parmDict['Pdabc'] = []      #dummy Pdabc
1786    parmDict.update(Inst2)      #put in real one if there
1787    if prevVaryList:
1788        varyList = prevVaryList[:]
1789    else:
1790        varyList = bakVary+insVary+peakVary
1791    fullvaryList = varyList[:]
1792    while True:
1793        begin = time.time()
1794        values =  np.array(Dict2Values(parmDict, varyList))
1795        Rvals = {}
1796        badVary = []
1797        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1798               args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1799        ncyc = int(result[2]['nfev']/2)
1800        runtime = time.time()-begin   
1801        chisq = np.sum(result[2]['fvec']**2)
1802        Values2Dict(parmDict, varyList, result[0])
1803        Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1804        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1805        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1806        if ncyc:
1807            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1808        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1809        sig = [0]*len(varyList)
1810        if len(varyList) == 0: break  # if nothing was refined
1811        try:
1812            sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1813            if np.any(np.isnan(sig)):
1814                print '*** Least squares aborted - some invalid esds possible ***'
1815            break                   #refinement succeeded - finish up!
1816        except ValueError:          #result[1] is None on singular matrix
1817            print '**** Refinement failed - singular matrix ****'
1818            Ipvt = result[2]['ipvt']
1819            for i,ipvt in enumerate(Ipvt):
1820                if not np.sum(result[2]['fjac'],axis=1)[i]:
1821                    print 'Removing parameter: ',varyList[ipvt-1]
1822                    badVary.append(varyList[ipvt-1])
1823                    del(varyList[ipvt-1])
1824                    break
1825            else: # nothing removed
1826                break
1827    if dlg: dlg.Destroy()
1828    sigDict = dict(zip(varyList,sig))
1829    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]
1830    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1831    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1832    GetBackgroundParms(parmDict,Background)
1833    if bakVary: BackgroundPrint(Background,sigDict)
1834    GetInstParms(parmDict,Inst,varyList,Peaks)
1835    if insVary: InstPrint(Inst,sigDict)
1836    GetPeaksParms(Inst,parmDict,Peaks,varyList)
1837    binsperFWHM = []
1838    for peak in Peaks:
1839        FWHM = getFWHM(peak[0],Inst)
1840        try:
1841            binsperFWHM.append(FWHM/cw[x.searchsorted(peak[0])])
1842        except IndexError:
1843            binsperFWHM.append(0.)
1844    if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
1845    if len(binsperFWHM):
1846        if min(binsperFWHM) < 1.:
1847            print '*** Warning: calculated peak widths are too narrow to refine profile coefficients ***'
1848            if 'T' in Inst['Type'][0]:
1849                print ' Manually increase sig-0, 1, or 2 in Instrument Parameters'
1850            else:
1851                print ' Manually increase W in Instrument Parameters'
1852        elif min(binsperFWHM) < 4.:
1853            print '*** Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement ***'
1854            print '*** recommended is 6-10; you have %.2f ***'%(min(binsperFWHM))
1855    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
1856   
1857def calcIncident(Iparm,xdata):
1858    'needs a doc string'
1859
1860    def IfunAdv(Iparm,xdata):
1861        Itype = Iparm['Itype']
1862        Icoef = Iparm['Icoeff']
1863        DYI = np.ones((12,xdata.shape[0]))
1864        YI = np.ones_like(xdata)*Icoef[0]
1865       
1866        x = xdata/1000.                 #expressions are in ms
1867        if Itype == 'Exponential':
1868            for i in [1,3,5,7,9]:
1869                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1870                YI += Icoef[i]*Eterm
1871                DYI[i] *= Eterm
1872                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
1873        elif 'Maxwell'in Itype:
1874            Eterm = np.exp(-Icoef[2]/x**2)
1875            DYI[1] = Eterm/x**5
1876            DYI[2] = -Icoef[1]*DYI[1]/x**2
1877            YI += (Icoef[1]*Eterm/x**5)
1878            if 'Exponential' in Itype:
1879                for i in range(3,11,2):
1880                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1881                    YI += Icoef[i]*Eterm
1882                    DYI[i] *= Eterm
1883                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
1884            else:   #Chebyschev
1885                T = (2./x)-1.
1886                Ccof = np.ones((12,xdata.shape[0]))
1887                Ccof[1] = T
1888                for i in range(2,12):
1889                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1890                for i in range(1,10):
1891                    YI += Ccof[i]*Icoef[i+2]
1892                    DYI[i+2] =Ccof[i]
1893        return YI,DYI
1894       
1895    Iesd = np.array(Iparm['Iesd'])
1896    Icovar = Iparm['Icovar']
1897    YI,DYI = IfunAdv(Iparm,xdata)
1898    YI = np.where(YI>0,YI,1.)
1899    WYI = np.zeros_like(xdata)
1900    vcov = np.zeros((12,12))
1901    k = 0
1902    for i in range(12):
1903        for j in range(i,12):
1904            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1905            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1906            k += 1
1907    M = np.inner(vcov,DYI.T)
1908    WYI = np.sum(M*DYI,axis=0)
1909    WYI = np.where(WYI>0.,WYI,0.)
1910    return YI,WYI
1911   
1912################################################################################
1913# Reflectometry calculations
1914################################################################################
1915
1916def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
1917    print 'fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0])
1918   
1919    class RandomDisplacementBounds(object):
1920        """random displacement with bounds"""
1921        def __init__(self, xmin, xmax, stepsize=0.5):
1922            self.xmin = xmin
1923            self.xmax = xmax
1924            self.stepsize = stepsize
1925   
1926        def __call__(self, x):
1927            """take a random step but ensure the new position is within the bounds"""
1928            while True:
1929                # this could be done in a much more clever way, but it will work for example purposes
1930                steps = self.xmax-self.xmin
1931                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
1932                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
1933                    break
1934            return xnew
1935   
1936    def GetModelParms():
1937        parmDict = {}
1938        varyList = []
1939        values = []
1940        bounds = []
1941        parmDict['dQ type'] = data['dQ type']
1942        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
1943        for parm in ['Scale','FltBack']:
1944            parmDict[parm] = data[parm][0]
1945            if data[parm][1]:
1946                varyList.append(parm)
1947                values.append(data[parm][0])
1948                bounds.append(Bounds[parm])
1949        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
1950        parmDict['nLayers'] = len(parmDict['Layer Seq'])
1951        for ilay,layer in enumerate(data['Layers']):
1952            name = layer['Name']
1953            cid = str(ilay)+';'
1954            parmDict[cid+'Name'] = name
1955            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
1956                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
1957                if layer.get(parm,[0.,False])[1]:
1958                    varyList.append(cid+parm)
1959                    value = layer[parm][0]
1960                    values.append(value)
1961                    if value:
1962                        bound = [value*Bfac,value/Bfac]
1963                    else:
1964                        bound = [0.,10.]
1965                    bounds.append(bound)
1966            if name not in ['vacuum','unit scatter']:
1967                parmDict[cid+'rho'] = Substances[name]['Scatt density']
1968                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
1969        return parmDict,varyList,values,bounds
1970   
1971    def SetModelParms():
1972        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
1973        if 'Scale' in varyList:
1974            data['Scale'][0] = parmDict['Scale']
1975            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
1976        print line
1977        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
1978        if 'FltBack' in varyList:
1979            data['FltBack'][0] = parmDict['FltBack']
1980            line += ' esd: %15.3g'%(sigDict['FltBack'])
1981        print line
1982        for ilay,layer in enumerate(data['Layers']):
1983            name = layer['Name']
1984            print ' Parameters for layer: %d %s'%(ilay,name)
1985            cid = str(ilay)+';'
1986            line = ' '
1987            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
1988            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
1989            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
1990                if parm in layer:
1991                    if parm == 'Rough':
1992                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
1993                    else:
1994                        layer[parm][0] = parmDict[cid+parm]
1995                    line += ' %s: %.3f'%(parm,layer[parm][0])
1996                    if cid+parm in varyList:
1997                        line += ' esd: %.3g'%(sigDict[cid+parm])
1998            print line
1999            print line2
2000   
2001    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2002        parmDict.update(zip(varyList,values))
2003        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2004        return M
2005   
2006    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2007        parmDict.update(zip(varyList,values))
2008        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2009        return np.sum(M**2)
2010   
2011    def getREFD(Q,Qsig,parmDict):
2012        Ic = np.ones_like(Q)*parmDict['FltBack']
2013        Scale = parmDict['Scale']
2014        Nlayers = parmDict['nLayers']
2015        Res = parmDict['Res']
2016        depth = np.zeros(Nlayers)
2017        rho = np.zeros(Nlayers)
2018        irho = np.zeros(Nlayers)
2019        sigma = np.zeros(Nlayers)
2020        for ilay,lay in enumerate(parmDict['Layer Seq']):
2021            cid = str(lay)+';'
2022            depth[ilay] = parmDict[cid+'Thick']
2023            sigma[ilay] = parmDict[cid+'Rough']
2024            if parmDict[cid+'Name'] == u'unit scatter':
2025                rho[ilay] = parmDict[cid+'DenMul']
2026                irho[ilay] = parmDict[cid+'iDenMul']
2027            elif 'vacuum' != parmDict[cid+'Name']:
2028                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
2029                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
2030            if cid+'Mag SLD' in parmDict:
2031                rho[ilay] += parmDict[cid+'Mag SLD']
2032        if parmDict['dQ type'] == 'None':
2033            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2034        elif 'const' in parmDict['dQ type']:
2035            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
2036        else:       #dQ/Q in data
2037            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
2038        Ic += AB*Scale
2039        return Ic
2040       
2041    def estimateT0(takestep):
2042        Mmax = -1.e-10
2043        Mmin = 1.e10
2044        for i in range(100):
2045            x0 = takestep(values)
2046            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2047            Mmin = min(M,Mmin)
2048            MMax = max(M,Mmax)
2049        return 1.5*(MMax-Mmin)
2050
2051    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2052    if data.get('2% weight'):
2053        wt = 1./(0.02*Io)**2
2054    Qmin = Limits[1][0]
2055    Qmax = Limits[1][1]
2056    wtFactor = ProfDict['wtFactor']
2057    Bfac = data['Toler']
2058    Ibeg = np.searchsorted(Q,Qmin)
2059    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
2060    Ic[:] = 0
2061    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
2062              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
2063    parmDict,varyList,values,bounds = GetModelParms()
2064    Msg = 'Failed to converge'
2065    if varyList:
2066        if data['Minimizer'] == 'LMLS': 
2067            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
2068                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2069            parmDict.update(zip(varyList,result[0]))
2070            chisq = np.sum(result[2]['fvec']**2)
2071            ncalc = result[2]['nfev']
2072            covM = result[1]
2073            newVals = result[0]
2074        elif data['Minimizer'] == 'Basin Hopping':
2075            xyrng = np.array(bounds).T
2076            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
2077            T0 = estimateT0(take_step)
2078            print ' Estimated temperature: %.3g'%(T0)
2079            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
2080                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
2081                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
2082            chisq = result.fun
2083            ncalc = result.nfev
2084            newVals = result.x
2085            covM = []
2086        elif data['Minimizer'] == 'MC/SA Anneal':
2087            xyrng = np.array(bounds).T
2088            result = G2mth.anneal(sumREFD, values, 
2089                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
2090                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
2091                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
2092                ranRange=0.20,autoRan=False,dlg=None)
2093            newVals = result[0]
2094            parmDict.update(zip(varyList,newVals))
2095            chisq = result[1]
2096            ncalc = result[3]
2097            covM = []
2098            print ' MC/SA final temperature: %.4g'%(result[2])
2099        elif data['Minimizer'] == 'L-BFGS-B':
2100            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
2101                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2102            parmDict.update(zip(varyList,result['x']))
2103            chisq = result.fun
2104            ncalc = result.nfev
2105            newVals = result.x
2106            covM = []
2107    else:   #nothing varied
2108        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2109        chisq = np.sum(M**2)
2110        ncalc = 0
2111        covM = []
2112        sig = []
2113        sigDict = {}
2114        result = []
2115    Rvals = {}
2116    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
2117    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
2118    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
2119    Ib[Ibeg:Ifin] = parmDict['FltBack']
2120    try:
2121        if not len(varyList):
2122            Msg += ' - nothing refined'
2123            raise ValueError
2124        Nans = np.isnan(newVals)
2125        if np.any(Nans):
2126            name = varyList[Nans.nonzero(True)[0]]
2127            Msg += ' Nan result for '+name+'!'
2128            raise ValueError
2129        Negs = np.less_equal(newVals,0.)
2130        if np.any(Negs):
2131            indx = Negs.nonzero()
2132            name = varyList[indx[0][0]]
2133            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
2134                Msg += ' negative coefficient for '+name+'!'
2135                raise ValueError
2136        if len(covM):
2137            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
2138            covMatrix = covM*Rvals['GOF']
2139        else:
2140            sig = np.zeros(len(varyList))
2141            covMatrix = []
2142        sigDict = dict(zip(varyList,sig))
2143        print ' Results of reflectometry data modelling fit:'
2144        print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
2145        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
2146        SetModelParms()
2147        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
2148    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
2149        print Msg
2150        return False,0,0,0,0,0,0,Msg
2151       
2152def makeSLDprofile(data,Substances):
2153   
2154    sq2 = np.sqrt(2.)
2155    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2156    Nlayers = len(laySeq)
2157    laySeq = np.array(laySeq,dtype=int)
2158    interfaces = np.zeros(Nlayers)
2159    rho = np.zeros(Nlayers)
2160    sigma = np.zeros(Nlayers)
2161    name = data['Layers'][0]['Name']
2162    thick = 0.
2163    for ilay,lay in enumerate(laySeq):
2164        layer = data['Layers'][lay]
2165        name = layer['Name']
2166        if 'Thick' in layer:
2167            thick += layer['Thick'][0]
2168            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
2169        if 'Rough' in layer:
2170            sigma[ilay] = max(0.001,layer['Rough'][0])
2171        if name != 'vacuum':
2172            if name == 'unit scatter':
2173                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
2174            else:
2175                rrho = Substances[name]['Scatt density']
2176                irho = Substances[name]['XImag density']
2177                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
2178        if 'Mag SLD' in layer:
2179            rho[ilay] += layer['Mag SLD'][0]
2180    name = data['Layers'][-1]['Name']
2181    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
2182    xr = np.flipud(x)
2183    interfaces[-1] = x[-1]
2184    y = np.ones_like(x)*rho[0]
2185    iBeg = 0
2186    for ilayer in range(Nlayers-1):
2187        delt = rho[ilayer+1]-rho[ilayer]
2188        iPos = np.searchsorted(x,interfaces[ilayer])
2189        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
2190        iBeg = iPos
2191    return x,xr,y   
2192
2193def REFDModelFxn(Profile,Inst,Limits,Substances,data):
2194   
2195    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2196    Qmin = Limits[1][0]
2197    Qmax = Limits[1][1]
2198    iBeg = np.searchsorted(Q,Qmin)
2199    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2200    Ib[:] = data['FltBack'][0]
2201    Ic[:] = 0
2202    Scale = data['Scale'][0]
2203    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2204    Nlayers = len(laySeq)
2205    depth = np.zeros(Nlayers)
2206    rho = np.zeros(Nlayers)
2207    irho = np.zeros(Nlayers)
2208    sigma = np.zeros(Nlayers)
2209    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
2210        layer = data['Layers'][lay]
2211        name = layer['Name']
2212        if 'Thick' in layer:    #skips first & last layers
2213            depth[ilay] = layer['Thick'][0]
2214        if 'Rough' in layer:    #skips first layer
2215            sigma[ilay] = layer['Rough'][0]
2216        if 'unit scatter' == name:
2217            rho[ilay] = layer['DenMul'][0]
2218            irho[ilay] = layer['iDenMul'][0]
2219        else:
2220            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2221            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
2222        if 'Mag SLD' in layer:
2223            rho[ilay] += layer['Mag SLD'][0]
2224    if data['dQ type'] == 'None':
2225        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2226    elif 'const' in data['dQ type']:
2227        res = data['Resolution'][0]/(100.*sateln2)
2228        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
2229    else:       #dQ/Q in data
2230        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
2231    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
2232
2233def abeles(kz, depth, rho, irho=0, sigma=0):
2234    """
2235    Optical matrix form of the reflectivity calculation.
2236    O.S. Heavens, Optical Properties of Thin Solid Films
2237   
2238    Reflectometry as a function of kz for a set of slabs.
2239
2240    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
2241        This is :math:`\\tfrac12 Q_z`.       
2242    :param depth:  float[m] (Ang).
2243        thickness of each layer.  The thickness of the incident medium
2244        and substrate are ignored.
2245    :param rho:  float[n,k] (1e-6/Ang^2)
2246        Real scattering length density for each layer for each kz
2247    :param irho:  float[n,k] (1e-6/Ang^2)
2248        Imaginary scattering length density for each layer for each kz
2249        Note: absorption cross section mu = 2 irho/lambda for neutrons
2250    :param sigma: float[m-1] (Ang)
2251        interfacial roughness.  This is the roughness between a layer
2252        and the previous layer. The sigma array should have m-1 entries.
2253
2254    Slabs are ordered with the surface SLD at index 0 and substrate at
2255    index -1, or reversed if kz < 0.
2256    """
2257    def calc(kz, depth, rho, irho, sigma):
2258        if len(kz) == 0: return kz
2259   
2260        # Complex index of refraction is relative to the incident medium.
2261        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
2262        # in place of kz^2, and ignoring rho_o
2263        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
2264        k = kz
2265   
2266        # According to Heavens, the initial matrix should be [ 1 F; F 1],
2267        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
2268        # multiply versus some coding convenience.
2269        B11 = 1
2270        B22 = 1
2271        B21 = 0
2272        B12 = 0
2273        for i in range(0, len(depth)-1):
2274            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
2275            F = (k - k_next) / (k + k_next)
2276            F *= np.exp(-2*k*k_next*sigma[i]**2)
2277            #print "==== layer",i
2278            #print "kz:", kz
2279            #print "k:", k
2280            #print "k_next:",k_next
2281            #print "F:",F
2282            #print "rho:",rho[:,i+1]
2283            #print "irho:",irho[:,i+1]
2284            #print "d:",depth[i],"sigma:",sigma[i]
2285            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
2286            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
2287            M21 = F*M11
2288            M12 = F*M22
2289            C1 = B11*M11 + B21*M12
2290            C2 = B11*M21 + B21*M22
2291            B11 = C1
2292            B21 = C2
2293            C1 = B12*M11 + B22*M12
2294            C2 = B12*M21 + B22*M22
2295            B12 = C1
2296            B22 = C2
2297            k = k_next
2298   
2299        r = B12/B11
2300        return np.absolute(r)**2
2301
2302    if np.isscalar(kz): kz = np.asarray([kz], 'd')
2303
2304    m = len(depth)
2305
2306    # Make everything into arrays
2307    depth = np.asarray(depth,'d')
2308    rho = np.asarray(rho,'d')
2309    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
2310    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
2311
2312    # Repeat rho,irho columns as needed
2313    if len(rho.shape) == 1:
2314        rho = rho[None,:]
2315        irho = irho[None,:]
2316
2317    return calc(kz, depth, rho, irho, sigma)
2318   
2319def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
2320    y = abeles(kz, depth, rho, irho, sigma)
2321    s = dq/2.
2322    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
2323    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
2324    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
2325    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
2326    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
2327    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
2328    y *= 0.137023
2329    return y
2330       
2331def makeRefdFFT(Limits,Profile):
2332    print 'make fft'
2333    Q,Io = Profile[:2]
2334    Qmin = Limits[1][0]
2335    Qmax = Limits[1][1]
2336    iBeg = np.searchsorted(Q,Qmin)
2337    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2338    Qf = np.linspace(0.,Q[iFin-1],5000)
2339    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
2340    If = QI(Qf)*Qf**4
2341    R = np.linspace(0.,5000.,5000)
2342    F = fft.rfft(If)
2343    return R,F
2344
2345   
2346################################################################################
2347# Stacking fault simulation codes
2348################################################################################
2349
2350def GetStackParms(Layers):
2351   
2352    Parms = []
2353#cell parms
2354    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
2355        Parms.append('cellA')
2356        Parms.append('cellC')
2357    else:
2358        Parms.append('cellA')
2359        Parms.append('cellB')
2360        Parms.append('cellC')
2361        if Layers['Laue'] != 'mmm':
2362            Parms.append('cellG')
2363#Transition parms
2364    for iY in range(len(Layers['Layers'])):
2365        for iX in range(len(Layers['Layers'])):
2366            Parms.append('TransP;%d;%d'%(iY,iX))
2367            Parms.append('TransX;%d;%d'%(iY,iX))
2368            Parms.append('TransY;%d;%d'%(iY,iX))
2369            Parms.append('TransZ;%d;%d'%(iY,iX))
2370    return Parms
2371
2372def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
2373    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
2374   
2375    :param dict Layers: dict with following items
2376
2377      ::
2378
2379       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
2380       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
2381       'Layers':[],'Stacking':[],'Transitions':[]}
2382       
2383    :param str ctrls: controls string to be written on DIFFaX controls.dif file
2384    :param float scale: scale factor
2385    :param dict background: background parameters
2386    :param list limits: min/max 2-theta to be calculated
2387    :param dict inst: instrument parameters dictionary
2388    :param list profile: powder pattern data
2389   
2390    Note that parameters all updated in place   
2391    '''
2392    import atmdata
2393    path = sys.path
2394    for name in path:
2395        if 'bin' in name:
2396            DIFFaX = name+'/DIFFaX.exe'
2397            print ' Execute ',DIFFaX
2398            break
2399    # make form factor file that DIFFaX wants - atom types are GSASII style
2400    sf = open('data.sfc','w')
2401    sf.write('GSASII special form factor file for DIFFaX\n\n')
2402    atTypes = Layers['AtInfo'].keys()
2403    if 'H' not in atTypes:
2404        atTypes.insert(0,'H')
2405    for atType in atTypes:
2406        if atType == 'H': 
2407            blen = -.3741
2408        else:
2409            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
2410        Adat = atmdata.XrayFF[atType]
2411        text = '%4s'%(atType.ljust(4))
2412        for i in range(4):
2413            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
2414        text += '%11.6f%11.6f'%(Adat['fc'],blen)
2415        text += '%3d\n'%(Adat['Z'])
2416        sf.write(text)
2417    sf.close()
2418    #make DIFFaX control.dif file - future use GUI to set some of these flags
2419    cf = open('control.dif','w')
2420    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2421        x0 = profile[0]
2422        iBeg = np.searchsorted(x0,limits[0])
2423        iFin = np.searchsorted(x0,limits[1])+1
2424        if iFin-iBeg > 20000:
2425            iFin = iBeg+20000
2426        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2427        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2428        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
2429    else:
2430        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2431        inst = {'Type':['XSC','XSC',]}
2432    cf.close()
2433    #make DIFFaX data file
2434    df = open('GSASII-DIFFaX.dat','w')
2435    df.write('INSTRUMENTAL\n')
2436    if 'X' in inst['Type'][0]:
2437        df.write('X-RAY\n')
2438    elif 'N' in inst['Type'][0]:
2439        df.write('NEUTRON\n')
2440    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
2441        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
2442        U = ateln2*inst['U'][1]/10000.
2443        V = ateln2*inst['V'][1]/10000.
2444        W = ateln2*inst['W'][1]/10000.
2445        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2446        HW = np.sqrt(np.mean(HWHM))
2447    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
2448        if 'Mean' in Layers['selInst']:
2449            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
2450        elif 'Gaussian' in Layers['selInst']:
2451            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
2452        else:
2453            df.write('None\n')
2454    else:
2455        df.write('0.10\nNone\n')
2456    df.write('STRUCTURAL\n')
2457    a,b,c = Layers['Cell'][1:4]
2458    gam = Layers['Cell'][6]
2459    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
2460    laue = Layers['Laue']
2461    if laue == '2/m(ab)':
2462        laue = '2/m(1)'
2463    elif laue == '2/m(c)':
2464        laue = '2/m(2)'
2465    if 'unknown' in Layers['Laue']:
2466        df.write('%s %.3f\n'%(laue,Layers['Toler']))
2467    else:   
2468        df.write('%s\n'%(laue))
2469    df.write('%d\n'%(len(Layers['Layers'])))
2470    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
2471        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
2472    layerNames = []
2473    for layer in Layers['Layers']:
2474        layerNames.append(layer['Name'])
2475    for il,layer in enumerate(Layers['Layers']):
2476        if layer['SameAs']:
2477            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
2478            continue
2479        df.write('LAYER %d\n'%(il+1))
2480        if '-1' in layer['Symm']:
2481            df.write('CENTROSYMMETRIC\n')
2482        else:
2483            df.write('NONE\n')
2484        for ia,atom in enumerate(layer['Atoms']):
2485            [name,atype,x,y,z,frac,Uiso] = atom
2486            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
2487                frac /= 2.
2488            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
2489    df.write('STACKING\n')
2490    df.write('%s\n'%(Layers['Stacking'][0]))
2491    if 'recursive' in Layers['Stacking'][0]:
2492        df.write('%s\n'%Layers['Stacking'][1])
2493    else:
2494        if 'list' in Layers['Stacking'][1]:
2495            Slen = len(Layers['Stacking'][2])
2496            iB = 0
2497            iF = 0
2498            while True:
2499                iF += 68
2500                if iF >= Slen:
2501                    break
2502                iF = min(iF,Slen)
2503                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
2504                iB = iF
2505        else:
2506            df.write('%s\n'%Layers['Stacking'][1])   
2507    df.write('TRANSITIONS\n')
2508    for iY in range(len(Layers['Layers'])):
2509        sumPx = 0.
2510        for iX in range(len(Layers['Layers'])):
2511            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
2512            p = round(p,3)
2513            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
2514            sumPx += p
2515        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
2516            print 'ERROR - Layer probabilities sum to ',sumPx,' DIFFaX will insist it = 1.0'
2517            df.close()
2518            os.remove('data.sfc')
2519            os.remove('control.dif')
2520            os.remove('GSASII-DIFFaX.dat')
2521            return       
2522    df.close()   
2523    time0 = time.time()
2524    try:
2525        subp.call(DIFFaX)
2526    except OSError:
2527        print ' DIFFax.exe is not available for this platform - under development'
2528    print ' DIFFaX time = %.2fs'%(time.time()-time0)
2529    if os.path.exists('GSASII-DIFFaX.spc'):
2530        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
2531        iFin = iBeg+Xpat.shape[1]
2532        bakType,backDict,backVary = SetBackgroundParms(background)
2533        backDict['Lam1'] = G2mth.getWave(inst)
2534        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2535        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
2536        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2537            rv = st.poisson(profile[3][iBeg:iFin])
2538            profile[1][iBeg:iFin] = rv.rvs()
2539            Z = np.ones_like(profile[3][iBeg:iFin])
2540            Z[1::2] *= -1
2541            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2542            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2543        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2544    #cleanup files..
2545        os.remove('GSASII-DIFFaX.spc')
2546    elif os.path.exists('GSASII-DIFFaX.sadp'):
2547        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
2548        Sadp = np.reshape(Sadp,(256,-1))
2549        Layers['Sadp']['Img'] = Sadp
2550        os.remove('GSASII-DIFFaX.sadp')
2551    os.remove('data.sfc')
2552    os.remove('control.dif')
2553    os.remove('GSASII-DIFFaX.dat')
2554   
2555def SetPWDRscan(inst,limits,profile):
2556   
2557    wave = G2mth.getMeanWave(inst)
2558    x0 = profile[0]
2559    iBeg = np.searchsorted(x0,limits[0])
2560    iFin = np.searchsorted(x0,limits[1])
2561    if iFin-iBeg > 20000:
2562        iFin = iBeg+20000
2563    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2564    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
2565    return iFin-iBeg
2566       
2567def SetStackingSF(Layers,debug):
2568# Load scattering factors into DIFFaX arrays
2569    import atmdata
2570    atTypes = Layers['AtInfo'].keys()
2571    aTypes = []
2572    for atype in atTypes:
2573        aTypes.append('%4s'%(atype.ljust(4)))
2574    SFdat = []
2575    for atType in atTypes:
2576        Adat = atmdata.XrayFF[atType]
2577        SF = np.zeros(9)
2578        SF[:8:2] = Adat['fa']
2579        SF[1:8:2] = Adat['fb']
2580        SF[8] = Adat['fc']
2581        SFdat.append(SF)
2582    SFdat = np.array(SFdat)
2583    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
2584   
2585def SetStackingClay(Layers,Type):
2586# Controls
2587    rand.seed()
2588    ranSeed = rand.randint(1,2**16-1)
2589    try:   
2590        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
2591            '6/m','6/mmm'].index(Layers['Laue'])+1
2592    except ValueError:  #for 'unknown'
2593        laueId = -1
2594    if 'SADP' in Type:
2595        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
2596        lmax = int(Layers['Sadp']['Lmax'])
2597    else:
2598        planeId = 0
2599        lmax = 0
2600# Sequences
2601    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
2602    try:
2603        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
2604    except ValueError:
2605        StkParm = -1
2606    if StkParm == 2:    #list
2607        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
2608        Nstk = len(StkSeq)
2609    else:
2610        Nstk = 1
2611        StkSeq = [0,]
2612    if StkParm == -1:
2613        StkParm = int(Layers['Stacking'][1])
2614    Wdth = Layers['Width'][0]
2615    mult = 1
2616    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
2617    LaueSym = Layers['Laue'].ljust(12)
2618    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
2619    return laueId,controls
2620   
2621def SetCellAtoms(Layers):
2622    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
2623# atoms in layers
2624    atTypes = Layers['AtInfo'].keys()
2625    AtomXOU = []
2626    AtomTp = []
2627    LayerSymm = []
2628    LayerNum = []
2629    layerNames = []
2630    Natm = 0
2631    Nuniq = 0
2632    for layer in Layers['Layers']:
2633        layerNames.append(layer['Name'])
2634    for il,layer in enumerate(Layers['Layers']):
2635        if layer['SameAs']:
2636            LayerNum.append(layerNames.index(layer['SameAs'])+1)
2637            continue
2638        else:
2639            LayerNum.append(il+1)
2640            Nuniq += 1
2641        if '-1' in layer['Symm']:
2642            LayerSymm.append(1)
2643        else:
2644            LayerSymm.append(0)
2645        for ia,atom in enumerate(layer['Atoms']):
2646            [name,atype,x,y,z,frac,Uiso] = atom
2647            Natm += 1
2648            AtomTp.append('%4s'%(atype.ljust(4)))
2649            Ta = atTypes.index(atype)+1
2650            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
2651    AtomXOU = np.array(AtomXOU)
2652    Nlayers = len(layerNames)
2653    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
2654    return Nlayers
2655   
2656def SetStackingTrans(Layers,Nlayers):
2657# Transitions
2658    TransX = []
2659    TransP = []
2660    for Ytrans in Layers['Transitions']:
2661        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
2662        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
2663    TransP = np.array(TransP,dtype='float').T
2664    TransX = np.array(TransX,dtype='float')
2665#    GSASIIpath.IPyBreak()
2666    pyx.pygettrans(Nlayers,TransP,TransX)
2667   
2668def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
2669# Scattering factors
2670    SetStackingSF(Layers,debug)
2671# Controls & sequences
2672    laueId,controls = SetStackingClay(Layers,'PWDR')
2673# cell & atoms
2674    Nlayers = SetCellAtoms(Layers)
2675    Volume = Layers['Cell'][7]   
2676# Transitions
2677    SetStackingTrans(Layers,Nlayers)
2678# PWDR scan
2679    Nsteps = SetPWDRscan(inst,limits,profile)
2680# result as Spec
2681    x0 = profile[0]
2682    profile[3] = np.zeros(len(profile[0]))
2683    profile[4] = np.zeros(len(profile[0]))
2684    profile[5] = np.zeros(len(profile[0]))
2685    iBeg = np.searchsorted(x0,limits[0])
2686    iFin = np.searchsorted(x0,limits[1])+1
2687    if iFin-iBeg > 20000:
2688        iFin = iBeg+20000
2689    Nspec = 20001       
2690    spec = np.zeros(Nspec,dtype='double')   
2691    time0 = time.time()
2692    pyx.pygetspc(controls,Nspec,spec)
2693    print ' GETSPC time = %.2fs'%(time.time()-time0)
2694    time0 = time.time()
2695    U = ateln2*inst['U'][1]/10000.
2696    V = ateln2*inst['V'][1]/10000.
2697    W = ateln2*inst['W'][1]/10000.
2698    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
2699    HW = np.sqrt(np.mean(HWHM))
2700    BrdSpec = np.zeros(Nsteps)
2701    if 'Mean' in Layers['selInst']:
2702        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
2703    elif 'Gaussian' in Layers['selInst']:
2704        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
2705    else:
2706        BrdSpec = spec[:Nsteps]
2707    BrdSpec /= Volume
2708    iFin = iBeg+Nsteps
2709    bakType,backDict,backVary = SetBackgroundParms(background)
2710    backDict['Lam1'] = G2mth.getWave(inst)
2711    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
2712    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
2713    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
2714        try:
2715            rv = st.poisson(profile[3][iBeg:iFin])
2716            profile[1][iBeg:iFin] = rv.rvs()
2717        except ValueError:
2718            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
2719        Z = np.ones_like(profile[3][iBeg:iFin])
2720        Z[1::2] *= -1
2721        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2722        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2723    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2724    print ' Broadening time = %.2fs'%(time.time()-time0)
2725   
2726def CalcStackingSADP(Layers,debug):
2727   
2728# Scattering factors
2729    SetStackingSF(Layers,debug)
2730# Controls & sequences
2731    laueId,controls = SetStackingClay(Layers,'SADP')
2732# cell & atoms
2733    Nlayers = SetCellAtoms(Layers)   
2734# Transitions
2735    SetStackingTrans(Layers,Nlayers)
2736# result as Sadp
2737    Nspec = 20001       
2738    spec = np.zeros(Nspec,dtype='double')   
2739    time0 = time.time()
2740    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
2741    Sapd = np.zeros((256,256))
2742    iB = 0
2743    for i in range(hkLim):
2744        iF = iB+Nblk
2745        p1 = 127+int(i*Incr)
2746        p2 = 128-int(i*Incr)
2747        if Nblk == 128:
2748            if i:
2749                Sapd[128:,p1] = spec[iB:iF]
2750                Sapd[:128,p1] = spec[iF:iB:-1]
2751            Sapd[128:,p2] = spec[iB:iF]
2752            Sapd[:128,p2] = spec[iF:iB:-1]
2753        else:
2754            if i:
2755                Sapd[:,p1] = spec[iB:iF]
2756            Sapd[:,p2] = spec[iB:iF]
2757        iB += Nblk
2758    Layers['Sadp']['Img'] = Sapd
2759    print ' GETSAD time = %.2fs'%(time.time()-time0)
2760#    GSASIIpath.IPyBreak()
2761   
2762#testing data
2763NeedTestData = True
2764def TestData():
2765    'needs a doc string'
2766#    global NeedTestData
2767    global bakType
2768    bakType = 'chebyschev'
2769    global xdata
2770    xdata = np.linspace(4.0,40.0,36000)
2771    global parmDict0
2772    parmDict0 = {
2773        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
2774        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
2775        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
2776        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
2777        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
2778        'Back0':5.384,'Back1':-0.015,'Back2':.004,
2779        }
2780    global parmDict1
2781    parmDict1 = {
2782        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
2783        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
2784        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
2785        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
2786        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
2787        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
2788        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
2789        'Back0':36.897,'Back1':-0.508,'Back2':.006,
2790        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2791        }
2792    global parmDict2
2793    parmDict2 = {
2794        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
2795        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
2796        'Back0':5.,'Back1':-0.02,'Back2':.004,
2797#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2798        }
2799    global varyList
2800    varyList = []
2801
2802def test0():
2803    if NeedTestData: TestData()
2804    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
2805    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
2806    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
2807    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
2808    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
2809    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
2810   
2811def test1():
2812    if NeedTestData: TestData()
2813    time0 = time.time()
2814    for i in range(100):
2815        getPeakProfile(parmDict1,xdata,varyList,bakType)
2816    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
2817   
2818def test2(name,delt):
2819    if NeedTestData: TestData()
2820    varyList = [name,]
2821    xdata = np.linspace(5.6,5.8,400)
2822    hplot = plotter.add('derivatives test for '+name).gca()
2823    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
2824    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2825    parmDict2[name] += delt
2826    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2827    hplot.plot(xdata,(y1-y0)/delt,'r+')
2828   
2829def test3(name,delt):
2830    if NeedTestData: TestData()
2831    names = ['pos','sig','gam','shl']
2832    idx = names.index(name)
2833    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
2834    xdata = np.linspace(5.6,5.8,800)
2835    dx = xdata[1]-xdata[0]
2836    hplot = plotter.add('derivatives test for '+name).gca()
2837    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
2838    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2839    myDict[name] += delt
2840    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2841    hplot.plot(xdata,(y1-y0)/delt,'r+')
2842
2843if __name__ == '__main__':
2844    import GSASIItestplot as plot
2845    global plotter
2846    plotter = plot.PlotNotebook()
2847#    test0()
2848#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
2849    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
2850        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
2851        test2(name,shft)
2852    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
2853        test3(name,shft)
2854    print "OK"
2855    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.