source: trunk/GSASIIpwd.py @ 4915

Last change on this file since 4915 was 4915, checked in by vondreele, 5 months ago

add some doc strings for profile functions in G2pwd
rename getPowderProfileDervMP as getPowderProfileDerv since it is more generally used than just multiprocessing

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 172.0 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-05-26 12:43:09 +0000 (Wed, 26 May 2021) $
10# $Author: vondreele $
11# $Revision: 4915 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 4915 2021-05-26 12:43:09Z 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: 4915 $"
37GSASIIpath.SetVersionNumber("$Revision: 4915 $")
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#### Using old gsas fortran routines for powder peak shapes & derivatives
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 /= np.sum(Df)
1135    return Df
1136
1137def getdFCJVoigt3(pos,sig,gam,shl,xdata):
1138    '''Compute analytic derivatives the Finger-Cox-Jepcoat modified Pseudo-Voigt
1139    function for a CW powder peak
1140    '''
1141    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1142    return Df,dFdp,dFds,dFdg,dFdsh
1143
1144def getPsVoigt(pos,sig,gam,xdata):
1145    '''Compute the simple Pseudo-Voigt function for a
1146    small angle Bragg peak in external Fortran routine
1147    '''
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    '''Compute the simple Pseudo-Voigt function derivatives for a
1155    small angle Bragg peak peak in external Fortran routine
1156    '''
1157   
1158    Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
1159    return Df,dFdp,dFds,dFdg
1160
1161def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
1162    '''Compute the double exponential Pseudo-Voigt convolution function for a
1163    neutron TOF & CW pink peak in external Fortran routine
1164    '''
1165   
1166    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1167    Df /= np.sum(Df)
1168    return Df 
1169   
1170def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
1171    '''Compute the double exponential Pseudo-Voigt convolution function derivatives for a
1172    neutron TOF & CW pink peak in external Fortran routine
1173    '''
1174   
1175    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1176    return Df,dFdp,dFda,dFdb,dFds,dFdg   
1177
1178def ellipseSize(H,Sij,GB):
1179    '''Implements r=1/sqrt(sum((1/S)*(q.v)^2) per note from Alexander Brady
1180    '''
1181   
1182    HX = np.inner(H.T,GB)
1183    lenHX = np.sqrt(np.sum(HX**2))
1184    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
1185    R = np.inner(HX/lenHX,Rsize)**2*Esize         #want column length for hkl in crystal
1186    lenR = 1./np.sqrt(np.sum(R))
1187    return lenR
1188
1189def ellipseSizeDerv(H,Sij,GB):
1190    '''Implements r=1/sqrt(sum((1/S)*(q.v)^2) derivative per note from Alexander Brady
1191    '''
1192   
1193    lenR = ellipseSize(H,Sij,GB)
1194    delt = 0.001
1195    dRdS = np.zeros(6)
1196    for i in range(6):
1197        Sij[i] -= delt
1198        lenM = ellipseSize(H,Sij,GB)
1199        Sij[i] += 2.*delt
1200        lenP = ellipseSize(H,Sij,GB)
1201        Sij[i] -= delt
1202        dRdS[i] = (lenP-lenM)/(2.*delt)
1203    return lenR,dRdS
1204
1205def getMustrain(HKL,G,SGData,muStrData):
1206    if muStrData[0] == 'isotropic':
1207        return np.ones(HKL.shape[1])*muStrData[1][0]
1208    elif muStrData[0] == 'uniaxial':
1209        H = np.array(HKL)
1210        P = np.array(muStrData[3])
1211        cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1212        Si = muStrData[1][0]
1213        Sa = muStrData[1][1]
1214        return Si*Sa/(np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1215    else:       #generalized - P.W. Stephens model
1216        H = np.array(HKL)
1217        rdsq = np.array([G2lat.calc_rDsq2(h,G) for h in H.T])
1218        Strms = np.array(G2spc.MustrainCoeff(H,SGData))
1219        Sum = np.sum(np.array(muStrData[4])[:,nxs]*Strms,axis=0)
1220        return np.sqrt(Sum)/rdsq
1221   
1222def getCrSize(HKL,G,GB,sizeData):
1223    if sizeData[0] == 'isotropic':
1224        return np.ones(HKL.shape[1])*sizeData[1][0]
1225    elif sizeData[0] == 'uniaxial':
1226        H = np.array(HKL)
1227        P = np.array(sizeData[3])
1228        cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1229        Si = sizeData[1][0]
1230        Sa = sizeData[1][1]
1231        return Si*Sa/(np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1232    else:
1233        Sij =[sizeData[4][i] for i in range(6)]
1234        H = np.array(HKL)
1235        return 1./np.array([ellipseSize(h,Sij,GB) for h in H.T])**2
1236
1237def getHKLpeak(dmin,SGData,A,Inst=None,nodup=False):
1238    '''
1239    Generates allowed by symmetry reflections with d >= dmin
1240    NB: GenHKLf & checkMagextc return True for extinct reflections
1241
1242    :param dmin:  minimum d-spacing
1243    :param SGData: space group data obtained from SpcGroup
1244    :param A: lattice parameter terms A1-A6
1245    :param Inst: instrument parameter info
1246    :returns: HKLs: np.array hkl, etc for allowed reflections
1247
1248    '''
1249    HKL = G2lat.GenHLaue(dmin,SGData,A)       
1250    HKLs = []
1251    ds = []
1252    for h,k,l,d in HKL:
1253        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1254        if ext and 'MagSpGrp' in SGData:
1255            ext = G2spc.checkMagextc([h,k,l],SGData)
1256        if not ext:
1257            if nodup and int(10000*d) in ds:
1258                continue
1259            ds.append(int(10000*d))
1260            if Inst == None:
1261                HKLs.append([h,k,l,d,0,-1])
1262            else:
1263                HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1264    return np.array(HKLs)
1265
1266def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1267    'needs a doc string'
1268    HKLs = []
1269    vec = np.array(Vec)
1270    vstar = np.sqrt(G2lat.calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1271    dvec = 1./(maxH*vstar+1./dmin)
1272    HKL = G2lat.GenHLaue(dvec,SGData,A)       
1273    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1274    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1275    ifMag = False
1276    if 'MagSpGrp' in SGData:
1277        ifMag = True
1278    for h,k,l,d in HKL:
1279        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1280        if not ext and d >= dmin:
1281            HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1282        for dH in SSdH:
1283            if dH:
1284                DH = SSdH[dH]
1285                H = [h+DH[0],k+DH[1],l+DH[2]]
1286                d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1287                if d >= dmin:
1288                    HKLM = np.array([h,k,l,dH])
1289                    if G2spc.checkSSextc(HKLM,SSGData) or ifMag:
1290                        HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])   
1291    return G2lat.sortHKLd(HKLs,True,True,True)
1292
1293peakInstPrmMode = True
1294'''Determines the mode used for peak fitting. When peakInstPrmMode=True peak
1295width parameters are computed from the instrument parameters (UVW,... or
1296alpha,... etc) unless the individual parameter is refined. This allows the
1297instrument parameters to be refined. When peakInstPrmMode=False, the instrument
1298parameters are not used and cannot be refined.
1299The default is peakFitMode=True.
1300'''
1301
1302def setPeakInstPrmMode(normal=True):
1303    '''Determines the mode used for peak fitting. If normal=True (default)
1304    peak width parameters are computed from the instrument parameters
1305    unless the individual parameter is refined. If normal=False,
1306    peak widths are used as supplied for each peak.
1307
1308    Note that normal=True unless this routine is called. Also,
1309    instrument parameters can only be refined with normal=True.
1310
1311    :param bool normal: setting to apply to global variable
1312      :data:`peakInstPrmMode`
1313    '''
1314    global peakInstPrmMode
1315    peakInstPrmMode = normal
1316
1317def getPeakProfile(dataType,parmDict,xdata,fixback,varyList,bakType):
1318    'Computes the profile for a powder pattern'
1319   
1320    yb = getBackground('',parmDict,bakType,dataType,xdata,fixback)[0]
1321    yc = np.zeros_like(yb)
1322    cw = np.diff(xdata)
1323    cw = np.append(cw,cw[-1])
1324    if 'C' in dataType:
1325        shl = max(parmDict['SH/L'],0.002)
1326        Ka2 = False
1327        if 'Lam1' in parmDict.keys():
1328            Ka2 = True
1329            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1330            kRatio = parmDict['I(L2)/I(L1)']
1331        iPeak = 0
1332        while True:
1333            try:
1334                pos = parmDict['pos'+str(iPeak)]
1335                tth = (pos-parmDict['Zero'])
1336                intens = parmDict['int'+str(iPeak)]
1337                sigName = 'sig'+str(iPeak)
1338                if sigName in varyList or not peakInstPrmMode:
1339                    sig = parmDict[sigName]
1340                else:
1341                    sig = G2mth.getCWsig(parmDict,tth)
1342                sig = max(sig,0.001)          #avoid neg sigma^2
1343                gamName = 'gam'+str(iPeak)
1344                if gamName in varyList or not peakInstPrmMode:
1345                    gam = parmDict[gamName]
1346                else:
1347                    gam = G2mth.getCWgam(parmDict,tth)
1348                gam = max(gam,0.001)             #avoid neg gamma
1349                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1350                iBeg = np.searchsorted(xdata,pos-fmin)
1351                iFin = np.searchsorted(xdata,pos+fmin)
1352                if not iBeg+iFin:       #peak below low limit
1353                    iPeak += 1
1354                    continue
1355                elif not iBeg-iFin:     #peak above high limit
1356                    return yb+yc
1357                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1358                if Ka2:
1359                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1360                    iBeg = np.searchsorted(xdata,pos2-fmin)
1361                    iFin = np.searchsorted(xdata,pos2+fmin)
1362                    if iBeg-iFin:
1363                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1364                iPeak += 1
1365            except KeyError:        #no more peaks to process
1366                return yb+yc
1367    elif 'B' in dataType:
1368        iPeak = 0
1369        dsp = 1.0 #for now - fix later
1370        while True:
1371            try:
1372                pos = parmDict['pos'+str(iPeak)]
1373                tth = (pos-parmDict['Zero'])
1374                intens = parmDict['int'+str(iPeak)]
1375                alpName = 'alp'+str(iPeak)
1376                if alpName in varyList or not peakInstPrmMode:
1377                    alp = parmDict[alpName]
1378                else:
1379                    alp = G2mth.getPinkalpha(parmDict,tth)
1380                alp = max(0.1,alp)
1381                betName = 'bet'+str(iPeak)
1382                if betName in varyList or not peakInstPrmMode:
1383                    bet = parmDict[betName]
1384                else:
1385                    bet = G2mth.getPinkbeta(parmDict,tth)
1386                bet = max(0.1,bet)
1387                sigName = 'sig'+str(iPeak)
1388                if sigName in varyList or not peakInstPrmMode:
1389                    sig = parmDict[sigName]
1390                else:
1391                    sig = G2mth.getCWsig(parmDict,tth)
1392                sig = max(sig,0.001)          #avoid neg sigma^2
1393                gamName = 'gam'+str(iPeak)
1394                if gamName in varyList or not peakInstPrmMode:
1395                    gam = parmDict[gamName]
1396                else:
1397                    gam = G2mth.getCWgam(parmDict,tth)
1398                gam = max(gam,0.001)             #avoid neg gamma
1399                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1400                iBeg = np.searchsorted(xdata,pos-fmin)
1401                iFin = np.searchsorted(xdata,pos+fmin)
1402                if not iBeg+iFin:       #peak below low limit
1403                    iPeak += 1
1404                    continue
1405                elif not iBeg-iFin:     #peak above high limit
1406                    return yb+yc
1407                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])
1408                iPeak += 1
1409            except KeyError:        #no more peaks to process
1410                return yb+yc       
1411    else:
1412        Pdabc = parmDict['Pdabc']
1413        difC = parmDict['difC']
1414        iPeak = 0
1415        while True:
1416            try:
1417                pos = parmDict['pos'+str(iPeak)]               
1418                tof = pos-parmDict['Zero']
1419                dsp = tof/difC
1420                intens = parmDict['int'+str(iPeak)]
1421                alpName = 'alp'+str(iPeak)
1422                if alpName in varyList or not peakInstPrmMode:
1423                    alp = parmDict[alpName]
1424                else:
1425                    if len(Pdabc):
1426                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1427                    else:
1428                        alp = G2mth.getTOFalpha(parmDict,dsp)
1429                alp = max(0.1,alp)
1430                betName = 'bet'+str(iPeak)
1431                if betName in varyList or not peakInstPrmMode:
1432                    bet = parmDict[betName]
1433                else:
1434                    if len(Pdabc):
1435                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1436                    else:
1437                        bet = G2mth.getTOFbeta(parmDict,dsp)
1438                bet = max(0.0001,bet)
1439                sigName = 'sig'+str(iPeak)
1440                if sigName in varyList or not peakInstPrmMode:
1441                    sig = parmDict[sigName]
1442                else:
1443                    sig = G2mth.getTOFsig(parmDict,dsp)
1444                gamName = 'gam'+str(iPeak)
1445                if gamName in varyList or not peakInstPrmMode:
1446                    gam = parmDict[gamName]
1447                else:
1448                    gam = G2mth.getTOFgamma(parmDict,dsp)
1449                gam = max(gam,0.001)             #avoid neg gamma
1450                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1451                iBeg = np.searchsorted(xdata,pos-fmin)
1452                iFin = np.searchsorted(xdata,pos+fmax)
1453                lenX = len(xdata)
1454                if not iBeg:
1455                    iFin = np.searchsorted(xdata,pos+fmax)
1456                elif iBeg == lenX:
1457                    iFin = iBeg
1458                else:
1459                    iFin = np.searchsorted(xdata,pos+fmax)
1460                if not iBeg+iFin:       #peak below low limit
1461                    iPeak += 1
1462                    continue
1463                elif not iBeg-iFin:     #peak above high limit
1464                    return yb+yc
1465                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1466                iPeak += 1
1467            except KeyError:        #no more peaks to process
1468                return yb+yc
1469           
1470def getPeakProfileDerv(dataType,parmDict,xdata,fixback,varyList,bakType):
1471    'needs a doc string'
1472# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1473    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1474    dMdb,dMddb,dMdpk,dMdfb = getBackgroundDerv('',parmDict,bakType,dataType,xdata,fixback)
1475    if 'Back;0' in varyList:            #background derivs are in front if present
1476        dMdv[0:len(dMdb)] = dMdb
1477    names = ['DebyeA','DebyeR','DebyeU']
1478    for name in varyList:
1479        if 'Debye' in name:
1480            parm,Id = name.split(';')
1481            ip = names.index(parm)
1482            dMdv[varyList.index(name)] = dMddb[3*int(Id)+ip]
1483    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1484    for name in varyList:
1485        if 'BkPk' in name:
1486            parm,Id = name.split(';')
1487            ip = names.index(parm)
1488            dMdv[varyList.index(name)] = dMdpk[4*int(Id)+ip]
1489    cw = np.diff(xdata)
1490    cw = np.append(cw,cw[-1])
1491    if 'C' in dataType:
1492        shl = max(parmDict['SH/L'],0.002)
1493        Ka2 = False
1494        if 'Lam1' in parmDict.keys():
1495            Ka2 = True
1496            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1497            kRatio = parmDict['I(L2)/I(L1)']
1498        iPeak = 0
1499        while True:
1500            try:
1501                pos = parmDict['pos'+str(iPeak)]
1502                tth = (pos-parmDict['Zero'])
1503                intens = parmDict['int'+str(iPeak)]
1504                sigName = 'sig'+str(iPeak)
1505                if sigName in varyList or not peakInstPrmMode:
1506                    sig = parmDict[sigName]
1507                    dsdU = dsdV = dsdW = 0
1508                else:
1509                    sig = G2mth.getCWsig(parmDict,tth)
1510                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1511                sig = max(sig,0.001)          #avoid neg sigma
1512                gamName = 'gam'+str(iPeak)
1513                if gamName in varyList or not peakInstPrmMode:
1514                    gam = parmDict[gamName]
1515                    dgdX = dgdY = dgdZ = 0
1516                else:
1517                    gam = G2mth.getCWgam(parmDict,tth)
1518                    dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1519                gam = max(gam,0.001)             #avoid neg gamma
1520                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1521                iBeg = np.searchsorted(xdata,pos-fmin)
1522                iFin = np.searchsorted(xdata,pos+fmin)
1523                if not iBeg+iFin:       #peak below low limit
1524                    iPeak += 1
1525                    continue
1526                elif not iBeg-iFin:     #peak above high limit
1527                    break
1528                dMdpk = np.zeros(shape=(6,len(xdata)))
1529                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1530                for i in range(1,5):
1531                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
1532                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
1533                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1534                if Ka2:
1535                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1536                    iBeg = np.searchsorted(xdata,pos2-fmin)
1537                    iFin = np.searchsorted(xdata,pos2+fmin)
1538                    if iBeg-iFin:
1539                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1540                        for i in range(1,5):
1541                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
1542                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
1543                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
1544                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1545                for parmName in ['pos','int','sig','gam']:
1546                    try:
1547                        idx = varyList.index(parmName+str(iPeak))
1548                        dMdv[idx] = dervDict[parmName]
1549                    except ValueError:
1550                        pass
1551                if 'U' in varyList:
1552                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1553                if 'V' in varyList:
1554                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1555                if 'W' in varyList:
1556                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1557                if 'X' in varyList:
1558                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1559                if 'Y' in varyList:
1560                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1561                if 'Z' in varyList:
1562                    dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1563                if 'SH/L' in varyList:
1564                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
1565                if 'I(L2)/I(L1)' in varyList:
1566                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1567                iPeak += 1
1568            except KeyError:        #no more peaks to process
1569                break
1570    elif 'B' in dataType:
1571        iPeak = 0
1572        while True:
1573            try:
1574                pos = parmDict['pos'+str(iPeak)]
1575                tth = (pos-parmDict['Zero'])
1576                intens = parmDict['int'+str(iPeak)]
1577                alpName = 'alp'+str(iPeak)
1578                if alpName in varyList or not peakInstPrmMode:
1579                    alp = parmDict[alpName]
1580                    dada0 = dada1 = 0.0
1581                else:
1582                    alp = G2mth.getPinkalpha(parmDict,tth)
1583                    dada0,dada1 = G2mth.getPinkalphaDeriv(tth)
1584                alp = max(0.0001,alp)
1585                betName = 'bet'+str(iPeak)
1586                if betName in varyList or not peakInstPrmMode:
1587                    bet = parmDict[betName]
1588                    dbdb0 = dbdb1 = 0.0
1589                else:
1590                    bet = G2mth.getPinkbeta(parmDict,tth)
1591                    dbdb0,dbdb1 = G2mth.getPinkbetaDeriv(tth)
1592                bet = max(0.0001,bet)
1593                sigName = 'sig'+str(iPeak)
1594                if sigName in varyList or not peakInstPrmMode:
1595                    sig = parmDict[sigName]
1596                    dsdU = dsdV = dsdW = 0
1597                else:
1598                    sig = G2mth.getCWsig(parmDict,tth)
1599                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1600                sig = max(sig,0.001)          #avoid neg sigma
1601                gamName = 'gam'+str(iPeak)
1602                if gamName in varyList or not peakInstPrmMode:
1603                    gam = parmDict[gamName]
1604                    dgdX = dgdY = dgdZ = 0
1605                else:
1606                    gam = G2mth.getCWgam(parmDict,tth)
1607                    dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1608                gam = max(gam,0.001)             #avoid neg gamma
1609                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig/1.e4,gam/100.)
1610                iBeg = np.searchsorted(xdata,pos-fmin)
1611                iFin = np.searchsorted(xdata,pos+fmin)
1612                if not iBeg+iFin:       #peak below low limit
1613                    iPeak += 1
1614                    continue
1615                elif not iBeg-iFin:     #peak above high limit
1616                    break
1617                dMdpk = np.zeros(shape=(7,len(xdata)))
1618                dMdipk = getdEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])
1619                for i in range(1,6):
1620                    dMdpk[i][iBeg:iFin] += cw[iBeg:iFin]*intens*dMdipk[i]
1621                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1622                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4]/1.e4,'gam':dMdpk[5]/100.}
1623                for parmName in ['pos','int','alp','bet','sig','gam']:
1624                    try:
1625                        idx = varyList.index(parmName+str(iPeak))
1626                        dMdv[idx] = dervDict[parmName]
1627                    except ValueError:
1628                        pass
1629                if 'U' in varyList:
1630                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1631                if 'V' in varyList:
1632                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1633                if 'W' in varyList:
1634                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1635                if 'X' in varyList:
1636                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1637                if 'Y' in varyList:
1638                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1639                if 'Z' in varyList:
1640                    dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1641                if 'alpha-0' in varyList:
1642                    dMdv[varyList.index('alpha-0')] += dada0*dervDict['alp']
1643                if 'alpha-1' in varyList:
1644                    dMdv[varyList.index('alpha-1')] += dada1*dervDict['alp']
1645                if 'beta-0' in varyList:
1646                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1647                if 'beta-1' in varyList:
1648                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1649                iPeak += 1
1650            except KeyError:        #no more peaks to process
1651                break       
1652    else:
1653        Pdabc = parmDict['Pdabc']
1654        difC = parmDict['difC']
1655        iPeak = 0
1656        while True:
1657            try:
1658                pos = parmDict['pos'+str(iPeak)]               
1659                tof = pos-parmDict['Zero']
1660                dsp = tof/difC
1661                intens = parmDict['int'+str(iPeak)]
1662                alpName = 'alp'+str(iPeak)
1663                if alpName in varyList or not peakInstPrmMode:
1664                    alp = parmDict[alpName]
1665                else:
1666                    if len(Pdabc):
1667                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1668                        dada0 = 0
1669                    else:
1670                        alp = G2mth.getTOFalpha(parmDict,dsp)
1671                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1672                betName = 'bet'+str(iPeak)
1673                if betName in varyList or not peakInstPrmMode:
1674                    bet = parmDict[betName]
1675                else:
1676                    if len(Pdabc):
1677                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1678                        dbdb0 = dbdb1 = dbdb2 = 0
1679                    else:
1680                        bet = G2mth.getTOFbeta(parmDict,dsp)
1681                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1682                sigName = 'sig'+str(iPeak)
1683                if sigName in varyList or not peakInstPrmMode:
1684                    sig = parmDict[sigName]
1685                    dsds0 = dsds1 = dsds2 = dsds3 = 0
1686                else:
1687                    sig = G2mth.getTOFsig(parmDict,dsp)
1688                    dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1689                gamName = 'gam'+str(iPeak)
1690                if gamName in varyList or not peakInstPrmMode:
1691                    gam = parmDict[gamName]
1692                    dsdX = dsdY = dsdZ = 0
1693                else:
1694                    gam = G2mth.getTOFgamma(parmDict,dsp)
1695                    dsdX,dsdY,dsdZ = G2mth.getTOFgammaDeriv(dsp)
1696                gam = max(gam,0.001)             #avoid neg gamma
1697                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1698                iBeg = np.searchsorted(xdata,pos-fmin)
1699                lenX = len(xdata)
1700                if not iBeg:
1701                    iFin = np.searchsorted(xdata,pos+fmax)
1702                elif iBeg == lenX:
1703                    iFin = iBeg
1704                else:
1705                    iFin = np.searchsorted(xdata,pos+fmax)
1706                if not iBeg+iFin:       #peak below low limit
1707                    iPeak += 1
1708                    continue
1709                elif not iBeg-iFin:     #peak above high limit
1710                    break
1711                dMdpk = np.zeros(shape=(7,len(xdata)))
1712                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1713                for i in range(1,6):
1714                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1715                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1716                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1717                for parmName in ['pos','int','alp','bet','sig','gam']:
1718                    try:
1719                        idx = varyList.index(parmName+str(iPeak))
1720                        dMdv[idx] = dervDict[parmName]
1721                    except ValueError:
1722                        pass
1723                if 'alpha' in varyList:
1724                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1725                if 'beta-0' in varyList:
1726                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1727                if 'beta-1' in varyList:
1728                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1729                if 'beta-q' in varyList:
1730                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1731                if 'sig-0' in varyList:
1732                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1733                if 'sig-1' in varyList:
1734                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1735                if 'sig-2' in varyList:
1736                    dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1737                if 'sig-q' in varyList:
1738                    dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1739                if 'X' in varyList:
1740                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1741                if 'Y' in varyList:
1742                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']
1743                if 'Z' in varyList:
1744                    dMdv[varyList.index('Z')] += dsdZ*dervDict['gam']
1745                iPeak += 1
1746            except KeyError:        #no more peaks to process
1747                break
1748    if 'BF mult' in varyList:
1749        dMdv[varyList.index('BF mult')] = fixback
1750       
1751    return dMdv
1752       
1753def Dict2Values(parmdict, varylist):
1754    '''Use before call to leastsq to setup list of values for the parameters
1755    in parmdict, as selected by key in varylist'''
1756    return [parmdict[key] for key in varylist] 
1757   
1758def Values2Dict(parmdict, varylist, values):
1759    ''' Use after call to leastsq to update the parameter dictionary with
1760    values corresponding to keys in varylist'''
1761    parmdict.update(zip(varylist,values))
1762   
1763def SetBackgroundParms(Background):
1764    'Loads background parameters into dicts/lists to create varylist & parmdict'
1765    if len(Background) == 1:            # fix up old backgrounds
1766        Background.append({'nDebye':0,'debyeTerms':[]})
1767    bakType,bakFlag = Background[0][:2]
1768    backVals = Background[0][3:]
1769    backNames = ['Back;'+str(i) for i in range(len(backVals))]
1770    Debye = Background[1]           #also has background peaks stuff
1771    backDict = dict(zip(backNames,backVals))
1772    backVary = []
1773    if bakFlag:
1774        backVary = backNames
1775
1776    backDict['nDebye'] = Debye['nDebye']
1777    debyeDict = {}
1778    debyeList = []
1779    for i in range(Debye['nDebye']):
1780        debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1781        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1782        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1783    debyeVary = []
1784    for item in debyeList:
1785        if item[1]:
1786            debyeVary.append(item[0])
1787    backDict.update(debyeDict)
1788    backVary += debyeVary
1789
1790    backDict['nPeaks'] = Debye['nPeaks']
1791    peaksDict = {}
1792    peaksList = []
1793    for i in range(Debye['nPeaks']):
1794        peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1795        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1796        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1797    peaksVary = []
1798    for item in peaksList:
1799        if item[1]:
1800            peaksVary.append(item[0])
1801    backDict.update(peaksDict)
1802    backVary += peaksVary
1803    if 'background PWDR' in Background[1]:
1804        backDict['Back File'] = Background[1]['background PWDR'][0]
1805        backDict['BF mult'] = Background[1]['background PWDR'][1]
1806        if len(Background[1]['background PWDR']) > 2:
1807            if Background[1]['background PWDR'][2]:
1808                backVary += ['BF mult',]
1809    return bakType,backDict,backVary
1810   
1811def DoCalibInst(IndexPeaks,Inst):
1812   
1813    def SetInstParms():
1814        dataType = Inst['Type'][0]
1815        insVary = []
1816        insNames = []
1817        insVals = []
1818        for parm in Inst:
1819            insNames.append(parm)
1820            insVals.append(Inst[parm][1])
1821            if parm in ['Lam','difC','difA','difB','Zero',]:
1822                if Inst[parm][2]:
1823                    insVary.append(parm)
1824        instDict = dict(zip(insNames,insVals))
1825        return dataType,instDict,insVary
1826       
1827    def GetInstParms(parmDict,Inst,varyList):
1828        for name in Inst:
1829            Inst[name][1] = parmDict[name]
1830       
1831    def InstPrint(Inst,sigDict):
1832        print ('Instrument Parameters:')
1833        if 'C' in Inst['Type'][0] or 'B' in Inst['Type'][0]:
1834            ptfmt = "%12.6f"
1835        else:
1836            ptfmt = "%12.3f"
1837        ptlbls = 'names :'
1838        ptstr =  'values:'
1839        sigstr = 'esds  :'
1840        for parm in Inst:
1841            if parm in  ['Lam','difC','difA','difB','Zero',]:
1842                ptlbls += "%s" % (parm.center(12))
1843                ptstr += ptfmt % (Inst[parm][1])
1844                if parm in sigDict:
1845                    sigstr += ptfmt % (sigDict[parm])
1846                else:
1847                    sigstr += 12*' '
1848        print (ptlbls)
1849        print (ptstr)
1850        print (sigstr)
1851       
1852    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1853        parmDict.update(zip(varyList,values))
1854        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1855
1856    peakPos = []
1857    peakDsp = []
1858    peakWt = []
1859    for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1860        if peak[2] and peak[3] and sig > 0.:
1861            peakPos.append(peak[0])
1862            peakDsp.append(peak[-1])    #d-calc
1863#            peakWt.append(peak[-1]**2/sig**2)   #weight by d**2
1864            peakWt.append(1./(sig*peak[-1]))   #
1865    peakPos = np.array(peakPos)
1866    peakDsp = np.array(peakDsp)
1867    peakWt = np.array(peakWt)
1868    dataType,insDict,insVary = SetInstParms()
1869    parmDict = {}
1870    parmDict.update(insDict)
1871    varyList = insVary
1872    if not len(varyList):
1873        G2fil.G2Print ('**** ERROR - nothing to refine! ****')
1874        return False
1875    while True:
1876        begin = time.time()
1877        values =  np.array(Dict2Values(parmDict, varyList))
1878        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1879            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1880        ncyc = int(result[2]['nfev']/2)
1881        runtime = time.time()-begin   
1882        chisq = np.sum(result[2]['fvec']**2)
1883        Values2Dict(parmDict, varyList, result[0])
1884        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
1885        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(peakPos),len(varyList)))
1886        G2fil.G2Print ('calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
1887        G2fil.G2Print ('chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF))
1888        try:
1889            sig = np.sqrt(np.diag(result[1])*GOF)
1890            if np.any(np.isnan(sig)):
1891                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***')
1892            break                   #refinement succeeded - finish up!
1893        except ValueError:          #result[1] is None on singular matrix
1894            G2fil.G2Print ('**** Refinement failed - singular matrix ****')
1895       
1896    sigDict = dict(zip(varyList,sig))
1897    GetInstParms(parmDict,Inst,varyList)
1898    InstPrint(Inst,sigDict)
1899    return True
1900           
1901def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],oneCycle=False,controls=None,wtFactor=1.0,dlg=None):
1902    '''Called to perform a peak fit, refining the selected items in the peak
1903    table as well as selected items in the background.
1904
1905    :param str FitPgm: type of fit to perform. At present this is ignored.
1906    :param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1907      four values followed by a refine flag where the values are: position, intensity,
1908      sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1909      tree entry, dict item "peaks"
1910    :param list Background: describes the background. List with two items.
1911      Item 0 specifies a background model and coefficients. Item 1 is a dict.
1912      From the Histogram/Background tree entry.
1913    :param list Limits: min and max x-value to use
1914    :param dict Inst: Instrument parameters
1915    :param dict Inst2: more Instrument parameters
1916    :param numpy.array data: a 5xn array. data[0] is the x-values,
1917      data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1918      calc, background and difference intensities, respectively.
1919    :param array fixback: fixed background array; same size as data[0-5]
1920    :param list prevVaryList: Used in sequential refinements to override the
1921      variable list. Defaults as an empty list.
1922    :param bool oneCycle: True if only one cycle of fitting should be performed
1923    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1924      and derivType = controls['deriv type']. If None default values are used.
1925    :param float wtFactor: weight multiplier; = 1.0 by default
1926    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1927      Defaults to None, which means no updates are done.
1928    '''
1929    def GetBackgroundParms(parmList,Background):
1930        iBak = 0
1931        while True:
1932            try:
1933                bakName = 'Back;'+str(iBak)
1934                Background[0][iBak+3] = parmList[bakName]
1935                iBak += 1
1936            except KeyError:
1937                break
1938        iDb = 0
1939        while True:
1940            names = ['DebyeA;','DebyeR;','DebyeU;']
1941            try:
1942                for i,name in enumerate(names):
1943                    val = parmList[name+str(iDb)]
1944                    Background[1]['debyeTerms'][iDb][2*i] = val
1945                iDb += 1
1946            except KeyError:
1947                break
1948        iDb = 0
1949        while True:
1950            names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1951            try:
1952                for i,name in enumerate(names):
1953                    val = parmList[name+str(iDb)]
1954                    Background[1]['peaksList'][iDb][2*i] = val
1955                iDb += 1
1956            except KeyError:
1957                break
1958        if 'BF mult' in parmList:
1959            Background[1]['background PWDR'][1] = parmList['BF mult']
1960               
1961    def BackgroundPrint(Background,sigDict):
1962        print ('Background coefficients for '+Background[0][0]+' function')
1963        ptfmt = "%12.5f"
1964        ptstr =  'value: '
1965        sigstr = 'esd  : '
1966        for i,back in enumerate(Background[0][3:]):
1967            ptstr += ptfmt % (back)
1968            if Background[0][1]:
1969                prm = 'Back;'+str(i)
1970                if prm in sigDict:
1971                    sigstr += ptfmt % (sigDict[prm])
1972                else:
1973                    sigstr += " "*12
1974            if len(ptstr) > 75:
1975                print (ptstr)
1976                if Background[0][1]: print (sigstr)
1977                ptstr =  'value: '
1978                sigstr = 'esd  : '
1979        if len(ptstr) > 8:
1980            print (ptstr)
1981            if Background[0][1]: print (sigstr)
1982
1983        if Background[1]['nDebye']:
1984            parms = ['DebyeA;','DebyeR;','DebyeU;']
1985            print ('Debye diffuse scattering coefficients')
1986            ptfmt = "%12.5f"
1987            print (' term       DebyeA       esd        DebyeR       esd        DebyeU        esd')
1988            for term in range(Background[1]['nDebye']):
1989                line = ' term %d'%(term)
1990                for ip,name in enumerate(parms):
1991                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1992                    if name+str(term) in sigDict:
1993                        line += ptfmt%(sigDict[name+str(term)])
1994                    else:
1995                        line += " "*12
1996                print (line)
1997        if Background[1]['nPeaks']:
1998            print ('Coefficients for Background Peaks')
1999            ptfmt = "%15.3f"
2000            for j,pl in enumerate(Background[1]['peaksList']):
2001                names =  'peak %3d:'%(j+1)
2002                ptstr =  'values  :'
2003                sigstr = 'esds    :'
2004                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
2005                    val = pl[2*i]
2006                    prm = lbl+";"+str(j)
2007                    names += '%15s'%(prm)
2008                    ptstr += ptfmt%(val)
2009                    if prm in sigDict:
2010                        sigstr += ptfmt%(sigDict[prm])
2011                    else:
2012                        sigstr += " "*15
2013                print (names)
2014                print (ptstr)
2015                print (sigstr)
2016        if 'BF mult' in sigDict:
2017            print('Background file mult: %.3f(%d)'%(Background[1]['background PWDR'][1],int(1000*sigDict['BF mult'])))
2018                           
2019    def SetInstParms(Inst):
2020        dataType = Inst['Type'][0]
2021        insVary = []
2022        insNames = []
2023        insVals = []
2024        for parm in Inst:
2025            insNames.append(parm)
2026            insVals.append(Inst[parm][1])
2027            if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
2028                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1'] and Inst[parm][2]:
2029                    insVary.append(parm)
2030        instDict = dict(zip(insNames,insVals))
2031#        instDict['X'] = max(instDict['X'],0.01)
2032#        instDict['Y'] = max(instDict['Y'],0.01)
2033        if 'SH/L' in instDict:
2034            instDict['SH/L'] = max(instDict['SH/L'],0.002)
2035        return dataType,instDict,insVary
2036       
2037    def GetInstParms(parmDict,Inst,varyList,Peaks):
2038        for name in Inst:
2039            Inst[name][1] = parmDict[name]
2040        iPeak = 0
2041        while True:
2042            try:
2043                sigName = 'sig'+str(iPeak)
2044                pos = parmDict['pos'+str(iPeak)]
2045                if sigName not in varyList and peakInstPrmMode:
2046                    if 'T' in Inst['Type'][0]:
2047                        dsp = G2lat.Pos2dsp(Inst,pos)
2048                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
2049                    else:
2050                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
2051                gamName = 'gam'+str(iPeak)
2052                if gamName not in varyList and peakInstPrmMode:
2053                    if 'T' in Inst['Type'][0]:
2054                        dsp = G2lat.Pos2dsp(Inst,pos)
2055                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
2056                    else:
2057                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
2058                iPeak += 1
2059            except KeyError:
2060                break
2061       
2062    def InstPrint(Inst,sigDict):
2063        print ('Instrument Parameters:')
2064        ptfmt = "%12.6f"
2065        ptlbls = 'names :'
2066        ptstr =  'values:'
2067        sigstr = 'esds  :'
2068        for parm in Inst:
2069            if parm in  ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
2070                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1']:
2071                ptlbls += "%s" % (parm.center(12))
2072                ptstr += ptfmt % (Inst[parm][1])
2073                if parm in sigDict:
2074                    sigstr += ptfmt % (sigDict[parm])
2075                else:
2076                    sigstr += 12*' '
2077        print (ptlbls)
2078        print (ptstr)
2079        print (sigstr)
2080
2081    def SetPeaksParms(dataType,Peaks):
2082        peakNames = []
2083        peakVary = []
2084        peakVals = []
2085        if 'C' in dataType:
2086            names = ['pos','int','sig','gam']
2087        else:   #'T' and 'B'
2088            names = ['pos','int','alp','bet','sig','gam']
2089        for i,peak in enumerate(Peaks):
2090            for j,name in enumerate(names):
2091                peakVals.append(peak[2*j])
2092                parName = name+str(i)
2093                peakNames.append(parName)
2094                if peak[2*j+1]:
2095                    peakVary.append(parName)
2096        return dict(zip(peakNames,peakVals)),peakVary
2097               
2098    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
2099        if 'C' in Inst['Type'][0]:
2100            names = ['pos','int','sig','gam']
2101        else:   #'T' & 'B'
2102            names = ['pos','int','alp','bet','sig','gam']
2103        for i,peak in enumerate(Peaks):
2104            pos = parmDict['pos'+str(i)]
2105            if 'difC' in Inst:
2106                dsp = pos/Inst['difC'][1]
2107            for j in range(len(names)):
2108                parName = names[j]+str(i)
2109                if parName in varyList or not peakInstPrmMode:
2110                    peak[2*j] = parmDict[parName]
2111                elif 'alp' in parName:
2112                    if 'T' in Inst['Type'][0]:
2113                        peak[2*j] = G2mth.getTOFalpha(parmDict,dsp)
2114                    else: #'B'
2115                        peak[2*j] = G2mth.getPinkalpha(parmDict,pos)
2116                elif 'bet' in parName:
2117                    if 'T' in Inst['Type'][0]:
2118                        peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
2119                    else:   #'B'
2120                        peak[2*j] = G2mth.getPinkbeta(parmDict,pos)
2121                elif 'sig' in parName:
2122                    if 'T' in Inst['Type'][0]:
2123                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
2124                    else:   #'C' & 'B'
2125                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
2126                elif 'gam' in parName:
2127                    if 'T' in Inst['Type'][0]:
2128                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
2129                    else:   #'C' & 'B'
2130                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
2131                       
2132    def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
2133        print ('Peak coefficients:')
2134        if 'C' in dataType:
2135            names = ['pos','int','sig','gam']
2136        else:   #'T' & 'B'
2137            names = ['pos','int','alp','bet','sig','gam']           
2138        head = 13*' '
2139        for name in names:
2140            if name in ['alp','bet']:
2141                head += name.center(8)+'esd'.center(8)
2142            else:
2143                head += name.center(10)+'esd'.center(10)
2144        head += 'bins'.center(8)
2145        print (head)
2146        if 'C' in dataType:
2147            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
2148        elif 'T' in dataType:
2149            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
2150        else: #'B'
2151            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'alp':"%8.2f",'bet':"%8.4f",'sig':"%10.3f",'gam':"%10.3f"}
2152        for i,peak in enumerate(Peaks):
2153            ptstr =  ':'
2154            for j in range(len(names)):
2155                name = names[j]
2156                parName = name+str(i)
2157                ptstr += ptfmt[name] % (parmDict[parName])
2158                if parName in varyList:
2159                    ptstr += ptfmt[name] % (sigDict[parName])
2160                else:
2161                    if name in ['alp','bet']:
2162                        ptstr += 8*' '
2163                    else:
2164                        ptstr += 10*' '
2165            ptstr += '%9.2f'%(ptsperFW[i])
2166            print ('%s'%(('Peak'+str(i+1)).center(8)),ptstr)
2167               
2168    def devPeakProfile(values,xdata,ydata,fixback, weights,dataType,parmdict,varylist,bakType,dlg):
2169        parmdict.update(zip(varylist,values))
2170        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,fixback,varylist,bakType)
2171           
2172    def errPeakProfile(values,xdata,ydata,fixback,weights,dataType,parmdict,varylist,bakType,dlg):       
2173        parmdict.update(zip(varylist,values))
2174        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,fixback,varylist,bakType)-ydata)
2175        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
2176        if dlg:
2177            dlg.Raise()
2178            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
2179            if not GoOn:
2180                return -M           #abort!!
2181        return M
2182
2183    # beginning of DoPeakFit
2184    if controls:
2185        Ftol = controls['min dM/M']
2186    else:
2187        Ftol = 0.0001
2188    if oneCycle:
2189        Ftol = 1.0
2190    x,y,w,yc,yb,yd = data   #these are numpy arrays - remove masks!
2191    if fixback is None:
2192        fixback = np.zeros_like(y)
2193    yc *= 0.                            #set calcd ones to zero
2194    yb *= 0.
2195    yd *= 0.
2196    cw = x[1:]-x[:-1]
2197    xBeg = np.searchsorted(x,Limits[0])
2198    xFin = np.searchsorted(x,Limits[1])+1
2199    bakType,bakDict,bakVary = SetBackgroundParms(Background)
2200    dataType,insDict,insVary = SetInstParms(Inst)
2201    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
2202    parmDict = {}
2203    parmDict.update(bakDict)
2204    parmDict.update(insDict)
2205    parmDict.update(peakDict)
2206    parmDict['Pdabc'] = []      #dummy Pdabc
2207    parmDict.update(Inst2)      #put in real one if there
2208    if prevVaryList:
2209        varyList = prevVaryList[:]
2210    else:
2211        varyList = bakVary+insVary+peakVary
2212    fullvaryList = varyList[:]
2213    if not peakInstPrmMode:
2214        for v in ('U','V','W','X','Y','Z','alpha','alpha-0','alpha-1',
2215            'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',):
2216            if v in varyList:
2217                raise Exception('Instrumental profile terms cannot be varied '+
2218                                    'after setPeakInstPrmMode(False) is used')
2219    while True:
2220        begin = time.time()
2221        values =  np.array(Dict2Values(parmDict, varyList))
2222        Rvals = {}
2223        badVary = []
2224        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
2225               args=(x[xBeg:xFin],y[xBeg:xFin],fixback[xBeg:xFin],wtFactor*w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
2226        ncyc = int(result[2]['nfev']/2)
2227        runtime = time.time()-begin   
2228        chisq = np.sum(result[2]['fvec']**2)
2229        Values2Dict(parmDict, varyList, result[0])
2230        Rvals['Rwp'] = np.sqrt(chisq/np.sum(wtFactor*w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
2231        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
2232        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList)))
2233        if ncyc:
2234            G2fil.G2Print ('fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
2235        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2236        sig = [0]*len(varyList)
2237        if len(varyList) == 0: break  # if nothing was refined
2238        try:
2239            sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
2240            if np.any(np.isnan(sig)):
2241                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***')
2242            break                   #refinement succeeded - finish up!
2243        except ValueError:          #result[1] is None on singular matrix
2244            G2fil.G2Print ('**** Refinement failed - singular matrix ****')
2245            Ipvt = result[2]['ipvt']
2246            for i,ipvt in enumerate(Ipvt):
2247                if not np.sum(result[2]['fjac'],axis=1)[i]:
2248                    G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
2249                    badVary.append(varyList[ipvt-1])
2250                    del(varyList[ipvt-1])
2251                    break
2252            else: # nothing removed
2253                break
2254    if dlg: dlg.Destroy()
2255    sigDict = dict(zip(varyList,sig))
2256    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin],fixback[xBeg:xFin])[0]
2257    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],fixback[xBeg:xFin],varyList,bakType)
2258    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
2259    GetBackgroundParms(parmDict,Background)
2260    if bakVary: BackgroundPrint(Background,sigDict)
2261    GetInstParms(parmDict,Inst,varyList,Peaks)
2262    if insVary: InstPrint(Inst,sigDict)
2263    GetPeaksParms(Inst,parmDict,Peaks,varyList)
2264    binsperFWHM = []
2265    for peak in Peaks:
2266        FWHM = getFWHM(peak[0],Inst)
2267        try:
2268            binsperFWHM.append(FWHM/cw[x.searchsorted(peak[0])])
2269        except IndexError:
2270            binsperFWHM.append(0.)
2271    if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
2272    if len(binsperFWHM):
2273        if min(binsperFWHM) < 1.:
2274            G2fil.G2Print ('*** Warning: calculated peak widths are too narrow to refine profile coefficients ***')
2275            if 'T' in Inst['Type'][0]:
2276                G2fil.G2Print (' Manually increase sig-0, 1, or 2 in Instrument Parameters')
2277            else:
2278                G2fil.G2Print (' Manually increase W in Instrument Parameters')
2279        elif min(binsperFWHM) < 4.:
2280            G2fil.G2Print ('*** Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement ***')
2281            G2fil.G2Print ('*** recommended is 6-10; you have %.2f ***'%(min(binsperFWHM)))
2282    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
2283   
2284def calcIncident(Iparm,xdata):
2285    'needs a doc string'
2286
2287    def IfunAdv(Iparm,xdata):
2288        Itype = Iparm['Itype']
2289        Icoef = Iparm['Icoeff']
2290        DYI = np.ones((12,xdata.shape[0]))
2291        YI = np.ones_like(xdata)*Icoef[0]
2292       
2293        x = xdata/1000.                 #expressions are in ms
2294        if Itype == 'Exponential':
2295            for i in [1,3,5,7,9]:
2296                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
2297                YI += Icoef[i]*Eterm
2298                DYI[i] *= Eterm
2299                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
2300        elif 'Maxwell'in Itype:
2301            Eterm = np.exp(-Icoef[2]/x**2)
2302            DYI[1] = Eterm/x**5
2303            DYI[2] = -Icoef[1]*DYI[1]/x**2
2304            YI += (Icoef[1]*Eterm/x**5)
2305            if 'Exponential' in Itype:
2306                for i in range(3,11,2):
2307                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
2308                    YI += Icoef[i]*Eterm
2309                    DYI[i] *= Eterm
2310                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
2311            else:   #Chebyschev
2312                T = (2./x)-1.
2313                Ccof = np.ones((12,xdata.shape[0]))
2314                Ccof[1] = T
2315                for i in range(2,12):
2316                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
2317                for i in range(1,10):
2318                    YI += Ccof[i]*Icoef[i+2]
2319                    DYI[i+2] =Ccof[i]
2320        return YI,DYI
2321       
2322    Iesd = np.array(Iparm['Iesd'])
2323    Icovar = Iparm['Icovar']
2324    YI,DYI = IfunAdv(Iparm,xdata)
2325    YI = np.where(YI>0,YI,1.)
2326    WYI = np.zeros_like(xdata)
2327    vcov = np.zeros((12,12))
2328    k = 0
2329    for i in range(12):
2330        for j in range(i,12):
2331            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
2332            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
2333            k += 1
2334    M = np.inner(vcov,DYI.T)
2335    WYI = np.sum(M*DYI,axis=0)
2336    WYI = np.where(WYI>0.,WYI,0.)
2337    return YI,WYI
2338
2339#### RMCutilities ################################################################################
2340def MakeInst(PWDdata,Name,Size,Mustrain,useSamBrd):
2341    inst = PWDdata['Instrument Parameters'][0]
2342    Xsb = 0.
2343    Ysb = 0.
2344    if 'T' in inst['Type'][1]:
2345        difC = inst['difC'][1]
2346        if useSamBrd[0]:
2347            if 'ellipsoidal' not in Size[0]:    #take the isotropic term only
2348                Xsb = 1.e-4*difC/Size[1][0]
2349        if useSamBrd[1]:
2350            if 'generalized' not in Mustrain[0]:    #take the isotropic term only
2351                Ysb = 1.e-6*difC*Mustrain[1][0]
2352        prms = ['Bank',
2353                'difC','difA','Zero','2-theta',
2354                'alpha','beta-0','beta-1',
2355                'sig-0','sig-1','sig-2',
2356                'Z','X','Y']
2357        fname = Name+'.inst'
2358        fl = open(fname,'w')
2359        fl.write('1\n')
2360        fl.write('%d\n'%int(inst[prms[0]][1]))
2361        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]))
2362        fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[5]][1],inst[prms[6]][1],inst[prms[7]][1]))
2363        fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[8]][1],inst[prms[9]][1],inst[prms[10]][1]))   
2364        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))
2365        fl.close()
2366    else:
2367        if useSamBrd[0]:
2368            wave = G2mth.getWave(inst)
2369            if 'ellipsoidal' not in Size[0]:    #take the isotropic term only
2370                Xsb = 1.8*wave/(np.pi*Size[1][0])
2371        if useSamBrd[1]:
2372            if 'generalized' not in Mustrain[0]:    #take the isotropic term only
2373                Ysb = 0.0180*Mustrain[1][0]/np.pi
2374        prms = ['Bank',
2375                'Lam','Zero','Polariz.',
2376                'U','V','W',
2377                'X','Y']
2378        fname = Name+'.inst'
2379        fl = open(fname,'w')
2380        fl.write('1\n')
2381        fl.write('%d\n'%int(inst[prms[0]][1]))
2382        fl.write('%10.5f%10.5f%10.4f%10d\n'%(inst[prms[1]][1],-100.*inst[prms[2]][1],inst[prms[3]][1],0))
2383        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[4]][1],inst[prms[5]][1],inst[prms[6]][1]))
2384        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[7]][1]+Xsb,inst[prms[8]][1]+Ysb,0.0))   
2385        fl.write('%10.3f%10.3f%10.3f\n'%(0.0,0.0,0.0))
2386        fl.close()
2387    return fname
2388   
2389def MakeBack(PWDdata,Name):
2390    Back = PWDdata['Background'][0]
2391    inst = PWDdata['Instrument Parameters'][0]
2392    if 'chebyschev-1' != Back[0]:
2393        return None
2394    Nback = Back[2]
2395    BackVals = Back[3:]
2396    fname = Name+'.back'
2397    fl = open(fname,'w')
2398    fl.write('%10d\n'%Nback)
2399    for val in BackVals:
2400        if 'T' in inst['Type'][1]:
2401            fl.write('%12.6g\n'%(float(val)))
2402        else:
2403            fl.write('%12.6g\n'%val)
2404    fl.close()
2405    return fname
2406
2407def findDup(Atoms):
2408    Dup = []
2409    Fracs = []
2410    for iat1,at1 in enumerate(Atoms):
2411        if any([at1[0] in dup for dup in Dup]):
2412            continue
2413        else:
2414            Dup.append([at1[0],])
2415            Fracs.append([at1[6],])
2416        for iat2,at2 in enumerate(Atoms[(iat1+1):]):
2417            if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
2418                Dup[-1] += [at2[0],]
2419                Fracs[-1]+= [at2[6],]
2420    return Dup,Fracs
2421
2422def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):   
2423   
2424    Meta = RMCPdict['metadata']
2425    Atseq = RMCPdict['atSeq']
2426    Supercell =  RMCPdict['SuperCell']
2427    generalData = Phase['General']
2428    Dups,Fracs = findDup(Phase['Atoms'])
2429    Sfracs = [np.cumsum(fracs) for fracs in Fracs]
2430    ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
2431    Sample = PWDdata['Sample Parameters']
2432    Meta['temperature'] = Sample['Temperature']
2433    Meta['pressure'] = Sample['Pressure']
2434    Cell = generalData['Cell'][1:7]
2435    Trans = np.eye(3)*np.array(Supercell)
2436    newPhase = copy.deepcopy(Phase)
2437    newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2438    newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans)
2439    GB = G2lat.cell2Gmat( newPhase['General']['Cell'][1:7])[0]
2440    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.
2441    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
2442    Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2443    Natm = np.count_nonzero(Natm-1)
2444    Atoms = newPhase['Atoms']
2445    reset = False
2446   
2447    if ifSfracs:
2448        Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2449        Natm = np.count_nonzero(Natm-1)
2450        Satoms = []
2451        for i in range(len(Atoms)//Natm):
2452            ind = i*Natm
2453            Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
2454        Natoms = []
2455        for satoms in Satoms:
2456            for idup,dup in enumerate(Dups):
2457                ldup = len(dup)
2458                natm = len(satoms)
2459                i = 0
2460                while i < natm:
2461                    if satoms[i][0] in dup:
2462                        atoms = satoms[i:i+ldup]
2463                        try:
2464                            atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2465                            Natoms.append(atom)
2466                        except IndexError:      #what about vacancies?
2467                            if 'Va' not in Atseq:
2468                                reset = True
2469                                Atseq.append('Va')
2470                                RMCPdict['aTypes']['Va'] = 0.0
2471                            atom = atoms[0]
2472                            atom[1] = 'Va'
2473                            Natoms.append(atom)
2474                        i += ldup
2475                    else:
2476                       i += 1
2477    else:
2478        Natoms = Atoms
2479   
2480    NAtype = np.zeros(len(Atseq))
2481    for atom in Natoms:
2482        NAtype[Atseq.index(atom[1])] += 1
2483    NAstr = ['%6d'%i for i in NAtype]
2484    Cell = newPhase['General']['Cell'][1:7]
2485    if os.path.exists(Name+'.his6f'):
2486        os.remove(Name+'.his6f')
2487    if os.path.exists(Name+'.neigh'):
2488        os.remove(Name+'.neigh')
2489    fname = Name+'.rmc6f'
2490    fl = open(fname,'w')
2491    fl.write('(Version 6f format configuration file)\n')
2492    for item in Meta:
2493        fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2494    fl.write('Atom types present:                 %s\n'%'    '.join(Atseq))
2495    fl.write('Number of each atom type:       %s\n'%''.join(NAstr))
2496    fl.write('Number of atoms:                %d\n'%len(Natoms))
2497    fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2498    fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2499            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2500    A,B = G2lat.cell2AB(Cell,True)
2501    fl.write('Lattice vectors (Ang):\n')   
2502    for i in [0,1,2]:
2503        fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2504    fl.write('Atoms (fractional coordinates):\n')
2505    nat = 0
2506    for atm in Atseq:
2507        for iat,atom in enumerate(Natoms):
2508            if atom[1] == atm:
2509                nat += 1
2510                atcode = Atcodes[iat].split(':')
2511                cell = [0,0,0]
2512                if '+' in atcode[1]:
2513                    cell = eval(atcode[1].split('+')[1])
2514                fl.write('%6d%4s  [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(       
2515                        nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2516    fl.close()
2517    return fname,reset
2518
2519def MakeBragg(PWDdata,Name,Phase):
2520    generalData = Phase['General']
2521    Vol = generalData['Cell'][7]
2522    Data = PWDdata['Data']
2523    Inst = PWDdata['Instrument Parameters'][0]
2524    Bank = int(Inst['Bank'][1])
2525    Sample = PWDdata['Sample Parameters']
2526    Scale = Sample['Scale'][0]
2527    if 'X' in Inst['Type'][1]:
2528        Scale *= 2.
2529    Limits = PWDdata['Limits'][1]
2530    Ibeg = np.searchsorted(Data[0],Limits[0])
2531    Ifin = np.searchsorted(Data[0],Limits[1])+1
2532    fname = Name+'.bragg'
2533    fl = open(fname,'w')
2534    fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol))
2535    if 'T' in Inst['Type'][0]:
2536        fl.write('%12s%12s\n'%('   TOF,ms','  I(obs)'))
2537        for i in range(Ibeg,Ifin-1):
2538            fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2539    else:
2540        fl.write('%12s%12s\n'%('   2-theta, deg','  I(obs)'))
2541        for i in range(Ibeg,Ifin-1):
2542            fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))       
2543    fl.close()
2544    return fname
2545
2546def MakeRMCPdat(PWDdata,Name,Phase,RMCPdict):
2547    Meta = RMCPdict['metadata']
2548    Times = RMCPdict['runTimes']
2549    Atseq = RMCPdict['atSeq']
2550    Atypes = RMCPdict['aTypes']
2551    atPairs = RMCPdict['Pairs']
2552    Files = RMCPdict['files']
2553    BraggWt = RMCPdict['histogram'][1]
2554    inst = PWDdata['Instrument Parameters'][0]
2555    try:
2556        refList = PWDdata['Reflection Lists'][Name]['RefList']
2557    except KeyError:
2558        return 'Error - missing reflection list; you must do Refine first'
2559    dMin = refList[-1][4]
2560    gsasType = 'xray2'
2561    if 'T' in inst['Type'][1]:
2562        gsasType = 'gsas3'
2563    elif 'X' in inst['Type'][1]:
2564        XFF = G2elem.GetFFtable(Atseq)
2565        Xfl = open(Name+'.xray','w')
2566        for atm in Atseq:
2567            fa = XFF[atm]['fa']
2568            fb = XFF[atm]['fb']
2569            fc = XFF[atm]['fc']
2570            Xfl.write('%2s  %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%(
2571                    atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc))
2572        Xfl.close()
2573    lenA = len(Atseq)
2574    Pairs = []
2575    for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]:
2576        Pairs += pair
2577    pairMin = [atPairs[pair]for pair in Pairs if pair in atPairs]
2578    maxMoves = [Atypes[atm] for atm in Atseq if atm in Atypes]
2579    fname = Name+'.dat'
2580    fl = open(fname,'w')
2581    fl.write(' %% Hand edit the following as needed\n')
2582    fl.write('TITLE :: '+Name+'\n')
2583    fl.write('MATERIAL :: '+Meta['material']+'\n')
2584    fl.write('PHASE :: '+Meta['phase']+'\n')
2585    fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n')
2586    fl.write('INVESTIGATOR :: '+Meta['owner']+'\n')
2587    minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin])
2588    minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin])
2589    maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin])
2590    fl.write('MINIMUM_DISTANCES ::   %s  Angstrom\n'%minHD)
2591    maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves])
2592    fl.write('MAXIMUM_MOVES ::   %s Angstrom\n'%maxMv)
2593    fl.write('R_SPACING ::  0.0200 Angstrom\n')
2594    fl.write('PRINT_PERIOD :: 100\n')
2595    fl.write('TIME_LIMIT ::     %.2f MINUTES\n'%Times[0])
2596    fl.write('SAVE_PERIOD ::    %.2f MINUTES\n'%Times[1])
2597    fl.write('\n')
2598    fl.write('ATOMS :: '+' '.join(Atseq)+'\n')
2599    fl.write('\n')
2600    fl.write('FLAGS ::\n')
2601    fl.write('  > NO_MOVEOUT\n')
2602    fl.write('  > NO_SAVE_CONFIGURATIONS\n')
2603    fl.write('  > NO_RESOLUTION_CONVOLUTION\n')
2604    fl.write('\n')
2605    fl.write('INPUT_CONFIGURATION_FORMAT ::  rmc6f\n')
2606    fl.write('SAVE_CONFIGURATION_FORMAT  ::  rmc6f\n')
2607    fl.write('IGNORE_HISTORY_FILE ::\n')
2608    fl.write('\n')
2609    fl.write('DISTANCE_WINDOW ::\n')
2610    fl.write('  > MNDIST :: %s\n'%minD)
2611    fl.write('  > MXDIST :: %s\n'%maxD)
2612    if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']):
2613        fl.write('\n')
2614        fl.write('POTENTIALS ::\n')
2615        fl.write('  > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.'])
2616        fl.write('  > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n')
2617        if len(RMCPdict['Potentials']['Stretch']):
2618            fl.write('  > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search'])
2619            for bond in RMCPdict['Potentials']['Stretch']:
2620                fl.write('  > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2]))       
2621        if len(RMCPdict['Potentials']['Angles']):
2622            fl.write('  > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
2623            for angle in RMCPdict['Potentials']['Angles']:
2624                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
2625                    (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
2626    if RMCPdict['useBVS']:
2627        fl.write('BVS ::\n')
2628        fl.write('  > ATOM :: '+' '.join(Atseq)+'\n')
2629        fl.write('  > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2630        oxid = []
2631        for val in RMCPdict['Oxid']:
2632            if len(val) == 3:
2633                oxid.append(val[0][1:])
2634            else:
2635                oxid.append(val[0][2:])
2636        fl.write('  > OXID :: %s\n'%' '.join(oxid))
2637        fl.write('  > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']]))
2638        fl.write('  > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']]))
2639        fl.write('  > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))       
2640        fl.write('  > SAVE :: 100000\n')
2641        fl.write('  > UPDATE :: 100000\n')
2642        if len(RMCPdict['Swap']):
2643            fl.write('\n')
2644            fl.write('SWAP_MULTI ::\n')
2645            for swap in RMCPdict['Swap']:
2646                try:
2647                    at1 = Atseq.index(swap[0])
2648                    at2 = Atseq.index(swap[1])
2649                except ValueError:
2650                    break
2651                fl.write('  > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2]))
2652       
2653    if len(RMCPdict['FxCN']):
2654        fl.write('FIXED_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['FxCN']))       
2655        for ifx,fxcn in enumerate(RMCPdict['FxCN']):
2656            try:
2657                at1 = Atseq.index(fxcn[0])
2658                at2 = Atseq.index(fxcn[1])
2659            except ValueError:
2660                break
2661            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]))
2662    if len(RMCPdict['AveCN']):
2663        fl.write('AVERAGE_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['AveCN']))
2664        for iav,avcn in enumerate(RMCPdict['AveCN']):
2665            try:
2666                at1 = Atseq.index(avcn[0])
2667                at2 = Atseq.index(avcn[1])
2668            except ValueError:
2669                break
2670            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]))
2671    for File in Files:
2672        if Files[File][0] and Files[File][0] != 'Select':
2673            if 'Xray' in File and 'F(Q)' in File:
2674                fqdata = open(Files[File][0],'r')
2675                lines = int(fqdata.readline()[:-1])
2676            fl.write('\n')
2677            fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_'))
2678            fl.write('  > FILENAME :: %s\n'%Files[File][0])
2679            fl.write('  > DATA_TYPE :: %s\n'%Files[File][2])
2680            fl.write('  > FIT_TYPE :: %s\n'%Files[File][2])
2681            if 'Xray' not in File:
2682                fl.write('  > START_POINT :: 1\n')
2683                fl.write('  > END_POINT :: 3000\n')
2684                fl.write('  > WEIGHT :: %.4f\n'%Files[File][1])
2685            fl.write('  > CONSTANT_OFFSET 0.000\n')
2686            fl.write('  > NO_FITTED_OFFSET\n')
2687            if RMCPdict['FitScale']:
2688                fl.write('  > FITTED_SCALE\n')
2689            else:
2690                fl.write('  > NO_FITTED_SCALE\n')
2691            if Files[File][3] !='RMC':
2692                fl.write('  > %s\n'%Files[File][3])
2693            if 'reciprocal' in File:
2694                fl.write('  > CONVOLVE ::\n')
2695                if 'Xray' in File:
2696                    fl.write('  > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines)
2697                    fl.write('  > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1]))
2698                    fl.write('  > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2))
2699                    fl.write('  > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,1./Files[File][1]))
2700    fl.write('\n')
2701    fl.write('BRAGG ::\n')
2702    fl.write('  > BRAGG_SHAPE :: %s\n'%gsasType)
2703    fl.write('  > RECALCUATE\n')
2704    fl.write('  > DMIN :: %.2f\n'%(dMin-0.02))
2705    fl.write('  > WEIGHT :: %10.3f\n'%BraggWt)
2706    fl.write('\n')
2707    fl.write('END  ::\n')
2708    fl.close()
2709    return fname
2710
2711# def FindBonds(Phase,RMCPdict):
2712#     generalData = Phase['General']
2713#     cx,ct,cs,cia = generalData['AtomPtrs']
2714#     atomData = Phase['Atoms']
2715#     Res = 'RMC'
2716#     if 'macro' in generalData['Type']:
2717#         Res = atomData[0][ct-3]
2718#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2719#     Pairs = RMCPdict['Pairs']   #dict!
2720#     BondList = []
2721#     notNames = []
2722#     for FrstName in AtDict:
2723#         nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
2724#         Atyp1 = AtDict[FrstName]
2725#         if 'Va' in Atyp1:
2726#             continue
2727#         for nbr in nbrs:
2728#             Atyp2 = AtDict[nbr[0]]
2729#             if 'Va' in Atyp2:
2730#                 continue
2731#             try:
2732#                 bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
2733#             except KeyError:
2734#                 bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
2735#             if any(bndData):
2736#                 if bndData[0] <= nbr[1] <= bndData[1]:
2737#                     bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
2738#                     revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
2739#                     if bondStr not in BondList and revbondStr not in BondList:
2740#                         BondList.append(bondStr)
2741#         notNames.append(FrstName)
2742#     return Res,BondList
2743
2744# def FindAngles(Phase,RMCPdict):
2745#     generalData = Phase['General']
2746#     Cell = generalData['Cell'][1:7]
2747#     Amat = G2lat.cell2AB(Cell)[0]
2748#     cx,ct,cs,cia = generalData['AtomPtrs']
2749#     atomData = Phase['Atoms']
2750#     AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
2751#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2752#     Angles = RMCPdict['Angles']
2753#     AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
2754#     AngleList = []
2755#     for MidName in AtDict:
2756#         nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
2757#         if len(nbrs) < 2: #need 2 neighbors to make an angle
2758#             continue
2759#         Atyp2 = AtDict[MidName]
2760#         for i,nbr1 in enumerate(nbrs):
2761#             Atyp1 = AtDict[nbr1[0]]
2762#             for j,nbr3 in enumerate(nbrs[i+1:]):
2763#                 Atyp3 = AtDict[nbr3[0]]
2764#                 IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
2765#                 try:
2766#                     angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
2767#                 except KeyError:
2768#                     try:
2769#                         angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
2770#                     except KeyError:
2771#                         continue
2772#                 XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
2773#                 calAngle = G2mth.getRestAngle(XYZ,Amat)
2774#                 if angData[0] <= calAngle <= angData[1]:
2775#                     angStr = str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n'
2776#                     revangStr = str((MidName,nbr3[0],nbr1[0])+tuple(angData))+',\n'
2777#                     if angStr not in AngleList and revangStr not in AngleList:
2778#                         AngleList.append(angStr)
2779#     return AngleList
2780
2781# def GetSqConvolution(XY,d):
2782
2783#     n = XY.shape[1]
2784#     snew = np.zeros(n)
2785#     dq = np.zeros(n)
2786#     sold = XY[1]
2787#     q = XY[0]
2788#     dq[1:] = np.diff(q)
2789#     dq[0] = dq[1]
2790   
2791#     for j in range(n):
2792#         for i in range(n):
2793#             b = abs(q[i]-q[j])
2794#             t = q[i]+q[j]
2795#             if j == i:
2796#                 snew[j] += q[i]*sold[i]*(d-np.sin(t*d)/t)*dq[i]
2797#             else:
2798#                 snew[j] += q[i]*sold[i]*(np.sin(b*d)/b-np.sin(t*d)/t)*dq[i]
2799#         snew[j] /= np.pi*q[j]
2800   
2801#     snew[0] = snew[1]
2802#     return snew
2803
2804# def GetMaxSphere(pdbName):
2805#     try:
2806#         pFil = open(pdbName,'r')
2807#     except FileNotFoundError:
2808#         return None
2809#     while True:
2810#         line = pFil.readline()
2811#         if 'Boundary' in line:
2812#             line = line.split()[3:]
2813#             G = np.array([float(item) for item in line])
2814#             G = np.reshape(G,(3,3))**2
2815#             G = nl.inv(G)
2816#             pFil.close()
2817#             break
2818#     dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
2819#     return min(dspaces)
2820   
2821def MakefullrmcRun(pName,Phase,RMCPdict):
2822    BondList = {}
2823    for k in RMCPdict['Pairs']:
2824        if RMCPdict['Pairs'][k][1]+RMCPdict['Pairs'][k][2]>0:
2825            BondList[k] = (RMCPdict['Pairs'][k][1],RMCPdict['Pairs'][k][2])
2826    AngleList = []
2827    for angle in RMCPdict['Angles']:
2828        if angle[3] == angle[4] or angle[5] >= angle[6] or angle[6] <= 0:
2829            continue
2830        for i in (0,1,2):
2831            angle[i] = angle[i].strip()
2832        AngleList.append(angle)
2833    rmin = RMCPdict['min Contact']
2834    cell = Phase['General']['Cell'][1:7]
2835    SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
2836    cx,ct,cs,cia = Phase['General']['AtomPtrs']
2837    atomsList = []
2838    for atom in Phase['Atoms']:
2839        el = ''.join([i for i in atom[ct] if i.isalpha()])
2840        atomsList.append([el] + atom[cx:cx+4])
2841    rname = pName+'-fullrmc.py'
2842    restart = '%s_restart.pdb'%pName
2843    Files = RMCPdict['files']
2844    rundata = ''
2845    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%rname
2846    rundata += '# created in '+__file__+" v"+filversion.split()[1]
2847    rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
2848    rundata += '''
2849# fullrmc imports (all that are potentially useful)
2850import os,glob
2851import time
2852#import matplotlib.pyplot as plt
2853import numpy as np
2854from fullrmc.Core import Collection
2855from fullrmc.Engine import Engine
2856import fullrmc.Constraints.PairDistributionConstraints as fPDF
2857from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint
2858from fullrmc.Constraints.DistanceConstraints import DistanceConstraint
2859from fullrmc.Constraints.BondConstraints import BondConstraint
2860from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
2861from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
2862from fullrmc.Generators.Swaps import SwapPositionsGenerator
2863
2864### When True, erases an existing enging to provide a fresh start
2865FRESH_START = {}
2866time0 = time.time()
2867'''.format(RMCPdict['ReStart'][0])
2868   
2869    rundata += '# setup structure\n'
2870    rundata += 'cell = ' + str(cell) + '\n'
2871    rundata += "SymOpList = "+str([i.lower() for i in SymOpList]) + '\n'
2872    rundata += 'atomList = ' + str(atomsList).replace('],','],\n  ') + '\n'
2873    rundata += 'supercell = ' + str(RMCPdict['SuperCell']) + '\n'
2874
2875    rundata += '\n# initialize engine\n'
2876    rundata += 'engineFileName = "%s.rmc"\n'%pName
2877
2878    rundata += '''\n# check Engine exists if so (and not FRESH_START) load it
2879# otherwise build it
2880ENGINE = Engine(path=None)
2881if not ENGINE.is_engine(engineFileName) or FRESH_START:
2882    ## create structure
2883    ENGINE = Engine(path=engineFileName, freshStart=True)
2884    ENGINE.build_crystal_set_pdb(symOps     = SymOpList,
2885                                 atoms      = atomList,
2886                                 unitcellBC = cell,
2887                                 supercell  = supercell)
2888    rmax = min( [ENGINE.boundaryConditions.get_a(), ENGINE.boundaryConditions.get_b(), ENGINE.boundaryConditions.get_c()] ) /2.
2889'''   
2890    import atmdata
2891    rundata += '# conversion factors (may be needed)\n'
2892    rundata += '    sumCiBi2 = 0.\n'
2893    for elem in Phase['General']['AtomTypes']:
2894        rundata += '    Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
2895        rundata += '    sumCiBi2 += (Ci*{})**2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
2896    rundata += '    rho0 = len(ENGINE.allNames)/ENGINE.volume\n'
2897    # settings that require a new Engine
2898    for File in Files:
2899        filDat = RMCPdict['files'][File]
2900        if not os.path.exists(filDat[0]): continue
2901        sfwt = 'neutronCohb'
2902        if 'Xray' in File:
2903            sfwt = 'atomicNumber'
2904        if 'G(r)' in File:
2905            rundata += '    GR = np.loadtxt("%s").T\n'%filDat[0]
2906            if filDat[3] == 0:
2907                rundata += '''    # read and xform G(r) as defined in RMCProfile
2908    # see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
2909                rundata += '    GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n'
2910                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
2911            elif filDat[3] == 1:
2912                rundata += '    # This is G(r) as defined in PDFFIT\n'
2913                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
2914            elif filDat[3] == 2:
2915                rundata += '    # This is g(r)\n'
2916                rundata += '    GofR = fPDF.PairCorrelationConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
2917            else:
2918                raise ValueError('Invalid G(r) type: '+str(filDat[3]))
2919            rundata += '    ENGINE.add_constraints([GofR])\n'
2920            rundata += '    GofR.set_limits((None, rmax))\n'
2921        elif '(Q)' in File:
2922            rundata += '    SOQ = np.loadtxt("%s").T\n'%filDat[0]
2923            if filDat[3] == 0:
2924                rundata += '    # Read & xform F(Q) as defined in RMCProfile to S(Q)-1\n'
2925                rundata += '    SOQ[1] *= 1 / sumCiBi2\n'
2926            elif filDat[3] == 1:
2927                rundata += '    # This is S(Q) as defined in PDFFIT\n'
2928                rundata += '    SOQ[1] -= 1\n'
2929            if filDat[4]:
2930                rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=rmax)\n'
2931            rundata += '    SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
2932            rundata += '    ENGINE.add_constraints([SofQ])\n'
2933        else:
2934            print('What is this?')
2935    rundata += '    ENGINE.add_constraints(DistanceConstraint(defaultLowerDistance={}))\n'.format(RMCPdict['min Contact'])
2936    if BondList:
2937        rundata += '''    B_CONSTRAINT   = BondConstraint()
2938    ENGINE.add_constraints(B_CONSTRAINT)
2939    B_CONSTRAINT.create_supercell_bonds(bondsDefinition=[
2940'''
2941        for pair in BondList:
2942            e1,e2 = pair.split('-')
2943            rundata += '            ("element","{}","{}",{},{}),\n'.format(
2944                                        e1.strip(),e2.strip(),*BondList[pair])
2945        rundata += '             ])\n'
2946    if AngleList:
2947        rundata += '''    A_CONSTRAINT   = BondsAngleConstraint()
2948    ENGINE.add_constraints(A_CONSTRAINT)
2949    A_CONSTRAINT.create_supercell_angles(anglesDefinition=[
2950'''
2951        for item in AngleList:
2952            rundata += ('            '+
2953               '("element","{1}","{0}","{2}",{5},{6},{5},{6},{3},{4}),\n'.format(*item))
2954        rundata += '             ])\n'
2955    rundata += '    for f in glob.glob("{}_*.log"): os.remove(f)\n'.format(pName)
2956    rundata += '''
2957    ENGINE.save()
2958else:
2959    ENGINE = ENGINE.load(path=engineFileName)
2960'''
2961    rundata += 'ENGINE.set_log_file("{}")\n'.format(pName)
2962    if RMCPdict['Swaps']:
2963        rundata += '\n#set up for site swaps\n'
2964        rundata += 'aN = ENGINE.allNames\n'
2965        rundata += 'SwapGen = {}\n'
2966        for swap in RMCPdict['Swaps']:
2967            rundata += 'SwapA = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[0]
2968            rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1]
2969            rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
2970        rundata += '    for swaps in SwapGen:\n'
2971        rundata += '        AB = swaps.split("-")\n'
2972        rundata += '        ENGINE.set_groups_as_atoms()\n'
2973        rundata += '        for g in ENGINE.groups:\n'
2974        rundata += '            if aN[g.indexes[0]]==AB[0]:\n'
2975        rundata += '                g.set_move_generator(SwapGen[swaps][0])\n'
2976        rundata += '            elif aN[g.indexes[0]]==AB[1]:\n'
2977        rundata += '                g.set_move_generator(SwapGen[swaps][1])\n'
2978        rundata += '            sProb = SwapGen[swaps][2]\n'
2979    rundata += '\n# set weights -- do this now so values can be changed without a restart\n'
2980    rundata += 'wtDict = {}\n'
2981    for File in Files:
2982        filDat = RMCPdict['files'][File]
2983        if not os.path.exists(filDat[0]): continue
2984        if 'Xray' in File:
2985            sfwt = 'atomicNumber'
2986        else:
2987            sfwt = 'neutronCohb'
2988        if 'G(r)' in File:
2989            typ = 'Pair'
2990        elif '(Q)' in File:
2991            typ = 'Struct'
2992        rundata += 'wtDict["{}-{}"] = {}\n'.format(typ,sfwt,filDat[1])
2993    rundata += 'for c in ENGINE.constraints:  # loop over predefined constraints\n'
2994    rundata += '    if type(c) is fPDF.PairDistributionConstraint:\n'
2995    rundata += '        c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
2996    rundata += '        c.set_limits((None,rmax))\n'
2997    if RMCPdict['FitScale']:
2998        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
2999    rundata += '    elif type(c) is ReducedStructureFactorConstraint:\n'
3000    rundata += '        c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
3001    if RMCPdict['FitScale']:
3002        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3003    # torsions difficult to implement, must be internal to cell & named with
3004    # fullrmc atom names
3005    # if len(RMCPdict['Torsions']):         # Torsions currently commented out in GUI
3006    #     rundata += 'for c in ENGINE.constraints:  # look for Dihedral Angle Constraints\n'
3007    #     rundata += '    if type(c) is DihedralAngleConstraint:\n'
3008    #     rundata += '        c.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
3009    #     rundata += '        c.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
3010    #     for torsion in RMCPdict['Torsions']:
3011    #         rundata += '    %s\n'%str(tuple(torsion))
3012    #     rundata += '        ]})\n'           
3013    rundata += '\n# setup runs for fullrmc\n'
3014
3015    rundata += 'steps = 10000\n'
3016    rundata += 'for _ in range({}):\n'.format(RMCPdict['Cycles'])
3017    rundata += '    ENGINE.set_groups_as_atoms()\n'
3018    rundata += '    expected = ENGINE.generated+steps\n'
3019   
3020    rundata += '    ENGINE.run(restartPdb="%s",numberOfSteps=steps, saveFrequency=1000)\n'%restart
3021    rundata += '    if ENGINE.generated != expected: break # run was stopped\n'
3022    rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
3023    rfile = open(rname,'w')
3024    rfile.writelines(rundata)
3025    rfile.close()
3026    return rname
3027   
3028def GetRMCBonds(general,RMCPdict,Atoms,bondList):
3029    bondDist = []
3030    Cell = general['Cell'][1:7]
3031    Supercell =  RMCPdict['SuperCell']
3032    Trans = np.eye(3)*np.array(Supercell)
3033    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3034    Amat,Bmat = G2lat.cell2AB(Cell)
3035    indices = (-1,0,1)
3036    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3037    for bonds in bondList:
3038        Oxyz = np.array(Atoms[bonds[0]][1:])
3039        Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])       
3040        Dx = np.array([Txyz-Oxyz+unit for unit in Units])
3041        Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
3042        for dx in Dx.T:
3043            bondDist.append(np.min(dx))
3044    return np.array(bondDist)
3045   
3046def GetRMCAngles(general,RMCPdict,Atoms,angleList):
3047    bondAngles = []
3048    Cell = general['Cell'][1:7]
3049    Supercell =  RMCPdict['SuperCell']
3050    Trans = np.eye(3)*np.array(Supercell)
3051    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3052    Amat,Bmat = G2lat.cell2AB(Cell)
3053    indices = (-1,0,1)
3054    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3055    for angle in angleList:
3056        Oxyz = np.array(Atoms[angle[0]][1:])
3057        TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])       
3058        TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])       
3059        DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
3060        DAx = np.sqrt(np.sum(DAxV**2,axis=2))
3061        DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
3062        DBx = np.sqrt(np.sum(DBxV**2,axis=2))
3063        iDAx = np.argmin(DAx,axis=0)
3064        iDBx = np.argmin(DBx,axis=0)
3065        for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
3066            DAv = DAxV[iA,i]/DAx[iA,i]
3067            DBv = DBxV[iB,i]/DBx[iB,i]
3068            bondAngles.append(npacosd(np.sum(DAv*DBv)))
3069    return np.array(bondAngles)
3070   
3071#### Reflectometry calculations ################################################################################
3072def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
3073    G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
3074   
3075    class RandomDisplacementBounds(object):
3076        """random displacement with bounds"""
3077        def __init__(self, xmin, xmax, stepsize=0.5):
3078            self.xmin = xmin
3079            self.xmax = xmax
3080            self.stepsize = stepsize
3081   
3082        def __call__(self, x):
3083            """take a random step but ensure the new position is within the bounds"""
3084            while True:
3085                # this could be done in a much more clever way, but it will work for example purposes
3086                steps = self.xmax-self.xmin
3087                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
3088                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
3089                    break
3090            return xnew
3091   
3092    def GetModelParms():
3093        parmDict = {}
3094        varyList = []
3095        values = []
3096        bounds = []
3097        parmDict['dQ type'] = data['dQ type']
3098        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
3099        for parm in ['Scale','FltBack']:
3100            parmDict[parm] = data[parm][0]
3101            if data[parm][1]:
3102                varyList.append(parm)
3103                values.append(data[parm][0])
3104                bounds.append(Bounds[parm])
3105        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
3106        parmDict['nLayers'] = len(parmDict['Layer Seq'])
3107        for ilay,layer in enumerate(data['Layers']):
3108            name = layer['Name']
3109            cid = str(ilay)+';'
3110            parmDict[cid+'Name'] = name
3111            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3112                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
3113                if layer.get(parm,[0.,False])[1]:
3114                    varyList.append(cid+parm)
3115                    value = layer[parm][0]
3116                    values.append(value)
3117                    if value:
3118                        bound = [value*Bfac,value/Bfac]
3119                    else:
3120                        bound = [0.,10.]
3121                    bounds.append(bound)
3122            if name not in ['vacuum','unit scatter']:
3123                parmDict[cid+'rho'] = Substances[name]['Scatt density']
3124                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
3125        return parmDict,varyList,values,bounds
3126   
3127    def SetModelParms():
3128        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
3129        if 'Scale' in varyList:
3130            data['Scale'][0] = parmDict['Scale']
3131            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
3132        G2fil.G2Print (line)
3133        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
3134        if 'FltBack' in varyList:
3135            data['FltBack'][0] = parmDict['FltBack']
3136            line += ' esd: %15.3g'%(sigDict['FltBack'])
3137        G2fil.G2Print (line)
3138        for ilay,layer in enumerate(data['Layers']):
3139            name = layer['Name']
3140            G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
3141            cid = str(ilay)+';'
3142            line = ' '
3143            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
3144            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
3145            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3146                if parm in layer:
3147                    if parm == 'Rough':
3148                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
3149                    else:
3150                        layer[parm][0] = parmDict[cid+parm]
3151                    line += ' %s: %.3f'%(parm,layer[parm][0])
3152                    if cid+parm in varyList:
3153                        line += ' esd: %.3g'%(sigDict[cid+parm])
3154            G2fil.G2Print (line)
3155            G2fil.G2Print (line2)
3156   
3157    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3158        parmDict.update(zip(varyList,values))
3159        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3160        return M
3161   
3162    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3163        parmDict.update(zip(varyList,values))
3164        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3165        return np.sum(M**2)
3166   
3167    def getREFD(Q,Qsig,parmDict):
3168        Ic = np.ones_like(Q)*parmDict['FltBack']
3169        Scale = parmDict['Scale']
3170        Nlayers = parmDict['nLayers']
3171        Res = parmDict['Res']
3172        depth = np.zeros(Nlayers)
3173        rho = np.zeros(Nlayers)
3174        irho = np.zeros(Nlayers)
3175        sigma = np.zeros(Nlayers)
3176        for ilay,lay in enumerate(parmDict['Layer Seq']):
3177            cid = str(lay)+';'
3178            depth[ilay] = parmDict[cid+'Thick']
3179            sigma[ilay] = parmDict[cid+'Rough']
3180            if parmDict[cid+'Name'] == u'unit scatter':
3181                rho[ilay] = parmDict[cid+'DenMul']
3182                irho[ilay] = parmDict[cid+'iDenMul']
3183            elif 'vacuum' != parmDict[cid+'Name']:
3184                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
3185                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
3186            if cid+'Mag SLD' in parmDict:
3187                rho[ilay] += parmDict[cid+'Mag SLD']
3188        if parmDict['dQ type'] == 'None':
3189            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
3190        elif 'const' in parmDict['dQ type']:
3191            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
3192        else:       #dQ/Q in data
3193            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
3194        Ic += AB*Scale
3195        return Ic
3196       
3197    def estimateT0(takestep):
3198        Mmax = -1.e-10
3199        Mmin = 1.e10
3200        for i in range(100):
3201            x0 = takestep(values)
3202            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3203            Mmin = min(M,Mmin)
3204            MMax = max(M,Mmax)
3205        return 1.5*(MMax-Mmin)
3206
3207    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3208    if data.get('2% weight'):
3209        wt = 1./(0.02*Io)**2
3210    Qmin = Limits[1][0]
3211    Qmax = Limits[1][1]
3212    wtFactor = ProfDict['wtFactor']
3213    Bfac = data['Toler']
3214    Ibeg = np.searchsorted(Q,Qmin)
3215    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
3216    Ic[:] = 0
3217    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
3218              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
3219    parmDict,varyList,values,bounds = GetModelParms()
3220    Msg = 'Failed to converge'
3221    if varyList:
3222        if data['Minimizer'] == 'LMLS': 
3223            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
3224                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3225            parmDict.update(zip(varyList,result[0]))
3226            chisq = np.sum(result[2]['fvec']**2)
3227            ncalc = result[2]['nfev']
3228            covM = result[1]
3229            newVals = result[0]
3230        elif data['Minimizer'] == 'Basin Hopping':
3231            xyrng = np.array(bounds).T
3232            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
3233            T0 = estimateT0(take_step)
3234            G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
3235            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
3236                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
3237                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
3238            chisq = result.fun
3239            ncalc = result.nfev
3240            newVals = result.x
3241            covM = []
3242        elif data['Minimizer'] == 'MC/SA Anneal':
3243            xyrng = np.array(bounds).T
3244            result = G2mth.anneal(sumREFD, values, 
3245                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
3246                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
3247                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
3248                ranRange=0.20,autoRan=False,dlg=None)
3249            newVals = result[0]
3250            parmDict.update(zip(varyList,newVals))
3251            chisq = result[1]
3252            ncalc = result[3]
3253            covM = []
3254            G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
3255        elif data['Minimizer'] == 'L-BFGS-B':
3256            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
3257                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3258            parmDict.update(zip(varyList,result['x']))
3259            chisq = result.fun
3260            ncalc = result.nfev
3261            newVals = result.x
3262            covM = []
3263    else:   #nothing varied
3264        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3265        chisq = np.sum(M**2)
3266        ncalc = 0
3267        covM = []
3268        sig = []
3269        sigDict = {}
3270        result = []
3271    Rvals = {}
3272    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
3273    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
3274    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
3275    Ib[Ibeg:Ifin] = parmDict['FltBack']
3276    try:
3277        if not len(varyList):
3278            Msg += ' - nothing refined'
3279            raise ValueError
3280        Nans = np.isnan(newVals)
3281        if np.any(Nans):
3282            name = varyList[Nans.nonzero(True)[0]]
3283            Msg += ' Nan result for '+name+'!'
3284            raise ValueError
3285        Negs = np.less_equal(newVals,0.)
3286        if np.any(Negs):
3287            indx = Negs.nonzero()
3288            name = varyList[indx[0][0]]
3289            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
3290                Msg += ' negative coefficient for '+name+'!'
3291                raise ValueError
3292        if len(covM):
3293            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
3294            covMatrix = covM*Rvals['GOF']
3295        else:
3296            sig = np.zeros(len(varyList))
3297            covMatrix = []
3298        sigDict = dict(zip(varyList,sig))
3299        G2fil.G2Print (' Results of reflectometry data modelling fit:')
3300        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
3301        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
3302        SetModelParms()
3303        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
3304    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
3305        G2fil.G2Print (Msg)
3306        return False,0,0,0,0,0,0,Msg
3307       
3308def makeSLDprofile(data,Substances):
3309   
3310    sq2 = np.sqrt(2.)
3311    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3312    Nlayers = len(laySeq)
3313    laySeq = np.array(laySeq,dtype=int)
3314    interfaces = np.zeros(Nlayers)
3315    rho = np.zeros(Nlayers)
3316    sigma = np.zeros(Nlayers)
3317    name = data['Layers'][0]['Name']
3318    thick = 0.
3319    for ilay,lay in enumerate(laySeq):
3320        layer = data['Layers'][lay]
3321        name = layer['Name']
3322        if 'Thick' in layer:
3323            thick += layer['Thick'][0]
3324            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
3325        if 'Rough' in layer:
3326            sigma[ilay] = max(0.001,layer['Rough'][0])
3327        if name != 'vacuum':
3328            if name == 'unit scatter':
3329                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
3330            else:
3331                rrho = Substances[name]['Scatt density']
3332                irho = Substances[name]['XImag density']
3333                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
3334        if 'Mag SLD' in layer:
3335            rho[ilay] += layer['Mag SLD'][0]
3336    name = data['Layers'][-1]['Name']
3337    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
3338    xr = np.flipud(x)
3339    interfaces[-1] = x[-1]
3340    y = np.ones_like(x)*rho[0]
3341    iBeg = 0
3342    for ilayer in range(Nlayers-1):
3343        delt = rho[ilayer+1]-rho[ilayer]
3344        iPos = np.searchsorted(x,interfaces[ilayer])
3345        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
3346        iBeg = iPos
3347    return x,xr,y   
3348
3349def REFDModelFxn(Profile,Inst,Limits,Substances,data):
3350   
3351    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3352    Qmin = Limits[1][0]
3353    Qmax = Limits[1][1]
3354    iBeg = np.searchsorted(Q,Qmin)
3355    iFin = np.searchsorted(Q,Qmax)+1    #include last point
3356    Ib[:] = data['FltBack'][0]
3357    Ic[:] = 0
3358    Scale = data['Scale'][0]
3359    if data['Layer Seq'] == []:
3360        return
3361    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3362    Nlayers = len(laySeq)
3363    depth = np.zeros(Nlayers)
3364    rho = np.zeros(Nlayers)
3365    irho = np.zeros(Nlayers)
3366    sigma = np.zeros(Nlayers)
3367    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
3368        layer = data['Layers'][lay]
3369        name = layer['Name']
3370        if 'Thick' in layer:    #skips first & last layers
3371            depth[ilay] = layer['Thick'][0]
3372        if 'Rough' in layer:    #skips first layer
3373            sigma[ilay] = layer['Rough'][0]
3374        if 'unit scatter' == name:
3375            rho[ilay] = layer['DenMul'][0]
3376            irho[ilay] = layer['iDenMul'][0]
3377        else:
3378            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
3379            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
3380        if 'Mag SLD' in layer:
3381            rho[ilay] += layer['Mag SLD'][0]
3382    if data['dQ type'] == 'None':
3383        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
3384    elif 'const' in data['dQ type']:
3385        res = data['Resolution'][0]/(100.*sateln2)
3386        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
3387    else:       #dQ/Q in data
3388        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
3389    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
3390
3391def abeles(kz, depth, rho, irho=0, sigma=0):
3392    """
3393    Optical matrix form of the reflectivity calculation.
3394    O.S. Heavens, Optical Properties of Thin Solid Films
3395   
3396    Reflectometry as a function of kz for a set of slabs.
3397
3398    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\\pi\\sin(\\theta)/\\lambda`.
3399        This is :math:`\\tfrac12 Q_z`.       
3400    :param depth:  float[m] (Ang).
3401        thickness of each layer.  The thickness of the incident medium
3402        and substrate are ignored.
3403    :param rho:  float[n,k] (1e-6/Ang^2)
3404        Real scattering length density for each layer for each kz
3405    :param irho:  float[n,k] (1e-6/Ang^2)
3406        Imaginary scattering length density for each layer for each kz
3407        Note: absorption cross section mu = 2 irho/lambda for neutrons
3408    :param sigma: float[m-1] (Ang)
3409        interfacial roughness.  This is the roughness between a layer
3410        and the previous layer. The sigma array should have m-1 entries.
3411
3412    Slabs are ordered with the surface SLD at index 0 and substrate at
3413    index -1, or reversed if kz < 0.
3414    """
3415    def calc(kz, depth, rho, irho, sigma):
3416        if len(kz) == 0: return kz
3417   
3418        # Complex index of refraction is relative to the incident medium.
3419        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
3420        # in place of kz^2, and ignoring rho_o
3421        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
3422        k = kz
3423   
3424        # According to Heavens, the initial matrix should be [ 1 F; F 1],
3425        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
3426        # multiply versus some coding convenience.
3427        B11 = 1
3428        B22 = 1
3429        B21 = 0
3430        B12 = 0
3431        for i in range(0, len(depth)-1):
3432            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
3433            F = (k - k_next) / (k + k_next)
3434            F *= np.exp(-2*k*k_next*sigma[i]**2)
3435            #print "==== layer",i
3436            #print "kz:", kz
3437            #print "k:", k
3438            #print "k_next:",k_next
3439            #print "F:",F
3440            #print "rho:",rho[:,i+1]
3441            #print "irho:",irho[:,i+1]
3442            #print "d:",depth[i],"sigma:",sigma[i]
3443            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
3444            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
3445            M21 = F*M11
3446            M12 = F*M22
3447            C1 = B11*M11 + B21*M12
3448            C2 = B11*M21 + B21*M22
3449            B11 = C1
3450            B21 = C2
3451            C1 = B12*M11 + B22*M12
3452            C2 = B12*M21 + B22*M22
3453            B12 = C1
3454            B22 = C2
3455            k = k_next
3456   
3457        r = B12/B11
3458        return np.absolute(r)**2
3459
3460    if np.isscalar(kz): kz = np.asarray([kz], 'd')
3461
3462    m = len(depth)
3463
3464    # Make everything into arrays
3465    depth = np.asarray(depth,'d')
3466    rho = np.asarray(rho,'d')
3467    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
3468    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
3469
3470    # Repeat rho,irho columns as needed
3471    if len(rho.shape) == 1:
3472        rho = rho[None,:]
3473        irho = irho[None,:]
3474
3475    return calc(kz, depth, rho, irho, sigma)
3476   
3477def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
3478    y = abeles(kz, depth, rho, irho, sigma)
3479    s = dq/2.
3480    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
3481    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
3482    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
3483    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
3484    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
3485    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
3486    y *= 0.137023
3487    return y
3488       
3489def makeRefdFFT(Limits,Profile):
3490    G2fil.G2Print ('make fft')
3491    Q,Io = Profile[:2]
3492    Qmin = Limits[1][0]
3493    Qmax = Limits[1][1]
3494    iBeg = np.searchsorted(Q,Qmin)
3495    iFin = np.searchsorted(Q,Qmax)+1    #include last point
3496    Qf = np.linspace(0.,Q[iFin-1],5000)
3497    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
3498    If = QI(Qf)*Qf**4
3499    R = np.linspace(0.,5000.,5000)
3500    F = fft.rfft(If)
3501    return R,F
3502
3503   
3504#### Stacking fault simulation codes ################################################################################
3505def GetStackParms(Layers):
3506   
3507    Parms = []
3508#cell parms
3509    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
3510        Parms.append('cellA')
3511        Parms.append('cellC')
3512    else:
3513        Parms.append('cellA')
3514        Parms.append('cellB')
3515        Parms.append('cellC')
3516        if Layers['Laue'] != 'mmm':
3517            Parms.append('cellG')
3518#Transition parms
3519    for iY in range(len(Layers['Layers'])):
3520        for iX in range(len(Layers['Layers'])):
3521            Parms.append('TransP;%d;%d'%(iY,iX))
3522            Parms.append('TransX;%d;%d'%(iY,iX))
3523            Parms.append('TransY;%d;%d'%(iY,iX))
3524            Parms.append('TransZ;%d;%d'%(iY,iX))
3525    return Parms
3526
3527def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
3528    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
3529   
3530    :param dict Layers: dict with following items
3531
3532      ::
3533
3534       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
3535       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
3536       'Layers':[],'Stacking':[],'Transitions':[]}
3537       
3538    :param str ctrls: controls string to be written on DIFFaX controls.dif file
3539    :param float scale: scale factor
3540    :param dict background: background parameters
3541    :param list limits: min/max 2-theta to be calculated
3542    :param dict inst: instrument parameters dictionary
3543    :param list profile: powder pattern data
3544   
3545    Note that parameters all updated in place   
3546    '''
3547    import atmdata
3548    path = sys.path
3549    for name in path:
3550        if 'bin' in name:
3551            DIFFaX = name+'/DIFFaX.exe'
3552            G2fil.G2Print (' Execute '+DIFFaX)
3553            break
3554    # make form factor file that DIFFaX wants - atom types are GSASII style
3555    sf = open('data.sfc','w')
3556    sf.write('GSASII special form factor file for DIFFaX\n\n')
3557    atTypes = list(Layers['AtInfo'].keys())
3558    if 'H' not in atTypes:
3559        atTypes.insert(0,'H')
3560    for atType in atTypes:
3561        if atType == 'H': 
3562            blen = -.3741
3563        else:
3564            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
3565        Adat = atmdata.XrayFF[atType]
3566        text = '%4s'%(atType.ljust(4))
3567        for i in range(4):
3568            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
3569        text += '%11.6f%11.6f'%(Adat['fc'],blen)
3570        text += '%3d\n'%(Adat['Z'])
3571        sf.write(text)
3572    sf.close()
3573    #make DIFFaX control.dif file - future use GUI to set some of these flags
3574    cf = open('control.dif','w')
3575    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
3576        x0 = profile[0]
3577        iBeg = np.searchsorted(x0,limits[0])
3578        iFin = np.searchsorted(x0,limits[1])+1
3579        if iFin-iBeg > 20000:
3580            iFin = iBeg+20000
3581        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3582        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3583        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
3584    else:
3585        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3586        inst = {'Type':['XSC','XSC',]}
3587    cf.close()
3588    #make DIFFaX data file
3589    df = open('GSASII-DIFFaX.dat','w')
3590    df.write('INSTRUMENTAL\n')
3591    if 'X' in inst['Type'][0]:
3592        df.write('X-RAY\n')
3593    elif 'N' in inst['Type'][0]:
3594        df.write('NEUTRON\n')
3595    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
3596        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
3597        U = ateln2*inst['U'][1]/10000.
3598        V = ateln2*inst['V'][1]/10000.
3599        W = ateln2*inst['W'][1]/10000.
3600        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
3601        HW = np.sqrt(np.mean(HWHM))
3602    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
3603        if 'Mean' in Layers['selInst']:
3604            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
3605        elif 'Gaussian' in Layers['selInst']:
3606            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
3607        else:
3608            df.write('None\n')
3609    else:
3610        df.write('0.10\nNone\n')
3611    df.write('STRUCTURAL\n')
3612    a,b,c = Layers['Cell'][1:4]
3613    gam = Layers['Cell'][6]
3614    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
3615    laue = Layers['Laue']
3616    if laue == '2/m(ab)':
3617        laue = '2/m(1)'
3618    elif laue == '2/m(c)':
3619        laue = '2/m(2)'
3620    if 'unknown' in Layers['Laue']:
3621        df.write('%s %.3f\n'%(laue,Layers['Toler']))
3622    else:   
3623        df.write('%s\n'%(laue))
3624    df.write('%d\n'%(len(Layers['Layers'])))
3625    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
3626        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
3627    layerNames = []
3628    for layer in Layers['Layers']:
3629        layerNames.append(layer['Name'])
3630    for il,layer in enumerate(Layers['Layers']):
3631        if layer['SameAs']:
3632            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
3633            continue
3634        df.write('LAYER %d\n'%(il+1))
3635        if '-1' in layer['Symm']:
3636            df.write('CENTROSYMMETRIC\n')
3637        else:
3638            df.write('NONE\n')
3639        for ia,atom in enumerate(layer['Atoms']):
3640            [name,atype,x,y,z,frac,Uiso] = atom
3641            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
3642                frac /= 2.
3643            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
3644    df.write('STACKING\n')
3645    df.write('%s\n'%(Layers['Stacking'][0]))
3646    if 'recursive' in Layers['Stacking'][0]:
3647        df.write('%s\n'%Layers['Stacking'][1])
3648    else:
3649        if 'list' in Layers['Stacking'][1]:
3650            Slen = len(Layers['Stacking'][2])
3651            iB = 0
3652            iF = 0
3653            while True:
3654                iF += 68
3655                if iF >= Slen:
3656                    break
3657                iF = min(iF,Slen)
3658                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
3659                iB = iF
3660        else:
3661            df.write('%s\n'%Layers['Stacking'][1])   
3662    df.write('TRANSITIONS\n')
3663    for iY in range(len(Layers['Layers'])):
3664        sumPx = 0.
3665        for iX in range(len(Layers['Layers'])):
3666            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
3667            p = round(p,3)
3668            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
3669            sumPx += p
3670        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
3671            G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
3672            df.close()
3673            os.remove('data.sfc')
3674            os.remove('control.dif')
3675            os.remove('GSASII-DIFFaX.dat')
3676            return       
3677    df.close()   
3678    time0 = time.time()
3679    try:
3680        subp.call(DIFFaX)
3681    except OSError:
3682        G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
3683    G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
3684    if os.path.exists('GSASII-DIFFaX.spc'):
3685        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
3686        iFin = iBeg+Xpat.shape[1]
3687        bakType,backDict,backVary = SetBackgroundParms(background)
3688        backDict['Lam1'] = G2mth.getWave(inst)
3689        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
3690        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
3691        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
3692            rv = st.poisson(profile[3][iBeg:iFin])
3693            profile[1][iBeg:iFin] = rv.rvs()
3694            Z = np.ones_like(profile[3][iBeg:iFin])
3695            Z[1::2] *= -1
3696            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3697            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3698        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3699    #cleanup files..
3700        os.remove('GSASII-DIFFaX.spc')
3701    elif os.path.exists('GSASII-DIFFaX.sadp'):
3702        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
3703        Sadp = np.reshape(Sadp,(256,-1))
3704        Layers['Sadp']['Img'] = Sadp
3705        os.remove('GSASII-DIFFaX.sadp')
3706    os.remove('data.sfc')
3707    os.remove('control.dif')
3708    os.remove('GSASII-DIFFaX.dat')
3709   
3710def SetPWDRscan(inst,limits,profile):
3711   
3712    wave = G2mth.getMeanWave(inst)
3713    x0 = profile[0]
3714    iBeg = np.searchsorted(x0,limits[0])
3715    iFin = np.searchsorted(x0,limits[1])
3716    if iFin-iBeg > 20000:
3717        iFin = iBeg+20000
3718    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3719    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
3720    return iFin-iBeg
3721       
3722def SetStackingSF(Layers,debug):
3723# Load scattering factors into DIFFaX arrays
3724    import atmdata
3725    atTypes = Layers['AtInfo'].keys()
3726    aTypes = []
3727    for atype in atTypes:
3728        aTypes.append('%4s'%(atype.ljust(4)))
3729    SFdat = []
3730    for atType in atTypes:
3731        Adat = atmdata.XrayFF[atType]
3732        SF = np.zeros(9)
3733        SF[:8:2] = Adat['fa']
3734        SF[1:8:2] = Adat['fb']
3735        SF[8] = Adat['fc']
3736        SFdat.append(SF)
3737    SFdat = np.array(SFdat)
3738    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
3739   
3740def SetStackingClay(Layers,Type):
3741# Controls
3742    rand.seed()
3743    ranSeed = rand.randint(1,2**16-1)
3744    try:   
3745        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
3746            '6/m','6/mmm'].index(Layers['Laue'])+1
3747    except ValueError:  #for 'unknown'
3748        laueId = -1
3749    if 'SADP' in Type:
3750        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
3751        lmax = int(Layers['Sadp']['Lmax'])
3752    else:
3753        planeId = 0
3754        lmax = 0
3755# Sequences
3756    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
3757    try:
3758        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
3759    except ValueError:
3760        StkParm = -1
3761    if StkParm == 2:    #list
3762        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
3763        Nstk = len(StkSeq)
3764    else:
3765        Nstk = 1
3766        StkSeq = [0,]
3767    if StkParm == -1:
3768        StkParm = int(Layers['Stacking'][1])
3769    Wdth = Layers['Width'][0]
3770    mult = 1
3771    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
3772    LaueSym = Layers['Laue'].ljust(12)
3773    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
3774    return laueId,controls
3775   
3776def SetCellAtoms(Layers):
3777    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
3778# atoms in layers
3779    atTypes = list(Layers['AtInfo'].keys())
3780    AtomXOU = []
3781    AtomTp = []
3782    LayerSymm = []
3783    LayerNum = []
3784    layerNames = []
3785    Natm = 0
3786    Nuniq = 0
3787    for layer in Layers['Layers']:
3788        layerNames.append(layer['Name'])
3789    for il,layer in enumerate(Layers['Layers']):
3790        if layer['SameAs']:
3791            LayerNum.append(layerNames.index(layer['SameAs'])+1)
3792            continue
3793        else:
3794            LayerNum.append(il+1)
3795            Nuniq += 1
3796        if '-1' in layer['Symm']:
3797            LayerSymm.append(1)
3798        else:
3799            LayerSymm.append(0)
3800        for ia,atom in enumerate(layer['Atoms']):
3801            [name,atype,x,y,z,frac,Uiso] = atom
3802            Natm += 1
3803            AtomTp.append('%4s'%(atype.ljust(4)))
3804            Ta = atTypes.index(atype)+1
3805            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
3806    AtomXOU = np.array(AtomXOU)
3807    Nlayers = len(layerNames)
3808    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
3809    return Nlayers
3810   
3811def SetStackingTrans(Layers,Nlayers):
3812# Transitions
3813    TransX = []
3814    TransP = []
3815    for Ytrans in Layers['Transitions']:
3816        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
3817        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
3818    TransP = np.array(TransP,dtype='float').T
3819    TransX = np.array(TransX,dtype='float')
3820#    GSASIIpath.IPyBreak()
3821    pyx.pygettrans(Nlayers,TransP,TransX)
3822   
3823def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
3824# Scattering factors
3825    SetStackingSF(Layers,debug)
3826# Controls & sequences
3827    laueId,controls = SetStackingClay(Layers,'PWDR')
3828# cell & atoms
3829    Nlayers = SetCellAtoms(Layers)
3830    Volume = Layers['Cell'][7]   
3831# Transitions
3832    SetStackingTrans(Layers,Nlayers)
3833# PWDR scan
3834    Nsteps = SetPWDRscan(inst,limits,profile)
3835# result as Spec
3836    x0 = profile[0]
3837    profile[3] = np.zeros(len(profile[0]))
3838    profile[4] = np.zeros(len(profile[0]))
3839    profile[5] = np.zeros(len(profile[0]))
3840    iBeg = np.searchsorted(x0,limits[0])
3841    iFin = np.searchsorted(x0,limits[1])+1
3842    if iFin-iBeg > 20000:
3843        iFin = iBeg+20000
3844    Nspec = 20001       
3845    spec = np.zeros(Nspec,dtype='double')   
3846    time0 = time.time()
3847    pyx.pygetspc(controls,Nspec,spec)
3848    G2fil.G2Print (' GETSPC time = %.2fs'%(time.time()-time0))
3849    time0 = time.time()
3850    U = ateln2*inst['U'][1]/10000.
3851    V = ateln2*inst['V'][1]/10000.
3852    W = ateln2*inst['W'][1]/10000.
3853    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
3854    HW = np.sqrt(np.mean(HWHM))
3855    BrdSpec = np.zeros(Nsteps)
3856    if 'Mean' in Layers['selInst']:
3857        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
3858    elif 'Gaussian' in Layers['selInst']:
3859        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
3860    else:
3861        BrdSpec = spec[:Nsteps]
3862    BrdSpec /= Volume
3863    iFin = iBeg+Nsteps
3864    bakType,backDict,backVary = SetBackgroundParms(background)
3865    backDict['Lam1'] = G2mth.getWave(inst)
3866    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
3867    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
3868    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
3869        try:
3870            rv = st.poisson(profile[3][iBeg:iFin])
3871            profile[1][iBeg:iFin] = rv.rvs()
3872        except ValueError:
3873            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
3874        Z = np.ones_like(profile[3][iBeg:iFin])
3875        Z[1::2] *= -1
3876        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3877        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3878    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3879    G2fil.G2Print (' Broadening time = %.2fs'%(time.time()-time0))
3880   
3881def CalcStackingSADP(Layers,debug):
3882   
3883# Scattering factors
3884    SetStackingSF(Layers,debug)
3885# Controls & sequences
3886    laueId,controls = SetStackingClay(Layers,'SADP')
3887# cell & atoms
3888    Nlayers = SetCellAtoms(Layers)   
3889# Transitions
3890    SetStackingTrans(Layers,Nlayers)
3891# result as Sadp
3892    Nspec = 20001       
3893    spec = np.zeros(Nspec,dtype='double')   
3894    time0 = time.time()
3895    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
3896    Sapd = np.zeros((256,256))
3897    iB = 0
3898    for i in range(hkLim):
3899        iF = iB+Nblk
3900        p1 = 127+int(i*Incr)
3901        p2 = 128-int(i*Incr)
3902        if Nblk == 128:
3903            if i:
3904                Sapd[128:,p1] = spec[iB:iF]
3905                Sapd[:128,p1] = spec[iF:iB:-1]
3906            Sapd[128:,p2] = spec[iB:iF]
3907            Sapd[:128,p2] = spec[iF:iB:-1]
3908        else:
3909            if i:
3910                Sapd[:,p1] = spec[iB:iF]
3911            Sapd[:,p2] = spec[iB:iF]
3912        iB += Nblk
3913    Layers['Sadp']['Img'] = Sapd
3914    G2fil.G2Print (' GETSAD time = %.2fs'%(time.time()-time0))
3915   
3916#### Maximum Entropy Method - Dysnomia ###############################################################################
3917def makePRFfile(data,MEMtype):
3918    ''' makes Dysnomia .prf control file from Dysnomia GUI controls
3919   
3920    :param dict data: GSAS-II phase data
3921    :param int MEMtype: 1 for neutron data with negative scattering lengths
3922                        0 otherwise
3923    :returns str: name of Dysnomia control file
3924    '''
3925
3926    generalData = data['General']
3927    pName = generalData['Name'].replace(' ','_')
3928    DysData = data['Dysnomia']
3929    prfName = pName+'.prf'
3930    prf = open(prfName,'w')
3931    prf.write('$PREFERENCES\n')
3932    prf.write(pName+'.mem\n') #or .fos?
3933    prf.write(pName+'.out\n')
3934    prf.write(pName+'.pgrid\n')
3935    prf.write(pName+'.fba\n')
3936    prf.write(pName+'_eps.raw\n')
3937    prf.write('%d\n'%MEMtype)
3938    if DysData['DenStart'] == 'uniform':
3939        prf.write('0\n')
3940    else:
3941        prf.write('1\n')
3942    if DysData['Optimize'] == 'ZSPA':
3943        prf.write('0\n')
3944    else:
3945        prf.write('1\n')
3946    prf.write('1\n')
3947    if DysData['Lagrange'][0] == 'user':
3948        prf.write('0\n')
3949    else:
3950        prf.write('1\n')
3951    prf.write('%.4f %d\n'%(DysData['Lagrange'][1],DysData['wt pwr']))
3952    prf.write('%.3f\n'%DysData['Lagrange'][2])
3953    prf.write('%.2f\n'%DysData['E_factor'])
3954    prf.write('1\n')
3955    prf.write('0\n')
3956    prf.write('%d\n'%DysData['Ncyc'])
3957    prf.write('1\n')
3958    prf.write('1 0 0 0 0 0 0 0\n')
3959    if DysData['prior'] == 'uniform':
3960        prf.write('0\n')
3961    else:
3962        prf.write('1\n')
3963    prf.close()
3964    return prfName
3965
3966def makeMEMfile(data,reflData,MEMtype,DYSNOMIA):
3967    ''' make Dysnomia .mem file of reflection data, etc.
3968
3969    :param dict data: GSAS-II phase data
3970    :param list reflData: GSAS-II reflection data
3971    :param int MEMtype: 1 for neutron data with negative scattering lengths
3972                        0 otherwise
3973    :param str DYSNOMIA: path to dysnomia.exe
3974    '''
3975   
3976    DysData = data['Dysnomia']
3977    generalData = data['General']
3978    cell = generalData['Cell'][1:7]
3979    A = G2lat.cell2A(cell)
3980    SGData = generalData['SGData']
3981    pName = generalData['Name'].replace(' ','_')
3982    memName = pName+'.mem'
3983    Map = generalData['Map']
3984    Type = Map['Type']
3985    UseList = Map['RefList']
3986    mem = open(memName,'w')
3987    mem.write('%s\n'%(generalData['Name']+' from '+UseList[0]))
3988    a,b,c,alp,bet,gam = cell
3989    mem.write('%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n'%(a,b,c,alp,bet,gam))
3990    mem.write('      0.0000000      0.0000000     -1    0    0    0     P\n')   #dummy PO stuff
3991    SGSym = generalData['SGData']['SpGrp']
3992    try:
3993        SGId = G2spc.spgbyNum.index(SGSym)
3994    except ValueError:
3995        return False
3996    org = 1
3997    if SGSym in G2spc.spg2origins:
3998        org = 2
3999    mapsize = Map['rho'].shape
4000    sumZ = 0.
4001    sumpos = 0.
4002    sumneg = 0.
4003    mem.write('%5d%5d%5d%5d%5d\n'%(SGId,org,mapsize[0],mapsize[1],mapsize[2]))
4004    for atm in generalData['NoAtoms']:
4005        Nat = generalData['NoAtoms'][atm]
4006        AtInfo = G2elem.GetAtomInfo(atm)
4007        sumZ += Nat*AtInfo['Z']
4008        isotope = generalData['Isotope'][atm]
4009        blen = generalData['Isotopes'][atm][isotope]['SL'][0]
4010        if blen < 0.:
4011            sumneg += blen*Nat
4012        else:
4013            sumpos += blen*Nat
4014    if 'X' in Type:
4015        mem.write('%10.2f  0.001\n'%sumZ)
4016    elif 'N' in Type and MEMtype:
4017        mem.write('%10.3f%10.3f 0.001\n'%(sumpos,sumneg))
4018    else:
4019        mem.write('%10.3f 0.001\n'%sumpos)
4020       
4021    dmin = DysData['MEMdmin']
4022    TOFlam = 2.0*dmin*npsind(80.0)
4023    refSet = G2lat.GenHLaue(dmin,SGData,A)      #list of h,k,l,d
4024    refDict = {'%d %d %d'%(ref[0],ref[1],ref[2]):ref for ref in refSet}
4025       
4026    refs = []
4027    prevpos = 0.
4028    for ref in reflData:
4029        if ref[3] < 0:
4030            continue
4031        if 'T' in Type:
4032            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,x,x,x,prfo = ref[:16]
4033            s = np.sqrt(max(sig,0.0001))   #var -> sig in deg
4034            FWHM = getgamFW(gam,s)
4035            if dsp < dmin:
4036                continue
4037            theta = npasind(TOFlam/(2.*dsp))
4038            FWHM *= nptand(theta)/pos
4039            pos = 2.*theta
4040        else:
4041            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,prfo = ref[:13]
4042            g = gam/100.    #centideg -> deg
4043            s = np.sqrt(max(sig,0.0001))/100.   #var -> sig in deg
4044            FWHM = getgamFW(g,s)
4045        delt = pos-prevpos
4046        refs.append([h,k,l,mult,pos,FWHM,Fobs,phase,delt])
4047        prevpos = pos
4048           
4049    ovlp = DysData['overlap']
4050    refs1 = []
4051    refs2 = []
4052    nref2 = 0
4053    iref = 0
4054    Nref = len(refs)
4055    start = False
4056    while iref < Nref-1:
4057        if refs[iref+1][-1] < ovlp*refs[iref][5]:
4058            if refs[iref][-1] > ovlp*refs[iref][5]:
4059                refs2.append([])
4060                start = True
4061            if nref2 == len(refs2):
4062                refs2.append([])
4063            refs2[nref2].append(refs[iref])
4064        else:
4065            if start:
4066                refs2[nref2].append(refs[iref])
4067                start = False
4068                nref2 += 1
4069            else:
4070                refs1.append(refs[iref])
4071        iref += 1
4072    if start:
4073        refs2[nref2].append(refs[iref])
4074    else:
4075        refs1.append(refs[iref])
4076   
4077    mem.write('%5d\n'%len(refs1))
4078    for ref in refs1:
4079        h,k,l = ref[:3]
4080        hkl = '%d %d %d'%(h,k,l)
4081        if hkl in refDict:
4082            del refDict[hkl]
4083        Fobs = np.sqrt(ref[6])
4084        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)))
4085    while True and nref2:
4086        if not len(refs2[-1]):
4087            del refs2[-1]
4088        else:
4089            break
4090    mem.write('%5d\n'%len(refs2))
4091    for iref2,ref2 in enumerate(refs2):
4092        mem.write('#%5d\n'%iref2)
4093        mem.write('%5d\n'%len(ref2))
4094        Gsum = 0.
4095        Msum = 0
4096        for ref in ref2:
4097            Gsum += ref[6]*ref[3]
4098            Msum += ref[3]
4099        G = np.sqrt(Gsum/Msum)
4100        h,k,l = ref2[0][:3]
4101        hkl = '%d %d %d'%(h,k,l)
4102        if hkl in refDict:
4103            del refDict[hkl]
4104        mem.write('%5d%5d%5d%10.3f%10.3f%5d\n'%(h,k,l,G,max(0.01*G,0.1),ref2[0][3]))
4105        for ref in ref2[1:]:
4106            h,k,l,m = ref[:4]
4107            mem.write('%5d%5d%5d%5d\n'%(h,k,l,m))
4108            hkl = '%d %d %d'%(h,k,l)
4109            if hkl in refDict:
4110                del refDict[hkl]
4111    if len(refDict):
4112        mem.write('%d\n'%len(refDict))
4113        for hkl in list(refDict.keys()):
4114            h,k,l = refDict[hkl][:3]
4115            mem.write('%5d%5d%5d\n'%(h,k,l))
4116    else:
4117        mem.write('0\n')
4118    mem.close()
4119    return True
4120
4121def MEMupdateReflData(prfName,data,reflData):
4122    ''' Update reflection data with new Fosq, phase result from Dysnomia
4123
4124    :param str prfName: phase.mem file name
4125    :param list reflData: GSAS-II reflection data
4126    '''
4127   
4128    generalData = data['General']
4129    Map = generalData['Map']
4130    Type = Map['Type']
4131    cell = generalData['Cell'][1:7]
4132    A = G2lat.cell2A(cell)
4133    reflDict = {}
4134    newRefs = []
4135    for iref,ref in enumerate(reflData):
4136        if ref[3] > 0:
4137            newRefs.append(ref)
4138            reflDict[hash('%5d%5d%5d'%(ref[0],ref[1],ref[2]))] = iref
4139    fbaName = os.path.splitext(prfName)[0]+'.fba'
4140    if os.path.isfile(fbaName):
4141        fba = open(fbaName,'r')
4142    else:
4143        return False
4144    fba.readline()
4145    Nref = int(fba.readline()[:-1])
4146    fbalines = fba.readlines()
4147    for line in fbalines[:Nref]:
4148        info = line.split()
4149        h = int(info[0])
4150        k = int(info[1])
4151        l = int(info[2])
4152        FoR = float(info[3])
4153        FoI = float(info[4])
4154        Fosq = FoR**2+FoI**2
4155        phase = npatan2d(FoI,FoR)
4156        try:
4157            refId = reflDict[hash('%5d%5d%5d'%(h,k,l))]
4158        except KeyError:    #added reflections at end skipped
4159            d = float(1/np.sqrt(G2lat.calc_rDsq([h,k,l],A)))
4160            if 'T' in Type:
4161                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])
4162            else:
4163                newRefs.append([h,k,l,-1,d,0.,0.01,1.0,Fosq,Fosq,phase,1.0,1.0,1.0,1.0])
4164            continue
4165        newRefs[refId][8] = Fosq
4166        newRefs[refId][10] = phase
4167    newRefs = np.array(newRefs)
4168    return True,newRefs
4169   
4170#### testing data
4171NeedTestData = True
4172def TestData():
4173    'needs a doc string'
4174#    global NeedTestData
4175    global bakType
4176    bakType = 'chebyschev'
4177    global xdata
4178    xdata = np.linspace(4.0,40.0,36000)
4179    global parmDict0
4180    parmDict0 = {
4181        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
4182        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
4183        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
4184        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
4185        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'Z':0.0,'SH/L':0.002,
4186        'Back0':5.384,'Back1':-0.015,'Back2':.004,
4187        }
4188    global parmDict1
4189    parmDict1 = {
4190        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
4191        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
4192        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
4193        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
4194        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
4195        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
4196        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'Z':0.0,'SH/L':0.002,
4197        'Back0':36.897,'Back1':-0.508,'Back2':.006,
4198        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
4199        }
4200    global parmDict2
4201    parmDict2 = {
4202        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
4203        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'Z':0.0,'SH/L':0.02,
4204        'Back0':5.,'Back1':-0.02,'Back2':.004,
4205#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
4206        }
4207    global varyList
4208    varyList = []
4209
4210def test0():
4211    if NeedTestData: TestData()
4212    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
4213    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
4214    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,np.zeros_like(xdata),varyList,bakType))
4215    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
4216    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
4217    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,np.zeros_like(xdata),varyList,bakType))
4218   
4219def test1():
4220    if NeedTestData: TestData()
4221    time0 = time.time()
4222    for i in range(100):
4223        getPeakProfile(parmDict1,xdata,np.zeros_like(xdata),varyList,bakType)
4224    G2fil.G2Print ('100+6*Ka1-2 peaks=1200 peaks %.2f'%time.time()-time0)
4225   
4226def test2(name,delt):
4227    if NeedTestData: TestData()
4228    varyList = [name,]
4229    xdata = np.linspace(5.6,5.8,400)
4230    hplot = plotter.add('derivatives test for '+name).gca()
4231    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,np.zeros_like(xdata),varyList,bakType)[0])
4232    y0 = getPeakProfile(parmDict2,xdata,np.zeros_like(xdata),varyList,bakType)
4233    parmDict2[name] += delt
4234    y1 = getPeakProfile(parmDict2,xdata,np.zeros_like(xdata),varyList,bakType)
4235    hplot.plot(xdata,(y1-y0)/delt,'r+')
4236   
4237def test3(name,delt):
4238    if NeedTestData: TestData()
4239    names = ['pos','sig','gam','shl']
4240    idx = names.index(name)
4241    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
4242    xdata = np.linspace(5.6,5.8,800)
4243    dx = xdata[1]-xdata[0]
4244    hplot = plotter.add('derivatives test for '+name).gca()
4245    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
4246    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
4247    myDict[name] += delt
4248    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
4249    hplot.plot(xdata,(y1-y0)/delt,'r+')
4250
4251if __name__ == '__main__':
4252    import GSASIItestplot as plot
4253    global plotter
4254    plotter = plot.PlotNotebook()
4255#    test0()
4256#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)']:
4257    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
4258        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['Z',0.0001],['SH/L',0.00005]]:
4259        test2(name,shft)
4260    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
4261        test3(name,shft)
4262    G2fil.G2Print ("OK")
4263    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.