source: trunk/GSASIIpwd.py @ 4800

Last change on this file since 4800 was 4800, checked in by toby, 9 months ago

finish up docs cleanups

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