source: trunk/GSASIIpwd.py @ 4830

Last change on this file since 4830 was 4830, checked in by vondreele, 8 months ago

show more places in Absorb output
improve image absorption calcs for cylinders now seems to work OK
reorganize Image controls GUI & add option for vertical samples
hide printing of LS cycles if Rrint=False in HessianLSQ
implement (& not use) Klein-Nishina PDF correction
Fix bug in if importer

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