source: trunk/GSASIIpwd.py @ 5029

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

upper limit on crystallite size = 10 (it was 4)
consolidate import subprocess as subp (remove the extra imports)
minor change to MakeInst? for bug in RMCProfile 6.7.9

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