source: trunk/GSASIIpwd.py @ 4328

Last change on this file since 4328 was 4328, checked in by vondreele, 20 months ago

Add FITTED_SCALE to RMCprofile options
fix issues with changing atom sequence in RMCProfile
remove read of x-ray G(R) - it isn't used by RMCProfile but is generated from F(Q)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 145.7 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: 2020-02-25 12:47:52 +0000 (Tue, 25 Feb 2020) $
10# $Author: vondreele $
11# $Revision: 4328 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 4328 2020-02-25 12:47:52Z vondreele $
14########### SVN repository information ###################
15from __future__ import division, print_function
16import sys
17import math
18import time
19import os
20import subprocess as subp
21import copy
22
23import numpy as np
24import numpy.linalg as nl
25import numpy.ma as ma
26import random as rand
27import numpy.fft as fft
28import scipy.interpolate as si
29import scipy.stats as st
30import scipy.optimize as so
31import scipy.special as sp
32
33import GSASIIpath
34GSASIIpath.SetVersionNumber("$Revision: 4328 $")
35import GSASIIlattice as G2lat
36import GSASIIspc as G2spc
37import GSASIIElem as G2elem
38import GSASIImath as G2mth
39try:
40    import pypowder as pyd
41except ImportError:
42    print ('pypowder is not available - profile calcs. not allowed')
43try:
44    import pydiffax as pyx
45except ImportError:
46    print ('pydiffax is not available for this platform')
47import GSASIIfiles as G2fil
48
49   
50# trig functions in degrees
51tand = lambda x: math.tan(x*math.pi/180.)
52atand = lambda x: 180.*math.atan(x)/math.pi
53atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
54cosd = lambda x: math.cos(x*math.pi/180.)
55acosd = lambda x: 180.*math.acos(x)/math.pi
56rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
57#numpy versions
58npsind = lambda x: np.sin(x*np.pi/180.)
59npasind = lambda x: 180.*np.arcsin(x)/math.pi
60npcosd = lambda x: np.cos(x*math.pi/180.)
61npacosd = lambda x: 180.*np.arccos(x)/math.pi
62nptand = lambda x: np.tan(x*math.pi/180.)
63npatand = lambda x: 180.*np.arctan(x)/np.pi
64npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
65npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave    #=d*
66npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)    #=2pi*d*
67ateln2 = 8.0*math.log(2.0)
68sateln2 = np.sqrt(ateln2)
69nxs = np.newaxis
70
71################################################################################
72#### Powder utilities
73################################################################################
74
75def PhaseWtSum(G2frame,histo):
76    '''
77    Calculate sum of phase mass*phase fraction for PWDR data (exclude magnetic phases)
78   
79    :param G2frame: GSASII main frame structure
80    :param str histo: histogram name
81    :returns: sum(scale*mass) for phases in histo
82    '''
83    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
84    wtSum = 0.0
85    for phase in Phases:
86        if Phases[phase]['General']['Type'] != 'magnetic':
87            if histo in Phases[phase]['Histograms']:
88                if not Phases[phase]['Histograms'][histo]['Use']: continue
89                mass = Phases[phase]['General']['Mass']
90                phFr = Phases[phase]['Histograms'][histo]['Scale'][0]
91                wtSum += mass*phFr
92    return wtSum
93   
94################################################################################
95#### GSASII pwdr & pdf calculation routines
96################################################################################
97       
98def Transmission(Geometry,Abs,Diam):
99    '''
100    Calculate sample transmission
101
102    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
103    :param float Abs: absorption coeff in cm-1
104    :param float Diam: sample thickness/diameter in mm
105    '''
106    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
107        MuR = Abs*Diam/20.0
108        if MuR <= 3.0:
109            T0 = 16/(3.*math.pi)
110            T1 = -0.045780
111            T2 = -0.02489
112            T3 = 0.003045
113            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
114            if T < -20.:
115                return 2.06e-9
116            else:
117                return math.exp(T)
118        else:
119            T1 = 1.433902
120            T2 = 0.013869+0.337894
121            T3 = 1.933433+1.163198
122            T4 = 0.044365-0.04259
123            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
124            return T/100.
125    elif 'plate' in Geometry:
126        MuR = Abs*Diam/10.
127        return math.exp(-MuR)
128    elif 'Bragg' in Geometry:
129        return 0.0
130       
131def SurfaceRough(SRA,SRB,Tth):
132    ''' Suortti (J. Appl. Cryst, 5,325-331, 1972) surface roughness correction
133    :param float SRA: Suortti surface roughness parameter
134    :param float SRB: Suortti surface roughness parameter
135    :param float Tth: 2-theta(deg) - can be numpy array
136   
137    '''
138    sth = npsind(Tth/2.)
139    T1 = np.exp(-SRB/sth)
140    T2 = SRA+(1.-SRA)*np.exp(-SRB)
141    return (SRA+(1.-SRA)*T1)/T2
142   
143def SurfaceRoughDerv(SRA,SRB,Tth):
144    ''' Suortti surface roughness correction derivatives
145    :param float SRA: Suortti surface roughness parameter (dimensionless)
146    :param float SRB: Suortti surface roughness parameter (dimensionless)
147    :param float Tth: 2-theta(deg) - can be numpy array
148    :return list: [dydSRA,dydSRB] derivatives to be used for intensity derivative
149    '''
150    sth = npsind(Tth/2.)
151    T1 = np.exp(-SRB/sth)
152    T2 = SRA+(1.-SRA)*np.exp(-SRB)
153    Trans = (SRA+(1.-SRA)*T1)/T2
154    dydSRA = ((1.-T1)*T2-(1.-np.exp(-SRB))*Trans)/T2**2
155    dydSRB = ((SRA-1.)*T1*T2/sth-Trans*(SRA-T2))/T2**2
156    return [dydSRA,dydSRB]
157
158def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
159    '''Calculate sample absorption
160    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
161    :param float MuR: absorption coeff * sample thickness/2 or radius
162    :param Tth: 2-theta scattering angle - can be numpy array
163    :param float Phi: flat plate tilt angle - future
164    :param float Psi: flat plate tilt axis - future
165    '''
166   
167    def muRunder3(MuR,Sth2):
168        T0 = 16.0/(3.*np.pi)
169        T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
170            0.109561*np.sqrt(Sth2)-26.04556
171        T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
172            1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
173        T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
174        Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
175        return np.exp(Trns)
176   
177    def muRover3(MuR,Sth2):
178        T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
179            10.02088*Sth2**3-3.36778*Sth2**4
180        T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
181            (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
182        T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
183            0.13576*np.sqrt(Sth2)+1.163198
184        T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
185        Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
186        return Trns/100.
187       
188    Sth2 = npsind(Tth/2.0)**2
189    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
190        if 'array' in str(type(MuR)):
191            MuRSTh2 = np.vstack((MuR,Sth2))
192            AbsCr = np.where(MuRSTh2[0]<=3.0,muRunder3(MuRSTh2[0],MuRSTh2[1]),muRover3(MuRSTh2[0],MuRSTh2[1]))
193            return AbsCr
194        else:
195            if MuR <= 3.0:
196                return muRunder3(MuR,Sth2)
197            else:
198                return muRover3(MuR,Sth2)
199    elif 'Bragg' in Geometry:
200        return 1.0
201    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
202        # and only defined for 2theta < 90
203        MuT = 2.*MuR
204        T1 = np.exp(-MuT)
205        T2 = np.exp(-MuT/npcosd(Tth))
206        Tb = MuT-MuT/npcosd(Tth)
207        return (T2-T1)/Tb
208    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
209        MuT = 2.*MuR
210        cth = npcosd(Tth/2.0)
211        return np.exp(-MuT/cth)/cth
212       
213def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
214    'needs a doc string'
215    dA = 0.001
216    AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
217    if MuR:
218        AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
219        return (AbsP-AbsM)/(2.0*dA)
220    else:
221        return (AbsP-1.)/dA
222       
223def Polarization(Pola,Tth,Azm=0.0):
224    """   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
225
226    :param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
227    :param Azm: azimuthal angle e.g. 0.0 in plane of polarization
228    :param Tth: 2-theta scattering angle - can be numpy array
229      which (if either) of these is "right"?
230    :return: (pola, dpdPola)
231      * pola = ((1-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \
232        (1-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
233      * dpdPola: derivative needed for least squares
234
235    """
236    cazm = npcosd(Azm)
237    sazm = npsind(Azm)
238    pola = ((1.0-Pola)*cazm**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
239        (1.0-Pola)*sazm**2+Pola*cazm**2
240    dpdPola = -npsind(Tth)**2*(sazm**2-cazm**2)
241    return pola,dpdPola
242   
243def Oblique(ObCoeff,Tth):
244    'currently assumes detector is normal to beam'
245    if ObCoeff:
246        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
247    else:
248        return 1.0
249               
250def Ruland(RulCoff,wave,Q,Compton):
251    'needs a doc string'
252    C = 2.9978e8
253    D = 1.5e-3
254    hmc = 0.024262734687
255    sinth2 = (Q*wave/(4.0*np.pi))**2
256    dlam = (wave**2)*Compton*Q/C
257    dlam_c = 2.0*hmc*sinth2-D*wave**2
258    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
259   
260def LorchWeight(Q):
261    'needs a doc string'
262    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
263           
264def GetAsfMean(ElList,Sthl2):
265    '''Calculate various scattering factor terms for PDF calcs
266
267    :param dict ElList: element dictionary contains scattering factor coefficients, etc.
268    :param np.array Sthl2: numpy array of sin theta/lambda squared values
269    :returns: mean(f^2), mean(f)^2, mean(compton)
270    '''
271    sumNoAtoms = 0.0
272    FF = np.zeros_like(Sthl2)
273    FF2 = np.zeros_like(Sthl2)
274    CF = np.zeros_like(Sthl2)
275    for El in ElList:
276        sumNoAtoms += ElList[El]['FormulaNo']
277    for El in ElList:
278        el = ElList[El]
279        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
280        cf = G2elem.ComptonFac(el,Sthl2)
281        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
282        FF2 += ff2*el['FormulaNo']/sumNoAtoms
283        CF += cf*el['FormulaNo']/sumNoAtoms
284    return FF2,FF**2,CF
285   
286def GetNumDensity(ElList,Vol):
287    'needs a doc string'
288    sumNoAtoms = 0.0
289    for El in ElList:
290        sumNoAtoms += ElList[El]['FormulaNo']
291    return sumNoAtoms/Vol
292           
293def CalcPDF(data,inst,limits,xydata):
294    '''Computes I(Q), S(Q) & G(r) from Sample, Bkg, etc. diffraction patterns loaded into
295    dict xydata; results are placed in xydata.
296    Calculation parameters are found in dicts data and inst and list limits.
297    The return value is at present an empty list.
298    '''
299    auxPlot = []
300    if 'T' in inst['Type'][0]:
301        Ibeg = 0
302        Ifin = len(xydata['Sample'][1][0])
303    else:
304        Ibeg = np.searchsorted(xydata['Sample'][1][0],limits[0])
305        Ifin = np.searchsorted(xydata['Sample'][1][0],limits[1])+1
306    #subtract backgrounds - if any & use PWDR limits
307    IofQ = copy.deepcopy(xydata['Sample'])
308    IofQ[1] = np.array(IofQ[1])[:,Ibeg:Ifin]
309    if data['Sample Bkg.']['Name']:
310        IofQ[1][1] += xydata['Sample Bkg.'][1][1][Ibeg:Ifin]*data['Sample Bkg.']['Mult']
311    if data['Container']['Name']:
312        xycontainer = xydata['Container'][1][1]*data['Container']['Mult']
313        if data['Container Bkg.']['Name']:
314            xycontainer += xydata['Container Bkg.'][1][1][Ibeg:Ifin]*data['Container Bkg.']['Mult']
315        IofQ[1][1] += xycontainer[Ibeg:Ifin]
316    data['IofQmin'] = IofQ[1][1][-1]
317    IofQ[1][1] -= data.get('Flat Bkg',0.)
318    #get element data & absorption coeff.
319    ElList = data['ElList']
320    Tth = IofQ[1][0]    #2-theta or TOF!
321    if 'X' in inst['Type'][0]:
322        Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
323        #Apply angle dependent corrections
324        MuR = Abs*data['Diam']/20.0
325        IofQ[1][1] /= Absorb(data['Geometry'],MuR,Tth)
326        IofQ[1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
327        if data['DetType'] == 'Image plate':
328            IofQ[1][1] *= Oblique(data['ObliqCoeff'],Tth)
329    elif 'T' in inst['Type'][0]:    #neutron TOF normalized data - needs wavelength dependent absorption
330        wave = 2.*G2lat.TOF2dsp(inst,IofQ[1][0])*npsind(inst['2-theta'][1]/2.)
331        Els = ElList.keys()
332        Isotope = {El:'Nat. abund.' for El in Els}
333        GD = {'AtomTypes':ElList,'Isotope':Isotope}
334        BLtables = G2elem.GetBLtable(GD)
335        FP,FPP = G2elem.BlenResTOF(Els,BLtables,wave)
336        Abs = np.zeros(len(wave))
337        for iel,El in enumerate(Els):
338            BL = BLtables[El][1]
339            SA = BL['SA']*wave/1.798197+4.0*np.pi*FPP[iel]**2 #+BL['SL'][1]?
340            SA *= ElList[El]['FormulaNo']/data['Form Vol']
341            Abs += SA
342        MuR = Abs*data['Diam']/2.
343        IofQ[1][1] /= Absorb(data['Geometry'],MuR,inst['2-theta'][1]*np.ones(len(wave)))       
344    XY = IofQ[1]   
345    #convert to Q
346#    nQpoints = len(XY[0])     #points for Q interpolation
347    nQpoints = 5000
348    if 'C' in inst['Type'][0]:
349        wave = G2mth.getWave(inst)
350        minQ = npT2q(Tth[0],wave)
351        maxQ = npT2q(Tth[-1],wave)   
352        Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
353        dq = Qpoints[1]-Qpoints[0]
354        XY[0] = npT2q(XY[0],wave)
355        Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0])    #interpolate I(Q)
356    elif 'T' in inst['Type'][0]:
357        difC = inst['difC'][1]
358        minQ = 2.*np.pi*difC/Tth[-1]
359        maxQ = 2.*np.pi*difC/Tth[0]
360        Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
361        dq = Qpoints[1]-Qpoints[0]
362        XY[0] = 2.*np.pi*difC/XY[0]
363        Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][-1])    #interpolate I(Q)
364    Qdata -= np.min(Qdata)*data['BackRatio']
365   
366    qLimits = data['QScaleLim']
367    maxQ = np.searchsorted(Qpoints,min(Qpoints[-1],qLimits[1]))+1
368    minQ = np.searchsorted(Qpoints,min(qLimits[0],0.90*Qpoints[-1]))
369    qLimits = [Qpoints[minQ],Qpoints[maxQ-1]]
370    newdata = []
371    if len(IofQ) < 3:
372        xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],'']
373    else:
374        xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],IofQ[2]]
375    for item in xydata['IofQ'][1]:
376        newdata.append(item[:maxQ])
377    xydata['IofQ'][1] = newdata
378   
379    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
380    if 'XC' in inst['Type'][0]:
381        FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
382    else: #TOF
383        CF = np.zeros(len(xydata['SofQ'][1][0]))
384        FFSq = np.ones(len(xydata['SofQ'][1][0]))
385        SqFF = np.ones(len(xydata['SofQ'][1][0]))
386    Q = xydata['SofQ'][1][0]
387#    auxPlot.append([Q,np.copy(CF),'CF-unCorr'])
388    if 'XC' in inst['Type'][0]:
389        ruland = Ruland(data['Ruland'],wave,Q,CF)
390#    auxPlot.append([Q,ruland,'Ruland'])     
391        CF *= ruland
392#    auxPlot.append([Q,CF,'CF-Corr'])
393    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
394    xydata['SofQ'][1][1] *= scale
395    if 'XC' in inst['Type'][0]:
396        xydata['SofQ'][1][1] -= CF
397    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
398    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
399    xydata['SofQ'][1][1] *= scale
400    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
401    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
402    if data['Lorch']:
403        xydata['FofQ'][1][1] *= LorchWeight(Q)   
404    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
405    nR = len(xydata['GofR'][1][1])
406    Rmax = GSASIIpath.GetConfigValue('PDF_Rmax',100.)
407    mul = int(round(2.*np.pi*nR/(Rmax*qLimits[1])))
408#    mul = int(round(2.*np.pi*nR/(data.get('Rmax',100.)*qLimits[1])))
409    xydata['GofR'][1][0] = 2.*np.pi*np.linspace(0,nR,nR,endpoint=True)/(mul*qLimits[1])
410    xydata['GofR'][1][1] = -dq*np.imag(fft.fft(xydata['FofQ'][1][1],mul*nR)[:nR])
411    if data.get('noRing',True):
412        xydata['GofR'][1][1] = np.where(xydata['GofR'][1][0]<0.5,0.,xydata['GofR'][1][1])
413    return auxPlot
414   
415def PDFPeakFit(peaks,data):
416    rs2pi = 1./np.sqrt(2*np.pi)
417   
418    def MakeParms(peaks):
419        varyList = []
420        parmDict = {'slope':peaks['Background'][1][1]}
421        if peaks['Background'][2]:
422            varyList.append('slope')
423        for i,peak in enumerate(peaks['Peaks']):
424            parmDict['PDFpos;'+str(i)] = peak[0]
425            parmDict['PDFmag;'+str(i)] = peak[1]
426            parmDict['PDFsig;'+str(i)] = peak[2]
427            if 'P' in peak[3]:
428                varyList.append('PDFpos;'+str(i))
429            if 'M' in peak[3]:
430                varyList.append('PDFmag;'+str(i))
431            if 'S' in peak[3]:
432                varyList.append('PDFsig;'+str(i))
433        return parmDict,varyList
434       
435    def SetParms(peaks,parmDict,varyList):
436        if 'slope' in varyList:
437            peaks['Background'][1][1] = parmDict['slope']
438        for i,peak in enumerate(peaks['Peaks']):
439            if 'PDFpos;'+str(i) in varyList:
440                peak[0] = parmDict['PDFpos;'+str(i)]
441            if 'PDFmag;'+str(i) in varyList:
442                peak[1] = parmDict['PDFmag;'+str(i)]
443            if 'PDFsig;'+str(i) in varyList:
444                peak[2] = parmDict['PDFsig;'+str(i)]
445       
446   
447    def CalcPDFpeaks(parmdict,Xdata):
448        Z = parmDict['slope']*Xdata
449        ipeak = 0
450        while True:
451            try:
452                pos = parmdict['PDFpos;'+str(ipeak)]
453                mag = parmdict['PDFmag;'+str(ipeak)]
454                wid = parmdict['PDFsig;'+str(ipeak)]
455                wid2 = 2.*wid**2
456                Z += mag*rs2pi*np.exp(-(Xdata-pos)**2/wid2)/wid
457                ipeak += 1
458            except KeyError:        #no more peaks to process
459                return Z
460               
461    def errPDFProfile(values,xdata,ydata,parmdict,varylist):       
462        parmdict.update(zip(varylist,values))
463        M = CalcPDFpeaks(parmdict,xdata)-ydata
464        return M
465           
466    newpeaks = copy.copy(peaks)
467    iBeg = np.searchsorted(data[1][0],newpeaks['Limits'][0])
468    iFin = np.searchsorted(data[1][0],newpeaks['Limits'][1])+1
469    X = data[1][0][iBeg:iFin]
470    Y = data[1][1][iBeg:iFin]
471    parmDict,varyList = MakeParms(peaks)
472    if not len(varyList):
473        G2fil.G2Print (' Nothing varied')
474        return newpeaks,None,None,None,None,None
475   
476    Rvals = {}
477    values =  np.array(Dict2Values(parmDict, varyList))
478    result = so.leastsq(errPDFProfile,values,full_output=True,ftol=0.0001,
479           args=(X,Y,parmDict,varyList))
480    chisq = np.sum(result[2]['fvec']**2)
481    Values2Dict(parmDict, varyList, result[0])
482    SetParms(peaks,parmDict,varyList)
483    Rvals['Rwp'] = np.sqrt(chisq/np.sum(Y**2))*100.      #to %
484    chisq = np.sum(result[2]['fvec']**2)/(len(X)-len(values))   #reduced chi^2 = M/(Nobs-Nvar)
485    sigList = list(np.sqrt(chisq*np.diag(result[1])))   
486    Z = CalcPDFpeaks(parmDict,X)
487    newpeaks['calc'] = [X,Z]
488    return newpeaks,result[0],varyList,sigList,parmDict,Rvals   
489   
490def MakeRDF(RDFcontrols,background,inst,pwddata):
491    import scipy.signal as signal
492    auxPlot = []
493    if 'C' in inst['Type'][0]:
494        Tth = pwddata[0]
495        wave = G2mth.getWave(inst)
496        minQ = npT2q(Tth[0],wave)
497        maxQ = npT2q(Tth[-1],wave)
498        powQ = npT2q(Tth,wave) 
499    elif 'T' in inst['Type'][0]:
500        TOF = pwddata[0]
501        difC = inst['difC'][1]
502        minQ = 2.*np.pi*difC/TOF[-1]
503        maxQ = 2.*np.pi*difC/TOF[0]
504        powQ = 2.*np.pi*difC/TOF
505    piDQ = np.pi/(maxQ-minQ)
506    Qpoints = np.linspace(minQ,maxQ,len(pwddata[0]),endpoint=True)
507    if RDFcontrols['UseObsCalc'] == 'obs-calc':
508        Qdata = si.griddata(powQ,pwddata[1]-pwddata[3],Qpoints,method=RDFcontrols['Smooth'],fill_value=0.)
509    elif RDFcontrols['UseObsCalc'] == 'obs-back':
510        Qdata = si.griddata(powQ,pwddata[1]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
511    elif RDFcontrols['UseObsCalc'] == 'calc-back':
512        Qdata = si.griddata(powQ,pwddata[3]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
513    Qdata *= np.sin((Qpoints-minQ)*piDQ)/piDQ
514    Qdata *= 0.5*np.sqrt(Qpoints)       #Qbin normalization
515#    GSASIIpath.IPyBreak()
516    dq = Qpoints[1]-Qpoints[0]
517    nR = len(Qdata)
518    R = 0.5*np.pi*np.linspace(0,nR,nR)/(4.*maxQ)
519    iFin = np.searchsorted(R,RDFcontrols['maxR'])+1
520    bBut,aBut = signal.butter(4,0.01)
521    Qsmooth = signal.filtfilt(bBut,aBut,Qdata)
522#    auxPlot.append([Qpoints,Qdata,'interpolate:'+RDFcontrols['Smooth']])
523#    auxPlot.append([Qpoints,Qsmooth,'interpolate:'+RDFcontrols['Smooth']])
524    DofR = dq*np.imag(fft.fft(Qsmooth,16*nR)[:nR])
525#    DofR = dq*np.imag(ft.fft(Qsmooth,16*nR)[:nR])
526    auxPlot.append([R[:iFin],DofR[:iFin],'D(R) for '+RDFcontrols['UseObsCalc']])   
527    return auxPlot
528
529# PDF optimization =============================================================
530def OptimizePDF(data,xydata,limits,inst,showFit=True,maxCycles=5):
531    import scipy.optimize as opt
532    numbDen = GetNumDensity(data['ElList'],data['Form Vol'])
533    Min,Init,Done = SetupPDFEval(data,xydata,limits,inst,numbDen)
534    xstart = Init()
535    bakMul = data['Sample Bkg.']['Mult']
536    if showFit:
537        rms = Min(xstart)
538        G2fil.G2Print('  Optimizing corrections to improve G(r) at low r')
539        if data['Sample Bkg.'].get('Refine',False):
540#            data['Flat Bkg'] = 0.
541            G2fil.G2Print('  start: Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f})'.format(
542                data['Ruland'],data['Sample Bkg.']['Mult'],rms))
543        else:
544            G2fil.G2Print('  start: Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} (RMS:{:.4f})'.format(
545                data['Flat Bkg'],data['BackRatio'],data['Ruland'],rms))
546    if data['Sample Bkg.'].get('Refine',False):
547        res = opt.minimize(Min,xstart,bounds=([0.01,1],[1.2*bakMul,0.8*bakMul]),
548                    method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001)
549    else:
550        res = opt.minimize(Min,xstart,bounds=([0,None],[0,1],[0.01,1]),
551                    method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001)
552    Done(res['x'])
553    if showFit:
554        if res['success']:
555            msg = 'Converged'
556        else:
557            msg = 'Not Converged'
558        if data['Sample Bkg.'].get('Refine',False):
559            G2fil.G2Print('  end:   Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f}) *** {} ***\n'.format(
560                data['Ruland'],data['Sample Bkg.']['Mult'],res['fun'],msg))
561        else:
562            G2fil.G2Print('  end:   Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f}) *** {} ***\n'.format(
563                data['Flat Bkg'],data['BackRatio'],data['Ruland'],res['fun'],msg))
564    return res
565
566def SetupPDFEval(data,xydata,limits,inst,numbDen):
567    Data = copy.deepcopy(data)
568    BkgMax = 1.
569    def EvalLowPDF(arg):
570        '''Objective routine -- evaluates the RMS deviations in G(r)
571        from -4(pi)*#density*r for for r<Rmin
572        arguments are ['Flat Bkg','BackRatio','Ruland'] scaled so that
573        the min & max values are between 0 and 1.
574        '''
575        if Data['Sample Bkg.'].get('Refine',False):
576            R,S = arg
577            Data['Sample Bkg.']['Mult'] = S
578        else:
579            F,B,R = arg
580            Data['Flat Bkg'] = F*BkgMax
581            Data['BackRatio'] = B
582        Data['Ruland'] = R/10.
583        CalcPDF(Data,inst,limits,xydata)
584        # test low r computation
585        g = xydata['GofR'][1][1]
586        r = xydata['GofR'][1][0]
587        g0 = g[r < Data['Rmin']] + 4*np.pi*r[r < Data['Rmin']]*numbDen
588        M = sum(g0**2)/len(g0)
589        return M
590    def GetCurrentVals():
591        '''Get the current ['Flat Bkg','BackRatio','Ruland'] with scaling
592        '''
593        if data['Sample Bkg.'].get('Refine',False):
594                return [max(10*data['Ruland'],.05),data['Sample']['Mult']]
595        try:
596            F = data['Flat Bkg']/BkgMax
597        except:
598            F = 0
599        return [F,data['BackRatio'],max(10*data['Ruland'],.05)]
600    def SetFinalVals(arg):
601        '''Set the 'Flat Bkg', 'BackRatio' & 'Ruland' values from the
602        scaled, refined values and plot corrected region of G(r)
603        '''
604        if data['Sample Bkg.'].get('Refine',False):
605            R,S = arg
606            data['Sample Bkg.']['Mult'] = S
607        else:
608            F,B,R = arg
609            data['Flat Bkg'] = F*BkgMax
610            data['BackRatio'] = B
611        data['Ruland'] = R/10.
612        CalcPDF(data,inst,limits,xydata)
613    EvalLowPDF(GetCurrentVals())
614    BkgMax = max(xydata['IofQ'][1][1])/50.
615    return EvalLowPDF,GetCurrentVals,SetFinalVals
616
617################################################################################       
618#### GSASII peak fitting routines: Finger, Cox & Jephcoat model       
619################################################################################
620
621def factorize(num):
622    ''' Provide prime number factors for integer num
623    :returns: dictionary of prime factors (keys) & power for each (data)
624    '''
625    factors = {}
626    orig = num
627
628    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
629    i, sqi = 2, 4
630    while sqi <= num:
631        while not num%i:
632            num /= i
633            factors[i] = factors.get(i, 0) + 1
634
635        sqi += 2*i + 1
636        i += 1
637
638    if num != 1 and num != orig:
639        factors[num] = factors.get(num, 0) + 1
640
641    if factors:
642        return factors
643    else:
644        return {num:1}          #a prime number!
645           
646def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
647    ''' Provide list of optimal data sizes for FFT calculations
648
649    :param int nmin: minimum data size >= 1
650    :param int nmax: maximum data size > nmin
651    :param int thresh: maximum prime factor allowed
652    :Returns: list of data sizes where the maximum prime factor is < thresh
653    ''' 
654    plist = []
655    nmin = max(1,nmin)
656    nmax = max(nmin+1,nmax)
657    for p in range(nmin,nmax):
658        if max(list(factorize(p).keys())) < thresh:
659            plist.append(p)
660    return plist
661
662np.seterr(divide='ignore')
663
664# Normal distribution
665
666# loc = mu, scale = std
667_norm_pdf_C = 1./math.sqrt(2*math.pi)
668class norm_gen(st.rv_continuous):
669    'needs a doc string'
670     
671    def pdf(self,x,*args,**kwds):
672        loc,scale=kwds['loc'],kwds['scale']
673        x = (x-loc)/scale
674        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
675       
676norm = norm_gen(name='norm',longname='A normal',extradoc="""
677
678Normal distribution
679
680The location (loc) keyword specifies the mean.
681The scale (scale) keyword specifies the standard deviation.
682
683normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
684""")
685
686## Cauchy
687
688# median = loc
689
690class cauchy_gen(st.rv_continuous):
691    'needs a doc string'
692
693    def pdf(self,x,*args,**kwds):
694        loc,scale=kwds['loc'],kwds['scale']
695        x = (x-loc)/scale
696        return 1.0/np.pi/(1.0+x*x) / scale
697       
698cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
699
700Cauchy distribution
701
702cauchy.pdf(x) = 1/(pi*(1+x**2))
703
704This is the t distribution with one degree of freedom.
705""")
706   
707   
708#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
709
710
711class fcjde_gen(st.rv_continuous):
712    """
713    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
714    Ref: J. Appl. Cryst. (1994) 27, 892-900.
715
716    :param x: array -1 to 1
717    :param t: 2-theta position of peak
718    :param s: sum(S/L,H/L); S: sample height, H: detector opening,
719      L: sample to detector opening distance
720    :param dx: 2-theta step size in deg
721
722    :returns: for fcj.pdf
723
724     * T = x*dx+t
725     * s = S/L+H/L
726     * if x < 0::
727
728        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|
729
730     * if x >= 0: fcj.pdf = 0   
731    """
732    def _pdf(self,x,t,s,dx):
733        T = dx*x+t
734        ax2 = abs(npcosd(T))
735        ax = ax2**2
736        bx = npcosd(t)**2
737        bx = np.where(ax>bx,bx,ax)
738        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
739        fx = np.where(fx > 0.,fx,0.0)
740        return fx
741             
742    def pdf(self,x,*args,**kwds):
743        loc=kwds['loc']
744        return self._pdf(x-loc,*args)
745       
746fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
747               
748def getWidthsCW(pos,sig,gam,shl):
749    '''Compute the peak widths used for computing the range of a peak
750    for constant wavelength data. On low-angle side, 50 FWHM are used,
751    on high-angle side 75 are used, low angle side extended for axial divergence
752    (for peaks above 90 deg, these are reversed.)
753    '''
754    widths = [np.sqrt(sig)/100.,gam/100.]
755    fwhm = 2.355*widths[0]+widths[1]
756    fmin = 50.*(fwhm+shl*abs(npcosd(pos)))
757    fmax = 75.0*fwhm
758    if pos > 90:
759        fmin,fmax = [fmax,fmin]         
760    return widths,fmin,fmax
761   
762def getWidthsTOF(pos,alp,bet,sig,gam):
763    '''Compute the peak widths used for computing the range of a peak
764    for constant wavelength data. 50 FWHM are used on both sides each
765    extended by exponential coeff.
766    '''
767    widths = [np.sqrt(sig),gam]
768    fwhm = 2.355*widths[0]+2.*widths[1]
769    fmin = 50.*fwhm*(1.+1./alp)   
770    fmax = 50.*fwhm*(1.+1./bet)
771    return widths,fmin,fmax
772   
773def getFWHM(pos,Inst):
774    '''Compute total FWHM from Thompson, Cox & Hastings (1987) , J. Appl. Cryst. 20, 79-83
775    via getgamFW(g,s).
776   
777    :param pos: float peak position in deg 2-theta or tof in musec
778    :param Inst: dict instrument parameters
779   
780    :returns float: total FWHM of pseudoVoigt in deg or musec
781    ''' 
782   
783    sig = lambda Th,U,V,W: np.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W))
784    sigTOF = lambda dsp,S0,S1,S2,Sq: np.sqrt(S0+S1*dsp**2+S2*dsp**4+Sq*dsp)
785    gam = lambda Th,X,Y,Z: Z+X/cosd(Th)+Y*tand(Th)
786    gamTOF = lambda dsp,X,Y,Z: Z+X*dsp+Y*dsp**2
787    alpTOF = lambda dsp,alp: alp/dsp
788    betTOF = lambda dsp,bet0,bet1,betq: bet0+bet1/dsp**4+betq/dsp**2
789    if 'C' in Inst['Type'][0]:
790        s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
791        g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
792        return getgamFW(g,s)/100.  #returns FWHM in deg
793    else:
794        dsp = pos/Inst['difC'][0]
795        alp = alpTOF(dsp,Inst['alpha'][0])
796        bet = betTOF(dsp,Inst['beta-0'][0],Inst['beta-1'][0],Inst['beta-q'][0])
797        s = sigTOF(dsp,Inst['sig-0'][1],Inst['sig-1'][1],Inst['sig-2'][1],Inst['sig-q'][1])
798        g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
799        return getgamFW(g,s)+np.log(2.0)*(alp+bet)/(alp*bet)
800   
801def getgamFW(g,s):
802    '''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
803    lambda fxn needs FWHM for both Gaussian & Lorentzian components
804   
805    :param g: float Lorentzian gamma = FWHM(L)
806    :param s: float Gaussian sig
807   
808    :returns float: total FWHM of pseudoVoigt
809    ''' 
810    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.)
811    return gamFW(2.35482*s,g)   #sqrt(8ln2)*sig = FWHM(G)
812               
813def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
814    '''Compute the Finger-Cox-Jepcoat modified Voigt function for a
815    CW powder peak by direct convolution. This version is not used.
816    '''
817    DX = xdata[1]-xdata[0]
818    widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
819    x = np.linspace(pos-fmin,pos+fmin,256)
820    dx = x[1]-x[0]
821    Norm = norm.pdf(x,loc=pos,scale=widths[0])
822    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
823    arg = [pos,shl/57.2958,dx,]
824    FCJ = fcjde.pdf(x,*arg,loc=pos)
825    if len(np.nonzero(FCJ)[0])>5:
826        z = np.column_stack([Norm,Cauchy,FCJ]).T
827        Z = fft.fft(z)
828        Df = fft.ifft(Z.prod(axis=0)).real
829    else:
830        z = np.column_stack([Norm,Cauchy]).T
831        Z = fft.fft(z)
832        Df = fft.fftshift(fft.ifft(Z.prod(axis=0))).real
833    Df /= np.sum(Df)
834    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
835    return intens*Df(xdata)*DX/dx
836
837def getBackground(pfx,parmDict,bakType,dataType,xdata,fixedBkg={}):
838    '''Computes the background from vars pulled from gpx file or tree.
839    '''
840    if 'T' in dataType:
841        q = 2.*np.pi*parmDict[pfx+'difC']/xdata
842    elif 'C' in dataType:
843        wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
844        q = npT2q(xdata,wave)
845    yb = np.zeros_like(xdata)
846    nBak = 0
847    cw = np.diff(xdata)
848    cw = np.append(cw,cw[-1])
849    sumBk = [0.,0.,0]
850    while True:
851        key = pfx+'Back;'+str(nBak)
852        if key in parmDict:
853            nBak += 1
854        else:
855            break
856#empirical functions
857    if bakType in ['chebyschev','cosine','chebyschev-1']:
858        dt = xdata[-1]-xdata[0]   
859        for iBak in range(nBak):
860            key = pfx+'Back;'+str(iBak)
861            if bakType == 'chebyschev':
862                ybi = parmDict[key]*(-1.+2.*(xdata-xdata[0])/dt)**iBak
863            elif bakType == 'chebyschev-1':
864                xpos = -1.+2.*(xdata-xdata[0])/dt
865                ybi = parmDict[key]*np.cos(iBak*np.arccos(xpos))
866            elif bakType == 'cosine':
867                ybi = parmDict[key]*npcosd(180.*xdata*iBak/xdata[-1])
868            yb += ybi
869        sumBk[0] = np.sum(yb)
870    elif bakType in ['Q^2 power series','Q^-2 power series']:
871        QT = 1.
872        yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
873        for iBak in range(nBak-1):
874            key = pfx+'Back;'+str(iBak+1)
875            if '-2' in bakType:
876                QT *= (iBak+1)*q**-2
877            else:
878                QT *= q**2/(iBak+1)
879            yb += QT*parmDict[key]
880        sumBk[0] = np.sum(yb)
881    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
882        if nBak == 1:
883            yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
884        elif nBak == 2:
885            dX = xdata[-1]-xdata[0]
886            T2 = (xdata-xdata[0])/dX
887            T1 = 1.0-T2
888            yb = parmDict[pfx+'Back;0']*T1+parmDict[pfx+'Back;1']*T2
889        else:
890            xnomask = ma.getdata(xdata)
891            xmin,xmax = xnomask[0],xnomask[-1]
892            if bakType == 'lin interpolate':
893                bakPos = np.linspace(xmin,xmax,nBak,True)
894            elif bakType == 'inv interpolate':
895                bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
896            elif bakType == 'log interpolate':
897                bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
898            bakPos[0] = xmin
899            bakPos[-1] = xmax
900            bakVals = np.zeros(nBak)
901            for i in range(nBak):
902                bakVals[i] = parmDict[pfx+'Back;'+str(i)]
903            bakInt = si.interp1d(bakPos,bakVals,'linear')
904            yb = bakInt(ma.getdata(xdata))
905        sumBk[0] = np.sum(yb)
906#Debye function       
907    if pfx+'difC' in parmDict:
908        ff = 1.
909    else:       
910        try:
911            wave = parmDict[pfx+'Lam']
912        except KeyError:
913            wave = parmDict[pfx+'Lam1']
914        SQ = (q/(4.*np.pi))**2
915        FF = G2elem.GetFormFactorCoeff('Si')[0]
916        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
917    iD = 0       
918    while True:
919        try:
920            dbA = parmDict[pfx+'DebyeA;'+str(iD)]
921            dbR = parmDict[pfx+'DebyeR;'+str(iD)]
922            dbU = parmDict[pfx+'DebyeU;'+str(iD)]
923            ybi = ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
924            yb += ybi
925            sumBk[1] += np.sum(ybi)
926            iD += 1       
927        except KeyError:
928            break
929#peaks
930    iD = 0
931    while True:
932        try:
933            pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
934            pkI = max(parmDict[pfx+'BkPkint;'+str(iD)],0.1)
935            pkS = max(parmDict[pfx+'BkPksig;'+str(iD)],1.)
936            pkG = max(parmDict[pfx+'BkPkgam;'+str(iD)],0.1)
937            if 'C' in dataType:
938                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
939            else: #'T'OF
940                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
941            iBeg = np.searchsorted(xdata,pkP-fmin)
942            iFin = np.searchsorted(xdata,pkP+fmax)
943            lenX = len(xdata)
944            if not iBeg:
945                iFin = np.searchsorted(xdata,pkP+fmax)
946            elif iBeg == lenX:
947                iFin = iBeg
948            else:
949                iFin = np.searchsorted(xdata,pkP+fmax)
950            if 'C' in dataType:
951                ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])
952                yb[iBeg:iFin] += ybi
953            else:   #'T'OF
954                ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
955                yb[iBeg:iFin] += ybi
956            sumBk[2] += np.sum(ybi)
957            iD += 1       
958        except KeyError:
959            break
960        except ValueError:
961            G2fil.G2Print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****')
962            break
963    # fixed background from file
964    if len(fixedBkg) >= 3:
965        mult = fixedBkg.get('_fixedMult',0.0)
966        if len(fixedBkg.get('_fixedValues',[])) != len(yb):
967            G2fil.G2Print('Lengths of backgrounds do not agree: yb={}, fixed={}'.format(
968                len(yb),len(fixedBkg.get('_fixedValues',[]))))
969        elif mult: 
970            yb -= mult*fixedBkg.get('_fixedValues',[]) # N.B. mult is negative
971            sumBk[0] = sum(yb)
972    return yb,sumBk
973   
974def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata):
975    'needs a doc string'
976    if 'T' in dataType:
977        q = 2.*np.pi*parmDict[hfx+'difC']/xdata
978    elif 'C' in dataType:
979        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
980        q = 2.*np.pi*npsind(xdata/2.)/wave
981    nBak = 0
982    while True:
983        key = hfx+'Back;'+str(nBak)
984        if key in parmDict:
985            nBak += 1
986        else:
987            break
988    dydb = np.zeros(shape=(nBak,len(xdata)))
989    dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
990    dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
991    cw = np.diff(xdata)
992    cw = np.append(cw,cw[-1])
993
994    if bakType in ['chebyschev','cosine','chebyschev-1']:
995        dt = xdata[-1]-xdata[0]   
996        for iBak in range(nBak):   
997            if bakType == 'chebyschev':
998                dydb[iBak] = (-1.+2.*(xdata-xdata[0])/dt)**iBak
999            elif bakType == 'chebyschev-1':
1000                xpos = -1.+2.*(xdata-xdata[0])/dt
1001                dydb[iBak] = np.cos(iBak*np.arccos(xpos))
1002            elif bakType == 'cosine':
1003                dydb[iBak] = npcosd(180.*xdata*iBak/xdata[-1])
1004    elif bakType in ['Q^2 power series','Q^-2 power series']:
1005        QT = 1.
1006        dydb[0] = np.ones_like(xdata)
1007        for iBak in range(nBak-1):
1008            if '-2' in bakType:
1009                QT *= (iBak+1)*q**-2
1010            else:
1011                QT *= q**2/(iBak+1)
1012            dydb[iBak+1] = QT
1013    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
1014        if nBak == 1:
1015            dydb[0] = np.ones_like(xdata)
1016        elif nBak == 2:
1017            dX = xdata[-1]-xdata[0]
1018            T2 = (xdata-xdata[0])/dX
1019            T1 = 1.0-T2
1020            dydb = [T1,T2]
1021        else:
1022            xnomask = ma.getdata(xdata)
1023            xmin,xmax = xnomask[0],xnomask[-1]
1024            if bakType == 'lin interpolate':
1025                bakPos = np.linspace(xmin,xmax,nBak,True)
1026            elif bakType == 'inv interpolate':
1027                bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
1028            elif bakType == 'log interpolate':
1029                bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
1030            bakPos[0] = xmin
1031            bakPos[-1] = xmax
1032            for i,pos in enumerate(bakPos):
1033                if i == 0:
1034                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
1035                elif i == len(bakPos)-1:
1036                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
1037                else:
1038                    dydb[i] = np.where(xdata>bakPos[i],
1039                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
1040                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
1041    if hfx+'difC' in parmDict:
1042        ff = 1.
1043    else:
1044        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1045        q = npT2q(xdata,wave)
1046        SQ = (q/(4*np.pi))**2
1047        FF = G2elem.GetFormFactorCoeff('Si')[0]
1048        ff = np.array(G2elem.ScatFac(FF,SQ)[0])*np.pi**2    #needs pi^2~10. for cw data (why?)
1049    iD = 0       
1050    while True:
1051        try:
1052            if hfx+'difC' in parmDict:
1053                q = 2*np.pi*parmDict[hfx+'difC']/xdata
1054            dbA = parmDict[hfx+'DebyeA;'+str(iD)]
1055            dbR = parmDict[hfx+'DebyeR;'+str(iD)]
1056            dbU = parmDict[hfx+'DebyeU;'+str(iD)]
1057            sqr = np.sin(q*dbR)/(q*dbR)
1058            cqr = np.cos(q*dbR)
1059            temp = np.exp(-dbU*q**2)
1060            dyddb[3*iD] = ff*sqr*temp
1061            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(dbR)
1062            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2
1063            iD += 1
1064        except KeyError:
1065            break
1066    iD = 0
1067    while True:
1068        try:
1069            pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
1070            pkI = max(parmDict[hfx+'BkPkint;'+str(iD)],0.1)
1071            pkS = max(parmDict[hfx+'BkPksig;'+str(iD)],1.0)
1072            pkG = max(parmDict[hfx+'BkPkgam;'+str(iD)],0.1)
1073            if 'C' in dataType:
1074                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
1075            else: #'T'OF
1076                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
1077            iBeg = np.searchsorted(xdata,pkP-fmin)
1078            iFin = np.searchsorted(xdata,pkP+fmax)
1079            lenX = len(xdata)
1080            if not iBeg:
1081                iFin = np.searchsorted(xdata,pkP+fmax)
1082            elif iBeg == lenX:
1083                iFin = iBeg
1084            else:
1085                iFin = np.searchsorted(xdata,pkP+fmax)
1086            if 'C' in dataType:
1087                Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
1088                dydpk[4*iD][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdp
1089                dydpk[4*iD+1][iBeg:iFin] += 100.*cw[iBeg:iFin]*Df
1090                dydpk[4*iD+2][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFds
1091                dydpk[4*iD+3][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdg
1092            else:   #'T'OF
1093                Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
1094                dydpk[4*iD][iBeg:iFin] += pkI*dFdp
1095                dydpk[4*iD+1][iBeg:iFin] += Df
1096                dydpk[4*iD+2][iBeg:iFin] += pkI*dFds
1097                dydpk[4*iD+3][iBeg:iFin] += pkI*dFdg
1098            iD += 1       
1099        except KeyError:
1100            break
1101        except ValueError:
1102            G2fil.G2Print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****')
1103            break       
1104    return dydb,dyddb,dydpk
1105
1106#use old fortran routine
1107def getFCJVoigt3(pos,sig,gam,shl,xdata):
1108    '''Compute the Finger-Cox-Jepcoat modified Pseudo-Voigt function for a
1109    CW powder peak in external Fortran routine
1110    '''
1111    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1112#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
1113    Df /= np.sum(Df)
1114    return Df
1115
1116def getdFCJVoigt3(pos,sig,gam,shl,xdata):
1117    '''Compute analytic derivatives the Finger-Cox-Jepcoat modified Pseudo-Voigt
1118    function for a CW powder peak
1119    '''
1120    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1121#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
1122    return Df,dFdp,dFds,dFdg,dFdsh
1123
1124def getPsVoigt(pos,sig,gam,xdata):
1125    'needs a doc string'
1126   
1127    Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
1128    Df /= np.sum(Df)
1129    return Df
1130
1131def getdPsVoigt(pos,sig,gam,xdata):
1132    'needs a doc string'
1133   
1134    Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
1135    return Df,dFdp,dFds,dFdg
1136
1137def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
1138    'needs a doc string'
1139    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1140    Df /= np.sum(Df)
1141    return Df 
1142   
1143def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
1144    'needs a doc string'
1145    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1146    return Df,dFdp,dFda,dFdb,dFds,dFdg   
1147
1148def ellipseSize(H,Sij,GB):
1149    'Implements r=1/sqrt(sum((1/S)*(q.v)^2) per note from Alexander Brady'
1150    HX = np.inner(H.T,GB)
1151    lenHX = np.sqrt(np.sum(HX**2))
1152    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
1153    R = np.inner(HX/lenHX,Rsize)**2*Esize         #want column length for hkl in crystal
1154    lenR = 1./np.sqrt(np.sum(R))
1155    return lenR
1156
1157def ellipseSizeDerv(H,Sij,GB):
1158    'needs a doc string'
1159    lenR = ellipseSize(H,Sij,GB)
1160    delt = 0.001
1161    dRdS = np.zeros(6)
1162    for i in range(6):
1163        Sij[i] -= delt
1164        lenM = ellipseSize(H,Sij,GB)
1165        Sij[i] += 2.*delt
1166        lenP = ellipseSize(H,Sij,GB)
1167        Sij[i] -= delt
1168        dRdS[i] = (lenP-lenM)/(2.*delt)
1169    return lenR,dRdS
1170
1171def getHKLpeak(dmin,SGData,A,Inst=None,nodup=False):
1172    '''
1173    Generates allowed by symmetry reflections with d >= dmin
1174    NB: GenHKLf & checkMagextc return True for extinct reflections
1175
1176    :param dmin:  minimum d-spacing
1177    :param SGData: space group data obtained from SpcGroup
1178    :param A: lattice parameter terms A1-A6
1179    :param Inst: instrument parameter info
1180    :returns: HKLs: np.array hkl, etc for allowed reflections
1181
1182    '''
1183    HKL = G2lat.GenHLaue(dmin,SGData,A)       
1184    HKLs = []
1185    ds = []
1186    for h,k,l,d in HKL:
1187        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1188        if ext and 'MagSpGrp' in SGData:
1189            ext = G2spc.checkMagextc([h,k,l],SGData)
1190        if not ext:
1191            if nodup and int(10000*d) in ds:
1192                continue
1193            ds.append(int(10000*d))
1194            if Inst == None:
1195                HKLs.append([h,k,l,d,0,-1])
1196            else:
1197                HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1198    return np.array(HKLs)
1199
1200def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1201    'needs a doc string'
1202    HKLs = []
1203    vec = np.array(Vec)
1204    vstar = np.sqrt(G2lat.calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1205    dvec = 1./(maxH*vstar+1./dmin)
1206    HKL = G2lat.GenHLaue(dvec,SGData,A)       
1207    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1208    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1209    ifMag = False
1210    if 'MagSpGrp' in SGData:
1211        ifMag = True
1212    for h,k,l,d in HKL:
1213        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1214        if not ext and d >= dmin:
1215            HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1216        for dH in SSdH:
1217            if dH:
1218                DH = SSdH[dH]
1219                H = [h+DH[0],k+DH[1],l+DH[2]]
1220                d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1221                if d >= dmin:
1222                    HKLM = np.array([h,k,l,dH])
1223                    if G2spc.checkSSextc(HKLM,SSGData) or ifMag:
1224                        HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])   
1225    return G2lat.sortHKLd(HKLs,True,True,True)
1226
1227def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
1228    'Computes the profile for a powder pattern'
1229   
1230    yb = getBackground('',parmDict,bakType,dataType,xdata)[0]
1231    yc = np.zeros_like(yb)
1232    cw = np.diff(xdata)
1233    cw = np.append(cw,cw[-1])
1234    if 'C' in dataType:
1235        shl = max(parmDict['SH/L'],0.002)
1236        Ka2 = False
1237        if 'Lam1' in parmDict.keys():
1238            Ka2 = True
1239            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1240            kRatio = parmDict['I(L2)/I(L1)']
1241        iPeak = 0
1242        while True:
1243            try:
1244                pos = parmDict['pos'+str(iPeak)]
1245                tth = (pos-parmDict['Zero'])
1246                intens = parmDict['int'+str(iPeak)]
1247                sigName = 'sig'+str(iPeak)
1248                if sigName in varyList:
1249                    sig = parmDict[sigName]
1250                else:
1251                    sig = G2mth.getCWsig(parmDict,tth)
1252                sig = max(sig,0.001)          #avoid neg sigma^2
1253                gamName = 'gam'+str(iPeak)
1254                if gamName in varyList:
1255                    gam = parmDict[gamName]
1256                else:
1257                    gam = G2mth.getCWgam(parmDict,tth)
1258                gam = max(gam,0.001)             #avoid neg gamma
1259                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1260                iBeg = np.searchsorted(xdata,pos-fmin)
1261                iFin = np.searchsorted(xdata,pos+fmin)
1262                if not iBeg+iFin:       #peak below low limit
1263                    iPeak += 1
1264                    continue
1265                elif not iBeg-iFin:     #peak above high limit
1266                    return yb+yc
1267                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1268                if Ka2:
1269                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1270                    iBeg = np.searchsorted(xdata,pos2-fmin)
1271                    iFin = np.searchsorted(xdata,pos2+fmin)
1272                    if iBeg-iFin:
1273                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1274                iPeak += 1
1275            except KeyError:        #no more peaks to process
1276                return yb+yc
1277    else:
1278        Pdabc = parmDict['Pdabc']
1279        difC = parmDict['difC']
1280        iPeak = 0
1281        while True:
1282            try:
1283                pos = parmDict['pos'+str(iPeak)]               
1284                tof = pos-parmDict['Zero']
1285                dsp = tof/difC
1286                intens = parmDict['int'+str(iPeak)]
1287                alpName = 'alp'+str(iPeak)
1288                if alpName in varyList:
1289                    alp = parmDict[alpName]
1290                else:
1291                    if len(Pdabc):
1292                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1293                    else:
1294                        alp = G2mth.getTOFalpha(parmDict,dsp)
1295                alp = max(0.1,alp)
1296                betName = 'bet'+str(iPeak)
1297                if betName in varyList:
1298                    bet = parmDict[betName]
1299                else:
1300                    if len(Pdabc):
1301                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1302                    else:
1303                        bet = G2mth.getTOFbeta(parmDict,dsp)
1304                bet = max(0.0001,bet)
1305                sigName = 'sig'+str(iPeak)
1306                if sigName in varyList:
1307                    sig = parmDict[sigName]
1308                else:
1309                    sig = G2mth.getTOFsig(parmDict,dsp)
1310                gamName = 'gam'+str(iPeak)
1311                if gamName in varyList:
1312                    gam = parmDict[gamName]
1313                else:
1314                    gam = G2mth.getTOFgamma(parmDict,dsp)
1315                gam = max(gam,0.001)             #avoid neg gamma
1316                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1317                iBeg = np.searchsorted(xdata,pos-fmin)
1318                iFin = np.searchsorted(xdata,pos+fmax)
1319                lenX = len(xdata)
1320                if not iBeg:
1321                    iFin = np.searchsorted(xdata,pos+fmax)
1322                elif iBeg == lenX:
1323                    iFin = iBeg
1324                else:
1325                    iFin = np.searchsorted(xdata,pos+fmax)
1326                if not iBeg+iFin:       #peak below low limit
1327                    iPeak += 1
1328                    continue
1329                elif not iBeg-iFin:     #peak above high limit
1330                    return yb+yc
1331                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1332                iPeak += 1
1333            except KeyError:        #no more peaks to process
1334                return yb+yc
1335           
1336def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
1337    'needs a doc string'
1338# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1339    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1340    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,dataType,xdata)
1341    if 'Back;0' in varyList:            #background derivs are in front if present
1342        dMdv[0:len(dMdb)] = dMdb
1343    names = ['DebyeA','DebyeR','DebyeU']
1344    for name in varyList:
1345        if 'Debye' in name:
1346            parm,Id = name.split(';')
1347            ip = names.index(parm)
1348            dMdv[varyList.index(name)] = dMddb[3*int(Id)+ip]
1349    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1350    for name in varyList:
1351        if 'BkPk' in name:
1352            parm,Id = name.split(';')
1353            ip = names.index(parm)
1354            dMdv[varyList.index(name)] = dMdpk[4*int(Id)+ip]
1355    cw = np.diff(xdata)
1356    cw = np.append(cw,cw[-1])
1357    if 'C' in dataType:
1358        shl = max(parmDict['SH/L'],0.002)
1359        Ka2 = False
1360        if 'Lam1' in parmDict.keys():
1361            Ka2 = True
1362            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1363            kRatio = parmDict['I(L2)/I(L1)']
1364        iPeak = 0
1365        while True:
1366            try:
1367                pos = parmDict['pos'+str(iPeak)]
1368                tth = (pos-parmDict['Zero'])
1369                intens = parmDict['int'+str(iPeak)]
1370                sigName = 'sig'+str(iPeak)
1371                if sigName in varyList:
1372                    sig = parmDict[sigName]
1373                    dsdU = dsdV = dsdW = 0
1374                else:
1375                    sig = G2mth.getCWsig(parmDict,tth)
1376                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1377                sig = max(sig,0.001)          #avoid neg sigma
1378                gamName = 'gam'+str(iPeak)
1379                if gamName in varyList:
1380                    gam = parmDict[gamName]
1381                    dgdX = dgdY = dgdZ = 0
1382                else:
1383                    gam = G2mth.getCWgam(parmDict,tth)
1384                    dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1385                gam = max(gam,0.001)             #avoid neg gamma
1386                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1387                iBeg = np.searchsorted(xdata,pos-fmin)
1388                iFin = np.searchsorted(xdata,pos+fmin)
1389                if not iBeg+iFin:       #peak below low limit
1390                    iPeak += 1
1391                    continue
1392                elif not iBeg-iFin:     #peak above high limit
1393                    break
1394                dMdpk = np.zeros(shape=(6,len(xdata)))
1395                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1396                for i in range(1,5):
1397                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
1398                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
1399                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1400                if Ka2:
1401                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1402                    iBeg = np.searchsorted(xdata,pos2-fmin)
1403                    iFin = np.searchsorted(xdata,pos2+fmin)
1404                    if iBeg-iFin:
1405                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1406                        for i in range(1,5):
1407                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
1408                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
1409                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
1410                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1411                for parmName in ['pos','int','sig','gam']:
1412                    try:
1413                        idx = varyList.index(parmName+str(iPeak))
1414                        dMdv[idx] = dervDict[parmName]
1415                    except ValueError:
1416                        pass
1417                if 'U' in varyList:
1418                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1419                if 'V' in varyList:
1420                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1421                if 'W' in varyList:
1422                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1423                if 'X' in varyList:
1424                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1425                if 'Y' in varyList:
1426                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1427                if 'Z' in varyList:
1428                    dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1429                if 'SH/L' in varyList:
1430                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
1431                if 'I(L2)/I(L1)' in varyList:
1432                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1433                iPeak += 1
1434            except KeyError:        #no more peaks to process
1435                break
1436    else:
1437        Pdabc = parmDict['Pdabc']
1438        difC = parmDict['difC']
1439        iPeak = 0
1440        while True:
1441            try:
1442                pos = parmDict['pos'+str(iPeak)]               
1443                tof = pos-parmDict['Zero']
1444                dsp = tof/difC
1445                intens = parmDict['int'+str(iPeak)]
1446                alpName = 'alp'+str(iPeak)
1447                if alpName in varyList:
1448                    alp = parmDict[alpName]
1449                else:
1450                    if len(Pdabc):
1451                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1452                        dada0 = 0
1453                    else:
1454                        alp = G2mth.getTOFalpha(parmDict,dsp)
1455                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1456                betName = 'bet'+str(iPeak)
1457                if betName in varyList:
1458                    bet = parmDict[betName]
1459                else:
1460                    if len(Pdabc):
1461                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1462                        dbdb0 = dbdb1 = dbdb2 = 0
1463                    else:
1464                        bet = G2mth.getTOFbeta(parmDict,dsp)
1465                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1466                sigName = 'sig'+str(iPeak)
1467                if sigName in varyList:
1468                    sig = parmDict[sigName]
1469                    dsds0 = dsds1 = dsds2 = dsds3 = 0
1470                else:
1471                    sig = G2mth.getTOFsig(parmDict,dsp)
1472                    dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1473                gamName = 'gam'+str(iPeak)
1474                if gamName in varyList:
1475                    gam = parmDict[gamName]
1476                    dsdX = dsdY = dsdZ = 0
1477                else:
1478                    gam = G2mth.getTOFgamma(parmDict,dsp)
1479                    dsdX,dsdY,dsdZ = G2mth.getTOFgammaDeriv(dsp)
1480                gam = max(gam,0.001)             #avoid neg gamma
1481                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1482                iBeg = np.searchsorted(xdata,pos-fmin)
1483                lenX = len(xdata)
1484                if not iBeg:
1485                    iFin = np.searchsorted(xdata,pos+fmax)
1486                elif iBeg == lenX:
1487                    iFin = iBeg
1488                else:
1489                    iFin = np.searchsorted(xdata,pos+fmax)
1490                if not iBeg+iFin:       #peak below low limit
1491                    iPeak += 1
1492                    continue
1493                elif not iBeg-iFin:     #peak above high limit
1494                    break
1495                dMdpk = np.zeros(shape=(7,len(xdata)))
1496                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1497                for i in range(1,6):
1498                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1499                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1500                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1501                for parmName in ['pos','int','alp','bet','sig','gam']:
1502                    try:
1503                        idx = varyList.index(parmName+str(iPeak))
1504                        dMdv[idx] = dervDict[parmName]
1505                    except ValueError:
1506                        pass
1507                if 'alpha' in varyList:
1508                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1509                if 'beta-0' in varyList:
1510                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1511                if 'beta-1' in varyList:
1512                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1513                if 'beta-q' in varyList:
1514                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1515                if 'sig-0' in varyList:
1516                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1517                if 'sig-1' in varyList:
1518                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1519                if 'sig-2' in varyList:
1520                    dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1521                if 'sig-q' in varyList:
1522                    dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1523                if 'X' in varyList:
1524                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1525                if 'Y' in varyList:
1526                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']
1527                if 'Z' in varyList:
1528                    dMdv[varyList.index('Z')] += dsdZ*dervDict['gam']
1529                iPeak += 1
1530            except KeyError:        #no more peaks to process
1531                break
1532    return dMdv
1533       
1534def Dict2Values(parmdict, varylist):
1535    '''Use before call to leastsq to setup list of values for the parameters
1536    in parmdict, as selected by key in varylist'''
1537    return [parmdict[key] for key in varylist] 
1538   
1539def Values2Dict(parmdict, varylist, values):
1540    ''' Use after call to leastsq to update the parameter dictionary with
1541    values corresponding to keys in varylist'''
1542    parmdict.update(zip(varylist,values))
1543   
1544def SetBackgroundParms(Background):
1545    'Loads background parameters into dicts/lists to create varylist & parmdict'
1546    if len(Background) == 1:            # fix up old backgrounds
1547        Background.append({'nDebye':0,'debyeTerms':[]})
1548    bakType,bakFlag = Background[0][:2]
1549    backVals = Background[0][3:]
1550    backNames = ['Back;'+str(i) for i in range(len(backVals))]
1551    Debye = Background[1]           #also has background peaks stuff
1552    backDict = dict(zip(backNames,backVals))
1553    backVary = []
1554    if bakFlag:
1555        backVary = backNames
1556
1557    backDict['nDebye'] = Debye['nDebye']
1558    debyeDict = {}
1559    debyeList = []
1560    for i in range(Debye['nDebye']):
1561        debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1562        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1563        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1564    debyeVary = []
1565    for item in debyeList:
1566        if item[1]:
1567            debyeVary.append(item[0])
1568    backDict.update(debyeDict)
1569    backVary += debyeVary
1570
1571    backDict['nPeaks'] = Debye['nPeaks']
1572    peaksDict = {}
1573    peaksList = []
1574    for i in range(Debye['nPeaks']):
1575        peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1576        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1577        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1578    peaksVary = []
1579    for item in peaksList:
1580        if item[1]:
1581            peaksVary.append(item[0])
1582    backDict.update(peaksDict)
1583    backVary += peaksVary
1584    return bakType,backDict,backVary
1585   
1586def DoCalibInst(IndexPeaks,Inst):
1587   
1588    def SetInstParms():
1589        dataType = Inst['Type'][0]
1590        insVary = []
1591        insNames = []
1592        insVals = []
1593        for parm in Inst:
1594            insNames.append(parm)
1595            insVals.append(Inst[parm][1])
1596            if parm in ['Lam','difC','difA','difB','Zero',]:
1597                if Inst[parm][2]:
1598                    insVary.append(parm)
1599        instDict = dict(zip(insNames,insVals))
1600        return dataType,instDict,insVary
1601       
1602    def GetInstParms(parmDict,Inst,varyList):
1603        for name in Inst:
1604            Inst[name][1] = parmDict[name]
1605       
1606    def InstPrint(Inst,sigDict):
1607        print ('Instrument Parameters:')
1608        if 'C' in Inst['Type'][0]:
1609            ptfmt = "%12.6f"
1610        else:
1611            ptfmt = "%12.3f"
1612        ptlbls = 'names :'
1613        ptstr =  'values:'
1614        sigstr = 'esds  :'
1615        for parm in Inst:
1616            if parm in  ['Lam','difC','difA','difB','Zero',]:
1617                ptlbls += "%s" % (parm.center(12))
1618                ptstr += ptfmt % (Inst[parm][1])
1619                if parm in sigDict:
1620                    sigstr += ptfmt % (sigDict[parm])
1621                else:
1622                    sigstr += 12*' '
1623        print (ptlbls)
1624        print (ptstr)
1625        print (sigstr)
1626       
1627    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1628        parmDict.update(zip(varyList,values))
1629        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1630
1631    peakPos = []
1632    peakDsp = []
1633    peakWt = []
1634    for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1635        if peak[2] and peak[3] and sig > 0.:
1636            peakPos.append(peak[0])
1637            peakDsp.append(peak[-1])    #d-calc
1638#            peakWt.append(peak[-1]**2/sig**2)   #weight by d**2
1639            peakWt.append(1./(sig*peak[-1]))   #
1640    peakPos = np.array(peakPos)
1641    peakDsp = np.array(peakDsp)
1642    peakWt = np.array(peakWt)
1643    dataType,insDict,insVary = SetInstParms()
1644    parmDict = {}
1645    parmDict.update(insDict)
1646    varyList = insVary
1647    if not len(varyList):
1648        G2fil.G2Print ('**** ERROR - nothing to refine! ****')
1649        return False
1650    while True:
1651        begin = time.time()
1652        values =  np.array(Dict2Values(parmDict, varyList))
1653        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1654            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1655        ncyc = int(result[2]['nfev']/2)
1656        runtime = time.time()-begin   
1657        chisq = np.sum(result[2]['fvec']**2)
1658        Values2Dict(parmDict, varyList, result[0])
1659        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
1660        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(peakPos),len(varyList)))
1661        G2fil.G2Print ('calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
1662        G2fil.G2Print ('chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF))
1663        try:
1664            sig = np.sqrt(np.diag(result[1])*GOF)
1665            if np.any(np.isnan(sig)):
1666                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***')
1667            break                   #refinement succeeded - finish up!
1668        except ValueError:          #result[1] is None on singular matrix
1669            G2fil.G2Print ('**** Refinement failed - singular matrix ****')
1670       
1671    sigDict = dict(zip(varyList,sig))
1672    GetInstParms(parmDict,Inst,varyList)
1673    InstPrint(Inst,sigDict)
1674    return True
1675           
1676def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
1677    '''Called to perform a peak fit, refining the selected items in the peak
1678    table as well as selected items in the background.
1679
1680    :param str FitPgm: type of fit to perform. At present this is ignored.
1681    :param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1682      four values followed by a refine flag where the values are: position, intensity,
1683      sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1684      tree entry, dict item "peaks"
1685    :param list Background: describes the background. List with two items.
1686      Item 0 specifies a background model and coefficients. Item 1 is a dict.
1687      From the Histogram/Background tree entry.
1688    :param list Limits: min and max x-value to use
1689    :param dict Inst: Instrument parameters
1690    :param dict Inst2: more Instrument parameters
1691    :param numpy.array data: a 5xn array. data[0] is the x-values,
1692      data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1693      calc, background and difference intensities, respectively.
1694    :param array fixback: fixed background values
1695    :param list prevVaryList: Used in sequential refinements to override the
1696      variable list. Defaults as an empty list.
1697    :param bool oneCycle: True if only one cycle of fitting should be performed
1698    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1699      and derivType = controls['deriv type']. If None default values are used.
1700    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1701      Defaults to None, which means no updates are done.
1702    '''
1703    def GetBackgroundParms(parmList,Background):
1704        iBak = 0
1705        while True:
1706            try:
1707                bakName = 'Back;'+str(iBak)
1708                Background[0][iBak+3] = parmList[bakName]
1709                iBak += 1
1710            except KeyError:
1711                break
1712        iDb = 0
1713        while True:
1714            names = ['DebyeA;','DebyeR;','DebyeU;']
1715            try:
1716                for i,name in enumerate(names):
1717                    val = parmList[name+str(iDb)]
1718                    Background[1]['debyeTerms'][iDb][2*i] = val
1719                iDb += 1
1720            except KeyError:
1721                break
1722        iDb = 0
1723        while True:
1724            names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1725            try:
1726                for i,name in enumerate(names):
1727                    val = parmList[name+str(iDb)]
1728                    Background[1]['peaksList'][iDb][2*i] = val
1729                iDb += 1
1730            except KeyError:
1731                break
1732               
1733    def BackgroundPrint(Background,sigDict):
1734        print ('Background coefficients for '+Background[0][0]+' function')
1735        ptfmt = "%12.5f"
1736        ptstr =  'value: '
1737        sigstr = 'esd  : '
1738        for i,back in enumerate(Background[0][3:]):
1739            ptstr += ptfmt % (back)
1740            if Background[0][1]:
1741                prm = 'Back;'+str(i)
1742                if prm in sigDict:
1743                    sigstr += ptfmt % (sigDict[prm])
1744                else:
1745                    sigstr += " "*12
1746            if len(ptstr) > 75:
1747                print (ptstr)
1748                if Background[0][1]: print (sigstr)
1749                ptstr =  'value: '
1750                sigstr = 'esd  : '
1751        if len(ptstr) > 8:
1752            print (ptstr)
1753            if Background[0][1]: print (sigstr)
1754
1755        if Background[1]['nDebye']:
1756            parms = ['DebyeA;','DebyeR;','DebyeU;']
1757            print ('Debye diffuse scattering coefficients')
1758            ptfmt = "%12.5f"
1759            print (' term       DebyeA       esd        DebyeR       esd        DebyeU        esd')
1760            for term in range(Background[1]['nDebye']):
1761                line = ' term %d'%(term)
1762                for ip,name in enumerate(parms):
1763                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1764                    if name+str(term) in sigDict:
1765                        line += ptfmt%(sigDict[name+str(term)])
1766                    else:
1767                        line += " "*12
1768                print (line)
1769        if Background[1]['nPeaks']:
1770            print ('Coefficients for Background Peaks')
1771            ptfmt = "%15.3f"
1772            for j,pl in enumerate(Background[1]['peaksList']):
1773                names =  'peak %3d:'%(j+1)
1774                ptstr =  'values  :'
1775                sigstr = 'esds    :'
1776                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1777                    val = pl[2*i]
1778                    prm = lbl+";"+str(j)
1779                    names += '%15s'%(prm)
1780                    ptstr += ptfmt%(val)
1781                    if prm in sigDict:
1782                        sigstr += ptfmt%(sigDict[prm])
1783                    else:
1784                        sigstr += " "*15
1785                print (names)
1786                print (ptstr)
1787                print (sigstr)
1788                           
1789    def SetInstParms(Inst):
1790        dataType = Inst['Type'][0]
1791        insVary = []
1792        insNames = []
1793        insVals = []
1794        for parm in Inst:
1795            insNames.append(parm)
1796            insVals.append(Inst[parm][1])
1797            if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
1798                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1799                    insVary.append(parm)
1800        instDict = dict(zip(insNames,insVals))
1801#        instDict['X'] = max(instDict['X'],0.01)
1802#        instDict['Y'] = max(instDict['Y'],0.01)
1803        if 'SH/L' in instDict:
1804            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1805        return dataType,instDict,insVary
1806       
1807    def GetInstParms(parmDict,Inst,varyList,Peaks):
1808        for name in Inst:
1809            Inst[name][1] = parmDict[name]
1810        iPeak = 0
1811        while True:
1812            try:
1813                sigName = 'sig'+str(iPeak)
1814                pos = parmDict['pos'+str(iPeak)]
1815                if sigName not in varyList:
1816                    if 'C' in Inst['Type'][0]:
1817                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1818                    else:
1819                        dsp = G2lat.Pos2dsp(Inst,pos)
1820                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1821                gamName = 'gam'+str(iPeak)
1822                if gamName not in varyList:
1823                    if 'C' in Inst['Type'][0]:
1824                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1825                    else:
1826                        dsp = G2lat.Pos2dsp(Inst,pos)
1827                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1828                iPeak += 1
1829            except KeyError:
1830                break
1831       
1832    def InstPrint(Inst,sigDict):
1833        print ('Instrument Parameters:')
1834        ptfmt = "%12.6f"
1835        ptlbls = 'names :'
1836        ptstr =  'values:'
1837        sigstr = 'esds  :'
1838        for parm in Inst:
1839            if parm in  ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
1840                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1841                ptlbls += "%s" % (parm.center(12))
1842                ptstr += ptfmt % (Inst[parm][1])
1843                if parm in sigDict:
1844                    sigstr += ptfmt % (sigDict[parm])
1845                else:
1846                    sigstr += 12*' '
1847        print (ptlbls)
1848        print (ptstr)
1849        print (sigstr)
1850
1851    def SetPeaksParms(dataType,Peaks):
1852        peakNames = []
1853        peakVary = []
1854        peakVals = []
1855        if 'C' in dataType:
1856            names = ['pos','int','sig','gam']
1857        else:
1858            names = ['pos','int','alp','bet','sig','gam']
1859        for i,peak in enumerate(Peaks):
1860            for j,name in enumerate(names):
1861                peakVals.append(peak[2*j])
1862                parName = name+str(i)
1863                peakNames.append(parName)
1864                if peak[2*j+1]:
1865                    peakVary.append(parName)
1866        return dict(zip(peakNames,peakVals)),peakVary
1867               
1868    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1869        if 'C' in Inst['Type'][0]:
1870            names = ['pos','int','sig','gam']
1871        else:   #'T'
1872            names = ['pos','int','alp','bet','sig','gam']
1873        for i,peak in enumerate(Peaks):
1874            pos = parmDict['pos'+str(i)]
1875            if 'difC' in Inst:
1876                dsp = pos/Inst['difC'][1]
1877            for j in range(len(names)):
1878                parName = names[j]+str(i)
1879                if parName in varyList:
1880                    peak[2*j] = parmDict[parName]
1881                elif 'alpha' in parName:
1882                    peak[2*j] = parmDict['alpha']/dsp
1883                elif 'beta' in parName:
1884                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1885                elif 'sig' in parName:
1886                    if 'C' in Inst['Type'][0]:
1887                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1888                    else:
1889                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1890                elif 'gam' in parName:
1891                    if 'C' in Inst['Type'][0]:
1892                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1893                    else:
1894                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1895                       
1896    def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
1897        print ('Peak coefficients:')
1898        if 'C' in dataType:
1899            names = ['pos','int','sig','gam']
1900        else:   #'T'
1901            names = ['pos','int','alp','bet','sig','gam']           
1902        head = 13*' '
1903        for name in names:
1904            if name in ['alp','bet']:
1905                head += name.center(8)+'esd'.center(8)
1906            else:
1907                head += name.center(10)+'esd'.center(10)
1908        head += 'bins'.center(8)
1909        print (head)
1910        if 'C' in dataType:
1911            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1912        else:
1913            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1914        for i,peak in enumerate(Peaks):
1915            ptstr =  ':'
1916            for j in range(len(names)):
1917                name = names[j]
1918                parName = name+str(i)
1919                ptstr += ptfmt[name] % (parmDict[parName])
1920                if parName in varyList:
1921                    ptstr += ptfmt[name] % (sigDict[parName])
1922                else:
1923                    if name in ['alp','bet']:
1924                        ptstr += 8*' '
1925                    else:
1926                        ptstr += 10*' '
1927            ptstr += '%9.2f'%(ptsperFW[i])
1928            print ('%s'%(('Peak'+str(i+1)).center(8)),ptstr)
1929               
1930    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1931        parmdict.update(zip(varylist,values))
1932        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1933           
1934    def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):       
1935        parmdict.update(zip(varylist,values))
1936        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1937        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1938        if dlg:
1939            dlg.Raise()
1940            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1941            if not GoOn:
1942                return -M           #abort!!
1943        return M
1944       
1945    if controls:
1946        Ftol = controls['min dM/M']
1947    else:
1948        Ftol = 0.0001
1949    if oneCycle:
1950        Ftol = 1.0
1951    x,y,w,yc,yb,yd = data   #these are numpy arrays - remove masks!
1952    if fixback is None:
1953        fixback = np.zeros_like(y)
1954    yc *= 0.                            #set calcd ones to zero
1955    yb *= 0.
1956    yd *= 0.
1957    cw = x[1:]-x[:-1]
1958    xBeg = np.searchsorted(x,Limits[0])
1959    xFin = np.searchsorted(x,Limits[1])+1
1960    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1961    dataType,insDict,insVary = SetInstParms(Inst)
1962    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1963    parmDict = {}
1964    parmDict.update(bakDict)
1965    parmDict.update(insDict)
1966    parmDict.update(peakDict)
1967    parmDict['Pdabc'] = []      #dummy Pdabc
1968    parmDict.update(Inst2)      #put in real one if there
1969    if prevVaryList:
1970        varyList = prevVaryList[:]
1971    else:
1972        varyList = bakVary+insVary+peakVary
1973    fullvaryList = varyList[:]
1974    while True:
1975        begin = time.time()
1976        values =  np.array(Dict2Values(parmDict, varyList))
1977        Rvals = {}
1978        badVary = []
1979        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1980               args=(x[xBeg:xFin],(y+fixback)[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1981        ncyc = int(result[2]['nfev']/2)
1982        runtime = time.time()-begin   
1983        chisq = np.sum(result[2]['fvec']**2)
1984        Values2Dict(parmDict, varyList, result[0])
1985        Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*(y+fixback)[xBeg:xFin]**2))*100.      #to %
1986        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1987        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList)))
1988        if ncyc:
1989            G2fil.G2Print ('fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
1990        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
1991        sig = [0]*len(varyList)
1992        if len(varyList) == 0: break  # if nothing was refined
1993        try:
1994            sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1995            if np.any(np.isnan(sig)):
1996                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***')
1997            break                   #refinement succeeded - finish up!
1998        except ValueError:          #result[1] is None on singular matrix
1999            G2fil.G2Print ('**** Refinement failed - singular matrix ****')
2000            Ipvt = result[2]['ipvt']
2001            for i,ipvt in enumerate(Ipvt):
2002                if not np.sum(result[2]['fjac'],axis=1)[i]:
2003                    G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
2004                    badVary.append(varyList[ipvt-1])
2005                    del(varyList[ipvt-1])
2006                    break
2007            else: # nothing removed
2008                break
2009    if dlg: dlg.Destroy()
2010    sigDict = dict(zip(varyList,sig))
2011    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]-fixback[xBeg:xFin]
2012    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)-fixback[xBeg:xFin]
2013    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
2014    GetBackgroundParms(parmDict,Background)
2015    if bakVary: BackgroundPrint(Background,sigDict)
2016    GetInstParms(parmDict,Inst,varyList,Peaks)
2017    if insVary: InstPrint(Inst,sigDict)
2018    GetPeaksParms(Inst,parmDict,Peaks,varyList)
2019    binsperFWHM = []
2020    for peak in Peaks:
2021        FWHM = getFWHM(peak[0],Inst)
2022        try:
2023            binsperFWHM.append(FWHM/cw[x.searchsorted(peak[0])])
2024        except IndexError:
2025            binsperFWHM.append(0.)
2026    if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
2027    if len(binsperFWHM):
2028        if min(binsperFWHM) < 1.:
2029            G2fil.G2Print ('*** Warning: calculated peak widths are too narrow to refine profile coefficients ***')
2030            if 'T' in Inst['Type'][0]:
2031                G2fil.G2Print (' Manually increase sig-0, 1, or 2 in Instrument Parameters')
2032            else:
2033                G2fil.G2Print (' Manually increase W in Instrument Parameters')
2034        elif min(binsperFWHM) < 4.:
2035            G2fil.G2Print ('*** Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement ***')
2036            G2fil.G2Print ('*** recommended is 6-10; you have %.2f ***'%(min(binsperFWHM)))
2037    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
2038   
2039def calcIncident(Iparm,xdata):
2040    'needs a doc string'
2041
2042    def IfunAdv(Iparm,xdata):
2043        Itype = Iparm['Itype']
2044        Icoef = Iparm['Icoeff']
2045        DYI = np.ones((12,xdata.shape[0]))
2046        YI = np.ones_like(xdata)*Icoef[0]
2047       
2048        x = xdata/1000.                 #expressions are in ms
2049        if Itype == 'Exponential':
2050            for i in [1,3,5,7,9]:
2051                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
2052                YI += Icoef[i]*Eterm
2053                DYI[i] *= Eterm
2054                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
2055        elif 'Maxwell'in Itype:
2056            Eterm = np.exp(-Icoef[2]/x**2)
2057            DYI[1] = Eterm/x**5
2058            DYI[2] = -Icoef[1]*DYI[1]/x**2
2059            YI += (Icoef[1]*Eterm/x**5)
2060            if 'Exponential' in Itype:
2061                for i in range(3,11,2):
2062                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
2063                    YI += Icoef[i]*Eterm
2064                    DYI[i] *= Eterm
2065                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
2066            else:   #Chebyschev
2067                T = (2./x)-1.
2068                Ccof = np.ones((12,xdata.shape[0]))
2069                Ccof[1] = T
2070                for i in range(2,12):
2071                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
2072                for i in range(1,10):
2073                    YI += Ccof[i]*Icoef[i+2]
2074                    DYI[i+2] =Ccof[i]
2075        return YI,DYI
2076       
2077    Iesd = np.array(Iparm['Iesd'])
2078    Icovar = Iparm['Icovar']
2079    YI,DYI = IfunAdv(Iparm,xdata)
2080    YI = np.where(YI>0,YI,1.)
2081    WYI = np.zeros_like(xdata)
2082    vcov = np.zeros((12,12))
2083    k = 0
2084    for i in range(12):
2085        for j in range(i,12):
2086            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
2087            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
2088            k += 1
2089    M = np.inner(vcov,DYI.T)
2090    WYI = np.sum(M*DYI,axis=0)
2091    WYI = np.where(WYI>0.,WYI,0.)
2092    return YI,WYI
2093
2094################################################################################
2095#### RMCutilities
2096################################################################################
2097   
2098def MakeInst(G2frame,Name,Phase,useSamBrd,PWId):
2099    PWDdata = G2frame.GetPWDRdatafromTree(PWId)
2100    histoName = G2frame.GPXtree.GetItemPyData(PWId)[2]
2101    Size = Phase['Histograms'][histoName]['Size']
2102    Mustrain = Phase['Histograms'][histoName]['Mustrain']
2103    inst = PWDdata['Instrument Parameters'][0]
2104    Xsb = 0.
2105    Ysb = 0.
2106    if 'T' in inst['Type'][1]:
2107        difC = inst['difC'][1]
2108        if useSamBrd[0]:
2109            if 'ellipsoidal' not in Size[0]:    #take the isotropic term only
2110                Xsb = 1.e-4*difC/Size[1][0]/2.
2111        if useSamBrd[1]:
2112            if 'generalized' not in Mustrain[0]:    #take the isotropic term only
2113                Ysb = 1.e-6*difC*Mustrain[1][0]/2.
2114        prms = ['Bank',
2115                'difC','difA','Zero','2-theta',
2116                'alpha','beta-0','beta-1',
2117                'sig-0','sig-1','sig-2',
2118                'Z','X','Y']
2119        fname = Name+'.inst'
2120        fl = open(fname,'w')
2121        fl.write('1\n')
2122        fl.write('%d\n'%int(inst[prms[0]][1]))
2123        fl.write('%19.11f%19.11f%19.11f%19.11f\n'%(inst[prms[1]][1],inst[prms[2]][1],inst[prms[3]][1],inst[prms[4]][1]))
2124        fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[5]][1],inst[prms[6]][1],inst[prms[7]][1]))
2125        fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[8]][1],inst[prms[9]][1],inst[prms[10]][1]))   
2126        fl.write('%12.6e%14.6e%14.6e%14.6e%14.6e\n'%(inst[prms[11]][1],inst[prms[12]][1]+Xsb,inst[prms[13]][1]+Ysb,0.0,0.0))
2127        fl.close()
2128    else:
2129        if useSamBrd[0]:
2130            wave = G2mth.getWave(inst)
2131            if 'ellipsoidal' not in Size[0]:    #take the isotropic term only
2132                Xsb = 1.8*wave/(np.pi*Size[1][0])
2133        if useSamBrd[1]:
2134            if 'generalized' not in Mustrain[0]:    #take the isotropic term only
2135                Ysb = 0.0180*Mustrain[1][0]/np.pi
2136        prms = ['Bank',
2137                'Lam','Zero','Polariz.',
2138                'U','V','W',
2139                'X','Y']
2140        fname = Name+'.inst'
2141        fl = open(fname,'w')
2142        fl.write('1\n')
2143        fl.write('%d\n'%int(inst[prms[0]][1]))
2144        fl.write('%10.5f%10.5f%10.4f%10d\n'%(inst[prms[1]][1],-100.*inst[prms[2]][1],inst[prms[3]][1],0))
2145        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[4]][1],inst[prms[5]][1],inst[prms[6]][1]))
2146        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[7]][1]+Xsb,inst[prms[8]][1]+Ysb,0.0))   
2147        fl.write('%10.3f%10.3f%10.3f\n'%(0.0,0.0,0.0))
2148        fl.close()
2149    return fname
2150   
2151def MakeBack(G2frame,Name,PWId):
2152    PWDdata = G2frame.GetPWDRdatafromTree(PWId)
2153    Back = PWDdata['Background'][0]
2154    inst = PWDdata['Instrument Parameters'][0]
2155    if 'chebyschev-1' != Back[0]:
2156        return None
2157    Nback = Back[2]
2158    BackVals = Back[3:]
2159    fname = Name+'.back'
2160    fl = open(fname,'w')
2161    fl.write('%10d\n'%Nback)
2162    for val in BackVals:
2163        if 'T' in inst['Type'][1]:
2164            fl.write('%12.6g\n'%(float(val)))
2165        else:
2166            fl.write('%12.6g\n'%val)
2167    fl.close()
2168    return fname
2169
2170def MakeRMC6f(G2frame,Name,Phase,RMCPdict,PWId):
2171   
2172    def findDup(Atoms):
2173        Dup = []
2174        Fracs = []
2175        for iat1,at1 in enumerate(Atoms):
2176            if any([at1[0] in dup for dup in Dup]):
2177                continue
2178            else:
2179                Dup.append([at1[0],])
2180                Fracs.append([at1[6],])
2181            for iat2,at2 in enumerate(Atoms[(iat1+1):]):
2182                if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
2183                    Dup[-1] += [at2[0],]
2184                    Fracs[-1]+= [at2[6],]
2185        return Dup,Fracs
2186   
2187    Meta = RMCPdict['metadata']
2188    Atseq = RMCPdict['atSeq']
2189    Supercell =  RMCPdict['SuperCell']
2190    PWDdata = G2frame.GetPWDRdatafromTree(PWId)
2191    generalData = Phase['General']
2192    Dups,Fracs = findDup(Phase['Atoms'])
2193    Sfracs = [np.cumsum(fracs) for fracs in Fracs]
2194    Sample = PWDdata['Sample Parameters']
2195    Meta['temperature'] = Sample['Temperature']
2196    Meta['pressure'] = Sample['Pressure']
2197    Cell = generalData['Cell'][1:7]
2198    Trans = np.eye(3)*np.array(Supercell)
2199    newPhase = copy.deepcopy(Phase)
2200    newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2201    newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans)
2202    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False)
2203    Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2204    Natm = np.count_nonzero(Natm-1)
2205    Atoms = newPhase['Atoms']
2206    Satoms = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms,5),4),3)
2207    Datoms = [[atom for atom in Satoms if atom[0] in dup] for dup in Dups]
2208    Natoms = []
2209    reset = False
2210    for idup,dup in enumerate(Dups):
2211        ldup = len(dup)
2212        datoms = Datoms[idup]
2213        natm = len(datoms)
2214        i = 0
2215        while i < natm:
2216            atoms = datoms[i:i+ldup]
2217            try:
2218                atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2219                Natoms.append(atom)
2220            except IndexError:      #what about vacancies?
2221                if 'Va' not in Atseq:
2222                    reset = True
2223                    Atseq.append('Va')
2224                    RMCPdict['aTypes']['Va'] = 0.0
2225                atom = atoms[0]
2226                atom[1] = 'Va'
2227                Natoms.append(atom)
2228            i += ldup
2229    NAtype = np.zeros(len(Atseq))
2230    for atom in Natoms:
2231        NAtype[Atseq.index(atom[1])] += 1
2232    NAstr = ['%6d'%i for i in NAtype]
2233    Cell = newPhase['General']['Cell'][1:7]
2234    if os.path.exists(Name+'.his6f'):
2235        os.remove(Name+'.his6f')
2236    if os.path.exists(Name+'.neigh'):
2237        os.remove(Name+'.neigh')
2238    fname = Name+'.rmc6f'
2239    fl = open(fname,'w')
2240    fl.write('(Version 6f format configuration file)\n')
2241    for item in Meta:
2242        fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2243    fl.write('Atom types present:                 %s\n'%'    '.join(Atseq))
2244    fl.write('Number of each atom type:       %s\n'%''.join(NAstr))
2245    fl.write('Number of atoms:                %d\n'%len(Natoms))
2246    fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2247    fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2248            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2249    A,B = G2lat.cell2AB(Cell,True)
2250    fl.write('Lattice vectors (Ang):\n')   
2251    for i in [0,1,2]:
2252        fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2253    fl.write('Atoms (fractional coordinates):\n')
2254    nat = 0
2255    for atm in Atseq:
2256        for iat,atom in enumerate(Natoms):
2257            if atom[1] == atm:
2258                nat += 1
2259                atcode = Atcodes[iat].split(':')
2260                cell = [0,0,0]
2261                if '+' in atcode[1]:
2262                    cell = eval(atcode[1].split('+')[1])
2263                fl.write('%6d%4s  [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(       
2264                        nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2265    fl.close()
2266    return fname,reset
2267
2268def MakeBragg(G2frame,Name,Phase,PWId):
2269    PWDdata = G2frame.GetPWDRdatafromTree(PWId)
2270    generalData = Phase['General']
2271    Vol = generalData['Cell'][7]
2272    Data = PWDdata['Data']
2273    Inst = PWDdata['Instrument Parameters'][0]
2274    Bank = int(Inst['Bank'][1])
2275    Sample = PWDdata['Sample Parameters']
2276    Scale = Sample['Scale'][0]
2277    if 'X' in Inst['Type'][0]:
2278        Scale *= 2.
2279    Limits = PWDdata['Limits'][1]
2280    Ibeg = np.searchsorted(Data[0],Limits[0])
2281    Ifin = np.searchsorted(Data[0],Limits[1])+1
2282    fname = Name+'.bragg'
2283    fl = open(fname,'w')
2284    fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol))
2285    if 'T' in Inst['Type'][0]:
2286        fl.write('%12s%12s\n'%('   TOF,ms','  I(obs)'))
2287        for i in range(Ibeg,Ifin-1):
2288            fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2289    else:
2290        fl.write('%12s%12s\n'%('   2-theta, deg','  I(obs)'))
2291        for i in range(Ibeg,Ifin-1):
2292            fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))       
2293    fl.close()
2294    return fname
2295
2296def MakeRMCPdat(G2frame,Name,Phase,RMCPdict,PWId):
2297    Meta = RMCPdict['metadata']
2298    Times = RMCPdict['runTimes']
2299    Atseq = RMCPdict['atSeq']
2300    Atypes = RMCPdict['aTypes']
2301    atPairs = RMCPdict['Pairs']
2302    Files = RMCPdict['files']
2303    BraggWt = RMCPdict['histogram'][1]
2304    PWDdata = G2frame.GetPWDRdatafromTree(PWId)
2305    inst = PWDdata['Instrument Parameters'][0]
2306    refList = PWDdata['Reflection Lists'][Name]['RefList']
2307    dMin = refList[-1][4]
2308    gsasType = 'xray2'
2309    if 'T' in inst['Type'][1]:
2310        gsasType = 'gsas3'
2311    elif 'X' in inst['Type'][1]:
2312        XFF = G2elem.GetFFtable(Atseq)
2313        Xfl = open(Name+'.xray','w')
2314        for atm in Atseq:
2315            fa = XFF[atm]['fa']
2316            fb = XFF[atm]['fb']
2317            fc = XFF[atm]['fc']
2318            Xfl.write('%2s  %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%(
2319                    atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc))
2320        Xfl.close()
2321    lenA = len(Atseq)
2322    Pairs = []
2323    for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]:
2324        Pairs += pair
2325    pairMin = [atPairs[pair]for pair in Pairs if pair in atPairs]
2326    maxMoves = [Atypes[atm] for atm in Atseq if atm in Atypes]
2327    fname = Name+'.dat'
2328    fl = open(fname,'w')
2329    fl.write(' %% Hand edit the following as needed\n')
2330    fl.write('TITLE :: '+Name+'\n')
2331    fl.write('MATERIAL :: '+Meta['material']+'\n')
2332    fl.write('PHASE :: '+Meta['phase']+'\n')
2333    fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n')
2334    fl.write('INVESTIGATOR :: '+Meta['owner']+'\n')
2335    minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin])
2336    minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin])
2337    maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin])
2338    fl.write('MINIMUM_DISTANCES ::   %s  Angstrom\n'%minHD)
2339    maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves])
2340    fl.write('MAXIMUM_MOVES ::   %s Angstrom\n'%maxMv)
2341    fl.write('R_SPACING ::  0.0200 Angstrom\n')
2342    fl.write('PRINT_PERIOD :: 100\n')
2343    fl.write('TIME_LIMIT ::     %.2f MINUTES\n'%Times[0])
2344    fl.write('SAVE_PERIOD ::    %.2f MINUTES\n'%Times[1])
2345    fl.write('\n')
2346    fl.write('ATOMS :: '+' '.join(Atseq)+'\n')
2347    fl.write('\n')
2348    fl.write('FLAGS ::\n')
2349    fl.write('  > NO_MOVEOUT\n')
2350    fl.write('  > NO_SAVE_CONFIGURATIONS\n')
2351    fl.write('  > NO_RESOLUTION_CONVOLUTION\n')
2352    fl.write('\n')
2353    fl.write('INPUT_CONFIGURATION_FORMAT ::  rmc6f\n')
2354    fl.write('SAVE_CONFIGURATION_FORMAT  ::  rmc6f\n')
2355    fl.write('IGNORE_HISTORY_FILE ::\n')
2356    fl.write('\n')
2357    fl.write('DISTANCE_WINDOW ::\n')
2358    fl.write('  > MNDIST :: %s\n'%minD)
2359    fl.write('  > MXDIST :: %s\n'%maxD)
2360    if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']):
2361        fl.write('\n')
2362        fl.write('POTENTIALS ::\n')
2363        fl.write('  > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.'])
2364        fl.write('  > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n')
2365        if len(RMCPdict['Potentials']['Stretch']):
2366            fl.write('  > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search'])
2367            for bond in RMCPdict['Potentials']['Stretch']:
2368                fl.write('  > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2]))       
2369        if len(RMCPdict['Potentials']['Angles']):
2370            fl.write('  > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
2371            for angle in RMCPdict['Potentials']['Angles']:
2372                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f Ang %.2f Ang\n'%
2373                    (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
2374    if RMCPdict['useBVS']:
2375        fl.write('BVS ::\n')
2376        fl.write('  > ATOM :: '+' '.join(Atseq)+'\n')
2377        fl.write('  > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2378        oxid = []
2379        for val in RMCPdict['Oxid']:
2380            if len(val) == 3:
2381                oxid.append(val[0][1:])
2382            else:
2383                oxid.append(val[0][2:])
2384        fl.write('  > OXID :: %s\n'%' '.join(oxid))
2385        fl.write('  > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']]))
2386        fl.write('  > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']]))
2387        fl.write('  > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][3] for bvs in RMCPdict['BVS']]))       
2388        fl.write('  > SAVE :: 100000\n')
2389        fl.write('  > UPDATE :: 100000\n')
2390        if len(RMCPdict['Swap']):
2391            fl.write('\n')
2392            fl.write('SWAP_MULTI ::\n')
2393            for swap in RMCPdict['Swap']:
2394                try:
2395                    at1 = Atseq.index(swap[0])
2396                    at2 = Atseq.index(swap[1])
2397                except ValueError:
2398                    break
2399                fl.write('  > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2]))
2400       
2401    if len(RMCPdict['FxCN']):
2402        fl.write('FIXED_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['FxCN']))       
2403        for ifx,fxcn in enumerate(RMCPdict['FxCN']):
2404            try:
2405                at1 = Atseq.index(fxcn[0])
2406                at2 = Atseq.index(fxcn[1])
2407            except ValueError:
2408                break
2409            fl.write('  > CSTR%d ::   %d %d %.2f %.2f %.2f %.2f %.6f\n'%(ifx+1,at1+1,at2+1,fxcn[2],fxcn[3],fxcn[4],fxcn[5],fxcn[6]))
2410    if len(RMCPdict['AveCN']):
2411        fl.write('AVERAGE_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['AveCN']))
2412        for iav,avcn in enumerate(RMCPdict['AveCN']):
2413            try:
2414                at1 = Atseq.index(avcn[0])
2415                at2 = Atseq.index(avcn[1])
2416            except ValueError:
2417                break
2418            fl.write('  > CAVSTR%d ::   %d %d %.2f %.2f %.2f %.6f\n'%(iav+1,at1+1,at2+1,avcn[2],avcn[3],avcn[4],avcn[5]))
2419    for File in Files:
2420        if Files[File][0]:
2421            if 'Xray' in File and 'F(Q)' in File:
2422                fqdata = open(Files[File][0],'r')
2423                lines = int(fqdata.readline()[:-1])
2424            fl.write('\n')
2425            fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_'))
2426            fl.write('  > FILENAME :: %s\n'%Files[File][0])
2427            fl.write('  > DATA_TYPE :: %s\n'%Files[File][2])
2428            fl.write('  > FIT_TYPE :: %s\n'%Files[File][2])
2429            if 'Xray' not in File:
2430                fl.write('  > START_POINT :: 1\n')
2431                fl.write('  > END_POINT :: 3000\n')
2432                fl.write('  > WEIGHT :: %.4f\n'%Files[File][1])
2433            fl.write('  > CONSTANT_OFFSET 0.000\n')
2434            fl.write('  > NO_FITTED_OFFSET\n')
2435            if RMCPdict['FitScale']:
2436                fl.write('  > FITTED_SCALE\n')
2437            else:
2438                fl.write('  > NO_FITTED_SCALE\n')
2439            if Files[File][3] !='RMC':
2440                fl.write('  > %s\n'%Files[File][3])
2441            if 'reciprocal' in File:
2442                fl.write('  > CONVOLVE ::\n')
2443                if 'Xray' in File:
2444                    fl.write('  > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines)
2445                    fl.write('  > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1]))
2446                    fl.write('  > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2))
2447                    fl.write('  > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,Files[File][1]))
2448    fl.write('\n')
2449    fl.write('BRAGG ::\n')
2450    fl.write('  > BRAGG_SHAPE :: %s\n'%gsasType)
2451    fl.write('  > RECALCUATE\n')
2452    fl.write('  > DMIN :: %.2f\n'%(dMin-0.02))
2453    fl.write('  > WEIGHT :: %10.3f\n'%BraggWt)
2454    fl.write('\n')
2455    fl.write('END  ::\n')
2456    fl.close()
2457    return fname   
2458
2459def MakePDB(G2frame,Name,Phase,Atseq,Supercell):
2460    generalData = Phase['General']
2461    Cell = generalData['Cell'][1:7]
2462    Trans = np.eye(3)*np.array(Supercell)
2463    newPhase = copy.deepcopy(Phase)
2464    newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2465    newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans.T)
2466    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False)
2467    Atoms = newPhase['Atoms']
2468    Cell = newPhase['General']['Cell'][1:7]
2469    A,B = G2lat. cell2AB(Cell)
2470    fname = Name+'.pdb'
2471    fl = open(fname,'w')
2472    fl.write('REMARK    this file is generated using GSASII\n')
2473    fl.write('CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n'%(
2474            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2475    fl.write('ORIGX1      1.000000  0.000000  0.000000        0.00000\n')
2476    fl.write('ORIGX2      0.000000  1.000000  0.000000        0.00000\n')
2477    fl.write('ORIGX3      0.000000  0.000000  1.000000        0.00000\n')
2478
2479    Natm = np.core.defchararray.count(np.array(Atcodes),'+')
2480    Natm = np.count_nonzero(Natm-1)
2481    nat = 0
2482    for atm in Atseq:
2483        for iat,atom in enumerate(Atoms):
2484            if atom[1] == atm:
2485                nat += 1
2486                XYZ = np.inner(A,np.array(atom[3:6])-0.5)    #shift origin to middle & make Cartesian
2487#ATOM      1 Ni   RMC     1     -22.113 -22.113 -22.113  1.00  0.00          ni                     
2488                fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %-2s\n'%(       
2489                        nat,atom[0],nat,XYZ[0],XYZ[1],XYZ[2],atom[1]))
2490    fl.close()
2491    return fname
2492   
2493################################################################################
2494#### Reflectometry calculations
2495################################################################################
2496
2497def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
2498    G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
2499   
2500    class RandomDisplacementBounds(object):
2501        """random displacement with bounds"""
2502        def __init__(self, xmin, xmax, stepsize=0.5):
2503            self.xmin = xmin
2504            self.xmax = xmax
2505            self.stepsize = stepsize
2506   
2507        def __call__(self, x):
2508            """take a random step but ensure the new position is within the bounds"""
2509            while True:
2510                # this could be done in a much more clever way, but it will work for example purposes
2511                steps = self.xmax-self.xmin
2512                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
2513                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
2514                    break
2515            return xnew
2516   
2517    def GetModelParms():
2518        parmDict = {}
2519        varyList = []
2520        values = []
2521        bounds = []
2522        parmDict['dQ type'] = data['dQ type']
2523        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
2524        for parm in ['Scale','FltBack']:
2525            parmDict[parm] = data[parm][0]
2526            if data[parm][1]:
2527                varyList.append(parm)
2528                values.append(data[parm][0])
2529                bounds.append(Bounds[parm])
2530        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
2531        parmDict['nLayers'] = len(parmDict['Layer Seq'])
2532        for ilay,layer in enumerate(data['Layers']):
2533            name = layer['Name']
2534            cid = str(ilay)+';'
2535            parmDict[cid+'Name'] = name
2536            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2537                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
2538                if layer.get(parm,[0.,False])[1]:
2539                    varyList.append(cid+parm)
2540                    value = layer[parm][0]
2541                    values.append(value)
2542                    if value:
2543                        bound = [value*Bfac,value/Bfac]
2544                    else:
2545                        bound = [0.,10.]
2546                    bounds.append(bound)
2547            if name not in ['vacuum','unit scatter']:
2548                parmDict[cid+'rho'] = Substances[name]['Scatt density']
2549                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
2550        return parmDict,varyList,values,bounds
2551   
2552    def SetModelParms():
2553        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
2554        if 'Scale' in varyList:
2555            data['Scale'][0] = parmDict['Scale']
2556            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
2557        G2fil.G2Print (line)
2558        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
2559        if 'FltBack' in varyList:
2560            data['FltBack'][0] = parmDict['FltBack']
2561            line += ' esd: %15.3g'%(sigDict['FltBack'])
2562        G2fil.G2Print (line)
2563        for ilay,layer in enumerate(data['Layers']):
2564            name = layer['Name']
2565            G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
2566            cid = str(ilay)+';'
2567            line = ' '
2568            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
2569            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
2570            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2571                if parm in layer:
2572                    if parm == 'Rough':
2573                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
2574                    else:
2575                        layer[parm][0] = parmDict[cid+parm]
2576                    line += ' %s: %.3f'%(parm,layer[parm][0])
2577                    if cid+parm in varyList:
2578                        line += ' esd: %.3g'%(sigDict[cid+parm])
2579            G2fil.G2Print (line)
2580            G2fil.G2Print (line2)
2581   
2582    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2583        parmDict.update(zip(varyList,values))
2584        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2585        return M
2586   
2587    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2588        parmDict.update(zip(varyList,values))
2589        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2590        return np.sum(M**2)
2591   
2592    def getREFD(Q,Qsig,parmDict):
2593        Ic = np.ones_like(Q)*parmDict['FltBack']
2594        Scale = parmDict['Scale']
2595        Nlayers = parmDict['nLayers']
2596        Res = parmDict['Res']
2597        depth = np.zeros(Nlayers)
2598        rho = np.zeros(Nlayers)
2599        irho = np.zeros(Nlayers)
2600        sigma = np.zeros(Nlayers)
2601        for ilay,lay in enumerate(parmDict['Layer Seq']):
2602            cid = str(lay)+';'
2603            depth[ilay] = parmDict[cid+'Thick']
2604            sigma[ilay] = parmDict[cid+'Rough']
2605            if parmDict[cid+'Name'] == u'unit scatter':
2606                rho[ilay] = parmDict[cid+'DenMul']
2607                irho[ilay] = parmDict[cid+'iDenMul']
2608            elif 'vacuum' != parmDict[cid+'Name']:
2609                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
2610                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
2611            if cid+'Mag SLD' in parmDict:
2612                rho[ilay] += parmDict[cid+'Mag SLD']
2613        if parmDict['dQ type'] == 'None':
2614            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2615        elif 'const' in parmDict['dQ type']:
2616            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
2617        else:       #dQ/Q in data
2618            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
2619        Ic += AB*Scale
2620        return Ic
2621       
2622    def estimateT0(takestep):
2623        Mmax = -1.e-10
2624        Mmin = 1.e10
2625        for i in range(100):
2626            x0 = takestep(values)
2627            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2628            Mmin = min(M,Mmin)
2629            MMax = max(M,Mmax)
2630        return 1.5*(MMax-Mmin)
2631
2632    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2633    if data.get('2% weight'):
2634        wt = 1./(0.02*Io)**2
2635    Qmin = Limits[1][0]
2636    Qmax = Limits[1][1]
2637    wtFactor = ProfDict['wtFactor']
2638    Bfac = data['Toler']
2639    Ibeg = np.searchsorted(Q,Qmin)
2640    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
2641    Ic[:] = 0
2642    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
2643              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
2644    parmDict,varyList,values,bounds = GetModelParms()
2645    Msg = 'Failed to converge'
2646    if varyList:
2647        if data['Minimizer'] == 'LMLS': 
2648            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
2649                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2650            parmDict.update(zip(varyList,result[0]))
2651            chisq = np.sum(result[2]['fvec']**2)
2652            ncalc = result[2]['nfev']
2653            covM = result[1]
2654            newVals = result[0]
2655        elif data['Minimizer'] == 'Basin Hopping':
2656            xyrng = np.array(bounds).T
2657            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
2658            T0 = estimateT0(take_step)
2659            G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
2660            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
2661                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
2662                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
2663            chisq = result.fun
2664            ncalc = result.nfev
2665            newVals = result.x
2666            covM = []
2667        elif data['Minimizer'] == 'MC/SA Anneal':
2668            xyrng = np.array(bounds).T
2669            result = G2mth.anneal(sumREFD, values, 
2670                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
2671                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
2672                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
2673                ranRange=0.20,autoRan=False,dlg=None)
2674            newVals = result[0]
2675            parmDict.update(zip(varyList,newVals))
2676            chisq = result[1]
2677            ncalc = result[3]
2678            covM = []
2679            G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
2680        elif data['Minimizer'] == 'L-BFGS-B':
2681            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
2682                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2683            parmDict.update(zip(varyList,result['x']))
2684            chisq = result.fun
2685            ncalc = result.nfev
2686            newVals = result.x
2687            covM = []
2688    else:   #nothing varied
2689        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2690        chisq = np.sum(M**2)
2691        ncalc = 0
2692        covM = []
2693        sig = []
2694        sigDict = {}
2695        result = []
2696    Rvals = {}
2697    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
2698    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
2699    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
2700    Ib[Ibeg:Ifin] = parmDict['FltBack']
2701    try:
2702        if not len(varyList):
2703            Msg += ' - nothing refined'
2704            raise ValueError
2705        Nans = np.isnan(newVals)
2706        if np.any(Nans):
2707            name = varyList[Nans.nonzero(True)[0]]
2708            Msg += ' Nan result for '+name+'!'
2709            raise ValueError
2710        Negs = np.less_equal(newVals,0.)
2711        if np.any(Negs):
2712            indx = Negs.nonzero()
2713            name = varyList[indx[0][0]]
2714            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
2715                Msg += ' negative coefficient for '+name+'!'
2716                raise ValueError
2717        if len(covM):
2718            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
2719            covMatrix = covM*Rvals['GOF']
2720        else:
2721            sig = np.zeros(len(varyList))
2722            covMatrix = []
2723        sigDict = dict(zip(varyList,sig))
2724        G2fil.G2Print (' Results of reflectometry data modelling fit:')
2725        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
2726        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2727        SetModelParms()
2728        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
2729    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
2730        G2fil.G2Print (Msg)
2731        return False,0,0,0,0,0,0,Msg
2732       
2733def makeSLDprofile(data,Substances):
2734   
2735    sq2 = np.sqrt(2.)
2736    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2737    Nlayers = len(laySeq)
2738    laySeq = np.array(laySeq,dtype=int)
2739    interfaces = np.zeros(Nlayers)
2740    rho = np.zeros(Nlayers)
2741    sigma = np.zeros(Nlayers)
2742    name = data['Layers'][0]['Name']
2743    thick = 0.
2744    for ilay,lay in enumerate(laySeq):
2745        layer = data['Layers'][lay]
2746        name = layer['Name']
2747        if 'Thick' in layer:
2748            thick += layer['Thick'][0]
2749            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
2750        if 'Rough' in layer:
2751            sigma[ilay] = max(0.001,layer['Rough'][0])
2752        if name != 'vacuum':
2753            if name == 'unit scatter':
2754                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
2755            else:
2756                rrho = Substances[name]['Scatt density']
2757                irho = Substances[name]['XImag density']
2758                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
2759        if 'Mag SLD' in layer:
2760            rho[ilay] += layer['Mag SLD'][0]
2761    name = data['Layers'][-1]['Name']
2762    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
2763    xr = np.flipud(x)
2764    interfaces[-1] = x[-1]
2765    y = np.ones_like(x)*rho[0]
2766    iBeg = 0
2767    for ilayer in range(Nlayers-1):
2768        delt = rho[ilayer+1]-rho[ilayer]
2769        iPos = np.searchsorted(x,interfaces[ilayer])
2770        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
2771        iBeg = iPos
2772    return x,xr,y   
2773
2774def REFDModelFxn(Profile,Inst,Limits,Substances,data):
2775   
2776    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2777    Qmin = Limits[1][0]
2778    Qmax = Limits[1][1]
2779    iBeg = np.searchsorted(Q,Qmin)
2780    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2781    Ib[:] = data['FltBack'][0]
2782    Ic[:] = 0
2783    Scale = data['Scale'][0]
2784    if data['Layer Seq'] == []:
2785        return
2786    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2787    Nlayers = len(laySeq)
2788    depth = np.zeros(Nlayers)
2789    rho = np.zeros(Nlayers)
2790    irho = np.zeros(Nlayers)
2791    sigma = np.zeros(Nlayers)
2792    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
2793        layer = data['Layers'][lay]
2794        name = layer['Name']
2795        if 'Thick' in layer:    #skips first & last layers
2796            depth[ilay] = layer['Thick'][0]
2797        if 'Rough' in layer:    #skips first layer
2798            sigma[ilay] = layer['Rough'][0]
2799        if 'unit scatter' == name:
2800            rho[ilay] = layer['DenMul'][0]
2801            irho[ilay] = layer['iDenMul'][0]
2802        else:
2803            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2804            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
2805        if 'Mag SLD' in layer:
2806            rho[ilay] += layer['Mag SLD'][0]
2807    if data['dQ type'] == 'None':
2808        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2809    elif 'const' in data['dQ type']:
2810        res = data['Resolution'][0]/(100.*sateln2)
2811        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
2812    else:       #dQ/Q in data
2813        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
2814    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
2815
2816def abeles(kz, depth, rho, irho=0, sigma=0):
2817    """
2818    Optical matrix form of the reflectivity calculation.
2819    O.S. Heavens, Optical Properties of Thin Solid Films
2820   
2821    Reflectometry as a function of kz for a set of slabs.
2822
2823    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
2824        This is :math:`\\tfrac12 Q_z`.       
2825    :param depth:  float[m] (Ang).
2826        thickness of each layer.  The thickness of the incident medium
2827        and substrate are ignored.
2828    :param rho:  float[n,k] (1e-6/Ang^2)
2829        Real scattering length density for each layer for each kz
2830    :param irho:  float[n,k] (1e-6/Ang^2)
2831        Imaginary scattering length density for each layer for each kz
2832        Note: absorption cross section mu = 2 irho/lambda for neutrons
2833    :param sigma: float[m-1] (Ang)
2834        interfacial roughness.  This is the roughness between a layer
2835        and the previous layer. The sigma array should have m-1 entries.
2836
2837    Slabs are ordered with the surface SLD at index 0 and substrate at
2838    index -1, or reversed if kz < 0.
2839    """
2840    def calc(kz, depth, rho, irho, sigma):
2841        if len(kz) == 0: return kz
2842   
2843        # Complex index of refraction is relative to the incident medium.
2844        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
2845        # in place of kz^2, and ignoring rho_o
2846        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
2847        k = kz
2848   
2849        # According to Heavens, the initial matrix should be [ 1 F; F 1],
2850        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
2851        # multiply versus some coding convenience.
2852        B11 = 1
2853        B22 = 1
2854        B21 = 0
2855        B12 = 0
2856        for i in range(0, len(depth)-1):
2857            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
2858            F = (k - k_next) / (k + k_next)
2859            F *= np.exp(-2*k*k_next*sigma[i]**2)
2860            #print "==== layer",i
2861            #print "kz:", kz
2862            #print "k:", k
2863            #print "k_next:",k_next
2864            #print "F:",F
2865            #print "rho:",rho[:,i+1]
2866            #print "irho:",irho[:,i+1]
2867            #print "d:",depth[i],"sigma:",sigma[i]
2868            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
2869            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
2870            M21 = F*M11
2871            M12 = F*M22
2872            C1 = B11*M11 + B21*M12
2873            C2 = B11*M21 + B21*M22
2874            B11 = C1
2875            B21 = C2
2876            C1 = B12*M11 + B22*M12
2877            C2 = B12*M21 + B22*M22
2878            B12 = C1
2879            B22 = C2
2880            k = k_next
2881   
2882        r = B12/B11
2883        return np.absolute(r)**2
2884
2885    if np.isscalar(kz): kz = np.asarray([kz], 'd')
2886
2887    m = len(depth)
2888
2889    # Make everything into arrays
2890    depth = np.asarray(depth,'d')
2891    rho = np.asarray(rho,'d')
2892    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
2893    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
2894
2895    # Repeat rho,irho columns as needed
2896    if len(rho.shape) == 1:
2897        rho = rho[None,:]
2898        irho = irho[None,:]
2899
2900    return calc(kz, depth, rho, irho, sigma)
2901   
2902def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
2903    y = abeles(kz, depth, rho, irho, sigma)
2904    s = dq/2.
2905    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
2906    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
2907    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
2908    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
2909    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
2910    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
2911    y *= 0.137023
2912    return y
2913       
2914def makeRefdFFT(Limits,Profile):
2915    G2fil.G2Print ('make fft')
2916    Q,Io = Profile[:2]
2917    Qmin = Limits[1][0]
2918    Qmax = Limits[1][1]
2919    iBeg = np.searchsorted(Q,Qmin)
2920    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2921    Qf = np.linspace(0.,Q[iFin-1],5000)
2922    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
2923    If = QI(Qf)*Qf**4
2924    R = np.linspace(0.,5000.,5000)
2925    F = fft.rfft(If)
2926    return R,F
2927
2928   
2929################################################################################
2930#### Stacking fault simulation codes
2931################################################################################
2932
2933def GetStackParms(Layers):
2934   
2935    Parms = []
2936#cell parms
2937    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
2938        Parms.append('cellA')
2939        Parms.append('cellC')
2940    else:
2941        Parms.append('cellA')
2942        Parms.append('cellB')
2943        Parms.append('cellC')
2944        if Layers['Laue'] != 'mmm':
2945            Parms.append('cellG')
2946#Transition parms
2947    for iY in range(len(Layers['Layers'])):
2948        for iX in range(len(Layers['Layers'])):
2949            Parms.append('TransP;%d;%d'%(iY,iX))
2950            Parms.append('TransX;%d;%d'%(iY,iX))
2951            Parms.append('TransY;%d;%d'%(iY,iX))
2952            Parms.append('TransZ;%d;%d'%(iY,iX))
2953    return Parms
2954
2955def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
2956    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
2957   
2958    :param dict Layers: dict with following items
2959
2960      ::
2961
2962       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
2963       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
2964       'Layers':[],'Stacking':[],'Transitions':[]}
2965       
2966    :param str ctrls: controls string to be written on DIFFaX controls.dif file
2967    :param float scale: scale factor
2968    :param dict background: background parameters
2969    :param list limits: min/max 2-theta to be calculated
2970    :param dict inst: instrument parameters dictionary
2971    :param list profile: powder pattern data
2972   
2973    Note that parameters all updated in place   
2974    '''
2975    import atmdata
2976    path = sys.path
2977    for name in path:
2978        if 'bin' in name:
2979            DIFFaX = name+'/DIFFaX.exe'
2980            G2fil.G2Print (' Execute '+DIFFaX)
2981            break
2982    # make form factor file that DIFFaX wants - atom types are GSASII style
2983    sf = open('data.sfc','w')
2984    sf.write('GSASII special form factor file for DIFFaX\n\n')
2985    atTypes = list(Layers['AtInfo'].keys())
2986    if 'H' not in atTypes:
2987        atTypes.insert(0,'H')
2988    for atType in atTypes:
2989        if atType == 'H': 
2990            blen = -.3741
2991        else:
2992            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
2993        Adat = atmdata.XrayFF[atType]
2994        text = '%4s'%(atType.ljust(4))
2995        for i in range(4):
2996            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
2997        text += '%11.6f%11.6f'%(Adat['fc'],blen)
2998        text += '%3d\n'%(Adat['Z'])
2999        sf.write(text)
3000    sf.close()
3001    #make DIFFaX control.dif file - future use GUI to set some of these flags
3002    cf = open('control.dif','w')
3003    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
3004        x0 = profile[0]
3005        iBeg = np.searchsorted(x0,limits[0])
3006        iFin = np.searchsorted(x0,limits[1])+1
3007        if iFin-iBeg > 20000:
3008            iFin = iBeg+20000
3009        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3010        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3011        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
3012    else:
3013        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3014        inst = {'Type':['XSC','XSC',]}
3015    cf.close()
3016    #make DIFFaX data file
3017    df = open('GSASII-DIFFaX.dat','w')
3018    df.write('INSTRUMENTAL\n')
3019    if 'X' in inst['Type'][0]:
3020        df.write('X-RAY\n')
3021    elif 'N' in inst['Type'][0]:
3022        df.write('NEUTRON\n')
3023    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
3024        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
3025        U = ateln2*inst['U'][1]/10000.
3026        V = ateln2*inst['V'][1]/10000.
3027        W = ateln2*inst['W'][1]/10000.
3028        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
3029        HW = np.sqrt(np.mean(HWHM))
3030    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
3031        if 'Mean' in Layers['selInst']:
3032            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
3033        elif 'Gaussian' in Layers['selInst']:
3034            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
3035        else:
3036            df.write('None\n')
3037    else:
3038        df.write('0.10\nNone\n')
3039    df.write('STRUCTURAL\n')
3040    a,b,c = Layers['Cell'][1:4]
3041    gam = Layers['Cell'][6]
3042    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
3043    laue = Layers['Laue']
3044    if laue == '2/m(ab)':
3045        laue = '2/m(1)'
3046    elif laue == '2/m(c)':
3047        laue = '2/m(2)'
3048    if 'unknown' in Layers['Laue']:
3049        df.write('%s %.3f\n'%(laue,Layers['Toler']))
3050    else:   
3051        df.write('%s\n'%(laue))
3052    df.write('%d\n'%(len(Layers['Layers'])))
3053    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
3054        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
3055    layerNames = []
3056    for layer in Layers['Layers']:
3057        layerNames.append(layer['Name'])
3058    for il,layer in enumerate(Layers['Layers']):
3059        if layer['SameAs']:
3060            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
3061            continue
3062        df.write('LAYER %d\n'%(il+1))
3063        if '-1' in layer['Symm']:
3064            df.write('CENTROSYMMETRIC\n')
3065        else:
3066            df.write('NONE\n')
3067        for ia,atom in enumerate(layer['Atoms']):
3068            [name,atype,x,y,z,frac,Uiso] = atom
3069            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
3070                frac /= 2.
3071            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
3072    df.write('STACKING\n')
3073    df.write('%s\n'%(Layers['Stacking'][0]))
3074    if 'recursive' in Layers['Stacking'][0]:
3075        df.write('%s\n'%Layers['Stacking'][1])
3076    else:
3077        if 'list' in Layers['Stacking'][1]:
3078            Slen = len(Layers['Stacking'][2])
3079            iB = 0
3080            iF = 0
3081            while True:
3082                iF += 68
3083                if iF >= Slen:
3084                    break
3085                iF = min(iF,Slen)
3086                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
3087                iB = iF
3088        else:
3089            df.write('%s\n'%Layers['Stacking'][1])   
3090    df.write('TRANSITIONS\n')
3091    for iY in range(len(Layers['Layers'])):
3092        sumPx = 0.
3093        for iX in range(len(Layers['Layers'])):
3094            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
3095            p = round(p,3)
3096            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
3097            sumPx += p
3098        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
3099            G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
3100            df.close()
3101            os.remove('data.sfc')
3102            os.remove('control.dif')
3103            os.remove('GSASII-DIFFaX.dat')
3104            return       
3105    df.close()   
3106    time0 = time.time()
3107    try:
3108        subp.call(DIFFaX)
3109    except OSError:
3110        G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
3111    G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
3112    if os.path.exists('GSASII-DIFFaX.spc'):
3113        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
3114        iFin = iBeg+Xpat.shape[1]
3115        bakType,backDict,backVary = SetBackgroundParms(background)
3116        backDict['Lam1'] = G2mth.getWave(inst)
3117        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
3118        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
3119        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
3120            rv = st.poisson(profile[3][iBeg:iFin])
3121            profile[1][iBeg:iFin] = rv.rvs()
3122            Z = np.ones_like(profile[3][iBeg:iFin])
3123            Z[1::2] *= -1
3124            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3125            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3126        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3127    #cleanup files..
3128        os.remove('GSASII-DIFFaX.spc')
3129    elif os.path.exists('GSASII-DIFFaX.sadp'):
3130        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
3131        Sadp = np.reshape(Sadp,(256,-1))
3132        Layers['Sadp']['Img'] = Sadp
3133        os.remove('GSASII-DIFFaX.sadp')
3134    os.remove('data.sfc')
3135    os.remove('control.dif')
3136    os.remove('GSASII-DIFFaX.dat')
3137   
3138def SetPWDRscan(inst,limits,profile):
3139   
3140    wave = G2mth.getMeanWave(inst)
3141    x0 = profile[0]
3142    iBeg = np.searchsorted(x0,limits[0])
3143    iFin = np.searchsorted(x0,limits[1])
3144    if iFin-iBeg > 20000:
3145        iFin = iBeg+20000
3146    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3147    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
3148    return iFin-iBeg
3149       
3150def SetStackingSF(Layers,debug):
3151# Load scattering factors into DIFFaX arrays
3152    import atmdata
3153    atTypes = Layers['AtInfo'].keys()
3154    aTypes = []
3155    for atype in atTypes:
3156        aTypes.append('%4s'%(atype.ljust(4)))
3157    SFdat = []
3158    for atType in atTypes:
3159        Adat = atmdata.XrayFF[atType]
3160        SF = np.zeros(9)
3161        SF[:8:2] = Adat['fa']
3162        SF[1:8:2] = Adat['fb']
3163        SF[8] = Adat['fc']
3164        SFdat.append(SF)
3165    SFdat = np.array(SFdat)
3166    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
3167   
3168def SetStackingClay(Layers,Type):
3169# Controls
3170    rand.seed()
3171    ranSeed = rand.randint(1,2**16-1)
3172    try:   
3173        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
3174            '6/m','6/mmm'].index(Layers['Laue'])+1
3175    except ValueError:  #for 'unknown'
3176        laueId = -1
3177    if 'SADP' in Type:
3178        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
3179        lmax = int(Layers['Sadp']['Lmax'])
3180    else:
3181        planeId = 0
3182        lmax = 0
3183# Sequences
3184    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
3185    try:
3186        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
3187    except ValueError:
3188        StkParm = -1
3189    if StkParm == 2:    #list
3190        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
3191        Nstk = len(StkSeq)
3192    else:
3193        Nstk = 1
3194        StkSeq = [0,]
3195    if StkParm == -1:
3196        StkParm = int(Layers['Stacking'][1])
3197    Wdth = Layers['Width'][0]
3198    mult = 1
3199    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
3200    LaueSym = Layers['Laue'].ljust(12)
3201    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
3202    return laueId,controls
3203   
3204def SetCellAtoms(Layers):
3205    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
3206# atoms in layers
3207    atTypes = list(Layers['AtInfo'].keys())
3208    AtomXOU = []
3209    AtomTp = []
3210    LayerSymm = []
3211    LayerNum = []
3212    layerNames = []
3213    Natm = 0
3214    Nuniq = 0
3215    for layer in Layers['Layers']:
3216        layerNames.append(layer['Name'])
3217    for il,layer in enumerate(Layers['Layers']):
3218        if layer['SameAs']:
3219            LayerNum.append(layerNames.index(layer['SameAs'])+1)
3220            continue
3221        else:
3222            LayerNum.append(il+1)
3223            Nuniq += 1
3224        if '-1' in layer['Symm']:
3225            LayerSymm.append(1)
3226        else:
3227            LayerSymm.append(0)
3228        for ia,atom in enumerate(layer['Atoms']):
3229            [name,atype,x,y,z,frac,Uiso] = atom
3230            Natm += 1
3231            AtomTp.append('%4s'%(atype.ljust(4)))
3232            Ta = atTypes.index(atype)+1
3233            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
3234    AtomXOU = np.array(AtomXOU)
3235    Nlayers = len(layerNames)
3236    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
3237    return Nlayers
3238   
3239def SetStackingTrans(Layers,Nlayers):
3240# Transitions
3241    TransX = []
3242    TransP = []
3243    for Ytrans in Layers['Transitions']:
3244        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
3245        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
3246    TransP = np.array(TransP,dtype='float').T
3247    TransX = np.array(TransX,dtype='float')
3248#    GSASIIpath.IPyBreak()
3249    pyx.pygettrans(Nlayers,TransP,TransX)
3250   
3251def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
3252# Scattering factors
3253    SetStackingSF(Layers,debug)
3254# Controls & sequences
3255    laueId,controls = SetStackingClay(Layers,'PWDR')
3256# cell & atoms
3257    Nlayers = SetCellAtoms(Layers)
3258    Volume = Layers['Cell'][7]   
3259# Transitions
3260    SetStackingTrans(Layers,Nlayers)
3261# PWDR scan
3262    Nsteps = SetPWDRscan(inst,limits,profile)
3263# result as Spec
3264    x0 = profile[0]
3265    profile[3] = np.zeros(len(profile[0]))
3266    profile[4] = np.zeros(len(profile[0]))
3267    profile[5] = np.zeros(len(profile[0]))
3268    iBeg = np.searchsorted(x0,limits[0])
3269    iFin = np.searchsorted(x0,limits[1])+1
3270    if iFin-iBeg > 20000:
3271        iFin = iBeg+20000
3272    Nspec = 20001       
3273    spec = np.zeros(Nspec,dtype='double')   
3274    time0 = time.time()
3275    pyx.pygetspc(controls,Nspec,spec)
3276    G2fil.G2Print (' GETSPC time = %.2fs'%(time.time()-time0))
3277    time0 = time.time()
3278    U = ateln2*inst['U'][1]/10000.
3279    V = ateln2*inst['V'][1]/10000.
3280    W = ateln2*inst['W'][1]/10000.
3281    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
3282    HW = np.sqrt(np.mean(HWHM))
3283    BrdSpec = np.zeros(Nsteps)
3284    if 'Mean' in Layers['selInst']:
3285        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
3286    elif 'Gaussian' in Layers['selInst']:
3287        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
3288    else:
3289        BrdSpec = spec[:Nsteps]
3290    BrdSpec /= Volume
3291    iFin = iBeg+Nsteps
3292    bakType,backDict,backVary = SetBackgroundParms(background)
3293    backDict['Lam1'] = G2mth.getWave(inst)
3294    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
3295    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
3296    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
3297        try:
3298            rv = st.poisson(profile[3][iBeg:iFin])
3299            profile[1][iBeg:iFin] = rv.rvs()
3300        except ValueError:
3301            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
3302        Z = np.ones_like(profile[3][iBeg:iFin])
3303        Z[1::2] *= -1
3304        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3305        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3306    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3307    G2fil.G2Print (' Broadening time = %.2fs'%(time.time()-time0))
3308   
3309def CalcStackingSADP(Layers,debug):
3310   
3311# Scattering factors
3312    SetStackingSF(Layers,debug)
3313# Controls & sequences
3314    laueId,controls = SetStackingClay(Layers,'SADP')
3315# cell & atoms
3316    Nlayers = SetCellAtoms(Layers)   
3317# Transitions
3318    SetStackingTrans(Layers,Nlayers)
3319# result as Sadp
3320    Nspec = 20001       
3321    spec = np.zeros(Nspec,dtype='double')   
3322    time0 = time.time()
3323    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
3324    Sapd = np.zeros((256,256))
3325    iB = 0
3326    for i in range(hkLim):
3327        iF = iB+Nblk
3328        p1 = 127+int(i*Incr)
3329        p2 = 128-int(i*Incr)
3330        if Nblk == 128:
3331            if i:
3332                Sapd[128:,p1] = spec[iB:iF]
3333                Sapd[:128,p1] = spec[iF:iB:-1]
3334            Sapd[128:,p2] = spec[iB:iF]
3335            Sapd[:128,p2] = spec[iF:iB:-1]
3336        else:
3337            if i:
3338                Sapd[:,p1] = spec[iB:iF]
3339            Sapd[:,p2] = spec[iB:iF]
3340        iB += Nblk
3341    Layers['Sadp']['Img'] = Sapd
3342    G2fil.G2Print (' GETSAD time = %.2fs'%(time.time()-time0))
3343   
3344###############################################################################
3345#### Maximum Entropy Method - Dysnomia
3346###############################################################################
3347   
3348def makePRFfile(data,MEMtype):
3349    ''' makes Dysnomia .prf control file from Dysnomia GUI controls
3350   
3351    :param dict data: GSAS-II phase data
3352    :param int MEMtype: 1 for neutron data with negative scattering lengths
3353                        0 otherwise
3354    :returns str: name of Dysnomia control file
3355    '''
3356
3357    generalData = data['General']
3358    pName = generalData['Name'].replace(' ','_')
3359    DysData = data['Dysnomia']
3360    prfName = pName+'.prf'
3361    prf = open(prfName,'w')
3362    prf.write('$PREFERENCES\n')
3363    prf.write(pName+'.mem\n') #or .fos?
3364    prf.write(pName+'.out\n')
3365    prf.write(pName+'.pgrid\n')
3366    prf.write(pName+'.fba\n')
3367    prf.write(pName+'_eps.raw\n')
3368    prf.write('%d\n'%MEMtype)
3369    if DysData['DenStart'] == 'uniform':
3370        prf.write('0\n')
3371    else:
3372        prf.write('1\n')
3373    if DysData['Optimize'] == 'ZSPA':
3374        prf.write('0\n')
3375    else:
3376        prf.write('1\n')
3377    prf.write('1\n')
3378    if DysData['Lagrange'][0] == 'user':
3379        prf.write('0\n')
3380    else:
3381        prf.write('1\n')
3382    prf.write('%.4f %d\n'%(DysData['Lagrange'][1],DysData['wt pwr']))
3383    prf.write('%.3f\n'%DysData['Lagrange'][2])
3384    prf.write('%.2f\n'%DysData['E_factor'])
3385    prf.write('1\n')
3386    prf.write('0\n')
3387    prf.write('%d\n'%DysData['Ncyc'])
3388    prf.write('1\n')
3389    prf.write('1 0 0 0 0 0 0 0\n')
3390    if DysData['prior'] == 'uniform':
3391        prf.write('0\n')
3392    else:
3393        prf.write('1\n')
3394    prf.close()
3395    return prfName
3396
3397def makeMEMfile(data,reflData,MEMtype,DYSNOMIA):
3398    ''' make Dysnomia .mem file of reflection data, etc.
3399
3400    :param dict data: GSAS-II phase data
3401    :param list reflData: GSAS-II reflection data
3402    :param int MEMtype: 1 for neutron data with negative scattering lengths
3403                        0 otherwise
3404    :param str DYSNOMIA: path to dysnomia.exe
3405    '''
3406   
3407    DysData = data['Dysnomia']
3408    generalData = data['General']
3409    cell = generalData['Cell'][1:7]
3410    A = G2lat.cell2A(cell)
3411    SGData = generalData['SGData']
3412    pName = generalData['Name'].replace(' ','_')
3413    memName = pName+'.mem'
3414    Map = generalData['Map']
3415    Type = Map['Type']
3416    UseList = Map['RefList']
3417    mem = open(memName,'w')
3418    mem.write('%s\n'%(generalData['Name']+' from '+UseList[0]))
3419    a,b,c,alp,bet,gam = cell
3420    mem.write('%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n'%(a,b,c,alp,bet,gam))
3421    mem.write('      0.0000000      0.0000000     -1    0    0    0     P\n')   #dummy PO stuff
3422    SGSym = generalData['SGData']['SpGrp']
3423    try:
3424        SGId = G2spc.spgbyNum.index(SGSym)
3425    except ValueError:
3426        return False
3427    org = 1
3428    if SGSym in G2spc.spg2origins:
3429        org = 2
3430    mapsize = Map['rho'].shape
3431    sumZ = 0.
3432    sumpos = 0.
3433    sumneg = 0.
3434    mem.write('%5d%5d%5d%5d%5d\n'%(SGId,org,mapsize[0],mapsize[1],mapsize[2]))
3435    for atm in generalData['NoAtoms']:
3436        Nat = generalData['NoAtoms'][atm]
3437        AtInfo = G2elem.GetAtomInfo(atm)
3438        sumZ += Nat*AtInfo['Z']
3439        isotope = generalData['Isotope'][atm]
3440        blen = generalData['Isotopes'][atm][isotope]['SL'][0]
3441        if blen < 0.:
3442            sumneg += blen*Nat
3443        else:
3444            sumpos += blen*Nat
3445    if 'X' in Type:
3446        mem.write('%10.2f  0.001\n'%sumZ)
3447    elif 'N' in Type and MEMtype:
3448        mem.write('%10.3f%10.3f 0.001\n'%(sumpos,sumneg))
3449    else:
3450        mem.write('%10.3f 0.001\n'%sumpos)
3451       
3452    dmin = DysData['MEMdmin']
3453    TOFlam = 2.0*dmin*npsind(80.0)
3454    refSet = G2lat.GenHLaue(dmin,SGData,A)      #list of h,k,l,d
3455    refDict = {'%d %d %d'%(ref[0],ref[1],ref[2]):ref for ref in refSet}
3456       
3457    refs = []
3458    prevpos = 0.
3459    for ref in reflData:
3460        if ref[3] < 0:
3461            continue
3462        if 'T' in Type:
3463            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,x,x,x,prfo = ref[:16]
3464            s = np.sqrt(max(sig,0.0001))   #var -> sig in deg
3465            FWHM = getgamFW(gam,s)
3466            if dsp < dmin:
3467                continue
3468            theta = npasind(TOFlam/(2.*dsp))
3469            FWHM *= nptand(theta)/pos
3470            pos = 2.*theta
3471        else:
3472            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,prfo = ref[:13]
3473            g = gam/100.    #centideg -> deg
3474            s = np.sqrt(max(sig,0.0001))/100.   #var -> sig in deg
3475            FWHM = getgamFW(g,s)
3476        delt = pos-prevpos
3477        refs.append([h,k,l,mult,pos,FWHM,Fobs,phase,delt])
3478        prevpos = pos
3479           
3480    ovlp = DysData['overlap']
3481    refs1 = []
3482    refs2 = []
3483    nref2 = 0
3484    iref = 0
3485    Nref = len(refs)
3486    start = False
3487    while iref < Nref-1:
3488        if refs[iref+1][-1] < ovlp*refs[iref][5]:
3489            if refs[iref][-1] > ovlp*refs[iref][5]:
3490                refs2.append([])
3491                start = True
3492            if nref2 == len(refs2):
3493                refs2.append([])
3494            refs2[nref2].append(refs[iref])
3495        else:
3496            if start:
3497                refs2[nref2].append(refs[iref])
3498                start = False
3499                nref2 += 1
3500            else:
3501                refs1.append(refs[iref])
3502        iref += 1
3503    if start:
3504        refs2[nref2].append(refs[iref])
3505    else:
3506        refs1.append(refs[iref])
3507   
3508    mem.write('%5d\n'%len(refs1))
3509    for ref in refs1:
3510        h,k,l = ref[:3]
3511        hkl = '%d %d %d'%(h,k,l)
3512        if hkl in refDict:
3513            del refDict[hkl]
3514        Fobs = np.sqrt(ref[6])
3515        mem.write('%5d%5d%5d%10.3f%10.3f%10.3f\n'%(h,k,l,Fobs*npcosd(ref[7]),Fobs*npsind(ref[7]),max(0.01*Fobs,0.1)))
3516    while True and nref2:
3517        if not len(refs2[-1]):
3518            del refs2[-1]
3519        else:
3520            break
3521    mem.write('%5d\n'%len(refs2))
3522    for iref2,ref2 in enumerate(refs2):
3523        mem.write('#%5d\n'%iref2)
3524        mem.write('%5d\n'%len(ref2))
3525        Gsum = 0.
3526        Msum = 0
3527        for ref in ref2:
3528            Gsum += ref[6]*ref[3]
3529            Msum += ref[3]
3530        G = np.sqrt(Gsum/Msum)
3531        h,k,l = ref2[0][:3]
3532        hkl = '%d %d %d'%(h,k,l)
3533        if hkl in refDict:
3534            del refDict[hkl]
3535        mem.write('%5d%5d%5d%10.3f%10.3f%5d\n'%(h,k,l,G,max(0.01*G,0.1),ref2[0][3]))
3536        for ref in ref2[1:]:
3537            h,k,l,m = ref[:4]
3538            mem.write('%5d%5d%5d%5d\n'%(h,k,l,m))
3539            hkl = '%d %d %d'%(h,k,l)
3540            if hkl in refDict:
3541                del refDict[hkl]
3542    if len(refDict):
3543        mem.write('%d\n'%len(refDict))
3544        for hkl in list(refDict.keys()):
3545            h,k,l = refDict[hkl][:3]
3546            mem.write('%5d%5d%5d\n'%(h,k,l))
3547    else:
3548        mem.write('0\n')
3549    mem.close()
3550    return True
3551
3552def MEMupdateReflData(prfName,data,reflData):
3553    ''' Update reflection data with new Fosq, phase result from Dysnomia
3554
3555    :param str prfName: phase.mem file name
3556    :param list reflData: GSAS-II reflection data
3557    '''
3558   
3559    generalData = data['General']
3560    Map = generalData['Map']
3561    Type = Map['Type']
3562    cell = generalData['Cell'][1:7]
3563    A = G2lat.cell2A(cell)
3564    reflDict = {}
3565    newRefs = []
3566    for iref,ref in enumerate(reflData):
3567        if ref[3] > 0:
3568            newRefs.append(ref)
3569            reflDict[hash('%5d%5d%5d'%(ref[0],ref[1],ref[2]))] = iref
3570    fbaName = os.path.splitext(prfName)[0]+'.fba'
3571    try: # patch for FileNotFoundError not in Python 2.7
3572        FileNotFoundError
3573    except NameError:
3574        FileNotFoundError = Exception
3575    try:
3576        fba = open(fbaName,'r')
3577    except FileNotFoundError:
3578        return False
3579    fba.readline()
3580    Nref = int(fba.readline()[:-1])
3581    fbalines = fba.readlines()
3582    for line in fbalines[:Nref]:
3583        info = line.split()
3584        h = int(info[0])
3585        k = int(info[1])
3586        l = int(info[2])
3587        FoR = float(info[3])
3588        FoI = float(info[4])
3589        Fosq = FoR**2+FoI**2
3590        phase = npatan2d(FoI,FoR)
3591        try:
3592            refId = reflDict[hash('%5d%5d%5d'%(h,k,l))]
3593        except KeyError:    #added reflections at end skipped
3594            d = float(1/np.sqrt(G2lat.calc_rDsq([h,k,l],A)))
3595            if 'T' in Type:
3596                newRefs.append([h,k,l,-1,d,0.,0.01,1.0,Fosq,Fosq,phase,1.0,1.0,1.0,1.0,1.0,1.0,1.0])
3597            else:
3598                newRefs.append([h,k,l,-1,d,0.,0.01,1.0,Fosq,Fosq,phase,1.0,1.0,1.0,1.0])
3599            continue
3600        newRefs[refId][8] = Fosq
3601        newRefs[refId][10] = phase
3602    newRefs = np.array(newRefs)
3603    return True,newRefs
3604   
3605#### testing data
3606NeedTestData = True
3607def TestData():
3608    'needs a doc string'
3609#    global NeedTestData
3610    global bakType
3611    bakType = 'chebyschev'
3612    global xdata
3613    xdata = np.linspace(4.0,40.0,36000)
3614    global parmDict0
3615    parmDict0 = {
3616        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
3617        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
3618        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
3619        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
3620        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'Z':0.0,'SH/L':0.002,
3621        'Back0':5.384,'Back1':-0.015,'Back2':.004,
3622        }
3623    global parmDict1
3624    parmDict1 = {
3625        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
3626        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
3627        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
3628        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
3629        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
3630        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
3631        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'Z':0.0,'SH/L':0.002,
3632        'Back0':36.897,'Back1':-0.508,'Back2':.006,
3633        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
3634        }
3635    global parmDict2
3636    parmDict2 = {
3637        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
3638        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'Z':0.0,'SH/L':0.02,
3639        'Back0':5.,'Back1':-0.02,'Back2':.004,
3640#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
3641        }
3642    global varyList
3643    varyList = []
3644
3645def test0():
3646    if NeedTestData: TestData()
3647    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
3648    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
3649    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
3650    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
3651    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
3652    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
3653   
3654def test1():
3655    if NeedTestData: TestData()
3656    time0 = time.time()
3657    for i in range(100):
3658        getPeakProfile(parmDict1,xdata,varyList,bakType)
3659    G2fil.G2Print ('100+6*Ka1-2 peaks=1200 peaks %.2f'%time.time()-time0)
3660   
3661def test2(name,delt):
3662    if NeedTestData: TestData()
3663    varyList = [name,]
3664    xdata = np.linspace(5.6,5.8,400)
3665    hplot = plotter.add('derivatives test for '+name).gca()
3666    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
3667    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
3668    parmDict2[name] += delt
3669    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
3670    hplot.plot(xdata,(y1-y0)/delt,'r+')
3671   
3672def test3(name,delt):
3673    if NeedTestData: TestData()
3674    names = ['pos','sig','gam','shl']
3675    idx = names.index(name)
3676    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
3677    xdata = np.linspace(5.6,5.8,800)
3678    dx = xdata[1]-xdata[0]
3679    hplot = plotter.add('derivatives test for '+name).gca()
3680    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
3681    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
3682    myDict[name] += delt
3683    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
3684    hplot.plot(xdata,(y1-y0)/delt,'r+')
3685
3686if __name__ == '__main__':
3687    import GSASIItestplot as plot
3688    global plotter
3689    plotter = plot.PlotNotebook()
3690#    test0()
3691#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)']:
3692    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
3693        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['Z',0.0001],['SH/L',0.00005]]:
3694        test2(name,shft)
3695    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
3696        test3(name,shft)
3697    G2fil.G2Print ("OK")
3698    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.