source: trunk/GSASIIpwd.py

Last change on this file was 5037, checked in by vondreele, 4 days ago

add use GPU option for running RMCProfile
remove scrolled sub panels in RMCProfile GUI - didn't work & interfered with normal use

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