source: trunk/GSASIIpwd.py @ 4821

Last change on this file since 4821 was 4821, checked in by toby, 8 months ago

new peakfit mode in scripting: allow peak width terms to remain unchanged (neither refined nor calc from inst. params.); not accesible from GUI.

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