source: trunk/GSASIIpwd.py @ 4951

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

working fullrmc implementation. Much more to do (see TODO in phsGUI), but good start. Needs external fullrmc Python from Bachir until 5.0 is released.

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