source: trunk/GSASIIpwd.py @ 4400

Last change on this file since 4400 was 4400, checked in by vondreele, 3 years ago

add rigid body save to make new pdb file with properly named atoms - needed for fullrmc
implement use of pdbparser (part of fullrmc) to prepare amorphous big box models from isolated molecules.

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