source: trunk/GSASIIpwd.py @ 4998

Last change on this file since 4998 was 4998, checked in by toby, 3 months ago

partial fix of magnetic constraints (still need to deal w/unvaried items; fullrmc minor stuff

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