source: trunk/GSASIIpwd.py @ 5021

Last change on this file since 5021 was 5021, checked in by vondreele, 3 months ago

use copy.deepcopy on element info in PDF stuff to avoid cross talk between PDFs
A currently commented out smoothing of I(Q) tried - seems not to matter
Allow PDF optimizer to handle negative fixed background & allow its input in the GUI

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 179.3 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASIIpwd: Powder calculations module*
5==============================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2021-09-03 20:19:48 +0000 (Fri, 03 Sep 2021) $
10# $Author: vondreele $
11# $Revision: 5021 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 5021 2021-09-03 20:19:48Z 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: 5021 $"
38GSASIIpath.SetVersionNumber("$Revision: 5021 $")
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.close()
2433    else:
2434        if useSamBrd[0]:
2435            wave = G2mth.getWave(inst)
2436            if 'ellipsoidal' not in Size[0]:    #take the isotropic term only
2437                Xsb = 1.8*wave/(np.pi*Size[1][0])
2438        if useSamBrd[1]:
2439            if 'generalized' not in Mustrain[0]:    #take the isotropic term only
2440                Ysb = 0.0180*Mustrain[1][0]/np.pi
2441        prms = ['Bank',
2442                'Lam','Zero','Polariz.',
2443                'U','V','W',
2444                'X','Y']
2445        fname = Name+'.inst'
2446        fl = open(fname,'w')
2447        fl.write('1\n')
2448        fl.write('%d\n'%int(inst[prms[0]][1]))
2449        fl.write('%10.5f%10.5f%10.4f%10d\n'%(inst[prms[1]][1],-100.*inst[prms[2]][1],inst[prms[3]][1],0))
2450        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[4]][1],inst[prms[5]][1],inst[prms[6]][1]))
2451        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[7]][1]+Xsb,inst[prms[8]][1]+Ysb,0.0))   
2452        fl.write('%10.3f%10.3f%10.3f\n'%(0.0,0.0,0.0))
2453        fl.close()
2454    return fname
2455   
2456def MakeBack(PWDdata,Name):
2457    Back = PWDdata['Background'][0]
2458    inst = PWDdata['Instrument Parameters'][0]
2459    if 'chebyschev-1' != Back[0]:
2460        return None
2461    Nback = Back[2]
2462    BackVals = Back[3:]
2463    fname = Name+'.back'
2464    fl = open(fname,'w')
2465    fl.write('%10d\n'%Nback)
2466    for val in BackVals:
2467        if 'T' in inst['Type'][1]:
2468            fl.write('%12.6g\n'%(float(val)))
2469        else:
2470            fl.write('%12.6g\n'%val)
2471    fl.close()
2472    return fname
2473
2474def findDup(Atoms):
2475    Dup = []
2476    Fracs = []
2477    for iat1,at1 in enumerate(Atoms):
2478        if any([at1[0] in dup for dup in Dup]):
2479            continue
2480        else:
2481            Dup.append([at1[0],])
2482            Fracs.append([at1[6],])
2483        for iat2,at2 in enumerate(Atoms[(iat1+1):]):
2484            if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
2485                Dup[-1] += [at2[0],]
2486                Fracs[-1]+= [at2[6],]
2487    return Dup,Fracs
2488
2489def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):   
2490   
2491    Meta = RMCPdict['metadata']
2492    Atseq = RMCPdict['atSeq']
2493    Supercell =  RMCPdict['SuperCell']
2494    generalData = Phase['General']
2495    Dups,Fracs = findDup(Phase['Atoms'])
2496    Sfracs = [np.cumsum(fracs) for fracs in Fracs]
2497    ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
2498    Sample = PWDdata['Sample Parameters']
2499    Meta['temperature'] = Sample['Temperature']
2500    Meta['pressure'] = Sample['Pressure']
2501    Cell = generalData['Cell'][1:7]
2502    Trans = np.eye(3)*np.array(Supercell)
2503    newPhase = copy.deepcopy(Phase)
2504    newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2505    newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans)
2506    GB = G2lat.cell2Gmat( newPhase['General']['Cell'][1:7])[0]
2507    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.
2508    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
2509    Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2510    Natm = np.count_nonzero(Natm-1)
2511    Atoms = newPhase['Atoms']
2512    reset = False
2513   
2514    if ifSfracs:
2515        Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2516        Natm = np.count_nonzero(Natm-1)
2517        Satoms = []
2518        for i in range(len(Atoms)//Natm):
2519            ind = i*Natm
2520            Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
2521        Natoms = []
2522        for satoms in Satoms:
2523            for idup,dup in enumerate(Dups):
2524                ldup = len(dup)
2525                natm = len(satoms)
2526                i = 0
2527                while i < natm:
2528                    if satoms[i][0] in dup:
2529                        atoms = satoms[i:i+ldup]
2530                        try:
2531                            atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2532                            Natoms.append(atom)
2533                        except IndexError:      #what about vacancies?
2534                            if 'Va' not in Atseq:
2535                                reset = True
2536                                Atseq.append('Va')
2537                                RMCPdict['aTypes']['Va'] = 0.0
2538                            atom = atoms[0]
2539                            atom[1] = 'Va'
2540                            Natoms.append(atom)
2541                        i += ldup
2542                    else:
2543                       i += 1
2544    else:
2545        Natoms = Atoms
2546   
2547    NAtype = np.zeros(len(Atseq))
2548    for atom in Natoms:
2549        NAtype[Atseq.index(atom[1])] += 1
2550    NAstr = ['%6d'%i for i in NAtype]
2551    Cell = newPhase['General']['Cell'][1:7]
2552    if os.path.exists(Name+'.his6f'):
2553        os.remove(Name+'.his6f')
2554    if os.path.exists(Name+'.neigh'):
2555        os.remove(Name+'.neigh')
2556    fname = Name+'.rmc6f'
2557    fl = open(fname,'w')
2558    fl.write('(Version 6f format configuration file)\n')
2559    for item in Meta:
2560        fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2561    fl.write('Atom types present:                 %s\n'%'    '.join(Atseq))
2562    fl.write('Number of each atom type:       %s\n'%''.join(NAstr))
2563    fl.write('Number of atoms:                %d\n'%len(Natoms))
2564    fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2565    fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2566            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2567    A,B = G2lat.cell2AB(Cell,True)
2568    fl.write('Lattice vectors (Ang):\n')   
2569    for i in [0,1,2]:
2570        fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2571    fl.write('Atoms (fractional coordinates):\n')
2572    nat = 0
2573    for atm in Atseq:
2574        for iat,atom in enumerate(Natoms):
2575            if atom[1] == atm:
2576                nat += 1
2577                atcode = Atcodes[iat].split(':')
2578                cell = [0,0,0]
2579                if '+' in atcode[1]:
2580                    cell = eval(atcode[1].split('+')[1])
2581                fl.write('%6d%4s  [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(       
2582                        nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2583    fl.close()
2584    return fname,reset
2585
2586def MakeBragg(PWDdata,Name,Phase):
2587    generalData = Phase['General']
2588    Vol = generalData['Cell'][7]
2589    Data = PWDdata['Data']
2590    Inst = PWDdata['Instrument Parameters'][0]
2591    Bank = int(Inst['Bank'][1])
2592    Sample = PWDdata['Sample Parameters']
2593    Scale = Sample['Scale'][0]
2594    if 'X' in Inst['Type'][1]:
2595        Scale *= 2.
2596    Limits = PWDdata['Limits'][1]
2597    Ibeg = np.searchsorted(Data[0],Limits[0])
2598    Ifin = np.searchsorted(Data[0],Limits[1])+1
2599    fname = Name+'.bragg'
2600    fl = open(fname,'w')
2601    fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol))
2602    if 'T' in Inst['Type'][0]:
2603        fl.write('%12s%12s\n'%('   TOF,ms','  I(obs)'))
2604        for i in range(Ibeg,Ifin-1):
2605            fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2606    else:
2607        fl.write('%12s%12s\n'%('   2-theta, deg','  I(obs)'))
2608        for i in range(Ibeg,Ifin-1):
2609            fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))       
2610    fl.close()
2611    return fname
2612
2613def MakeRMCPdat(PWDdata,Name,Phase,RMCPdict):
2614    Meta = RMCPdict['metadata']
2615    Times = RMCPdict['runTimes']
2616    Atseq = RMCPdict['atSeq']
2617    Atypes = RMCPdict['aTypes']
2618    atPairs = RMCPdict['Pairs']
2619    Files = RMCPdict['files']
2620    BraggWt = RMCPdict['histogram'][1]
2621    inst = PWDdata['Instrument Parameters'][0]
2622    try:
2623        refList = PWDdata['Reflection Lists'][Name]['RefList']
2624    except KeyError:
2625        return 'Error - missing reflection list; you must do Refine first'
2626    dMin = refList[-1][4]
2627    gsasType = 'xray2'
2628    if 'T' in inst['Type'][1]:
2629        gsasType = 'gsas3'
2630    elif 'X' in inst['Type'][1]:
2631        XFF = G2elem.GetFFtable(Atseq)
2632        Xfl = open(Name+'.xray','w')
2633        for atm in Atseq:
2634            fa = XFF[atm]['fa']
2635            fb = XFF[atm]['fb']
2636            fc = XFF[atm]['fc']
2637            Xfl.write('%2s  %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%(
2638                    atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc))
2639        Xfl.close()
2640    lenA = len(Atseq)
2641    Pairs = []
2642    for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]:
2643        Pairs += pair
2644    pairMin = [atPairs[pair]for pair in Pairs if pair in atPairs]
2645    maxMoves = [Atypes[atm] for atm in Atseq if atm in Atypes]
2646    fname = Name+'.dat'
2647    fl = open(fname,'w')
2648    fl.write(' %% Hand edit the following as needed\n')
2649    fl.write('TITLE :: '+Name+'\n')
2650    fl.write('MATERIAL :: '+Meta['material']+'\n')
2651    fl.write('PHASE :: '+Meta['phase']+'\n')
2652    fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n')
2653    fl.write('INVESTIGATOR :: '+Meta['owner']+'\n')
2654    minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin])
2655    minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin])
2656    maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin])
2657    fl.write('MINIMUM_DISTANCES ::   %s  Angstrom\n'%minHD)
2658    maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves])
2659    fl.write('MAXIMUM_MOVES ::   %s Angstrom\n'%maxMv)
2660    fl.write('R_SPACING ::  0.0200 Angstrom\n')
2661    fl.write('PRINT_PERIOD :: 100\n')
2662    fl.write('TIME_LIMIT ::     %.2f MINUTES\n'%Times[0])
2663    fl.write('SAVE_PERIOD ::    %.2f MINUTES\n'%Times[1])
2664    fl.write('\n')
2665    fl.write('ATOMS :: '+' '.join(Atseq)+'\n')
2666    fl.write('\n')
2667    fl.write('FLAGS ::\n')
2668    fl.write('  > NO_MOVEOUT\n')
2669    fl.write('  > NO_SAVE_CONFIGURATIONS\n')
2670    fl.write('  > NO_RESOLUTION_CONVOLUTION\n')
2671    fl.write('\n')
2672    fl.write('INPUT_CONFIGURATION_FORMAT ::  rmc6f\n')
2673    fl.write('SAVE_CONFIGURATION_FORMAT  ::  rmc6f\n')
2674    fl.write('IGNORE_HISTORY_FILE ::\n')
2675    fl.write('\n')
2676    fl.write('DISTANCE_WINDOW ::\n')
2677    fl.write('  > MNDIST :: %s\n'%minD)
2678    fl.write('  > MXDIST :: %s\n'%maxD)
2679    if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']):
2680        fl.write('\n')
2681        fl.write('POTENTIALS ::\n')
2682        fl.write('  > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.'])
2683        fl.write('  > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n')
2684        if len(RMCPdict['Potentials']['Stretch']):
2685            fl.write('  > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search'])
2686            for bond in RMCPdict['Potentials']['Stretch']:
2687                fl.write('  > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2]))       
2688        if len(RMCPdict['Potentials']['Angles']):
2689            fl.write('  > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
2690            for angle in RMCPdict['Potentials']['Angles']:
2691                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
2692                    (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
2693    if RMCPdict['useBVS']:
2694        fl.write('BVS ::\n')
2695        fl.write('  > ATOM :: '+' '.join(Atseq)+'\n')
2696        fl.write('  > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2697        oxid = []
2698        for val in RMCPdict['Oxid']:
2699            if len(val) == 3:
2700                oxid.append(val[0][1:])
2701            else:
2702                oxid.append(val[0][2:])
2703        fl.write('  > OXID :: %s\n'%' '.join(oxid))
2704        fl.write('  > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']]))
2705        fl.write('  > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']]))
2706        fl.write('  > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))       
2707        fl.write('  > SAVE :: 100000\n')
2708        fl.write('  > UPDATE :: 100000\n')
2709        if len(RMCPdict['Swap']):
2710            fl.write('\n')
2711            fl.write('SWAP_MULTI ::\n')
2712            for swap in RMCPdict['Swap']:
2713                try:
2714                    at1 = Atseq.index(swap[0])
2715                    at2 = Atseq.index(swap[1])
2716                except ValueError:
2717                    break
2718                fl.write('  > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2]))
2719       
2720    if len(RMCPdict['FxCN']):
2721        fl.write('FIXED_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['FxCN']))       
2722        for ifx,fxcn in enumerate(RMCPdict['FxCN']):
2723            try:
2724                at1 = Atseq.index(fxcn[0])
2725                at2 = Atseq.index(fxcn[1])
2726            except ValueError:
2727                break
2728            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]))
2729    if len(RMCPdict['AveCN']):
2730        fl.write('AVERAGE_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['AveCN']))
2731        for iav,avcn in enumerate(RMCPdict['AveCN']):
2732            try:
2733                at1 = Atseq.index(avcn[0])
2734                at2 = Atseq.index(avcn[1])
2735            except ValueError:
2736                break
2737            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]))
2738    for File in Files:
2739        if Files[File][0] and Files[File][0] != 'Select':
2740            if 'Xray' in File and 'F(Q)' in File:
2741                fqdata = open(Files[File][0],'r')
2742                lines = int(fqdata.readline()[:-1])
2743            fl.write('\n')
2744            fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_'))
2745            fl.write('  > FILENAME :: %s\n'%Files[File][0])
2746            fl.write('  > DATA_TYPE :: %s\n'%Files[File][2])
2747            fl.write('  > FIT_TYPE :: %s\n'%Files[File][2])
2748            if 'Xray' not in File:
2749                fl.write('  > START_POINT :: 1\n')
2750                fl.write('  > END_POINT :: 3000\n')
2751                fl.write('  > WEIGHT :: %.4f\n'%Files[File][1])
2752            fl.write('  > CONSTANT_OFFSET 0.000\n')
2753            fl.write('  > NO_FITTED_OFFSET\n')
2754            if RMCPdict['FitScale']:
2755                fl.write('  > FITTED_SCALE\n')
2756            else:
2757                fl.write('  > NO_FITTED_SCALE\n')
2758            if Files[File][3] !='RMC':
2759                fl.write('  > %s\n'%Files[File][3])
2760            if 'reciprocal' in File:
2761                fl.write('  > CONVOLVE ::\n')
2762                if 'Xray' in File:
2763                    fl.write('  > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines)
2764                    fl.write('  > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1]))
2765                    fl.write('  > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2))
2766                    fl.write('  > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,1./Files[File][1]))
2767    fl.write('\n')
2768    fl.write('BRAGG ::\n')
2769    fl.write('  > BRAGG_SHAPE :: %s\n'%gsasType)
2770    fl.write('  > RECALCUATE\n')
2771    fl.write('  > DMIN :: %.2f\n'%(dMin-0.02))
2772    fl.write('  > WEIGHT :: %10.3f\n'%BraggWt)
2773    fl.write('\n')
2774    fl.write('END  ::\n')
2775    fl.close()
2776    return fname
2777
2778# def FindBonds(Phase,RMCPdict):
2779#     generalData = Phase['General']
2780#     cx,ct,cs,cia = generalData['AtomPtrs']
2781#     atomData = Phase['Atoms']
2782#     Res = 'RMC'
2783#     if 'macro' in generalData['Type']:
2784#         Res = atomData[0][ct-3]
2785#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2786#     Pairs = RMCPdict['Pairs']   #dict!
2787#     BondList = []
2788#     notNames = []
2789#     for FrstName in AtDict:
2790#         nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
2791#         Atyp1 = AtDict[FrstName]
2792#         if 'Va' in Atyp1:
2793#             continue
2794#         for nbr in nbrs:
2795#             Atyp2 = AtDict[nbr[0]]
2796#             if 'Va' in Atyp2:
2797#                 continue
2798#             try:
2799#                 bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
2800#             except KeyError:
2801#                 bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
2802#             if any(bndData):
2803#                 if bndData[0] <= nbr[1] <= bndData[1]:
2804#                     bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
2805#                     revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
2806#                     if bondStr not in BondList and revbondStr not in BondList:
2807#                         BondList.append(bondStr)
2808#         notNames.append(FrstName)
2809#     return Res,BondList
2810
2811# def FindAngles(Phase,RMCPdict):
2812#     generalData = Phase['General']
2813#     Cell = generalData['Cell'][1:7]
2814#     Amat = G2lat.cell2AB(Cell)[0]
2815#     cx,ct,cs,cia = generalData['AtomPtrs']
2816#     atomData = Phase['Atoms']
2817#     AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
2818#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2819#     Angles = RMCPdict['Angles']
2820#     AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
2821#     AngleList = []
2822#     for MidName in AtDict:
2823#         nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
2824#         if len(nbrs) < 2: #need 2 neighbors to make an angle
2825#             continue
2826#         Atyp2 = AtDict[MidName]
2827#         for i,nbr1 in enumerate(nbrs):
2828#             Atyp1 = AtDict[nbr1[0]]
2829#             for j,nbr3 in enumerate(nbrs[i+1:]):
2830#                 Atyp3 = AtDict[nbr3[0]]
2831#                 IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
2832#                 try:
2833#                     angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
2834#                 except KeyError:
2835#                     try:
2836#                         angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
2837#                     except KeyError:
2838#                         continue
2839#                 XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
2840#                 calAngle = G2mth.getRestAngle(XYZ,Amat)
2841#                 if angData[0] <= calAngle <= angData[1]:
2842#                     angStr = str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n'
2843#                     revangStr = str((MidName,nbr3[0],nbr1[0])+tuple(angData))+',\n'
2844#                     if angStr not in AngleList and revangStr not in AngleList:
2845#                         AngleList.append(angStr)
2846#     return AngleList
2847
2848# def GetSqConvolution(XY,d):
2849
2850#     n = XY.shape[1]
2851#     snew = np.zeros(n)
2852#     dq = np.zeros(n)
2853#     sold = XY[1]
2854#     q = XY[0]
2855#     dq[1:] = np.diff(q)
2856#     dq[0] = dq[1]
2857   
2858#     for j in range(n):
2859#         for i in range(n):
2860#             b = abs(q[i]-q[j])
2861#             t = q[i]+q[j]
2862#             if j == i:
2863#                 snew[j] += q[i]*sold[i]*(d-np.sin(t*d)/t)*dq[i]
2864#             else:
2865#                 snew[j] += q[i]*sold[i]*(np.sin(b*d)/b-np.sin(t*d)/t)*dq[i]
2866#         snew[j] /= np.pi*q[j]
2867   
2868#     snew[0] = snew[1]
2869#     return snew
2870
2871# def GetMaxSphere(pdbName):
2872#     try:
2873#         pFil = open(pdbName,'r')
2874#     except FileNotFoundError:
2875#         return None
2876#     while True:
2877#         line = pFil.readline()
2878#         if 'Boundary' in line:
2879#             line = line.split()[3:]
2880#             G = np.array([float(item) for item in line])
2881#             G = np.reshape(G,(3,3))**2
2882#             G = nl.inv(G)
2883#             pFil.close()
2884#             break
2885#     dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
2886#     return min(dspaces)
2887
2888def findfullrmc():
2889    '''Find where fullrmc is installed. Tries the following:
2890   
2891         1. Returns the Config var 'fullrmc_exec', if defined. No check
2892            is done that the interpreter has fullrmc
2893         2. The current Python interpreter if fullrmc can be imported
2894            and fullrmc is version 5+
2895         3. The path is checked for a fullrmc image as named by Bachir
2896
2897    :returns: the full path to a python executable that is assumed to
2898      have fullrmc installed or None, if it was not found.
2899    '''
2900    is_exe = lambda fpath: os.path.isfile(fpath) and os.access(fpath, os.X_OK)
2901    if GSASIIpath.GetConfigValue('fullrmc_exec') is not None and is_exe(
2902            GSASIIpath.GetConfigValue('fullrmc_exec')):
2903        return GSASIIpath.GetConfigValue('fullrmc_exec')
2904    try:
2905        import fullrmc
2906        if int(fullrmc.__version__.split('.')[0]) >= 5:
2907            return sys.executable
2908    except:
2909        pass
2910    pathlist = os.environ["PATH"].split(os.pathsep)
2911    for p in (GSASIIpath.path2GSAS2,GSASIIpath.binaryPath,os.getcwd(),
2912                  os.path.split(sys.executable)[0]):
2913        if p not in pathlist: pathlist.insert(0,p)
2914    import glob
2915    for p in pathlist:
2916        if sys.platform == "darwin":
2917            lookfor = "fullrmc*macOS*i386-64bit"
2918        elif sys.platform == "win32":
2919            lookfor = "fullrmc*.exe"
2920        else:
2921            lookfor = "fullrmc*"
2922        fl = glob.glob(lookfor)
2923        if len(fl) > 0:
2924            return os.path.abspath(sorted(fl)[0])
2925       
2926def MakefullrmcRun(pName,Phase,RMCPdict):
2927    '''Creates a script to run fullrmc. Returns the name of the file that was
2928    created.
2929    '''
2930    BondList = {}
2931    for k in RMCPdict['Pairs']:
2932        if RMCPdict['Pairs'][k][1]+RMCPdict['Pairs'][k][2]>0:
2933            BondList[k] = (RMCPdict['Pairs'][k][1],RMCPdict['Pairs'][k][2])
2934    AngleList = []
2935    for angle in RMCPdict['Angles']:
2936        if angle[3] == angle[4] or angle[5] >= angle[6] or angle[6] <= 0:
2937            continue
2938        for i in (0,1,2):
2939            angle[i] = angle[i].strip()
2940        AngleList.append(angle)
2941    # rmin = RMCPdict['min Contact']
2942    cell = Phase['General']['Cell'][1:7]
2943    SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
2944    cx,ct,cs,cia = Phase['General']['AtomPtrs']
2945    atomsList = []
2946    for atom in Phase['Atoms']:
2947        el = ''.join([i for i in atom[ct] if i.isalpha()])
2948        atomsList.append([el] + atom[cx:cx+4])
2949    projDir,projName = os.path.split(os.path.abspath(pName))
2950    scrname = pName+'-fullrmc.py'
2951    restart = '%s_restart.pdb'%pName
2952    Files = RMCPdict['files']
2953    rundata = ''
2954    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%scrname
2955    rundata += '# created in '+__file__+" v"+filversion.split()[1]
2956    rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
2957    rundata += '''
2958# fullrmc imports (all that are potentially useful)
2959import os,glob
2960import time
2961import pickle
2962import numpy as np
2963from fullrmc.Core import Collection
2964from fullrmc.Engine import Engine
2965import fullrmc.Constraints.PairDistributionConstraints as fPDF
2966from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint
2967from fullrmc.Constraints.DistanceConstraints import DistanceConstraint
2968from fullrmc.Constraints.BondConstraints import BondConstraint
2969from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
2970from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
2971from fullrmc.Generators.Swaps import SwapPositionsGenerator
2972# utility routines
2973def writeHeader(ENGINE,statFP):
2974    'header for stats file'
2975    statFP.write('generated-steps, total-error, ')
2976    for c in ENGINE.constraints:
2977        statFP.write(c.constraintName)
2978        statFP.write(', ')
2979    statFP.write('\\n')
2980    statFP.flush()
2981   
2982def writeCurrentStatus(ENGINE,statFP,plotF):
2983    'line in stats file & current constraint plots'
2984    statFP.write(str(ENGINE.generated))
2985    statFP.write(', ')
2986    statFP.write(str(ENGINE.totalStandardError))
2987    statFP.write(', ')
2988    for c in ENGINE.constraints:
2989        statFP.write(str(c.standardError))
2990        statFP.write(', ')
2991    statFP.write('\\n')
2992    statFP.flush()
2993    mpl.use('agg')
2994    fp = open(plotF,'wb')
2995    for c in ENGINE.constraints:
2996        p = c.plot(show=False)
2997        p[0].canvas.draw()
2998        image = p[0].canvas.buffer_rgba()
2999        pickle.dump(c.constraintName,fp)
3000        pickle.dump(np.array(image),fp)
3001    fp.close()
3002
3003def calcRmax(ENGINE):
3004    'from Bachir, works for non-othorhombic'
3005    a,b,c = ENGINE.basisVectors
3006    lens = []
3007    ts    = np.linalg.norm(np.cross(a,b))/2
3008    lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(b)] )
3009    ts = np.linalg.norm(np.cross(b,c))/2
3010    lens.extend( [ts/np.linalg.norm(b), ts/np.linalg.norm(c)] )
3011    ts = np.linalg.norm(np.cross(a,c))/2
3012    lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(c)] )
3013    return min(lens)
3014'''
3015    rundata += '''
3016### When True, erases an existing engine to provide a fresh start
3017FRESH_START = {:}
3018dirName = "{:}"
3019prefix = "{:}"
3020project = prefix + "-fullrmc"
3021time0 = time.time()
3022'''.format(RMCPdict['ReStart'][0],projDir,projName)
3023   
3024    rundata += '# setup structure\n'
3025    rundata += 'cell = ' + str(cell) + '\n'
3026    rundata += "SymOpList = "+str([i.lower() for i in SymOpList]) + '\n'
3027    rundata += 'atomList = ' + str(atomsList).replace('],','],\n  ') + '\n'
3028    rundata += 'supercell = ' + str(RMCPdict['SuperCell']) + '\n'
3029
3030    rundata += '\n# initialize engine\n'
3031    rundata += '''
3032engineFileName = os.path.join(dirName, project + '.rmc')
3033projectStats = os.path.join(dirName, project + '.stats')
3034projectPlots = os.path.join(dirName, project + '.plots')
3035pdbFile = os.path.join(dirName, project + '_restart.pdb')
3036# check Engine exists if so (and not FRESH_START) load it
3037# otherwise build it
3038ENGINE = Engine(path=None)
3039if not ENGINE.is_engine(engineFileName) or FRESH_START:
3040    ## create structure
3041    ENGINE = Engine(path=engineFileName, freshStart=True)
3042    ENGINE.build_crystal_set_pdb(symOps     = SymOpList,
3043                                 atoms      = atomList,
3044                                 unitcellBC = cell,
3045                                 supercell  = supercell)
3046'''   
3047    import atmdata
3048    # rundata += '    # conversion factors (may be needed)\n'
3049    # rundata += '    sumCiBi2 = 0.\n'
3050    # for elem in Phase['General']['AtomTypes']:
3051    #     rundata += '    Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
3052    #     rundata += '    sumCiBi2 += (Ci*{})**2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
3053    rundata += '    rho0 = len(ENGINE.allNames)/ENGINE.volume\n'
3054    # settings that require a new Engine
3055    for File in Files:
3056        filDat = RMCPdict['files'][File]
3057        if not os.path.exists(filDat[0]): continue
3058        sfwt = 'neutronCohb'
3059        if 'Xray' in File:
3060            sfwt = 'atomicNumber'
3061        if 'G(r)' in File:
3062            rundata += '    GR = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3063            if filDat[3] == 0:
3064                #rundata += '''    # read and xform G(r) as defined in RMCProfile
3065    # see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
3066                #rundata += '    GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n'
3067                #rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3068                rundata += '    # G(r) as defined in RMCProfile\n'
3069                rundata += '    GofR = fullrmc.Constraints.RadialDistributionConstraints.RadialDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3070            elif filDat[3] == 1:
3071                rundata += '    # This is G(r) as defined in PDFFIT\n'
3072                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3073            elif filDat[3] == 2:
3074                rundata += '    # This is g(r)\n'
3075                rundata += '    GofR = fPDF.PairCorrelationConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3076            else:
3077                raise ValueError('Invalid G(r) type: '+str(filDat[3]))
3078            rundata += '    ENGINE.add_constraints([GofR])\n'
3079            rundata += '    GofR.set_limits((None, calcRmax(ENGINE)))\n'
3080        elif '(Q)' in File:
3081            rundata += '    SOQ = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3082            if filDat[3] == 0:
3083                rundata += '    # F(Q) as defined in RMCProfile\n'
3084                #rundata += '    SOQ[1] *= 1 / sumCiBi2\n'
3085                if filDat[4]:
3086                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3087                rundata += '    SofQ = fullrmc.Constraints.StructureFactorConstraints.NormalizedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3088            elif filDat[3] == 1:
3089                rundata += '    # S(Q) as defined in PDFFIT\n'
3090                rundata += '    SOQ[1] -= 1\n'
3091                if filDat[4]:
3092                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3093                rundata += '    SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3094            else:
3095                raise ValueError('Invalid S(Q) type: '+str(filDat[3]))
3096            rundata += '    ENGINE.add_constraints([SofQ])\n'
3097        else:
3098            print('What is this?')
3099    rundata += '    ENGINE.add_constraints(DistanceConstraint(defaultLowerDistance={}))\n'.format(RMCPdict['min Contact'])
3100    if BondList:
3101        rundata += '''    B_CONSTRAINT   = BondConstraint()
3102    ENGINE.add_constraints(B_CONSTRAINT)
3103    B_CONSTRAINT.create_supercell_bonds(bondsDefinition=[
3104'''
3105        for pair in BondList:
3106            e1,e2 = pair.split('-')
3107            d1,d2 = BondList[pair]
3108            if d1 == 0: continue
3109            if d2 == 0:
3110                print('Ignoring min distance without maximum')
3111                #rundata += '            ("element","{}","{}",{}),\n'.format(
3112                #                        e1.strip(),e2.strip(),d1)
3113            else:
3114                rundata += '            ("element","{}","{}",{},{}),\n'.format(
3115                                        e1.strip(),e2.strip(),d1,d2)
3116        rundata += '             ])\n'
3117    if AngleList:
3118        rundata += '''    A_CONSTRAINT   = BondsAngleConstraint()
3119    ENGINE.add_constraints(A_CONSTRAINT)
3120    A_CONSTRAINT.create_supercell_angles(anglesDefinition=[
3121'''
3122        for item in AngleList:
3123            rundata += ('            '+
3124               '("element","{1}","{0}","{2}",{5},{6},{5},{6},{3},{4}),\n'.format(*item))
3125        rundata += '             ])\n'
3126    rundata += '''
3127    for f in glob.glob(os.path.join(dirName,prefix+"_*.log")): os.remove(f)
3128    ENGINE.save()
3129else:
3130    ENGINE = ENGINE.load(path=engineFileName)
3131
3132ENGINE.set_log_file(os.path.join(dirName,prefix))
3133'''
3134    if RMCPdict['Swaps']:
3135        rundata += '\n#set up for site swaps\n'
3136        rundata += 'aN = ENGINE.allNames\n'
3137        rundata += 'SwapGen = {}\n'
3138        for swap in RMCPdict['Swaps']:
3139            rundata += 'SwapA = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[0]
3140            rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1]
3141            rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
3142        rundata += '    for swaps in SwapGen:\n'
3143        rundata += '        AB = swaps.split("-")\n'
3144        rundata += '        ENGINE.set_groups_as_atoms()\n'
3145        rundata += '        for g in ENGINE.groups:\n'
3146        rundata += '            if aN[g.indexes[0]]==AB[0]:\n'
3147        rundata += '                g.set_move_generator(SwapGen[swaps][0])\n'
3148        rundata += '            elif aN[g.indexes[0]]==AB[1]:\n'
3149        rundata += '                g.set_move_generator(SwapGen[swaps][1])\n'
3150        rundata += '            sProb = SwapGen[swaps][2]\n'
3151    rundata += '\n# set weights -- do this now so values can be changed without a restart\n'
3152    # rundata += 'wtDict = {}\n'
3153    # for File in Files:
3154    #     filDat = RMCPdict['files'][File]
3155    #     if not os.path.exists(filDat[0]): continue
3156    #     if 'Xray' in File:
3157    #         sfwt = 'atomicNumber'
3158    #     else:
3159    #         sfwt = 'neutronCohb'
3160    #     if 'G(r)' in File:
3161    #         typ = 'Pair'
3162    #     elif '(Q)' in File:
3163    #         typ = 'Struct'
3164    #     rundata += 'wtDict["{}-{}"] = {}\n'.format(typ,sfwt,filDat[1])
3165    rundata += 'for c in ENGINE.constraints:  # loop over predefined constraints\n'
3166    rundata += '    if type(c) is fPDF.PairDistributionConstraint:\n'
3167    # rundata += '        c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
3168    rundata += '        c.set_limits((None,calcRmax(ENGINE)))\n'
3169    if RMCPdict['FitScale']:
3170        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3171    # rundata += '        c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
3172    if RMCPdict['FitScale']:
3173        rundata += '    elif type(c) is ReducedStructureFactorConstraint:\n'
3174        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3175    # torsions difficult to implement, must be internal to cell & named with
3176    # fullrmc atom names
3177    # if len(RMCPdict['Torsions']):         # Torsions currently commented out in GUI
3178    #     rundata += 'for c in ENGINE.constraints:  # look for Dihedral Angle Constraints\n'
3179    #     rundata += '    if type(c) is DihedralAngleConstraint:\n'
3180    #     rundata += '        c.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
3181    #     rundata += '        c.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
3182    #     for torsion in RMCPdict['Torsions']:
3183    #         rundata += '    %s\n'%str(tuple(torsion))
3184    #     rundata += '        ]})\n'           
3185    rundata += '''
3186if FRESH_START:
3187    # initialize engine with one step to get starting config energetics
3188    ENGINE.run(restartPdb=pdbFile,numberOfSteps=1, saveFrequency=1)
3189    statFP = open(projectStats,'w')
3190    writeHeader(ENGINE,statFP)
3191    writeCurrentStatus(ENGINE,statFP,projectPlots)
3192else:
3193    statFP = open(projectStats,'a')
3194
3195# setup runs for fullrmc
3196'''
3197    rundata += 'steps = {}\n'.format(RMCPdict['Steps/cycle'])
3198    rundata += 'for _ in range({}):\n'.format(RMCPdict['Cycles'])
3199    rundata += '    ENGINE.set_groups_as_atoms()\n'
3200    rundata += '    expected = ENGINE.generated+steps\n'
3201   
3202    rundata += '    ENGINE.run(restartPdb=pdbFile,numberOfSteps=steps, saveFrequency=steps)\n'
3203    rundata += '    writeCurrentStatus(ENGINE,statFP,projectPlots)\n'
3204    rundata += '    if ENGINE.generated != expected: break # run was stopped\n'
3205    rundata += 'statFP.close()\n'
3206    rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
3207    rfile = open(scrname,'w')
3208    rfile.writelines(rundata)
3209    rfile.close()
3210    return scrname
3211   
3212def GetRMCBonds(general,RMCPdict,Atoms,bondList):
3213    bondDist = []
3214    Cell = general['Cell'][1:7]
3215    Supercell =  RMCPdict['SuperCell']
3216    Trans = np.eye(3)*np.array(Supercell)
3217    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3218    Amat,Bmat = G2lat.cell2AB(Cell)
3219    indices = (-1,0,1)
3220    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3221    for bonds in bondList:
3222        Oxyz = np.array(Atoms[bonds[0]][1:])
3223        Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])       
3224        Dx = np.array([Txyz-Oxyz+unit for unit in Units])
3225        Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
3226        for dx in Dx.T:
3227            bondDist.append(np.min(dx))
3228    return np.array(bondDist)
3229   
3230def GetRMCAngles(general,RMCPdict,Atoms,angleList):
3231    bondAngles = []
3232    Cell = general['Cell'][1:7]
3233    Supercell =  RMCPdict['SuperCell']
3234    Trans = np.eye(3)*np.array(Supercell)
3235    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3236    Amat,Bmat = G2lat.cell2AB(Cell)
3237    indices = (-1,0,1)
3238    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3239    for angle in angleList:
3240        Oxyz = np.array(Atoms[angle[0]][1:])
3241        TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])       
3242        TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])       
3243        DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
3244        DAx = np.sqrt(np.sum(DAxV**2,axis=2))
3245        DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
3246        DBx = np.sqrt(np.sum(DBxV**2,axis=2))
3247        iDAx = np.argmin(DAx,axis=0)
3248        iDBx = np.argmin(DBx,axis=0)
3249        for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
3250            DAv = DAxV[iA,i]/DAx[iA,i]
3251            DBv = DBxV[iB,i]/DBx[iB,i]
3252            bondAngles.append(npacosd(np.sum(DAv*DBv)))
3253    return np.array(bondAngles)
3254   
3255#### Reflectometry calculations ################################################################################
3256def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
3257    G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
3258   
3259    class RandomDisplacementBounds(object):
3260        """random displacement with bounds"""
3261        def __init__(self, xmin, xmax, stepsize=0.5):
3262            self.xmin = xmin
3263            self.xmax = xmax
3264            self.stepsize = stepsize
3265   
3266        def __call__(self, x):
3267            """take a random step but ensure the new position is within the bounds"""
3268            while True:
3269                # this could be done in a much more clever way, but it will work for example purposes
3270                steps = self.xmax-self.xmin
3271                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
3272                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
3273                    break
3274            return xnew
3275   
3276    def GetModelParms():
3277        parmDict = {}
3278        varyList = []
3279        values = []
3280        bounds = []
3281        parmDict['dQ type'] = data['dQ type']
3282        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
3283        for parm in ['Scale','FltBack']:
3284            parmDict[parm] = data[parm][0]
3285            if data[parm][1]:
3286                varyList.append(parm)
3287                values.append(data[parm][0])
3288                bounds.append(Bounds[parm])
3289        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
3290        parmDict['nLayers'] = len(parmDict['Layer Seq'])
3291        for ilay,layer in enumerate(data['Layers']):
3292            name = layer['Name']
3293            cid = str(ilay)+';'
3294            parmDict[cid+'Name'] = name
3295            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3296                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
3297                if layer.get(parm,[0.,False])[1]:
3298                    varyList.append(cid+parm)
3299                    value = layer[parm][0]
3300                    values.append(value)
3301                    if value:
3302                        bound = [value*Bfac,value/Bfac]
3303                    else:
3304                        bound = [0.,10.]
3305                    bounds.append(bound)
3306            if name not in ['vacuum','unit scatter']:
3307                parmDict[cid+'rho'] = Substances[name]['Scatt density']
3308                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
3309        return parmDict,varyList,values,bounds
3310   
3311    def SetModelParms():
3312        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
3313        if 'Scale' in varyList:
3314            data['Scale'][0] = parmDict['Scale']
3315            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
3316        G2fil.G2Print (line)
3317        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
3318        if 'FltBack' in varyList:
3319            data['FltBack'][0] = parmDict['FltBack']
3320            line += ' esd: %15.3g'%(sigDict['FltBack'])
3321        G2fil.G2Print (line)
3322        for ilay,layer in enumerate(data['Layers']):
3323            name = layer['Name']
3324            G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
3325            cid = str(ilay)+';'
3326            line = ' '
3327            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
3328            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
3329            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3330                if parm in layer:
3331                    if parm == 'Rough':
3332                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
3333                    else:
3334                        layer[parm][0] = parmDict[cid+parm]
3335                    line += ' %s: %.3f'%(parm,layer[parm][0])
3336                    if cid+parm in varyList:
3337                        line += ' esd: %.3g'%(sigDict[cid+parm])
3338            G2fil.G2Print (line)
3339            G2fil.G2Print (line2)
3340   
3341    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3342        parmDict.update(zip(varyList,values))
3343        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3344        return M
3345   
3346    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3347        parmDict.update(zip(varyList,values))
3348        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3349        return np.sum(M**2)
3350   
3351    def getREFD(Q,Qsig,parmDict):
3352        Ic = np.ones_like(Q)*parmDict['FltBack']
3353        Scale = parmDict['Scale']
3354        Nlayers = parmDict['nLayers']
3355        Res = parmDict['Res']
3356        depth = np.zeros(Nlayers)
3357        rho = np.zeros(Nlayers)
3358        irho = np.zeros(Nlayers)
3359        sigma = np.zeros(Nlayers)
3360        for ilay,lay in enumerate(parmDict['Layer Seq']):
3361            cid = str(lay)+';'
3362            depth[ilay] = parmDict[cid+'Thick']
3363            sigma[ilay] = parmDict[cid+'Rough']
3364            if parmDict[cid+'Name'] == u'unit scatter':
3365                rho[ilay] = parmDict[cid+'DenMul']
3366                irho[ilay] = parmDict[cid+'iDenMul']
3367            elif 'vacuum' != parmDict[cid+'Name']:
3368                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
3369                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
3370            if cid+'Mag SLD' in parmDict:
3371                rho[ilay] += parmDict[cid+'Mag SLD']
3372        if parmDict['dQ type'] == 'None':
3373            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
3374        elif 'const' in parmDict['dQ type']:
3375            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
3376        else:       #dQ/Q in data
3377            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
3378        Ic += AB*Scale
3379        return Ic
3380       
3381    def estimateT0(takestep):
3382        Mmax = -1.e-10
3383        Mmin = 1.e10
3384        for i in range(100):
3385            x0 = takestep(values)
3386            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3387            Mmin = min(M,Mmin)
3388            MMax = max(M,Mmax)
3389        return 1.5*(MMax-Mmin)
3390
3391    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3392    if data.get('2% weight'):
3393        wt = 1./(0.02*Io)**2
3394    Qmin = Limits[1][0]
3395    Qmax = Limits[1][1]
3396    wtFactor = ProfDict['wtFactor']
3397    Bfac = data['Toler']
3398    Ibeg = np.searchsorted(Q,Qmin)
3399    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
3400    Ic[:] = 0
3401    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
3402              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
3403    parmDict,varyList,values,bounds = GetModelParms()
3404    Msg = 'Failed to converge'
3405    if varyList:
3406        if data['Minimizer'] == 'LMLS': 
3407            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
3408                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3409            parmDict.update(zip(varyList,result[0]))
3410            chisq = np.sum(result[2]['fvec']**2)
3411            ncalc = result[2]['nfev']
3412            covM = result[1]
3413            newVals = result[0]
3414        elif data['Minimizer'] == 'Basin Hopping':
3415            xyrng = np.array(bounds).T
3416            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
3417            T0 = estimateT0(take_step)
3418            G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
3419            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
3420                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
3421                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
3422            chisq = result.fun
3423            ncalc = result.nfev
3424            newVals = result.x
3425            covM = []
3426        elif data['Minimizer'] == 'MC/SA Anneal':
3427            xyrng = np.array(bounds).T
3428            result = G2mth.anneal(sumREFD, values, 
3429                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
3430                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
3431                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
3432                ranRange=0.20,autoRan=False,dlg=None)
3433            newVals = result[0]
3434            parmDict.update(zip(varyList,newVals))
3435            chisq = result[1]
3436            ncalc = result[3]
3437            covM = []
3438            G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
3439        elif data['Minimizer'] == 'L-BFGS-B':
3440            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
3441                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3442            parmDict.update(zip(varyList,result['x']))
3443            chisq = result.fun
3444            ncalc = result.nfev
3445            newVals = result.x
3446            covM = []
3447    else:   #nothing varied
3448        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3449        chisq = np.sum(M**2)
3450        ncalc = 0
3451        covM = []
3452        sig = []
3453        sigDict = {}
3454        result = []
3455    Rvals = {}
3456    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
3457    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
3458    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
3459    Ib[Ibeg:Ifin] = parmDict['FltBack']
3460    try:
3461        if not len(varyList):
3462            Msg += ' - nothing refined'
3463            raise ValueError
3464        Nans = np.isnan(newVals)
3465        if np.any(Nans):
3466            name = varyList[Nans.nonzero(True)[0]]
3467            Msg += ' Nan result for '+name+'!'
3468            raise ValueError
3469        Negs = np.less_equal(newVals,0.)
3470        if np.any(Negs):
3471            indx = Negs.nonzero()
3472            name = varyList[indx[0][0]]
3473            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
3474                Msg += ' negative coefficient for '+name+'!'
3475                raise ValueError
3476        if len(covM):
3477            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
3478            covMatrix = covM*Rvals['GOF']
3479        else:
3480            sig = np.zeros(len(varyList))
3481            covMatrix = []
3482        sigDict = dict(zip(varyList,sig))
3483        G2fil.G2Print (' Results of reflectometry data modelling fit:')
3484        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
3485        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
3486        SetModelParms()
3487        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
3488    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
3489        G2fil.G2Print (Msg)
3490        return False,0,0,0,0,0,0,Msg
3491       
3492def makeSLDprofile(data,Substances):
3493   
3494    sq2 = np.sqrt(2.)
3495    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3496    Nlayers = len(laySeq)
3497    laySeq = np.array(laySeq,dtype=int)
3498    interfaces = np.zeros(Nlayers)
3499    rho = np.zeros(Nlayers)
3500    sigma = np.zeros(Nlayers)
3501    name = data['Layers'][0]['Name']
3502    thick = 0.
3503    for ilay,lay in enumerate(laySeq):
3504        layer = data['Layers'][lay]
3505        name = layer['Name']
3506        if 'Thick' in layer:
3507            thick += layer['Thick'][0]
3508            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
3509        if 'Rough' in layer:
3510            sigma[ilay] = max(0.001,layer['Rough'][0])
3511        if name != 'vacuum':
3512            if name == 'unit scatter':
3513                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
3514            else:
3515                rrho = Substances[name]['Scatt density']
3516                irho = Substances[name]['XImag density']
3517                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
3518        if 'Mag SLD' in layer:
3519            rho[ilay] += layer['Mag SLD'][0]
3520    name = data['Layers'][-1]['Name']
3521    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
3522    xr = np.flipud(x)
3523    interfaces[-1] = x[-1]
3524    y = np.ones_like(x)*rho[0]
3525    iBeg = 0
3526    for ilayer in range(Nlayers-1):
3527        delt = rho[ilayer+1]-rho[ilayer]
3528        iPos = np.searchsorted(x,interfaces[ilayer])
3529        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
3530        iBeg = iPos
3531    return x,xr,y   
3532
3533def REFDModelFxn(Profile,Inst,Limits,Substances,data):
3534   
3535    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3536    Qmin = Limits[1][0]
3537    Qmax = Limits[1][1]
3538    iBeg = np.searchsorted(Q,Qmin)
3539    iFin = np.searchsorted(Q,Qmax)+1    #include last point
3540    Ib[:] = data['FltBack'][0]
3541    Ic[:] = 0
3542    Scale = data['Scale'][0]
3543    if data['Layer Seq'] == []:
3544        return
3545    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3546    Nlayers = len(laySeq)
3547    depth = np.zeros(Nlayers)
3548    rho = np.zeros(Nlayers)
3549    irho = np.zeros(Nlayers)
3550    sigma = np.zeros(Nlayers)
3551    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
3552        layer = data['Layers'][lay]
3553        name = layer['Name']
3554        if 'Thick' in layer:    #skips first & last layers
3555            depth[ilay] = layer['Thick'][0]
3556        if 'Rough' in layer:    #skips first layer
3557            sigma[ilay] = layer['Rough'][0]
3558        if 'unit scatter' == name:
3559            rho[ilay] = layer['DenMul'][0]
3560            irho[ilay] = layer['iDenMul'][0]
3561        else:
3562            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
3563            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
3564        if 'Mag SLD' in layer:
3565            rho[ilay] += layer['Mag SLD'][0]
3566    if data['dQ type'] == 'None':
3567        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
3568    elif 'const' in data['dQ type']:
3569        res = data['Resolution'][0]/(100.*sateln2)
3570        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
3571    else:       #dQ/Q in data
3572        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
3573    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
3574
3575def abeles(kz, depth, rho, irho=0, sigma=0):
3576    """
3577    Optical matrix form of the reflectivity calculation.
3578    O.S. Heavens, Optical Properties of Thin Solid Films
3579   
3580    Reflectometry as a function of kz for a set of slabs.
3581
3582    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\\pi\\sin(\\theta)/\\lambda`.
3583        This is :math:`\\tfrac12 Q_z`.       
3584    :param depth:  float[m] (Ang).
3585        thickness of each layer.  The thickness of the incident medium
3586        and substrate are ignored.
3587    :param rho:  float[n,k] (1e-6/Ang^2)
3588        Real scattering length density for each layer for each kz
3589    :param irho:  float[n,k] (1e-6/Ang^2)
3590        Imaginary scattering length density for each layer for each kz
3591        Note: absorption cross section mu = 2 irho/lambda for neutrons
3592    :param sigma: float[m-1] (Ang)
3593        interfacial roughness.  This is the roughness between a layer
3594        and the previous layer. The sigma array should have m-1 entries.
3595
3596    Slabs are ordered with the surface SLD at index 0 and substrate at
3597    index -1, or reversed if kz < 0.
3598    """
3599    def calc(kz, depth, rho, irho, sigma):
3600        if len(kz) == 0: return kz
3601   
3602        # Complex index of refraction is relative to the incident medium.
3603        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
3604        # in place of kz^2, and ignoring rho_o
3605        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
3606        k = kz
3607   
3608        # According to Heavens, the initial matrix should be [ 1 F; F 1],
3609        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
3610        # multiply versus some coding convenience.
3611        B11 = 1
3612        B22 = 1
3613        B21 = 0
3614        B12 = 0
3615        for i in range(0, len(depth)-1):
3616            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
3617            F = (k - k_next) / (k + k_next)
3618            F *= np.exp(-2*k*k_next*sigma[i]**2)
3619            #print "==== layer",i
3620            #print "kz:", kz
3621            #print "k:", k
3622            #print "k_next:",k_next
3623            #print "F:",F
3624            #print "rho:",rho[:,i+1]
3625            #print "irho:",irho[:,i+1]
3626            #print "d:",depth[i],"sigma:",sigma[i]
3627            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
3628            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
3629            M21 = F*M11
3630            M12 = F*M22
3631            C1 = B11*M11 + B21*M12
3632            C2 = B11*M21 + B21*M22
3633            B11 = C1
3634            B21 = C2
3635            C1 = B12*M11 + B22*M12
3636            C2 = B12*M21 + B22*M22
3637            B12 = C1
3638            B22 = C2
3639            k = k_next
3640   
3641        r = B12/B11
3642        return np.absolute(r)**2
3643
3644    if np.isscalar(kz): kz = np.asarray([kz], 'd')
3645
3646    m = len(depth)
3647
3648    # Make everything into arrays
3649    depth = np.asarray(depth,'d')
3650    rho = np.asarray(rho,'d')
3651    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
3652    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
3653
3654    # Repeat rho,irho columns as needed
3655    if len(rho.shape) == 1:
3656        rho = rho[None,:]
3657        irho = irho[None,:]
3658
3659    return calc(kz, depth, rho, irho, sigma)
3660   
3661def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
3662    y = abeles(kz, depth, rho, irho, sigma)
3663    s = dq/2.
3664    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
3665    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
3666    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
3667    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
3668    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
3669    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
3670    y *= 0.137023
3671    return y
3672       
3673def makeRefdFFT(Limits,Profile):
3674    G2fil.G2Print ('make fft')
3675    Q,Io = Profile[:2]
3676    Qmin = Limits[1][0]
3677    Qmax = Limits[1][1]
3678    iBeg = np.searchsorted(Q,Qmin)
3679    iFin = np.searchsorted(Q,Qmax)+1    #include last point
3680    Qf = np.linspace(0.,Q[iFin-1],5000)
3681    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
3682    If = QI(Qf)*Qf**4
3683    R = np.linspace(0.,5000.,5000)
3684    F = fft.rfft(If)
3685    return R,F
3686
3687   
3688#### Stacking fault simulation codes ################################################################################
3689def GetStackParms(Layers):
3690   
3691    Parms = []
3692#cell parms
3693    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
3694        Parms.append('cellA')
3695        Parms.append('cellC')
3696    else:
3697        Parms.append('cellA')
3698        Parms.append('cellB')
3699        Parms.append('cellC')
3700        if Layers['Laue'] != 'mmm':
3701            Parms.append('cellG')
3702#Transition parms
3703    for iY in range(len(Layers['Layers'])):
3704        for iX in range(len(Layers['Layers'])):
3705            Parms.append('TransP;%d;%d'%(iY,iX))
3706            Parms.append('TransX;%d;%d'%(iY,iX))
3707            Parms.append('TransY;%d;%d'%(iY,iX))
3708            Parms.append('TransZ;%d;%d'%(iY,iX))
3709    return Parms
3710
3711def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
3712    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
3713   
3714    :param dict Layers: dict with following items
3715
3716      ::
3717
3718       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
3719       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
3720       'Layers':[],'Stacking':[],'Transitions':[]}
3721       
3722    :param str ctrls: controls string to be written on DIFFaX controls.dif file
3723    :param float scale: scale factor
3724    :param dict background: background parameters
3725    :param list limits: min/max 2-theta to be calculated
3726    :param dict inst: instrument parameters dictionary
3727    :param list profile: powder pattern data
3728   
3729    Note that parameters all updated in place   
3730    '''
3731    import atmdata
3732    path = sys.path
3733    for name in path:
3734        if 'bin' in name:
3735            DIFFaX = name+'/DIFFaX.exe'
3736            G2fil.G2Print (' Execute '+DIFFaX)
3737            break
3738    # make form factor file that DIFFaX wants - atom types are GSASII style
3739    sf = open('data.sfc','w')
3740    sf.write('GSASII special form factor file for DIFFaX\n\n')
3741    atTypes = list(Layers['AtInfo'].keys())
3742    if 'H' not in atTypes:
3743        atTypes.insert(0,'H')
3744    for atType in atTypes:
3745        if atType == 'H': 
3746            blen = -.3741
3747        else:
3748            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
3749        Adat = atmdata.XrayFF[atType]
3750        text = '%4s'%(atType.ljust(4))
3751        for i in range(4):
3752            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
3753        text += '%11.6f%11.6f'%(Adat['fc'],blen)
3754        text += '%3d\n'%(Adat['Z'])
3755        sf.write(text)
3756    sf.close()
3757    #make DIFFaX control.dif file - future use GUI to set some of these flags
3758    cf = open('control.dif','w')
3759    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
3760        x0 = profile[0]
3761        iBeg = np.searchsorted(x0,limits[0])
3762        iFin = np.searchsorted(x0,limits[1])+1
3763        if iFin-iBeg > 20000:
3764            iFin = iBeg+20000
3765        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3766        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3767        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
3768    else:
3769        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3770        inst = {'Type':['XSC','XSC',]}
3771    cf.close()
3772    #make DIFFaX data file
3773    df = open('GSASII-DIFFaX.dat','w')
3774    df.write('INSTRUMENTAL\n')
3775    if 'X' in inst['Type'][0]:
3776        df.write('X-RAY\n')
3777    elif 'N' in inst['Type'][0]:
3778        df.write('NEUTRON\n')
3779    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
3780        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
3781        U = ateln2*inst['U'][1]/10000.
3782        V = ateln2*inst['V'][1]/10000.
3783        W = ateln2*inst['W'][1]/10000.
3784        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
3785        HW = np.sqrt(np.mean(HWHM))
3786    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
3787        if 'Mean' in Layers['selInst']:
3788            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
3789        elif 'Gaussian' in Layers['selInst']:
3790            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
3791        else:
3792            df.write('None\n')
3793    else:
3794        df.write('0.10\nNone\n')
3795    df.write('STRUCTURAL\n')
3796    a,b,c = Layers['Cell'][1:4]
3797    gam = Layers['Cell'][6]
3798    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
3799    laue = Layers['Laue']
3800    if laue == '2/m(ab)':
3801        laue = '2/m(1)'
3802    elif laue == '2/m(c)':
3803        laue = '2/m(2)'
3804    if 'unknown' in Layers['Laue']:
3805        df.write('%s %.3f\n'%(laue,Layers['Toler']))
3806    else:   
3807        df.write('%s\n'%(laue))
3808    df.write('%d\n'%(len(Layers['Layers'])))
3809    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
3810        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
3811    layerNames = []
3812    for layer in Layers['Layers']:
3813        layerNames.append(layer['Name'])
3814    for il,layer in enumerate(Layers['Layers']):
3815        if layer['SameAs']:
3816            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
3817            continue
3818        df.write('LAYER %d\n'%(il+1))
3819        if '-1' in layer['Symm']:
3820            df.write('CENTROSYMMETRIC\n')
3821        else:
3822            df.write('NONE\n')
3823        for ia,atom in enumerate(layer['Atoms']):
3824            [name,atype,x,y,z,frac,Uiso] = atom
3825            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
3826                frac /= 2.
3827            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
3828    df.write('STACKING\n')
3829    df.write('%s\n'%(Layers['Stacking'][0]))
3830    if 'recursive' in Layers['Stacking'][0]:
3831        df.write('%s\n'%Layers['Stacking'][1])
3832    else:
3833        if 'list' in Layers['Stacking'][1]:
3834            Slen = len(Layers['Stacking'][2])
3835            iB = 0
3836            iF = 0
3837            while True:
3838                iF += 68
3839                if iF >= Slen:
3840                    break
3841                iF = min(iF,Slen)
3842                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
3843                iB = iF
3844        else:
3845            df.write('%s\n'%Layers['Stacking'][1])   
3846    df.write('TRANSITIONS\n')
3847    for iY in range(len(Layers['Layers'])):
3848        sumPx = 0.
3849        for iX in range(len(Layers['Layers'])):
3850            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
3851            p = round(p,3)
3852            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
3853            sumPx += p
3854        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
3855            G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
3856            df.close()
3857            os.remove('data.sfc')
3858            os.remove('control.dif')
3859            os.remove('GSASII-DIFFaX.dat')
3860            return       
3861    df.close()   
3862    time0 = time.time()
3863    try:
3864        subp.call(DIFFaX)
3865    except OSError:
3866        G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
3867    G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
3868    if os.path.exists('GSASII-DIFFaX.spc'):
3869        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
3870        iFin = iBeg+Xpat.shape[1]
3871        bakType,backDict,backVary = SetBackgroundParms(background)
3872        backDict['Lam1'] = G2mth.getWave(inst)
3873        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
3874        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
3875        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
3876            rv = st.poisson(profile[3][iBeg:iFin])
3877            profile[1][iBeg:iFin] = rv.rvs()
3878            Z = np.ones_like(profile[3][iBeg:iFin])
3879            Z[1::2] *= -1
3880            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3881            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3882        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3883    #cleanup files..
3884        os.remove('GSASII-DIFFaX.spc')
3885    elif os.path.exists('GSASII-DIFFaX.sadp'):
3886        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
3887        Sadp = np.reshape(Sadp,(256,-1))
3888        Layers['Sadp']['Img'] = Sadp
3889        os.remove('GSASII-DIFFaX.sadp')
3890    os.remove('data.sfc')
3891    os.remove('control.dif')
3892    os.remove('GSASII-DIFFaX.dat')
3893   
3894def SetPWDRscan(inst,limits,profile):
3895   
3896    wave = G2mth.getMeanWave(inst)
3897    x0 = profile[0]
3898    iBeg = np.searchsorted(x0,limits[0])
3899    iFin = np.searchsorted(x0,limits[1])
3900    if iFin-iBeg > 20000:
3901        iFin = iBeg+20000
3902    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3903    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
3904    return iFin-iBeg
3905       
3906def SetStackingSF(Layers,debug):
3907# Load scattering factors into DIFFaX arrays
3908    import atmdata
3909    atTypes = Layers['AtInfo'].keys()
3910    aTypes = []
3911    for atype in atTypes:
3912        aTypes.append('%4s'%(atype.ljust(4)))
3913    SFdat = []
3914    for atType in atTypes:
3915        Adat = atmdata.XrayFF[atType]
3916        SF = np.zeros(9)
3917        SF[:8:2] = Adat['fa']
3918        SF[1:8:2] = Adat['fb']
3919        SF[8] = Adat['fc']
3920        SFdat.append(SF)
3921    SFdat = np.array(SFdat)
3922    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
3923   
3924def SetStackingClay(Layers,Type):
3925# Controls
3926    rand.seed()
3927    ranSeed = rand.randint(1,2**16-1)
3928    try:   
3929        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
3930            '6/m','6/mmm'].index(Layers['Laue'])+1
3931    except ValueError:  #for 'unknown'
3932        laueId = -1
3933    if 'SADP' in Type:
3934        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
3935        lmax = int(Layers['Sadp']['Lmax'])
3936    else:
3937        planeId = 0
3938        lmax = 0
3939# Sequences
3940    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
3941    try:
3942        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
3943    except ValueError:
3944        StkParm = -1
3945    if StkParm == 2:    #list
3946        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
3947        Nstk = len(StkSeq)
3948    else:
3949        Nstk = 1
3950        StkSeq = [0,]
3951    if StkParm == -1:
3952        StkParm = int(Layers['Stacking'][1])
3953    Wdth = Layers['Width'][0]
3954    mult = 1
3955    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
3956    LaueSym = Layers['Laue'].ljust(12)
3957    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
3958    return laueId,controls
3959   
3960def SetCellAtoms(Layers):
3961    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
3962# atoms in layers
3963    atTypes = list(Layers['AtInfo'].keys())
3964    AtomXOU = []
3965    AtomTp = []
3966    LayerSymm = []
3967    LayerNum = []
3968    layerNames = []
3969    Natm = 0
3970    Nuniq = 0
3971    for layer in Layers['Layers']:
3972        layerNames.append(layer['Name'])
3973    for il,layer in enumerate(Layers['Layers']):
3974        if layer['SameAs']:
3975            LayerNum.append(layerNames.index(layer['SameAs'])+1)
3976            continue
3977        else:
3978            LayerNum.append(il+1)
3979            Nuniq += 1
3980        if '-1' in layer['Symm']:
3981            LayerSymm.append(1)
3982        else:
3983            LayerSymm.append(0)
3984        for ia,atom in enumerate(layer['Atoms']):
3985            [name,atype,x,y,z,frac,Uiso] = atom
3986            Natm += 1
3987            AtomTp.append('%4s'%(atype.ljust(4)))
3988            Ta = atTypes.index(atype)+1
3989            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
3990    AtomXOU = np.array(AtomXOU)
3991    Nlayers = len(layerNames)
3992    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
3993    return Nlayers
3994   
3995def SetStackingTrans(Layers,Nlayers):
3996# Transitions
3997    TransX = []
3998    TransP = []
3999    for Ytrans in Layers['Transitions']:
4000        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
4001        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
4002    TransP = np.array(TransP,dtype='float').T
4003    TransX = np.array(TransX,dtype='float')
4004#    GSASIIpath.IPyBreak()
4005    pyx.pygettrans(Nlayers,TransP,TransX)
4006   
4007def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
4008# Scattering factors
4009    SetStackingSF(Layers,debug)
4010# Controls & sequences
4011    laueId,controls = SetStackingClay(Layers,'PWDR')
4012# cell & atoms
4013    Nlayers = SetCellAtoms(Layers)
4014    Volume = Layers['Cell'][7]   
4015# Transitions
4016    SetStackingTrans(Layers,Nlayers)
4017# PWDR scan
4018    Nsteps = SetPWDRscan(inst,limits,profile)
4019# result as Spec
4020    x0 = profile[0]
4021    profile[3] = np.zeros(len(profile[0]))
4022    profile[4] = np.zeros(len(profile[0]))
4023    profile[5] = np.zeros(len(profile[0]))
4024    iBeg = np.searchsorted(x0,limits[0])
4025    iFin = np.searchsorted(x0,limits[1])+1
4026    if iFin-iBeg > 20000:
4027        iFin = iBeg+20000
4028    Nspec = 20001       
4029    spec = np.zeros(Nspec,dtype='double')   
4030    time0 = time.time()
4031    pyx.pygetspc(controls,Nspec,spec)
4032    G2fil.G2Print (' GETSPC time = %.2fs'%(time.time()-time0))
4033    time0 = time.time()
4034    U = ateln2*inst['U'][1]/10000.
4035    V = ateln2*inst['V'][1]/10000.
4036    W = ateln2*inst['W'][1]/10000.
4037    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
4038    HW = np.sqrt(np.mean(HWHM))
4039    BrdSpec = np.zeros(Nsteps)
4040    if 'Mean' in Layers['selInst']:
4041        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
4042    elif 'Gaussian' in Layers['selInst']:
4043        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
4044    else:
4045        BrdSpec = spec[:Nsteps]
4046    BrdSpec /= Volume
4047    iFin = iBeg+Nsteps
4048    bakType,backDict,backVary = SetBackgroundParms(background)
4049    backDict['Lam1'] = G2mth.getWave(inst)
4050    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
4051    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
4052    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
4053        try:
4054            rv = st.poisson(profile[3][iBeg:iFin])
4055            profile[1][iBeg:iFin] = rv.rvs()
4056        except ValueError:
4057            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
4058        Z = np.ones_like(profile[3][iBeg:iFin])
4059        Z[1::2] *= -1
4060        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
4061        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
4062    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
4063    G2fil.G2Print (' Broadening time = %.2fs'%(time.time()-time0))
4064   
4065def CalcStackingSADP(Layers,debug):
4066   
4067# Scattering factors
4068    SetStackingSF(Layers,debug)
4069# Controls & sequences
4070    laueId,controls = SetStackingClay(Layers,'SADP')
4071# cell & atoms
4072    Nlayers = SetCellAtoms(Layers)   
4073# Transitions
4074    SetStackingTrans(Layers,Nlayers)
4075# result as Sadp
4076    Nspec = 20001       
4077    spec = np.zeros(Nspec,dtype='double')   
4078    time0 = time.time()
4079    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
4080    Sapd = np.zeros((256,256))
4081    iB = 0
4082    for i in range(hkLim):
4083        iF = iB+Nblk
4084        p1 = 127+int(i*Incr)
4085        p2 = 128-int(i*Incr)
4086        if Nblk == 128:
4087            if i:
4088                Sapd[128:,p1] = spec[iB:iF]
4089                Sapd[:128,p1] = spec[iF:iB:-1]
4090            Sapd[128:,p2] = spec[iB:iF]
4091            Sapd[:128,p2] = spec[iF:iB:-1]
4092        else:
4093            if i:
4094                Sapd[:,p1] = spec[iB:iF]
4095            Sapd[:,p2] = spec[iB:iF]
4096        iB += Nblk
4097    Layers['Sadp']['Img'] = Sapd
4098    G2fil.G2Print (' GETSAD time = %.2fs'%(time.time()-time0))
4099   
4100#### Maximum Entropy Method - Dysnomia ###############################################################################
4101def makePRFfile(data,MEMtype):
4102    ''' makes Dysnomia .prf control file from Dysnomia GUI controls
4103   
4104    :param dict data: GSAS-II phase data
4105    :param int MEMtype: 1 for neutron data with negative scattering lengths
4106                        0 otherwise
4107    :returns str: name of Dysnomia control file
4108    '''
4109
4110    generalData = data['General']
4111    pName = generalData['Name'].replace(' ','_')
4112    DysData = data['Dysnomia']
4113    prfName = pName+'.prf'
4114    prf = open(prfName,'w')
4115    prf.write('$PREFERENCES\n')
4116    prf.write(pName+'.mem\n') #or .fos?
4117    prf.write(pName+'.out\n')
4118    prf.write(pName+'.pgrid\n')
4119    prf.write(pName+'.fba\n')
4120    prf.write(pName+'_eps.raw\n')
4121    prf.write('%d\n'%MEMtype)
4122    if DysData['DenStart'] == 'uniform':
4123        prf.write('0\n')
4124    else:
4125        prf.write('1\n')
4126    if DysData['Optimize'] == 'ZSPA':
4127        prf.write('0\n')
4128    else:
4129        prf.write('1\n')
4130    prf.write('1\n')
4131    if DysData['Lagrange'][0] == 'user':
4132        prf.write('0\n')
4133    else:
4134        prf.write('1\n')
4135    prf.write('%.4f %d\n'%(DysData['Lagrange'][1],DysData['wt pwr']))
4136    prf.write('%.3f\n'%DysData['Lagrange'][2])
4137    prf.write('%.2f\n'%DysData['E_factor'])
4138    prf.write('1\n')
4139    prf.write('0\n')
4140    prf.write('%d\n'%DysData['Ncyc'])
4141    prf.write('1\n')
4142    prf.write('1 0 0 0 0 0 0 0\n')
4143    if DysData['prior'] == 'uniform':
4144        prf.write('0\n')
4145    else:
4146        prf.write('1\n')
4147    prf.close()
4148    return prfName
4149
4150def makeMEMfile(data,reflData,MEMtype,DYSNOMIA):
4151    ''' make Dysnomia .mem file of reflection data, etc.
4152
4153    :param dict data: GSAS-II phase data
4154    :param list reflData: GSAS-II reflection data
4155    :param int MEMtype: 1 for neutron data with negative scattering lengths
4156                        0 otherwise
4157    :param str DYSNOMIA: path to dysnomia.exe
4158    '''
4159   
4160    DysData = data['Dysnomia']
4161    generalData = data['General']
4162    cell = generalData['Cell'][1:7]
4163    A = G2lat.cell2A(cell)
4164    SGData = generalData['SGData']
4165    pName = generalData['Name'].replace(' ','_')
4166    memName = pName+'.mem'
4167    Map = generalData['Map']
4168    Type = Map['Type']
4169    UseList = Map['RefList']
4170    mem = open(memName,'w')
4171    mem.write('%s\n'%(generalData['Name']+' from '+UseList[0]))
4172    a,b,c,alp,bet,gam = cell
4173    mem.write('%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n'%(a,b,c,alp,bet,gam))
4174    mem.write('      0.0000000      0.0000000     -1    0    0    0     P\n')   #dummy PO stuff
4175    SGSym = generalData['SGData']['SpGrp']
4176    try:
4177        SGId = G2spc.spgbyNum.index(SGSym)
4178    except ValueError:
4179        return False
4180    org = 1
4181    if SGSym in G2spc.spg2origins:
4182        org = 2
4183    mapsize = Map['rho'].shape
4184    sumZ = 0.
4185    sumpos = 0.
4186    sumneg = 0.
4187    mem.write('%5d%5d%5d%5d%5d\n'%(SGId,org,mapsize[0],mapsize[1],mapsize[2]))
4188    for atm in generalData['NoAtoms']:
4189        Nat = generalData['NoAtoms'][atm]
4190        AtInfo = G2elem.GetAtomInfo(atm)
4191        sumZ += Nat*AtInfo['Z']
4192        isotope = generalData['Isotope'][atm]
4193        blen = generalData['Isotopes'][atm][isotope]['SL'][0]
4194        if blen < 0.:
4195            sumneg += blen*Nat
4196        else:
4197            sumpos += blen*Nat
4198    if 'X' in Type:
4199        mem.write('%10.2f  0.001\n'%sumZ)
4200    elif 'N' in Type and MEMtype:
4201        mem.write('%10.3f%10.3f 0.001\n'%(sumpos,sumneg))
4202    else:
4203        mem.write('%10.3f 0.001\n'%sumpos)
4204       
4205    dmin = DysData['MEMdmin']
4206    TOFlam = 2.0*dmin*npsind(80.0)
4207    refSet = G2lat.GenHLaue(dmin,SGData,A)      #list of h,k,l,d
4208    refDict = {'%d %d %d'%(ref[0],ref[1],ref[2]):ref for ref in refSet}
4209       
4210    refs = []
4211    prevpos = 0.
4212    for ref in reflData:
4213        if ref[3] < 0:
4214            continue
4215        if 'T' in Type:
4216            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,x,x,x,prfo = ref[:16]
4217            s = np.sqrt(max(sig,0.0001))   #var -> sig in deg
4218            FWHM = getgamFW(gam,s)
4219            if dsp < dmin:
4220                continue
4221            theta = npasind(TOFlam/(2.*dsp))
4222            FWHM *= nptand(theta)/pos
4223            pos = 2.*theta
4224        else:
4225            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,prfo = ref[:13]
4226            g = gam/100.    #centideg -> deg
4227            s = np.sqrt(max(sig,0.0001))/100.   #var -> sig in deg
4228            FWHM = getgamFW(g,s)
4229        delt = pos-prevpos
4230        refs.append([h,k,l,mult,pos,FWHM,Fobs,phase,delt])
4231        prevpos = pos
4232           
4233    ovlp = DysData['overlap']
4234    refs1 = []
4235    refs2 = []
4236    nref2 = 0
4237    iref = 0
4238    Nref = len(refs)
4239    start = False
4240    while iref < Nref-1:
4241        if refs[iref+1][-1] < ovlp*refs[iref][5]:
4242            if refs[iref][-1] > ovlp*refs[iref][5]:
4243                refs2.append([])
4244                start = True
4245            if nref2 == len(refs2):
4246                refs2.append([])
4247            refs2[nref2].append(refs[iref])
4248        else:
4249            if start:
4250                refs2[nref2].append(refs[iref])
4251                start = False
4252                nref2 += 1
4253            else:
4254                refs1.append(refs[iref])
4255        iref += 1
4256    if start:
4257        refs2[nref2].append(refs[iref])
4258    else:
4259        refs1.append(refs[iref])
4260   
4261    mem.write('%5d\n'%len(refs1))
4262    for ref in refs1:
4263        h,k,l = ref[:3]
4264        hkl = '%d %d %d'%(h,k,l)
4265        if hkl in refDict:
4266            del refDict[hkl]
4267        Fobs = np.sqrt(ref[6])
4268        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)))
4269    while True and nref2:
4270        if not len(refs2[-1]):
4271            del refs2[-1]
4272        else:
4273            break
4274    mem.write('%5d\n'%len(refs2))
4275    for iref2,ref2 in enumerate(refs2):
4276        mem.write('#%5d\n'%iref2)
4277        mem.write('%5d\n'%len(ref2))
4278        Gsum = 0.
4279        Msum = 0
4280        for ref in ref2:
4281            Gsum += ref[6]*ref[3]
4282            Msum += ref[3]
4283        G = np.sqrt(Gsum/Msum)
4284        h,k,l = ref2[0][:3]
4285        hkl = '%d %d %d'%(h,k,l)
4286        if hkl in refDict:
4287            del refDict[hkl]
4288        mem.write('%5d%5d%5d%10.3f%10.3f%5d\n'%(h,k,l,G,max(0.01*G,0.1),ref2[0][3]))
4289        for ref in ref2[1:]:
4290            h,k,l,m = ref[:4]
4291            mem.write('%5d%5d%5d%5d\n'%(h,k,l,m))
4292            hkl = '%d %d %d'%(h,k,l)
4293            if hkl in refDict:
4294                del refDict[hkl]
4295    if len(refDict):
4296        mem.write('%d\n'%len(refDict))
4297        for hkl in list(refDict.keys()):
4298            h,k,l = refDict[hkl][:3]
4299            mem.write('%5d%5d%5d\n'%(h,k,l))
4300    else:
4301        mem.write('0\n')
4302    mem.close()
4303    return True
4304
4305def MEMupdateReflData(prfName,data,reflData):
4306    ''' Update reflection data with new Fosq, phase result from Dysnomia
4307
4308    :param str prfName: phase.mem file name
4309    :param list reflData: GSAS-II reflection data
4310    '''
4311   
4312    generalData = data['General']
4313    Map = generalData['Map']
4314    Type = Map['Type']
4315    cell = generalData['Cell'][1:7]
4316    A = G2lat.cell2A(cell)
4317    reflDict = {}
4318    newRefs = []
4319    for iref,ref in enumerate(reflData):
4320        if ref[3] > 0:
4321            newRefs.append(ref)
4322            reflDict[hash('%5d%5d%5d'%(ref[0],ref[1],ref[2]))] = iref
4323    fbaName = os.path.splitext(prfName)[0]+'.fba'
4324    if os.path.isfile(fbaName):
4325        fba = open(fbaName,'r')
4326    else:
4327        return False
4328    fba.readline()
4329    Nref = int(fba.readline()[:-1])
4330    fbalines = fba.readlines()
4331    for line in fbalines[:Nref]:
4332        info = line.split()
4333        h = int(info[0])
4334        k = int(info[1])
4335        l = int(info[2])
4336        FoR = float(info[3])
4337        FoI = float(info[4])
4338        Fosq = FoR**2+FoI**2
4339        phase = npatan2d(FoI,FoR)
4340        try:
4341            refId = reflDict[hash('%5d%5d%5d'%(h,k,l))]
4342        except KeyError:    #added reflections at end skipped
4343            d = float(1/np.sqrt(G2lat.calc_rDsq([h,k,l],A)))
4344            if 'T' in Type:
4345                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])
4346            else:
4347                newRefs.append([h,k,l,-1,d,0.,0.01,1.0,Fosq,Fosq,phase,1.0,1.0,1.0,1.0])
4348            continue
4349        newRefs[refId][8] = Fosq
4350        newRefs[refId][10] = phase
4351    newRefs = np.array(newRefs)
4352    return True,newRefs
4353   
4354#### testing data
4355NeedTestData = True
4356def TestData():
4357    'needs a doc string'
4358#    global NeedTestData
4359    global bakType
4360    bakType = 'chebyschev'
4361    global xdata
4362    xdata = np.linspace(4.0,40.0,36000)
4363    global parmDict0
4364    parmDict0 = {
4365        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
4366        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
4367        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
4368        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
4369        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'Z':0.0,'SH/L':0.002,
4370        'Back0':5.384,'Back1':-0.015,'Back2':.004,
4371        }
4372    global parmDict1
4373    parmDict1 = {
4374        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0