source: trunk/GSASIIpwd.py @ 5041

Last change on this file since 5041 was 5041, checked in by vondreele, 9 months ago

begin to add PDFfit to RMC stuff

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