source: trunk/GSASIIpwd.py @ 4953

Last change on this file since 4953 was 4953, checked in by toby, 4 months ago

fix fullrmc rmax bug

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