source: trunk/GSASIIpwd.py @ 4919

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

change powder peak shape functions to include effect of varying step size.
Impacts scale factors & peak intensities.
New notification system invoked to let users know of change.

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