source: trunk/GSASIIpwd.py @ 4389

Last change on this file since 4389 was 4389, checked in by vondreele, 19 months ago

adding fullrmc View & Execute methods. Still to do Setup.
Replace SetText? with SetItemLabel? for the Refine menu item . SetText? is obsolete
remove Load/Save? menu items in RMC tab - won't be used
put in citation info for RMCProfile & fullrmc
remove G2frame from RMCProfile Make routines in G2pwdr. It's a wx object & shouldn't be in there.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 149.3 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-04-02 18:44:42 +0000 (Thu, 02 Apr 2020) $
10# $Author: vondreele $
11# $Revision: 4389 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 4389 2020-04-02 18:44:42Z 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: 4389 $")
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(PWDdata,Name,Size,Mustrain,useSamBrd):
2110    inst = PWDdata['Instrument Parameters'][0]
2111    Xsb = 0.
2112    Ysb = 0.
2113    if 'T' in inst['Type'][1]:
2114        difC = inst['difC'][1]
2115        if useSamBrd[0]:
2116            if 'ellipsoidal' not in Size[0]:    #take the isotropic term only
2117                Xsb = 1.e-4*difC/Size[1][0]
2118        if useSamBrd[1]:
2119            if 'generalized' not in Mustrain[0]:    #take the isotropic term only
2120                Ysb = 1.e-6*difC*Mustrain[1][0]
2121        prms = ['Bank',
2122                'difC','difA','Zero','2-theta',
2123                'alpha','beta-0','beta-1',
2124                'sig-0','sig-1','sig-2',
2125                'Z','X','Y']
2126        fname = Name+'.inst'
2127        fl = open(fname,'w')
2128        fl.write('1\n')
2129        fl.write('%d\n'%int(inst[prms[0]][1]))
2130        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]))
2131        fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[5]][1],inst[prms[6]][1],inst[prms[7]][1]))
2132        fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[8]][1],inst[prms[9]][1],inst[prms[10]][1]))   
2133        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))
2134        fl.close()
2135    else:
2136        if useSamBrd[0]:
2137            wave = G2mth.getWave(inst)
2138            if 'ellipsoidal' not in Size[0]:    #take the isotropic term only
2139                Xsb = 1.8*wave/(np.pi*Size[1][0])
2140        if useSamBrd[1]:
2141            if 'generalized' not in Mustrain[0]:    #take the isotropic term only
2142                Ysb = 0.0180*Mustrain[1][0]/np.pi
2143        prms = ['Bank',
2144                'Lam','Zero','Polariz.',
2145                'U','V','W',
2146                'X','Y']
2147        fname = Name+'.inst'
2148        fl = open(fname,'w')
2149        fl.write('1\n')
2150        fl.write('%d\n'%int(inst[prms[0]][1]))
2151        fl.write('%10.5f%10.5f%10.4f%10d\n'%(inst[prms[1]][1],-100.*inst[prms[2]][1],inst[prms[3]][1],0))
2152        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[4]][1],inst[prms[5]][1],inst[prms[6]][1]))
2153        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[7]][1]+Xsb,inst[prms[8]][1]+Ysb,0.0))   
2154        fl.write('%10.3f%10.3f%10.3f\n'%(0.0,0.0,0.0))
2155        fl.close()
2156    return fname
2157   
2158def MakeBack(PWDdata,Name):
2159    Back = PWDdata['Background'][0]
2160    inst = PWDdata['Instrument Parameters'][0]
2161    if 'chebyschev-1' != Back[0]:
2162        return None
2163    Nback = Back[2]
2164    BackVals = Back[3:]
2165    fname = Name+'.back'
2166    fl = open(fname,'w')
2167    fl.write('%10d\n'%Nback)
2168    for val in BackVals:
2169        if 'T' in inst['Type'][1]:
2170            fl.write('%12.6g\n'%(float(val)))
2171        else:
2172            fl.write('%12.6g\n'%val)
2173    fl.close()
2174    return fname
2175
2176def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):
2177   
2178    def findDup(Atoms):
2179        Dup = []
2180        Fracs = []
2181        for iat1,at1 in enumerate(Atoms):
2182            if any([at1[0] in dup for dup in Dup]):
2183                continue
2184            else:
2185                Dup.append([at1[0],])
2186                Fracs.append([at1[6],])
2187            for iat2,at2 in enumerate(Atoms[(iat1+1):]):
2188                if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
2189                    Dup[-1] += [at2[0],]
2190                    Fracs[-1]+= [at2[6],]
2191        return Dup,Fracs
2192   
2193    Meta = RMCPdict['metadata']
2194    Atseq = RMCPdict['atSeq']
2195    Supercell =  RMCPdict['SuperCell']
2196    generalData = Phase['General']
2197    Dups,Fracs = findDup(Phase['Atoms'])
2198    Sfracs = [np.cumsum(fracs) for fracs in Fracs]
2199    Sample = PWDdata['Sample Parameters']
2200    Meta['temperature'] = Sample['Temperature']
2201    Meta['pressure'] = Sample['Pressure']
2202    Cell = generalData['Cell'][1:7]
2203    Trans = np.eye(3)*np.array(Supercell)
2204    newPhase = copy.deepcopy(Phase)
2205    newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2206    newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans)
2207    GB = G2lat.cell2Gmat( newPhase['General']['Cell'][1:7])[0]
2208    RMCPdict['Rmax'] = np.min(np.sqrt(np.array([1./G2lat.calc_rDsq2(H,GB) for H in [[1,0,0],[0,1,0],[0,0,1]]])))/2.
2209    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False)
2210    Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2211    Natm = np.count_nonzero(Natm-1)
2212    Atoms = newPhase['Atoms']
2213    Satoms = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms,5),4),3)
2214    Datoms = [[atom for atom in Satoms if atom[0] in dup] for dup in Dups]
2215    Natoms = []
2216    reset = False
2217    for idup,dup in enumerate(Dups):
2218        ldup = len(dup)
2219        datoms = Datoms[idup]
2220        natm = len(datoms)
2221        i = 0
2222        while i < natm:
2223            atoms = datoms[i:i+ldup]
2224            try:
2225                atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2226                Natoms.append(atom)
2227            except IndexError:      #what about vacancies?
2228                if 'Va' not in Atseq:
2229                    reset = True
2230                    Atseq.append('Va')
2231                    RMCPdict['aTypes']['Va'] = 0.0
2232                atom = atoms[0]
2233                atom[1] = 'Va'
2234                Natoms.append(atom)
2235            i += ldup
2236    NAtype = np.zeros(len(Atseq))
2237    for atom in Natoms:
2238        NAtype[Atseq.index(atom[1])] += 1
2239    NAstr = ['%6d'%i for i in NAtype]
2240    Cell = newPhase['General']['Cell'][1:7]
2241    if os.path.exists(Name+'.his6f'):
2242        os.remove(Name+'.his6f')
2243    if os.path.exists(Name+'.neigh'):
2244        os.remove(Name+'.neigh')
2245    fname = Name+'.rmc6f'
2246    fl = open(fname,'w')
2247    fl.write('(Version 6f format configuration file)\n')
2248    for item in Meta:
2249        fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2250    fl.write('Atom types present:                 %s\n'%'    '.join(Atseq))
2251    fl.write('Number of each atom type:       %s\n'%''.join(NAstr))
2252    fl.write('Number of atoms:                %d\n'%len(Natoms))
2253    fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2254    fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2255            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2256    A,B = G2lat.cell2AB(Cell,True)
2257    fl.write('Lattice vectors (Ang):\n')   
2258    for i in [0,1,2]:
2259        fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2260    fl.write('Atoms (fractional coordinates):\n')
2261    nat = 0
2262    for atm in Atseq:
2263        for iat,atom in enumerate(Natoms):
2264            if atom[1] == atm:
2265                nat += 1
2266                atcode = Atcodes[iat].split(':')
2267                cell = [0,0,0]
2268                if '+' in atcode[1]:
2269                    cell = eval(atcode[1].split('+')[1])
2270                fl.write('%6d%4s  [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(       
2271                        nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2272    fl.close()
2273    return fname,reset
2274
2275def MakeBragg(PWDdata,Name,Phase):
2276    generalData = Phase['General']
2277    Vol = generalData['Cell'][7]
2278    Data = PWDdata['Data']
2279    Inst = PWDdata['Instrument Parameters'][0]
2280    Bank = int(Inst['Bank'][1])
2281    Sample = PWDdata['Sample Parameters']
2282    Scale = Sample['Scale'][0]
2283    if 'X' in Inst['Type'][1]:
2284        Scale *= 2.
2285    Limits = PWDdata['Limits'][1]
2286    Ibeg = np.searchsorted(Data[0],Limits[0])
2287    Ifin = np.searchsorted(Data[0],Limits[1])+1
2288    fname = Name+'.bragg'
2289    fl = open(fname,'w')
2290    fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol))
2291    if 'T' in Inst['Type'][0]:
2292        fl.write('%12s%12s\n'%('   TOF,ms','  I(obs)'))
2293        for i in range(Ibeg,Ifin-1):
2294            fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2295    else:
2296        fl.write('%12s%12s\n'%('   2-theta, deg','  I(obs)'))
2297        for i in range(Ibeg,Ifin-1):
2298            fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))       
2299    fl.close()
2300    return fname
2301
2302def MakeRMCPdat(PWDdata,Name,Phase,RMCPdict):
2303    Meta = RMCPdict['metadata']
2304    Times = RMCPdict['runTimes']
2305    Atseq = RMCPdict['atSeq']
2306    Atypes = RMCPdict['aTypes']
2307    atPairs = RMCPdict['Pairs']
2308    Files = RMCPdict['files']
2309    BraggWt = RMCPdict['histogram'][1]
2310    inst = PWDdata['Instrument Parameters'][0]
2311    refList = PWDdata['Reflection Lists'][Name]['RefList']
2312    dMin = refList[-1][4]
2313    gsasType = 'xray2'
2314    if 'T' in inst['Type'][1]:
2315        gsasType = 'gsas3'
2316    elif 'X' in inst['Type'][1]:
2317        XFF = G2elem.GetFFtable(Atseq)
2318        Xfl = open(Name+'.xray','w')
2319        for atm in Atseq:
2320            fa = XFF[atm]['fa']
2321            fb = XFF[atm]['fb']
2322            fc = XFF[atm]['fc']
2323            Xfl.write('%2s  %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%(
2324                    atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc))
2325        Xfl.close()
2326    lenA = len(Atseq)
2327    Pairs = []
2328    for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]:
2329        Pairs += pair
2330    pairMin = [atPairs[pair]for pair in Pairs if pair in atPairs]
2331    maxMoves = [Atypes[atm] for atm in Atseq if atm in Atypes]
2332    fname = Name+'.dat'
2333    fl = open(fname,'w')
2334    fl.write(' %% Hand edit the following as needed\n')
2335    fl.write('TITLE :: '+Name+'\n')
2336    fl.write('MATERIAL :: '+Meta['material']+'\n')
2337    fl.write('PHASE :: '+Meta['phase']+'\n')
2338    fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n')
2339    fl.write('INVESTIGATOR :: '+Meta['owner']+'\n')
2340    minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin])
2341    minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin])
2342    maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin])
2343    fl.write('MINIMUM_DISTANCES ::   %s  Angstrom\n'%minHD)
2344    maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves])
2345    fl.write('MAXIMUM_MOVES ::   %s Angstrom\n'%maxMv)
2346    fl.write('R_SPACING ::  0.0200 Angstrom\n')
2347    fl.write('PRINT_PERIOD :: 100\n')
2348    fl.write('TIME_LIMIT ::     %.2f MINUTES\n'%Times[0])
2349    fl.write('SAVE_PERIOD ::    %.2f MINUTES\n'%Times[1])
2350    fl.write('\n')
2351    fl.write('ATOMS :: '+' '.join(Atseq)+'\n')
2352    fl.write('\n')
2353    fl.write('FLAGS ::\n')
2354    fl.write('  > NO_MOVEOUT\n')
2355    fl.write('  > NO_SAVE_CONFIGURATIONS\n')
2356    fl.write('  > NO_RESOLUTION_CONVOLUTION\n')
2357    fl.write('\n')
2358    fl.write('INPUT_CONFIGURATION_FORMAT ::  rmc6f\n')
2359    fl.write('SAVE_CONFIGURATION_FORMAT  ::  rmc6f\n')
2360    fl.write('IGNORE_HISTORY_FILE ::\n')
2361    fl.write('\n')
2362    fl.write('DISTANCE_WINDOW ::\n')
2363    fl.write('  > MNDIST :: %s\n'%minD)
2364    fl.write('  > MXDIST :: %s\n'%maxD)
2365    if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']):
2366        fl.write('\n')
2367        fl.write('POTENTIALS ::\n')
2368        fl.write('  > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.'])
2369        fl.write('  > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n')
2370        if len(RMCPdict['Potentials']['Stretch']):
2371            fl.write('  > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search'])
2372            for bond in RMCPdict['Potentials']['Stretch']:
2373                fl.write('  > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2]))       
2374        if len(RMCPdict['Potentials']['Angles']):
2375            fl.write('  > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
2376            for angle in RMCPdict['Potentials']['Angles']:
2377                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
2378                    (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
2379    if RMCPdict['useBVS']:
2380        fl.write('BVS ::\n')
2381        fl.write('  > ATOM :: '+' '.join(Atseq)+'\n')
2382        fl.write('  > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2383        oxid = []
2384        for val in RMCPdict['Oxid']:
2385            if len(val) == 3:
2386                oxid.append(val[0][1:])
2387            else:
2388                oxid.append(val[0][2:])
2389        fl.write('  > OXID :: %s\n'%' '.join(oxid))
2390        fl.write('  > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']]))
2391        fl.write('  > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']]))
2392        fl.write('  > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))       
2393        fl.write('  > SAVE :: 100000\n')
2394        fl.write('  > UPDATE :: 100000\n')
2395        if len(RMCPdict['Swap']):
2396            fl.write('\n')
2397            fl.write('SWAP_MULTI ::\n')
2398            for swap in RMCPdict['Swap']:
2399                try:
2400                    at1 = Atseq.index(swap[0])
2401                    at2 = Atseq.index(swap[1])
2402                except ValueError:
2403                    break
2404                fl.write('  > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2]))
2405       
2406    if len(RMCPdict['FxCN']):
2407        fl.write('FIXED_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['FxCN']))       
2408        for ifx,fxcn in enumerate(RMCPdict['FxCN']):
2409            try:
2410                at1 = Atseq.index(fxcn[0])
2411                at2 = Atseq.index(fxcn[1])
2412            except ValueError:
2413                break
2414            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]))
2415    if len(RMCPdict['AveCN']):
2416        fl.write('AVERAGE_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['AveCN']))
2417        for iav,avcn in enumerate(RMCPdict['AveCN']):
2418            try:
2419                at1 = Atseq.index(avcn[0])
2420                at2 = Atseq.index(avcn[1])
2421            except ValueError:
2422                break
2423            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]))
2424    for File in Files:
2425        if Files[File][0]:
2426            if 'Xray' in File and 'F(Q)' in File:
2427                fqdata = open(Files[File][0],'r')
2428                lines = int(fqdata.readline()[:-1])
2429            fl.write('\n')
2430            fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_'))
2431            fl.write('  > FILENAME :: %s\n'%Files[File][0])
2432            fl.write('  > DATA_TYPE :: %s\n'%Files[File][2])
2433            fl.write('  > FIT_TYPE :: %s\n'%Files[File][2])
2434            if 'Xray' not in File:
2435                fl.write('  > START_POINT :: 1\n')
2436                fl.write('  > END_POINT :: 3000\n')
2437                fl.write('  > WEIGHT :: %.4f\n'%Files[File][1])
2438            fl.write('  > CONSTANT_OFFSET 0.000\n')
2439            fl.write('  > NO_FITTED_OFFSET\n')
2440            if RMCPdict['FitScale']:
2441                fl.write('  > FITTED_SCALE\n')
2442            else:
2443                fl.write('  > NO_FITTED_SCALE\n')
2444            if Files[File][3] !='RMC':
2445                fl.write('  > %s\n'%Files[File][3])
2446            if 'reciprocal' in File:
2447                fl.write('  > CONVOLVE ::\n')
2448                if 'Xray' in File:
2449                    fl.write('  > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines)
2450                    fl.write('  > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1]))
2451                    fl.write('  > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2))
2452                    fl.write('  > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,1./Files[File][1]))
2453    fl.write('\n')
2454    fl.write('BRAGG ::\n')
2455    fl.write('  > BRAGG_SHAPE :: %s\n'%gsasType)
2456    fl.write('  > RECALCUATE\n')
2457    fl.write('  > DMIN :: %.2f\n'%(dMin-0.02))
2458    fl.write('  > WEIGHT :: %10.3f\n'%BraggWt)
2459    fl.write('\n')
2460    fl.write('END  ::\n')
2461    fl.close()
2462    return fname
2463
2464def MakefullrmcRun(Name,Phase,RMCPdict):
2465    print(' nothing happened yet')
2466    rname = Name+'-run.py'
2467    rundata = ''
2468    rundata += '#### fullrmc $s file; edit by hand if you so choose #####\n'&rname
2469    rundata += '''
2470import numpy as np
2471
2472# fullrmc library imports
2473from fullrmc.Globals import LOGGER, FLOAT_TYPE
2474from fullrmc.Engine import Engine
2475from fullrmc.Constraints.PairDistributionConstraints import PairDistributionConstraint
2476from fullrmc.Constraints.PairCorrelationConstraints import PairCorrelationConstraint
2477from fullrmc.Constraints.DistanceConstraints import InterMolecularDistanceConstraint
2478from fullrmc.Constraints.BondConstraints import BondConstraint
2479from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
2480from fullrmc.Constraints.ImproperAngleConstraints import ImproperAngleConstraint
2481from fullrmc.Core.Collection import convert_Gr_to_gr
2482from fullrmc.Core.MoveGenerator import MoveGeneratorCollector
2483from fullrmc.Core.GroupSelector import RecursiveGroupSelector
2484from fullrmc.Selectors.RandomSelectors import RandomSelector
2485from fullrmc.Selectors.OrderedSelectors import DefinedOrderSelector
2486from fullrmc.Generators.Translations import TranslationGenerator, TranslationAlongSymmetryAxisGenerator
2487from fullrmc.Generators.Rotations import RotationGenerator, RotationAboutSymmetryAxisGenerator
2488from fullrmc.Generators.Agitations import DistanceAgitationGenerator, AngleAgitationGenerator
2489
2490# engine file names - replace name with phase name
2491engineFileName = "name_engine.rmc"
2492expFileName    = "name_pdf.exp"
2493pdbFileName    = "nme.pdb"
2494freshStart     = False      #make TRUE for a restart
2495
2496            '''
2497    rfile = open(rname,'w')
2498    rfile.writelines(rundata)
2499    rfile.close()
2500   
2501
2502def MakePDB(Name,Phase,RMCPdict):
2503    return None
2504    # generalData = Phase['General']
2505    # Cell = generalData['Cell'][1:7]
2506    # Trans = np.eye(3)*np.array(Supercell)
2507    # newPhase = copy.deepcopy(Phase)
2508    # newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2509    # newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans.T)
2510    # newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False)
2511    # Atoms = newPhase['Atoms']
2512    # Cell = newPhase['General']['Cell'][1:7]
2513    # A,B = G2lat. cell2AB(Cell)
2514    # fname = Name+'.pdb'
2515    # fl = open(fname,'w')
2516    # fl.write('REMARK    this file is generated using GSASII\n')
2517    # fl.write('CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n'%(
2518    #         Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2519    # fl.write('ORIGX1      1.000000  0.000000  0.000000        0.00000\n')
2520    # fl.write('ORIGX2      0.000000  1.000000  0.000000        0.00000\n')
2521    # fl.write('ORIGX3      0.000000  0.000000  1.000000        0.00000\n')
2522
2523    # Natm = np.core.defchararray.count(np.array(Atcodes),'+')
2524    # Natm = np.count_nonzero(Natm-1)
2525    # nat = 0
2526    # for atm in Atseq:
2527    #     for iat,atom in enumerate(Atoms):
2528    #         if atom[1] == atm:
2529    #             nat += 1
2530    #             XYZ = np.inner(A,np.array(atom[3:6])-0.5)    #shift origin to middle & make Cartesian
2531    #             fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %-2s\n'%(       
2532    #                     nat,atom[0],nat,XYZ[0],XYZ[1],XYZ[2],atom[1]))
2533    # fl.close()
2534    # return fname
2535
2536def GetRMCBonds(general,RMCPdict,Atoms,bondList):
2537    bondDist = []
2538    Cell = general['Cell'][1:7]
2539    Supercell =  RMCPdict['SuperCell']
2540    Trans = np.eye(3)*np.array(Supercell)
2541    Cell = G2lat.TransformCell(Cell,Trans)[:6]
2542    Amat,Bmat = G2lat.cell2AB(Cell)
2543    indices = (-1,0,1)
2544    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
2545    for bonds in bondList:
2546        Oxyz = np.array(Atoms[bonds[0]][1:])
2547        Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])       
2548        Dx = np.array([Txyz-Oxyz+unit for unit in Units])
2549        Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
2550        for dx in Dx.T:
2551            bondDist.append(np.min(dx))
2552    return np.array(bondDist)
2553   
2554def GetRMCAngles(general,RMCPdict,Atoms,angleList):
2555    bondAngles = []
2556    Cell = general['Cell'][1:7]
2557    Supercell =  RMCPdict['SuperCell']
2558    Trans = np.eye(3)*np.array(Supercell)
2559    Cell = G2lat.TransformCell(Cell,Trans)[:6]
2560    Amat,Bmat = G2lat.cell2AB(Cell)
2561    indices = (-1,0,1)
2562    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
2563    for angle in angleList:
2564        Oxyz = np.array(Atoms[angle[0]][1:])
2565        TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])       
2566        TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])       
2567        DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
2568        DAx = np.sqrt(np.sum(DAxV**2,axis=2))
2569        DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
2570        DBx = np.sqrt(np.sum(DBxV**2,axis=2))
2571        iDAx = np.argmin(DAx,axis=0)
2572        iDBx = np.argmin(DBx,axis=0)
2573        for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
2574            DAv = DAxV[iA,i]/DAx[iA,i]
2575            DBv = DBxV[iB,i]/DBx[iB,i]
2576            bondAngles.append(npacosd(np.sum(DAv*DBv)))
2577    return np.array(bondAngles)
2578   
2579################################################################################
2580#### Reflectometry calculations
2581################################################################################
2582
2583def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
2584    G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
2585   
2586    class RandomDisplacementBounds(object):
2587        """random displacement with bounds"""
2588        def __init__(self, xmin, xmax, stepsize=0.5):
2589            self.xmin = xmin
2590            self.xmax = xmax
2591            self.stepsize = stepsize
2592   
2593        def __call__(self, x):
2594            """take a random step but ensure the new position is within the bounds"""
2595            while True:
2596                # this could be done in a much more clever way, but it will work for example purposes
2597                steps = self.xmax-self.xmin
2598                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
2599                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
2600                    break
2601            return xnew
2602   
2603    def GetModelParms():
2604        parmDict = {}
2605        varyList = []
2606        values = []
2607        bounds = []
2608        parmDict['dQ type'] = data['dQ type']
2609        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
2610        for parm in ['Scale','FltBack']:
2611            parmDict[parm] = data[parm][0]
2612            if data[parm][1]:
2613                varyList.append(parm)
2614                values.append(data[parm][0])
2615                bounds.append(Bounds[parm])
2616        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
2617        parmDict['nLayers'] = len(parmDict['Layer Seq'])
2618        for ilay,layer in enumerate(data['Layers']):
2619            name = layer['Name']
2620            cid = str(ilay)+';'
2621            parmDict[cid+'Name'] = name
2622            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2623                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
2624                if layer.get(parm,[0.,False])[1]:
2625                    varyList.append(cid+parm)
2626                    value = layer[parm][0]
2627                    values.append(value)
2628                    if value:
2629                        bound = [value*Bfac,value/Bfac]
2630                    else:
2631                        bound = [0.,10.]
2632                    bounds.append(bound)
2633            if name not in ['vacuum','unit scatter']:
2634                parmDict[cid+'rho'] = Substances[name]['Scatt density']
2635                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
2636        return parmDict,varyList,values,bounds
2637   
2638    def SetModelParms():
2639        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
2640        if 'Scale' in varyList:
2641            data['Scale'][0] = parmDict['Scale']
2642            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
2643        G2fil.G2Print (line)
2644        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
2645        if 'FltBack' in varyList:
2646            data['FltBack'][0] = parmDict['FltBack']
2647            line += ' esd: %15.3g'%(sigDict['FltBack'])
2648        G2fil.G2Print (line)
2649        for ilay,layer in enumerate(data['Layers']):
2650            name = layer['Name']
2651            G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
2652            cid = str(ilay)+';'
2653            line = ' '
2654            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
2655            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
2656            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2657                if parm in layer:
2658                    if parm == 'Rough':
2659                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
2660                    else:
2661                        layer[parm][0] = parmDict[cid+parm]
2662                    line += ' %s: %.3f'%(parm,layer[parm][0])
2663                    if cid+parm in varyList:
2664                        line += ' esd: %.3g'%(sigDict[cid+parm])
2665            G2fil.G2Print (line)
2666            G2fil.G2Print (line2)
2667   
2668    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2669        parmDict.update(zip(varyList,values))
2670        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2671        return M
2672   
2673    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2674        parmDict.update(zip(varyList,values))
2675        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2676        return np.sum(M**2)
2677   
2678    def getREFD(Q,Qsig,parmDict):
2679        Ic = np.ones_like(Q)*parmDict['FltBack']
2680        Scale = parmDict['Scale']
2681        Nlayers = parmDict['nLayers']
2682        Res = parmDict['Res']
2683        depth = np.zeros(Nlayers)
2684        rho = np.zeros(Nlayers)
2685        irho = np.zeros(Nlayers)
2686        sigma = np.zeros(Nlayers)
2687        for ilay,lay in enumerate(parmDict['Layer Seq']):
2688            cid = str(lay)+';'
2689            depth[ilay] = parmDict[cid+'Thick']
2690            sigma[ilay] = parmDict[cid+'Rough']
2691            if parmDict[cid+'Name'] == u'unit scatter':
2692                rho[ilay] = parmDict[cid+'DenMul']
2693                irho[ilay] = parmDict[cid+'iDenMul']
2694            elif 'vacuum' != parmDict[cid+'Name']:
2695                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
2696                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
2697            if cid+'Mag SLD' in parmDict:
2698                rho[ilay] += parmDict[cid+'Mag SLD']
2699        if parmDict['dQ type'] == 'None':
2700            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2701        elif 'const' in parmDict['dQ type']:
2702            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
2703        else:       #dQ/Q in data
2704            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
2705        Ic += AB*Scale
2706        return Ic
2707       
2708    def estimateT0(takestep):
2709        Mmax = -1.e-10
2710        Mmin = 1.e10
2711        for i in range(100):
2712            x0 = takestep(values)
2713            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2714            Mmin = min(M,Mmin)
2715            MMax = max(M,Mmax)
2716        return 1.5*(MMax-Mmin)
2717
2718    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2719    if data.get('2% weight'):
2720        wt = 1./(0.02*Io)**2
2721    Qmin = Limits[1][0]
2722    Qmax = Limits[1][1]
2723    wtFactor = ProfDict['wtFactor']
2724    Bfac = data['Toler']
2725    Ibeg = np.searchsorted(Q,Qmin)
2726    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
2727    Ic[:] = 0
2728    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
2729              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
2730    parmDict,varyList,values,bounds = GetModelParms()
2731    Msg = 'Failed to converge'
2732    if varyList:
2733        if data['Minimizer'] == 'LMLS': 
2734            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
2735                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2736            parmDict.update(zip(varyList,result[0]))
2737            chisq = np.sum(result[2]['fvec']**2)
2738            ncalc = result[2]['nfev']
2739            covM = result[1]
2740            newVals = result[0]
2741        elif data['Minimizer'] == 'Basin Hopping':
2742            xyrng = np.array(bounds).T
2743            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
2744            T0 = estimateT0(take_step)
2745            G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
2746            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
2747                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
2748                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
2749            chisq = result.fun
2750            ncalc = result.nfev
2751            newVals = result.x
2752            covM = []
2753        elif data['Minimizer'] == 'MC/SA Anneal':
2754            xyrng = np.array(bounds).T
2755            result = G2mth.anneal(sumREFD, values, 
2756                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
2757                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
2758                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
2759                ranRange=0.20,autoRan=False,dlg=None)
2760            newVals = result[0]
2761            parmDict.update(zip(varyList,newVals))
2762            chisq = result[1]
2763            ncalc = result[3]
2764            covM = []
2765            G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
2766        elif data['Minimizer'] == 'L-BFGS-B':
2767            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
2768                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2769            parmDict.update(zip(varyList,result['x']))
2770            chisq = result.fun
2771            ncalc = result.nfev
2772            newVals = result.x
2773            covM = []
2774    else:   #nothing varied
2775        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2776        chisq = np.sum(M**2)
2777        ncalc = 0
2778        covM = []
2779        sig = []
2780        sigDict = {}
2781        result = []
2782    Rvals = {}
2783    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
2784    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
2785    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
2786    Ib[Ibeg:Ifin] = parmDict['FltBack']
2787    try:
2788        if not len(varyList):
2789            Msg += ' - nothing refined'
2790            raise ValueError
2791        Nans = np.isnan(newVals)
2792        if np.any(Nans):
2793            name = varyList[Nans.nonzero(True)[0]]
2794            Msg += ' Nan result for '+name+'!'
2795            raise ValueError
2796        Negs = np.less_equal(newVals,0.)
2797        if np.any(Negs):
2798            indx = Negs.nonzero()
2799            name = varyList[indx[0][0]]
2800            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
2801                Msg += ' negative coefficient for '+name+'!'
2802                raise ValueError
2803        if len(covM):
2804            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
2805            covMatrix = covM*Rvals['GOF']
2806        else:
2807            sig = np.zeros(len(varyList))
2808            covMatrix = []
2809        sigDict = dict(zip(varyList,sig))
2810        G2fil.G2Print (' Results of reflectometry data modelling fit:')
2811        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
2812        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2813        SetModelParms()
2814        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
2815    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
2816        G2fil.G2Print (Msg)
2817        return False,0,0,0,0,0,0,Msg
2818       
2819def makeSLDprofile(data,Substances):
2820   
2821    sq2 = np.sqrt(2.)
2822    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2823    Nlayers = len(laySeq)
2824    laySeq = np.array(laySeq,dtype=int)
2825    interfaces = np.zeros(Nlayers)
2826    rho = np.zeros(Nlayers)
2827    sigma = np.zeros(Nlayers)
2828    name = data['Layers'][0]['Name']
2829    thick = 0.
2830    for ilay,lay in enumerate(laySeq):
2831        layer = data['Layers'][lay]
2832        name = layer['Name']
2833        if 'Thick' in layer:
2834            thick += layer['Thick'][0]
2835            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
2836        if 'Rough' in layer:
2837            sigma[ilay] = max(0.001,layer['Rough'][0])
2838        if name != 'vacuum':
2839            if name == 'unit scatter':
2840                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
2841            else:
2842                rrho = Substances[name]['Scatt density']
2843                irho = Substances[name]['XImag density']
2844                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
2845        if 'Mag SLD' in layer:
2846            rho[ilay] += layer['Mag SLD'][0]
2847    name = data['Layers'][-1]['Name']
2848    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
2849    xr = np.flipud(x)
2850    interfaces[-1] = x[-1]
2851    y = np.ones_like(x)*rho[0]
2852    iBeg = 0
2853    for ilayer in range(Nlayers-1):
2854        delt = rho[ilayer+1]-rho[ilayer]
2855        iPos = np.searchsorted(x,interfaces[ilayer])
2856        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
2857        iBeg = iPos
2858    return x,xr,y   
2859
2860def REFDModelFxn(Profile,Inst,Limits,Substances,data):
2861   
2862    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2863    Qmin = Limits[1][0]
2864    Qmax = Limits[1][1]
2865    iBeg = np.searchsorted(Q,Qmin)
2866    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2867    Ib[:] = data['FltBack'][0]
2868    Ic[:] = 0
2869    Scale = data['Scale'][0]
2870    if data['Layer Seq'] == []:
2871        return
2872    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2873    Nlayers = len(laySeq)
2874    depth = np.zeros(Nlayers)
2875    rho = np.zeros(Nlayers)
2876    irho = np.zeros(Nlayers)
2877    sigma = np.zeros(Nlayers)
2878    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
2879        layer = data['Layers'][lay]
2880        name = layer['Name']
2881        if 'Thick' in layer:    #skips first & last layers
2882            depth[ilay] = layer['Thick'][0]
2883        if 'Rough' in layer:    #skips first layer
2884            sigma[ilay] = layer['Rough'][0]
2885        if 'unit scatter' == name:
2886            rho[ilay] = layer['DenMul'][0]
2887            irho[ilay] = layer['iDenMul'][0]
2888        else:
2889            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2890            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
2891        if 'Mag SLD' in layer:
2892            rho[ilay] += layer['Mag SLD'][0]
2893    if data['dQ type'] == 'None':
2894        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2895    elif 'const' in data['dQ type']:
2896        res = data['Resolution'][0]/(100.*sateln2)
2897        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
2898    else:       #dQ/Q in data
2899        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
2900    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
2901
2902def abeles(kz, depth, rho, irho=0, sigma=0):
2903    """
2904    Optical matrix form of the reflectivity calculation.
2905    O.S. Heavens, Optical Properties of Thin Solid Films
2906   
2907    Reflectometry as a function of kz for a set of slabs.
2908
2909    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
2910        This is :math:`\\tfrac12 Q_z`.       
2911    :param depth:  float[m] (Ang).
2912        thickness of each layer.  The thickness of the incident medium
2913        and substrate are ignored.
2914    :param rho:  float[n,k] (1e-6/Ang^2)
2915        Real scattering length density for each layer for each kz
2916    :param irho:  float[n,k] (1e-6/Ang^2)
2917        Imaginary scattering length density for each layer for each kz
2918        Note: absorption cross section mu = 2 irho/lambda for neutrons
2919    :param sigma: float[m-1] (Ang)
2920        interfacial roughness.  This is the roughness between a layer
2921        and the previous layer. The sigma array should have m-1 entries.
2922
2923    Slabs are ordered with the surface SLD at index 0 and substrate at
2924    index -1, or reversed if kz < 0.
2925    """
2926    def calc(kz, depth, rho, irho, sigma):
2927        if len(kz) == 0: return kz
2928   
2929        # Complex index of refraction is relative to the incident medium.
2930        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
2931        # in place of kz^2, and ignoring rho_o
2932        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
2933        k = kz
2934   
2935        # According to Heavens, the initial matrix should be [ 1 F; F 1],
2936        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
2937        # multiply versus some coding convenience.
2938        B11 = 1
2939        B22 = 1
2940        B21 = 0
2941        B12 = 0
2942        for i in range(0, len(depth)-1):
2943            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
2944            F = (k - k_next) / (k + k_next)
2945            F *= np.exp(-2*k*k_next*sigma[i]**2)
2946            #print "==== layer",i
2947            #print "kz:", kz
2948            #print "k:", k
2949            #print "k_next:",k_next
2950            #print "F:",F
2951            #print "rho:",rho[:,i+1]
2952            #print "irho:",irho[:,i+1]
2953            #print "d:",depth[i],"sigma:",sigma[i]
2954            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
2955            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
2956            M21 = F*M11
2957            M12 = F*M22
2958            C1 = B11*M11 + B21*M12
2959            C2 = B11*M21 + B21*M22
2960            B11 = C1
2961            B21 = C2
2962            C1 = B12*M11 + B22*M12
2963            C2 = B12*M21 + B22*M22
2964            B12 = C1
2965            B22 = C2
2966            k = k_next
2967   
2968        r = B12/B11
2969        return np.absolute(r)**2
2970
2971    if np.isscalar(kz): kz = np.asarray([kz], 'd')
2972
2973    m = len(depth)
2974
2975    # Make everything into arrays
2976    depth = np.asarray(depth,'d')
2977    rho = np.asarray(rho,'d')
2978    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
2979    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
2980
2981    # Repeat rho,irho columns as needed
2982    if len(rho.shape) == 1:
2983        rho = rho[None,:]
2984        irho = irho[None,:]
2985
2986    return calc(kz, depth, rho, irho, sigma)
2987   
2988def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
2989    y = abeles(kz, depth, rho, irho, sigma)
2990    s = dq/2.
2991    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
2992    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
2993    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
2994    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
2995    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
2996    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
2997    y *= 0.137023
2998    return y
2999       
3000def makeRefdFFT(Limits,Profile):
3001    G2fil.G2Print ('make fft')
3002    Q,Io = Profile[:2]
3003    Qmin = Limits[1][0]
3004    Qmax = Limits[1][1]
3005    iBeg = np.searchsorted(Q,Qmin)
3006    iFin = np.searchsorted(Q,Qmax)+1    #include last point
3007    Qf = np.linspace(0.,Q[iFin-1],5000)
3008    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
3009    If = QI(Qf)*Qf**4
3010    R = np.linspace(0.,5000.,5000)
3011    F = fft.rfft(If)
3012    return R,F
3013
3014   
3015################################################################################
3016#### Stacking fault simulation codes
3017################################################################################
3018
3019def GetStackParms(Layers):
3020   
3021    Parms = []
3022#cell parms
3023    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
3024        Parms.append('cellA')
3025        Parms.append('cellC')
3026    else:
3027        Parms.append('cellA')
3028        Parms.append('cellB')
3029        Parms.append('cellC')
3030        if Layers['Laue'] != 'mmm':
3031            Parms.append('cellG')
3032#Transition parms
3033    for iY in range(len(Layers['Layers'])):
3034        for iX in range(len(Layers['Layers'])):
3035            Parms.append('TransP;%d;%d'%(iY,iX))
3036            Parms.append('TransX;%d;%d'%(iY,iX))
3037            Parms.append('TransY;%d;%d'%(iY,iX))
3038            Parms.append('TransZ;%d;%d'%(iY,iX))
3039    return Parms
3040
3041def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
3042    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
3043   
3044    :param dict Layers: dict with following items
3045
3046      ::
3047
3048       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
3049       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
3050       'Layers':[],'Stacking':[],'Transitions':[]}
3051       
3052    :param str ctrls: controls string to be written on DIFFaX controls.dif file
3053    :param float scale: scale factor
3054    :param dict background: background parameters
3055    :param list limits: min/max 2-theta to be calculated
3056    :param dict inst: instrument parameters dictionary
3057    :param list profile: powder pattern data
3058   
3059    Note that parameters all updated in place   
3060    '''
3061    import atmdata
3062    path = sys.path
3063    for name in path:
3064        if 'bin' in name:
3065            DIFFaX = name+'/DIFFaX.exe'
3066            G2fil.G2Print (' Execute '+DIFFaX)
3067            break
3068    # make form factor file that DIFFaX wants - atom types are GSASII style
3069    sf = open('data.sfc','w')
3070    sf.write('GSASII special form factor file for DIFFaX\n\n')
3071    atTypes = list(Layers['AtInfo'].keys())
3072    if 'H' not in atTypes:
3073        atTypes.insert(0,'H')
3074    for atType in atTypes:
3075        if atType == 'H': 
3076            blen = -.3741
3077        else:
3078            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
3079        Adat = atmdata.XrayFF[atType]
3080        text = '%4s'%(atType.ljust(4))
3081        for i in range(4):
3082            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
3083        text += '%11.6f%11.6f'%(Adat['fc'],blen)
3084        text += '%3d\n'%(Adat['Z'])
3085        sf.write(text)
3086    sf.close()
3087    #make DIFFaX control.dif file - future use GUI to set some of these flags
3088    cf = open('control.dif','w')
3089    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
3090        x0 = profile[0]
3091        iBeg = np.searchsorted(x0,limits[0])
3092        iFin = np.searchsorted(x0,limits[1])+1
3093        if iFin-iBeg > 20000:
3094            iFin = iBeg+20000
3095        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3096        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3097        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
3098    else:
3099        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3100        inst = {'Type':['XSC','XSC',]}
3101    cf.close()
3102    #make DIFFaX data file
3103    df = open('GSASII-DIFFaX.dat','w')
3104    df.write('INSTRUMENTAL\n')
3105    if 'X' in inst['Type'][0]:
3106        df.write('X-RAY\n')
3107    elif 'N' in inst['Type'][0]:
3108        df.write('NEUTRON\n')
3109    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
3110        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
3111        U = ateln2*inst['U'][1]/10000.
3112        V = ateln2*inst['V'][1]/10000.
3113        W = ateln2*inst['W'][1]/10000.
3114        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
3115        HW = np.sqrt(np.mean(HWHM))
3116    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
3117        if 'Mean' in Layers['selInst']:
3118            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
3119        elif 'Gaussian' in Layers['selInst']:
3120            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
3121        else:
3122            df.write('None\n')
3123    else:
3124        df.write('0.10\nNone\n')
3125    df.write('STRUCTURAL\n')
3126    a,b,c = Layers['Cell'][1:4]
3127    gam = Layers['Cell'][6]
3128    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
3129    laue = Layers['Laue']
3130    if laue == '2/m(ab)':
3131        laue = '2/m(1)'
3132    elif laue == '2/m(c)':
3133        laue = '2/m(2)'
3134    if 'unknown' in Layers['Laue']:
3135        df.write('%s %.3f\n'%(laue,Layers['Toler']))
3136    else:   
3137        df.write('%s\n'%(laue))
3138    df.write('%d\n'%(len(Layers['Layers'])))
3139    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
3140        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
3141    layerNames = []
3142    for layer in Layers['Layers']:
3143        layerNames.append(layer['Name'])
3144    for il,layer in enumerate(Layers['Layers']):
3145        if layer['SameAs']:
3146            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
3147            continue
3148        df.write('LAYER %d\n'%(il+1))
3149        if '-1' in layer['Symm']:
3150            df.write('CENTROSYMMETRIC\n')
3151        else:
3152            df.write('NONE\n')
3153        for ia,atom in enumerate(layer['Atoms']):
3154            [name,atype,x,y,z,frac,Uiso] = atom
3155            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
3156                frac /= 2.
3157            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
3158    df.write('STACKING\n')
3159    df.write('%s\n'%(Layers['Stacking'][0]))
3160    if 'recursive' in Layers['Stacking'][0]:
3161        df.write('%s\n'%Layers['Stacking'][1])
3162    else:
3163        if 'list' in Layers['Stacking'][1]:
3164            Slen = len(Layers['Stacking'][2])
3165            iB = 0
3166            iF = 0
3167            while True:
3168                iF += 68
3169                if iF >= Slen:
3170                    break
3171                iF = min(iF,Slen)
3172                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
3173                iB = iF
3174        else:
3175            df.write('%s\n'%Layers['Stacking'][1])   
3176    df.write('TRANSITIONS\n')
3177    for iY in range(len(Layers['Layers'])):
3178        sumPx = 0.
3179        for iX in range(len(Layers['Layers'])):
3180            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
3181            p = round(p,3)
3182            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
3183            sumPx += p
3184        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
3185            G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
3186            df.close()
3187            os.remove('data.sfc')
3188            os.remove('control.dif')
3189            os.remove('GSASII-DIFFaX.dat')
3190            return       
3191    df.close()   
3192    time0 = time.time()
3193    try:
3194        subp.call(DIFFaX)
3195    except OSError:
3196        G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
3197    G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
3198    if os.path.exists('GSASII-DIFFaX.spc'):
3199        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
3200        iFin = iBeg+Xpat.shape[1]
3201        bakType,backDict,backVary = SetBackgroundParms(background)
3202        backDict['Lam1'] = G2mth.getWave(inst)
3203        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
3204        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
3205        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
3206            rv = st.poisson(profile[3][iBeg:iFin])
3207            profile[1][iBeg:iFin] = rv.rvs()
3208            Z = np.ones_like(profile[3][iBeg:iFin])
3209            Z[1::2] *= -1
3210            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3211            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3212        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3213    #cleanup files..
3214        os.remove('GSASII-DIFFaX.spc')
3215    elif os.path.exists('GSASII-DIFFaX.sadp'):
3216        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
3217        Sadp = np.reshape(Sadp,(256,-1))
3218        Layers['Sadp']['Img'] = Sadp
3219        os.remove('GSASII-DIFFaX.sadp')
3220    os.remove('data.sfc')
3221    os.remove('control.dif')
3222    os.remove('GSASII-DIFFaX.dat')
3223   
3224def SetPWDRscan(inst,limits,profile):
3225   
3226    wave = G2mth.getMeanWave(inst)
3227    x0 = profile[0]
3228    iBeg = np.searchsorted(x0,limits[0])
3229    iFin = np.searchsorted(x0,limits[1])
3230    if iFin-iBeg > 20000:
3231        iFin = iBeg+20000
3232    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3233    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
3234    return iFin-iBeg
3235       
3236def SetStackingSF(Layers,debug):
3237# Load scattering factors into DIFFaX arrays
3238    import atmdata
3239    atTypes = Layers['AtInfo'].keys()
3240    aTypes = []
3241    for atype in atTypes:
3242        aTypes.append('%4s'%(atype.ljust(4)))
3243    SFdat = []
3244    for atType in atTypes:
3245        Adat = atmdata.XrayFF[atType]
3246        SF = np.zeros(9)
3247        SF[:8:2] = Adat['fa']
3248        SF[1:8:2] = Adat['fb']
3249        SF[8] = Adat['fc']
3250        SFdat.append(SF)
3251    SFdat = np.array(SFdat)
3252    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
3253   
3254def SetStackingClay(Layers,Type):
3255# Controls
3256    rand.seed()
3257    ranSeed = rand.randint(1,2**16-1)
3258    try:   
3259        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
3260            '6/m','6/mmm'].index(Layers['Laue'])+1
3261    except ValueError:  #for 'unknown'
3262        laueId = -1
3263    if 'SADP' in Type:
3264        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
3265        lmax = int(Layers['Sadp']['Lmax'])
3266    else:
3267        planeId = 0
3268        lmax = 0
3269# Sequences
3270    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
3271    try:
3272        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
3273    except ValueError:
3274        StkParm = -1
3275    if StkParm == 2:    #list
3276        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
3277        Nstk = len(StkSeq)
3278    else:
3279        Nstk = 1
3280        StkSeq = [0,]
3281    if StkParm == -1:
3282        StkParm = int(Layers['Stacking'][1])
3283    Wdth = Layers['Width'][0]
3284    mult = 1
3285    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
3286    LaueSym = Layers['Laue'].ljust(12)
3287    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
3288    return laueId,controls
3289   
3290def SetCellAtoms(Layers):
3291    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
3292# atoms in layers
3293    atTypes = list(Layers['AtInfo'].keys())
3294    AtomXOU = []
3295    AtomTp = []
3296    LayerSymm = []
3297    LayerNum = []
3298    layerNames = []
3299    Natm = 0
3300    Nuniq = 0
3301    for layer in Layers['Layers']:
3302        layerNames.append(layer['Name'])
3303    for il,layer in enumerate(Layers['Layers']):
3304        if layer['SameAs']:
3305            LayerNum.append(layerNames.index(layer['SameAs'])+1)
3306            continue
3307        else:
3308            LayerNum.append(il+1)
3309            Nuniq += 1
3310        if '-1' in layer['Symm']:
3311            LayerSymm.append(1)
3312        else:
3313            LayerSymm.append(0)
3314        for ia,atom in enumerate(layer['Atoms']):
3315            [name,atype,x,y,z,frac,Uiso] = atom
3316            Natm += 1
3317            AtomTp.append('%4s'%(atype.ljust(4)))
3318            Ta = atTypes.index(atype)+1
3319            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
3320    AtomXOU = np.array(AtomXOU)
3321    Nlayers = len(layerNames)
3322    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
3323    return Nlayers
3324   
3325def SetStackingTrans(Layers,Nlayers):
3326# Transitions
3327    TransX = []
3328    TransP = []
3329    for Ytrans in Layers['Transitions']:
3330        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
3331        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
3332    TransP = np.array(TransP,dtype='float').T
3333    TransX = np.array(TransX,dtype='float')
3334#    GSASIIpath.IPyBreak()
3335    pyx.pygettrans(Nlayers,TransP,TransX)
3336   
3337def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
3338# Scattering factors
3339    SetStackingSF(Layers,debug)
3340# Controls & sequences
3341    laueId,controls = SetStackingClay(Layers,'PWDR')
3342# cell & atoms
3343    Nlayers = SetCellAtoms(Layers)
3344    Volume = Layers['Cell'][7]   
3345# Transitions
3346    SetStackingTrans(Layers,Nlayers)
3347# PWDR scan
3348    Nsteps = SetPWDRscan(inst,limits,profile)
3349# result as Spec
3350    x0 = profile[0]
3351    profile[3] = np.zeros(len(profile[0]))
3352    profile[4] = np.zeros(len(profile[0]))
3353    profile[5] = np.zeros(len(profile[0]))
3354    iBeg = np.searchsorted(x0,limits[0])
3355    iFin = np.searchsorted(x0,limits[1])+1
3356    if iFin-iBeg > 20000:
3357        iFin = iBeg+20000
3358    Nspec = 20001       
3359    spec = np.zeros(Nspec,dtype='double')   
3360    time0 = time.time()
3361    pyx.pygetspc(controls,Nspec,spec)
3362    G2fil.G2Print (' GETSPC time = %.2fs'%(time.time()-time0))
3363    time0 = time.time()
3364    U = ateln2*inst['U'][1]/10000.
3365    V = ateln2*inst['V'][1]/10000.
3366    W = ateln2*inst['W'][1]/10000.
3367    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
3368    HW = np.sqrt(np.mean(HWHM))
3369    BrdSpec = np.zeros(Nsteps)
3370    if 'Mean' in Layers['selInst']:
3371        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
3372    elif 'Gaussian' in Layers['selInst']:
3373        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
3374    else:
3375        BrdSpec = spec[:Nsteps]
3376    BrdSpec /= Volume
3377    iFin = iBeg+Nsteps
3378    bakType,backDict,backVary = SetBackgroundParms(background)
3379    backDict['Lam1'] = G2mth.getWave(inst)
3380    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
3381    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
3382    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
3383        try:
3384            rv = st.poisson(profile[3][iBeg:iFin])
3385            profile[1][iBeg:iFin] = rv.rvs()
3386        except ValueError:
3387            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
3388        Z = np.ones_like(profile[3][iBeg:iFin])
3389        Z[1::2] *= -1
3390        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3391        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3392    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3393    G2fil.G2Print (' Broadening time = %.2fs'%(time.time()-time0))
3394   
3395def CalcStackingSADP(Layers,debug):
3396   
3397# Scattering factors
3398    SetStackingSF(Layers,debug)
3399# Controls & sequences
3400    laueId,controls = SetStackingClay(Layers,'SADP')
3401# cell & atoms
3402    Nlayers = SetCellAtoms(Layers)   
3403# Transitions
3404    SetStackingTrans(Layers,Nlayers)
3405# result as Sadp
3406    Nspec = 20001       
3407    spec = np.zeros(Nspec,dtype='double')   
3408    time0 = time.time()
3409    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
3410    Sapd = np.zeros((256,256))
3411    iB = 0
3412    for i in range(hkLim):
3413        iF = iB+Nblk
3414        p1 = 127+int(i*Incr)
3415        p2 = 128-int(i*Incr)
3416        if Nblk == 128:
3417            if i:
3418                Sapd[128:,p1] = spec[iB:iF]
3419                Sapd[:128,p1] = spec[iF:iB:-1]
3420            Sapd[128:,p2] = spec[iB:iF]
3421            Sapd[:128,p2] = spec[iF:iB:-1]
3422        else:
3423            if i:
3424                Sapd[:,p1] = spec[iB:iF]
3425            Sapd[:,p2] = spec[iB:iF]
3426        iB += Nblk
3427    Layers['Sadp']['Img'] = Sapd
3428    G2fil.G2Print (' GETSAD time = %.2fs'%(time.time()-time0))
3429   
3430###############################################################################
3431#### Maximum Entropy Method - Dysnomia
3432###############################################################################
3433   
3434def makePRFfile(data,MEMtype):
3435    ''' makes Dysnomia .prf control file from Dysnomia GUI controls
3436   
3437    :param dict data: GSAS-II phase data
3438    :param int MEMtype: 1 for neutron data with negative scattering lengths
3439                        0 otherwise
3440    :returns str: name of Dysnomia control file
3441    '''
3442
3443    generalData = data['General']
3444    pName = generalData['Name'].replace(' ','_')
3445    DysData = data['Dysnomia']
3446    prfName = pName+'.prf'
3447    prf = open(prfName,'w')
3448    prf.write('$PREFERENCES\n')
3449    prf.write(pName+'.mem\n') #or .fos?
3450    prf.write(pName+'.out\n')
3451    prf.write(pName+'.pgrid\n')
3452    prf.write(pName+'.fba\n')
3453    prf.write(pName+'_eps.raw\n')
3454    prf.write('%d\n'%MEMtype)
3455    if DysData['DenStart'] == 'uniform':
3456        prf.write('0\n')
3457    else:
3458        prf.write('1\n')
3459    if DysData['Optimize'] == 'ZSPA':
3460        prf.write('0\n')
3461    else:
3462        prf.write('1\n')
3463    prf.write('1\n')
3464    if DysData['Lagrange'][0] == 'user':
3465        prf.write('0\n')
3466    else:
3467        prf.write('1\n')
3468    prf.write('%.4f %d\n'%(DysData['Lagrange'][1],DysData['wt pwr']))
3469    prf.write('%.3f\n'%DysData['Lagrange'][2])
3470    prf.write('%.2f\n'%DysData['E_factor'])
3471    prf.write('1\n')
3472    prf.write('0\n')
3473    prf.write('%d\n'%DysData['Ncyc'])
3474    prf.write('1\n')
3475    prf.write('1 0 0 0 0 0 0 0\n')
3476    if DysData['prior'] == 'uniform':
3477        prf.write('0\n')
3478    else:
3479        prf.write('1\n')
3480    prf.close()
3481    return prfName
3482
3483def makeMEMfile(data,reflData,MEMtype,DYSNOMIA):
3484    ''' make Dysnomia .mem file of reflection data, etc.
3485
3486    :param dict data: GSAS-II phase data
3487    :param list reflData: GSAS-II reflection data
3488    :param int MEMtype: 1 for neutron data with negative scattering lengths
3489                        0 otherwise
3490    :param str DYSNOMIA: path to dysnomia.exe
3491    '''
3492   
3493    DysData = data['Dysnomia']
3494    generalData = data['General']
3495    cell = generalData['Cell'][1:7]
3496    A = G2lat.cell2A(cell)
3497    SGData = generalData['SGData']
3498    pName = generalData['Name'].replace(' ','_')
3499    memName = pName+'.mem'
3500    Map = generalData['Map']
3501    Type = Map['Type']
3502    UseList = Map['RefList']
3503    mem = open(memName,'w')
3504    mem.write('%s\n'%(generalData['Name']+' from '+UseList[0]))
3505    a,b,c,alp,bet,gam = cell
3506    mem.write('%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n'%(a,b,c,alp,bet,gam))
3507    mem.write('      0.0000000      0.0000000     -1    0    0    0     P\n')   #dummy PO stuff
3508    SGSym = generalData['SGData']['SpGrp']
3509    try:
3510        SGId = G2spc.spgbyNum.index(SGSym)
3511    except ValueError:
3512        return False
3513    org = 1
3514    if SGSym in G2spc.spg2origins:
3515        org = 2
3516    mapsize = Map['rho'].shape
3517    sumZ = 0.
3518    sumpos = 0.
3519    sumneg = 0.
3520    mem.write('%5d%5d%5d%5d%5d\n'%(SGId,org,mapsize[0],mapsize[1],mapsize[2]))
3521    for atm in generalData['NoAtoms']:
3522        Nat = generalData['NoAtoms'][atm]
3523        AtInfo = G2elem.GetAtomInfo(atm)
3524        sumZ += Nat*AtInfo['Z']
3525        isotope = generalData['Isotope'][atm]
3526        blen = generalData['Isotopes'][atm][isotope]['SL'][0]
3527        if blen < 0.:
3528            sumneg += blen*Nat
3529        else:
3530            sumpos += blen*Nat
3531    if 'X' in Type:
3532        mem.write('%10.2f  0.001\n'%sumZ)
3533    elif 'N' in Type and MEMtype:
3534        mem.write('%10.3f%10.3f 0.001\n'%(sumpos,sumneg))
3535    else:
3536        mem.write('%10.3f 0.001\n'%sumpos)
3537       
3538    dmin = DysData['MEMdmin']
3539    TOFlam = 2.0*dmin*npsind(80.0)
3540    refSet = G2lat.GenHLaue(dmin,SGData,A)      #list of h,k,l,d
3541    refDict = {'%d %d %d'%(ref[0],ref[1],ref[2]):ref for ref in refSet}
3542       
3543    refs = []
3544    prevpos = 0.
3545    for ref in reflData:
3546        if ref[3] < 0:
3547            continue
3548        if 'T' in Type:
3549            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,x,x,x,prfo = ref[:16]
3550            s = np.sqrt(max(sig,0.0001))   #var -> sig in deg
3551            FWHM = getgamFW(gam,s)
3552            if dsp < dmin:
3553                continue
3554            theta = npasind(TOFlam/(2.*dsp))
3555            FWHM *= nptand(theta)/pos
3556            pos = 2.*theta
3557        else:
3558            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,prfo = ref[:13]
3559            g = gam/100.    #centideg -> deg
3560            s = np.sqrt(max(sig,0.0001))/100.   #var -> sig in deg
3561            FWHM = getgamFW(g,s)
3562        delt = pos-prevpos
3563        refs.append([h,k,l,mult,pos,FWHM,Fobs,phase,delt])
3564        prevpos = pos
3565           
3566    ovlp = DysData['overlap']
3567    refs1 = []
3568    refs2 = []
3569    nref2 = 0
3570    iref = 0
3571    Nref = len(refs)
3572    start = False
3573    while iref < Nref-1:
3574        if refs[iref+1][-1] < ovlp*refs[iref][5]:
3575            if refs[iref][-1] > ovlp*refs[iref][5]:
3576                refs2.append([])
3577                start = True
3578            if nref2 == len(refs2):
3579                refs2.append([])
3580            refs2[nref2].append(refs[iref])
3581        else:
3582            if start:
3583                refs2[nref2].append(refs[iref])
3584                start = False
3585                nref2 += 1
3586            else:
3587                refs1.append(refs[iref])
3588        iref += 1
3589    if start:
3590        refs2[nref2].append(refs[iref])
3591    else:
3592        refs1.append(refs[iref])
3593   
3594    mem.write('%5d\n'%len(refs1))
3595    for ref in refs1:
3596        h,k,l = ref[:3]
3597        hkl = '%d %d %d'%(h,k,l)
3598        if hkl in refDict:
3599            del refDict[hkl]
3600        Fobs = np.sqrt(ref[6])
3601        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)))
3602    while True and nref2:
3603        if not len(refs2[-1]):
3604            del refs2[-1]
3605        else:
3606            break
3607    mem.write('%5d\n'%len(refs2))
3608    for iref2,ref2 in enumerate(refs2):
3609        mem.write('#%5d\n'%iref2)
3610        mem.write('%5d\n'%len(ref2))
3611        Gsum = 0.
3612        Msum = 0
3613        for ref in ref2:
3614            Gsum += ref[6]*ref[3]
3615            Msum += ref[3]
3616        G = np.sqrt(Gsum/Msum)
3617        h,k,l = ref2[0][:3]
3618        hkl = '%d %d %d'%(h,k,l)
3619        if hkl in refDict:
3620            del refDict[hkl]
3621        mem.write('%5d%5d%5d%10.3f%10.3f%5d\n'%(h,k,l,G,max(0.01*G,0.1),ref2[0][3]))
3622        for ref in ref2[1:]:
3623            h,k,l,m = ref[:4]
3624            mem.write('%5d%5d%5d%5d\n'%(h,k,l,m))
3625            hkl = '%d %d %d'%(h,k,l)
3626            if hkl in refDict:
3627                del refDict[hkl]
3628    if len(refDict):
3629        mem.write('%d\n'%len(refDict))
3630        for hkl in list(refDict.keys()):
3631            h,k,l = refDict[hkl][:3]
3632            mem.write('%5d%5d%5d\n'%(h,k,l))
3633    else:
3634        mem.write('0\n')
3635    mem.close()
3636    return True
3637
3638def MEMupdateReflData(prfName,data,reflData):
3639    ''' Update reflection data with new Fosq, phase result from Dysnomia
3640
3641    :param str prfName: phase.mem file name
3642    :param list reflData: GSAS-II reflection data
3643    '''
3644   
3645    generalData = data['General']
3646    Map = generalData['Map']
3647    Type = Map['Type']
3648    cell = generalData['Cell'][1:7]
3649    A = G2lat.cell2A(cell)
3650    reflDict = {}
3651    newRefs = []
3652    for iref,ref in enumerate(reflData):
3653        if ref[3] > 0:
3654            newRefs.append(ref)
3655            reflDict[hash('%5d%5d%5d'%(ref[0],ref[1],ref[2]))] = iref
3656    fbaName = os.path.splitext(prfName)[0]+'.fba'
3657    try: # patch for FileNotFoundError not in Python 2.7
3658        FileNotFoundError
3659    except NameError:
3660        FileNotFoundError = Exception
3661    try:
3662        fba = open(fbaName,'r')
3663    except FileNotFoundError:
3664        return False
3665    fba.readline()
3666    Nref = int(fba.readline()[:-1])
3667    fbalines = fba.readlines()
3668    for line in fbalines[:Nref]:
3669        info = line.split()
3670        h = int(info[0])
3671        k = int(info[1])
3672        l = int(info[2])
3673        FoR = float(info[3])
3674        FoI = float(info[4])
3675        Fosq = FoR**2+FoI**2
3676        phase = npatan2d(FoI,FoR)
3677        try:
3678            refId = reflDict[hash('%5d%5d%5d'%(h,k,l))]
3679        except KeyError:    #added reflections at end skipped
3680            d = float(1/np.sqrt(G2lat.calc_rDsq([h,k,l],A)))
3681            if 'T' in Type:
3682                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])
3683            else:
3684                newRefs.append([h,k,l,-1,d,0.,0.01,1.0,Fosq,Fosq,phase,1.0,1.0,1.0,1.0])
3685            continue
3686        newRefs[refId][8] = Fosq
3687        newRefs[refId][10] = phase
3688    newRefs = np.array(newRefs)
3689    return True,newRefs
3690   
3691#### testing data
3692NeedTestData = True
3693def TestData():
3694    'needs a doc string'
3695#    global NeedTestData
3696    global bakType
3697    bakType = 'chebyschev'
3698    global xdata
3699    xdata = np.linspace(4.0,40.0,36000)
3700    global parmDict0
3701    parmDict0 = {
3702        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
3703        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
3704        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
3705        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
3706        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'Z':0.0,'SH/L':0.002,
3707        'Back0':5.384,'Back1':-0.015,'Back2':.004,
3708        }
3709    global parmDict1
3710    parmDict1 = {
3711        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
3712        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
3713        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
3714        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
3715        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
3716        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
3717        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'Z':0.0,'SH/L':0.002,
3718        'Back0':36.897,'Back1':-0.508,'Back2':.006,
3719        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
3720        }
3721    global parmDict2
3722    parmDict2 = {
3723        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
3724        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'Z':0.0,'SH/L':0.02,
3725        'Back0':5.,'Back1':-0.02,'Back2':.004,
3726#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
3727        }
3728    global varyList
3729    varyList = []
3730
3731def test0():
3732    if NeedTestData: TestData()
3733    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
3734    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
3735    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
3736    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
3737    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
3738    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
3739   
3740def test1():
3741    if NeedTestData: TestData()
3742    time0 = time.time()
3743    for i in range(100):
3744        getPeakProfile(parmDict1,xdata,varyList,bakType)
3745    G2fil.G2Print ('100+6*Ka1-2 peaks=1200 peaks %.2f'%time.time()-time0)
3746   
3747def test2(name,delt):
3748    if NeedTestData: TestData()
3749    varyList = [name,]
3750    xdata = np.linspace(5.6,5.8,400)
3751    hplot = plotter.add('derivatives test for '+name).gca()
3752    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
3753    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
3754    parmDict2[name] += delt
3755    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
3756    hplot.plot(xdata,(y1-y0)/delt,'r+')
3757   
3758def test3(name,delt):
3759    if NeedTestData: TestData()
3760    names = ['pos','sig','gam','shl']
3761    idx = names.index(name)
3762    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
3763    xdata = np.linspace(5.6,5.8,800)
3764    dx = xdata[1]-xdata[0]
3765    hplot = plotter.add('derivatives test for '+name).gca()
3766    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
3767    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
3768    myDict[name] += delt
3769    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
3770    hplot.plot(xdata,(y1-y0)/delt,'r+')
3771
3772if __name__ == '__main__':
3773    import GSASIItestplot as plot
3774    global plotter
3775    plotter = plot.PlotNotebook()
3776#    test0()
3777#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)']:
3778    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
3779        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['Z',0.0001],['SH/L',0.00005]]:
3780        test2(name,shft)
3781    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
3782        test3(name,shft)
3783    G2fil.G2Print ("OK")
3784    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.