source: trunk/GSASIIpwd.py @ 4342

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

correct scale factor for xray .bragg file
try 1/w for xray real space & w for xray reciprocal space weights - may hve to revisit

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