source: trunk/GSASIIpwd.py @ 4952

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

switch to native fullrmc input modes rather than convert files

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