source: trunk/GSASIIpwd.py @ 4347

Last change on this file since 4347 was 4347, checked in by vondreele, 20 months ago

restrict G(R) plots to Rmax - sphere enclosed by big box model in RMCProfile

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