source: trunk/GSASIIpwd.py @ 4334

Last change on this file since 4334 was 4334, checked in by vondreele, 3 years ago

make new f(q) & g(r) files suitable for RMCProfile. The g(r) == PDF.
Fix 2X issue in xray scale factor for .bragg files

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 146.0 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-01 07:19:50 +0000 (Sun, 01 Mar 2020) $
10# $Author: vondreele $
11# $Revision: 4334 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 4334 2020-03-01 07:19:50Z 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: 4334 $")
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]/2.
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]/2.
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]+Xsb,inst[prms[13]][1]+Ysb,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    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False)
2214    Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2215    Natm = np.count_nonzero(Natm-1)
2216    Atoms = newPhase['Atoms']
2217    Satoms = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms,5),4),3)
2218    Datoms = [[atom for atom in Satoms if atom[0] in dup] for dup in Dups]
2219    Natoms = []
2220    reset = False
2221    for idup,dup in enumerate(Dups):
2222        ldup = len(dup)
2223        datoms = Datoms[idup]
2224        natm = len(datoms)
2225        i = 0
2226        while i < natm:
2227            atoms = datoms[i:i+ldup]
2228            try:
2229                atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2230                Natoms.append(atom)
2231            except IndexError:      #what about vacancies?
2232                if 'Va' not in Atseq:
2233                    reset = True
2234                    Atseq.append('Va')
2235                    RMCPdict['aTypes']['Va'] = 0.0
2236                atom = atoms[0]
2237                atom[1] = 'Va'
2238                Natoms.append(atom)
2239            i += ldup
2240    NAtype = np.zeros(len(Atseq))
2241    for atom in Natoms:
2242        NAtype[Atseq.index(atom[1])] += 1
2243    NAstr = ['%6d'%i for i in NAtype]
2244    Cell = newPhase['General']['Cell'][1:7]
2245    if os.path.exists(Name+'.his6f'):
2246        os.remove(Name+'.his6f')
2247    if os.path.exists(Name+'.neigh'):
2248        os.remove(Name+'.neigh')
2249    fname = Name+'.rmc6f'
2250    fl = open(fname,'w')
2251    fl.write('(Version 6f format configuration file)\n')
2252    for item in Meta:
2253        fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2254    fl.write('Atom types present:                 %s\n'%'    '.join(Atseq))
2255    fl.write('Number of each atom type:       %s\n'%''.join(NAstr))
2256    fl.write('Number of atoms:                %d\n'%len(Natoms))
2257    fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2258    fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2259            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2260    A,B = G2lat.cell2AB(Cell,True)
2261    fl.write('Lattice vectors (Ang):\n')   
2262    for i in [0,1,2]:
2263        fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2264    fl.write('Atoms (fractional coordinates):\n')
2265    nat = 0
2266    for atm in Atseq:
2267        for iat,atom in enumerate(Natoms):
2268            if atom[1] == atm:
2269                nat += 1
2270                atcode = Atcodes[iat].split(':')
2271                cell = [0,0,0]
2272                if '+' in atcode[1]:
2273                    cell = eval(atcode[1].split('+')[1])
2274                fl.write('%6d%4s  [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(       
2275                        nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2276    fl.close()
2277    return fname,reset
2278
2279def MakeBragg(G2frame,Name,Phase,PWId):
2280    PWDdata = G2frame.GetPWDRdatafromTree(PWId)
2281    generalData = Phase['General']
2282    Vol = generalData['Cell'][7]
2283    Data = PWDdata['Data']
2284    Inst = PWDdata['Instrument Parameters'][0]
2285    Bank = int(Inst['Bank'][1])
2286    Sample = PWDdata['Sample Parameters']
2287    Scale = Sample['Scale'][0]
2288    Limits = PWDdata['Limits'][1]
2289    Ibeg = np.searchsorted(Data[0],Limits[0])
2290    Ifin = np.searchsorted(Data[0],Limits[1])+1
2291    fname = Name+'.bragg'
2292    fl = open(fname,'w')
2293    fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol))
2294    if 'T' in Inst['Type'][0]:
2295        fl.write('%12s%12s\n'%('   TOF,ms','  I(obs)'))
2296        for i in range(Ibeg,Ifin-1):
2297            fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2298    else:
2299        fl.write('%12s%12s\n'%('   2-theta, deg','  I(obs)'))
2300        for i in range(Ibeg,Ifin-1):
2301            fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))       
2302    fl.close()
2303    return fname
2304
2305def MakeRMCPdat(G2frame,Name,Phase,RMCPdict,PWId):
2306    Meta = RMCPdict['metadata']
2307    Times = RMCPdict['runTimes']
2308    Atseq = RMCPdict['atSeq']
2309    Atypes = RMCPdict['aTypes']
2310    atPairs = RMCPdict['Pairs']
2311    Files = RMCPdict['files']
2312    BraggWt = RMCPdict['histogram'][1]
2313    PWDdata = G2frame.GetPWDRdatafromTree(PWId)
2314    inst = PWDdata['Instrument Parameters'][0]
2315    refList = PWDdata['Reflection Lists'][Name]['RefList']
2316    dMin = refList[-1][4]
2317    gsasType = 'xray2'
2318    if 'T' in inst['Type'][1]:
2319        gsasType = 'gsas3'
2320    elif 'X' in inst['Type'][1]:
2321        XFF = G2elem.GetFFtable(Atseq)
2322        Xfl = open(Name+'.xray','w')
2323        for atm in Atseq:
2324            fa = XFF[atm]['fa']
2325            fb = XFF[atm]['fb']
2326            fc = XFF[atm]['fc']
2327            Xfl.write('%2s  %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%(
2328                    atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc))
2329        Xfl.close()
2330    lenA = len(Atseq)
2331    Pairs = []
2332    for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]:
2333        Pairs += pair
2334    pairMin = [atPairs[pair]for pair in Pairs if pair in atPairs]
2335    maxMoves = [Atypes[atm] for atm in Atseq if atm in Atypes]
2336    fname = Name+'.dat'
2337    fl = open(fname,'w')
2338    fl.write(' %% Hand edit the following as needed\n')
2339    fl.write('TITLE :: '+Name+'\n')
2340    fl.write('MATERIAL :: '+Meta['material']+'\n')
2341    fl.write('PHASE :: '+Meta['phase']+'\n')
2342    fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n')
2343    fl.write('INVESTIGATOR :: '+Meta['owner']+'\n')
2344    minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin])
2345    minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin])
2346    maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin])
2347    fl.write('MINIMUM_DISTANCES ::   %s  Angstrom\n'%minHD)
2348    maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves])
2349    fl.write('MAXIMUM_MOVES ::   %s Angstrom\n'%maxMv)
2350    fl.write('R_SPACING ::  0.0200 Angstrom\n')
2351    fl.write('PRINT_PERIOD :: 100\n')
2352    fl.write('TIME_LIMIT ::     %.2f MINUTES\n'%Times[0])
2353    fl.write('SAVE_PERIOD ::    %.2f MINUTES\n'%Times[1])
2354    fl.write('\n')
2355    fl.write('ATOMS :: '+' '.join(Atseq)+'\n')
2356    fl.write('\n')
2357    fl.write('FLAGS ::\n')
2358    fl.write('  > NO_MOVEOUT\n')
2359    fl.write('  > NO_SAVE_CONFIGURATIONS\n')
2360    fl.write('  > NO_RESOLUTION_CONVOLUTION\n')
2361    fl.write('\n')
2362    fl.write('INPUT_CONFIGURATION_FORMAT ::  rmc6f\n')
2363    fl.write('SAVE_CONFIGURATION_FORMAT  ::  rmc6f\n')
2364    fl.write('IGNORE_HISTORY_FILE ::\n')
2365    fl.write('\n')
2366    fl.write('DISTANCE_WINDOW ::\n')
2367    fl.write('  > MNDIST :: %s\n'%minD)
2368    fl.write('  > MXDIST :: %s\n'%maxD)
2369    if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']):
2370        fl.write('\n')
2371        fl.write('POTENTIALS ::\n')
2372        fl.write('  > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.'])
2373        fl.write('  > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n')
2374        if len(RMCPdict['Potentials']['Stretch']):
2375            fl.write('  > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search'])
2376            for bond in RMCPdict['Potentials']['Stretch']:
2377                fl.write('  > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2]))       
2378        if len(RMCPdict['Potentials']['Angles']):
2379            fl.write('  > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
2380            for angle in RMCPdict['Potentials']['Angles']:
2381                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f Ang %.2f Ang\n'%
2382                    (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
2383    if RMCPdict['useBVS']:
2384        fl.write('BVS ::\n')
2385        fl.write('  > ATOM :: '+' '.join(Atseq)+'\n')
2386        fl.write('  > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2387        oxid = []
2388        for val in RMCPdict['Oxid']:
2389            if len(val) == 3:
2390                oxid.append(val[0][1:])
2391            else:
2392                oxid.append(val[0][2:])
2393        fl.write('  > OXID :: %s\n'%' '.join(oxid))
2394        fl.write('  > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']]))
2395        fl.write('  > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']]))
2396        fl.write('  > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][3] for bvs in RMCPdict['BVS']]))       
2397        fl.write('  > SAVE :: 100000\n')
2398        fl.write('  > UPDATE :: 100000\n')
2399        if len(RMCPdict['Swap']):
2400            fl.write('\n')
2401            fl.write('SWAP_MULTI ::\n')
2402            for swap in RMCPdict['Swap']:
2403                try:
2404                    at1 = Atseq.index(swap[0])
2405                    at2 = Atseq.index(swap[1])
2406                except ValueError:
2407                    break
2408                fl.write('  > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2]))
2409       
2410    if len(RMCPdict['FxCN']):
2411        fl.write('FIXED_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['FxCN']))       
2412        for ifx,fxcn in enumerate(RMCPdict['FxCN']):
2413            try:
2414                at1 = Atseq.index(fxcn[0])
2415                at2 = Atseq.index(fxcn[1])
2416            except ValueError:
2417                break
2418            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]))
2419    if len(RMCPdict['AveCN']):
2420        fl.write('AVERAGE_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['AveCN']))
2421        for iav,avcn in enumerate(RMCPdict['AveCN']):
2422            try:
2423                at1 = Atseq.index(avcn[0])
2424                at2 = Atseq.index(avcn[1])
2425            except ValueError:
2426                break
2427            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]))
2428    for File in Files:
2429        if Files[File][0]:
2430            if 'Xray' in File and 'F(Q)' in File:
2431                fqdata = open(Files[File][0],'r')
2432                lines = int(fqdata.readline()[:-1])
2433            fl.write('\n')
2434            fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_'))
2435            fl.write('  > FILENAME :: %s\n'%Files[File][0])
2436            fl.write('  > DATA_TYPE :: %s\n'%Files[File][2])
2437            fl.write('  > FIT_TYPE :: %s\n'%Files[File][2])
2438            if 'Xray' not in File:
2439                fl.write('  > START_POINT :: 1\n')
2440                fl.write('  > END_POINT :: 3000\n')
2441                fl.write('  > WEIGHT :: %.4f\n'%Files[File][1])
2442            fl.write('  > CONSTANT_OFFSET 0.000\n')
2443            fl.write('  > NO_FITTED_OFFSET\n')
2444            if RMCPdict['FitScale']:
2445                fl.write('  > FITTED_SCALE\n')
2446            else:
2447                fl.write('  > NO_FITTED_SCALE\n')
2448            if Files[File][3] !='RMC':
2449                fl.write('  > %s\n'%Files[File][3])
2450            if 'reciprocal' in File:
2451                fl.write('  > CONVOLVE ::\n')
2452                if 'Xray' in File:
2453                    fl.write('  > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines)
2454                    fl.write('  > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1]))
2455                    fl.write('  > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2))
2456                    fl.write('  > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,Files[File][1]))
2457    fl.write('\n')
2458    fl.write('BRAGG ::\n')
2459    fl.write('  > BRAGG_SHAPE :: %s\n'%gsasType)
2460    fl.write('  > RECALCUATE\n')
2461    fl.write('  > DMIN :: %.2f\n'%(dMin-0.02))
2462    fl.write('  > WEIGHT :: %10.3f\n'%BraggWt)
2463    fl.write('\n')
2464    fl.write('END  ::\n')
2465    fl.close()
2466    return fname   
2467
2468def MakePDB(G2frame,Name,Phase,Atseq,Supercell):
2469    generalData = Phase['General']
2470    Cell = generalData['Cell'][1:7]
2471    Trans = np.eye(3)*np.array(Supercell)
2472    newPhase = copy.deepcopy(Phase)
2473    newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2474    newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans.T)
2475    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False)
2476    Atoms = newPhase['Atoms']
2477    Cell = newPhase['General']['Cell'][1:7]
2478    A,B = G2lat. cell2AB(Cell)
2479    fname = Name+'.pdb'
2480    fl = open(fname,'w')
2481    fl.write('REMARK    this file is generated using GSASII\n')
2482    fl.write('CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n'%(
2483            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2484    fl.write('ORIGX1      1.000000  0.000000  0.000000        0.00000\n')
2485    fl.write('ORIGX2      0.000000  1.000000  0.000000        0.00000\n')
2486    fl.write('ORIGX3      0.000000  0.000000  1.000000        0.00000\n')
2487
2488    Natm = np.core.defchararray.count(np.array(Atcodes),'+')
2489    Natm = np.count_nonzero(Natm-1)
2490    nat = 0
2491    for atm in Atseq:
2492        for iat,atom in enumerate(Atoms):
2493            if atom[1] == atm:
2494                nat += 1
2495                XYZ = np.inner(A,np.array(atom[3:6])-0.5)    #shift origin to middle & make Cartesian
2496#ATOM      1 Ni   RMC     1     -22.113 -22.113 -22.113  1.00  0.00          ni                     
2497                fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %-2s\n'%(       
2498                        nat,atom[0],nat,XYZ[0],XYZ[1],XYZ[2],atom[1]))
2499    fl.close()
2500    return fname
2501   
2502################################################################################
2503#### Reflectometry calculations
2504################################################################################
2505
2506def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
2507    G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
2508   
2509    class RandomDisplacementBounds(object):
2510        """random displacement with bounds"""
2511        def __init__(self, xmin, xmax, stepsize=0.5):
2512            self.xmin = xmin
2513            self.xmax = xmax
2514            self.stepsize = stepsize
2515   
2516        def __call__(self, x):
2517            """take a random step but ensure the new position is within the bounds"""
2518            while True:
2519                # this could be done in a much more clever way, but it will work for example purposes
2520                steps = self.xmax-self.xmin
2521                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
2522                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
2523                    break
2524            return xnew
2525   
2526    def GetModelParms():
2527        parmDict = {}
2528        varyList = []
2529        values = []
2530        bounds = []
2531        parmDict['dQ type'] = data['dQ type']
2532        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
2533        for parm in ['Scale','FltBack']:
2534            parmDict[parm] = data[parm][0]
2535            if data[parm][1]:
2536                varyList.append(parm)
2537                values.append(data[parm][0])
2538                bounds.append(Bounds[parm])
2539        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
2540        parmDict['nLayers'] = len(parmDict['Layer Seq'])
2541        for ilay,layer in enumerate(data['Layers']):
2542            name = layer['Name']
2543            cid = str(ilay)+';'
2544            parmDict[cid+'Name'] = name
2545            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2546                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
2547                if layer.get(parm,[0.,False])[1]:
2548                    varyList.append(cid+parm)
2549                    value = layer[parm][0]
2550                    values.append(value)
2551                    if value:
2552                        bound = [value*Bfac,value/Bfac]
2553                    else:
2554                        bound = [0.,10.]
2555                    bounds.append(bound)
2556            if name not in ['vacuum','unit scatter']:
2557                parmDict[cid+'rho'] = Substances[name]['Scatt density']
2558                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
2559        return parmDict,varyList,values,bounds
2560   
2561    def SetModelParms():
2562        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
2563        if 'Scale' in varyList:
2564            data['Scale'][0] = parmDict['Scale']
2565            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
2566        G2fil.G2Print (line)
2567        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
2568        if 'FltBack' in varyList:
2569            data['FltBack'][0] = parmDict['FltBack']
2570            line += ' esd: %15.3g'%(sigDict['FltBack'])
2571        G2fil.G2Print (line)
2572        for ilay,layer in enumerate(data['Layers']):
2573            name = layer['Name']
2574            G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
2575            cid = str(ilay)+';'
2576            line = ' '
2577            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
2578            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
2579            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2580                if parm in layer:
2581                    if parm == 'Rough':
2582                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
2583                    else:
2584                        layer[parm][0] = parmDict[cid+parm]
2585                    line += ' %s: %.3f'%(parm,layer[parm][0])
2586                    if cid+parm in varyList:
2587                        line += ' esd: %.3g'%(sigDict[cid+parm])
2588            G2fil.G2Print (line)
2589            G2fil.G2Print (line2)
2590   
2591    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2592        parmDict.update(zip(varyList,values))
2593        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2594        return M
2595   
2596    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2597        parmDict.update(zip(varyList,values))
2598        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2599        return np.sum(M**2)
2600   
2601    def getREFD(Q,Qsig,parmDict):
2602        Ic = np.ones_like(Q)*parmDict['FltBack']
2603        Scale = parmDict['Scale']
2604        Nlayers = parmDict['nLayers']
2605        Res = parmDict['Res']
2606        depth = np.zeros(Nlayers)
2607        rho = np.zeros(Nlayers)
2608        irho = np.zeros(Nlayers)
2609        sigma = np.zeros(Nlayers)
2610        for ilay,lay in enumerate(parmDict['Layer Seq']):
2611            cid = str(lay)+';'
2612            depth[ilay] = parmDict[cid+'Thick']
2613            sigma[ilay] = parmDict[cid+'Rough']
2614            if parmDict[cid+'Name'] == u'unit scatter':
2615                rho[ilay] = parmDict[cid+'DenMul']
2616                irho[ilay] = parmDict[cid+'iDenMul']
2617            elif 'vacuum' != parmDict[cid+'Name']:
2618                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
2619                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
2620            if cid+'Mag SLD' in parmDict:
2621                rho[ilay] += parmDict[cid+'Mag SLD']
2622        if parmDict['dQ type'] == 'None':
2623            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2624        elif 'const' in parmDict['dQ type']:
2625            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
2626        else:       #dQ/Q in data
2627            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
2628        Ic += AB*Scale
2629        return Ic
2630       
2631    def estimateT0(takestep):
2632        Mmax = -1.e-10
2633        Mmin = 1.e10
2634        for i in range(100):
2635            x0 = takestep(values)
2636            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2637            Mmin = min(M,Mmin)
2638            MMax = max(M,Mmax)
2639        return 1.5*(MMax-Mmin)
2640
2641    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2642    if data.get('2% weight'):
2643        wt = 1./(0.02*Io)**2
2644    Qmin = Limits[1][0]
2645    Qmax = Limits[1][1]
2646    wtFactor = ProfDict['wtFactor']
2647    Bfac = data['Toler']
2648    Ibeg = np.searchsorted(Q,Qmin)
2649    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
2650    Ic[:] = 0
2651    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
2652              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
2653    parmDict,varyList,values,bounds = GetModelParms()
2654    Msg = 'Failed to converge'
2655    if varyList:
2656        if data['Minimizer'] == 'LMLS': 
2657            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
2658                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2659            parmDict.update(zip(varyList,result[0]))
2660            chisq = np.sum(result[2]['fvec']**2)
2661            ncalc = result[2]['nfev']
2662            covM = result[1]
2663            newVals = result[0]
2664        elif data['Minimizer'] == 'Basin Hopping':
2665            xyrng = np.array(bounds).T
2666            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
2667            T0 = estimateT0(take_step)
2668            G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
2669            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
2670                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
2671                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
2672            chisq = result.fun
2673            ncalc = result.nfev
2674            newVals = result.x
2675            covM = []
2676        elif data['Minimizer'] == 'MC/SA Anneal':
2677            xyrng = np.array(bounds).T
2678            result = G2mth.anneal(sumREFD, values, 
2679                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
2680                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
2681                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
2682                ranRange=0.20,autoRan=False,dlg=None)
2683            newVals = result[0]
2684            parmDict.update(zip(varyList,newVals))
2685            chisq = result[1]
2686            ncalc = result[3]
2687            covM = []
2688            G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
2689        elif data['Minimizer'] == 'L-BFGS-B':
2690            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
2691                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2692            parmDict.update(zip(varyList,result['x']))
2693            chisq = result.fun
2694            ncalc = result.nfev
2695            newVals = result.x
2696            covM = []
2697    else:   #nothing varied
2698        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2699        chisq = np.sum(M**2)
2700        ncalc = 0
2701        covM = []
2702        sig = []
2703        sigDict = {}
2704        result = []
2705    Rvals = {}
2706    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
2707    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
2708    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
2709    Ib[Ibeg:Ifin] = parmDict['FltBack']
2710    try:
2711        if not len(varyList):
2712            Msg += ' - nothing refined'
2713            raise ValueError
2714        Nans = np.isnan(newVals)
2715        if np.any(Nans):
2716            name = varyList[Nans.nonzero(True)[0]]
2717            Msg += ' Nan result for '+name+'!'
2718            raise ValueError
2719        Negs = np.less_equal(newVals,0.)
2720        if np.any(Negs):
2721            indx = Negs.nonzero()
2722            name = varyList[indx[0][0]]
2723            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
2724                Msg += ' negative coefficient for '+name+'!'
2725                raise ValueError
2726        if len(covM):
2727            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
2728            covMatrix = covM*Rvals['GOF']
2729        else:
2730            sig = np.zeros(len(varyList))
2731            covMatrix = []
2732        sigDict = dict(zip(varyList,sig))
2733        G2fil.G2Print (' Results of reflectometry data modelling fit:')
2734        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
2735        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2736        SetModelParms()
2737        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
2738    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
2739        G2fil.G2Print (Msg)
2740        return False,0,0,0,0,0,0,Msg
2741       
2742def makeSLDprofile(data,Substances):
2743   
2744    sq2 = np.sqrt(2.)
2745    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2746    Nlayers = len(laySeq)
2747    laySeq = np.array(laySeq,dtype=int)
2748    interfaces = np.zeros(Nlayers)
2749    rho = np.zeros(Nlayers)
2750    sigma = np.zeros(Nlayers)
2751    name = data['Layers'][0]['Name']
2752    thick = 0.
2753    for ilay,lay in enumerate(laySeq):
2754        layer = data['Layers'][lay]
2755        name = layer['Name']
2756        if 'Thick' in layer:
2757            thick += layer['Thick'][0]
2758            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
2759        if 'Rough' in layer:
2760            sigma[ilay] = max(0.001,layer['Rough'][0])
2761        if name != 'vacuum':
2762            if name == 'unit scatter':
2763                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
2764            else:
2765                rrho = Substances[name]['Scatt density']
2766                irho = Substances[name]['XImag density']
2767                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
2768        if 'Mag SLD' in layer:
2769            rho[ilay] += layer['Mag SLD'][0]
2770    name = data['Layers'][-1]['Name']
2771    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
2772    xr = np.flipud(x)
2773    interfaces[-1] = x[-1]
2774    y = np.ones_like(x)*rho[0]
2775    iBeg = 0
2776    for ilayer in range(Nlayers-1):
2777        delt = rho[ilayer+1]-rho[ilayer]
2778        iPos = np.searchsorted(x,interfaces[ilayer])
2779        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
2780        iBeg = iPos
2781    return x,xr,y   
2782
2783def REFDModelFxn(Profile,Inst,Limits,Substances,data):
2784   
2785    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2786    Qmin = Limits[1][0]
2787    Qmax = Limits[1][1]
2788    iBeg = np.searchsorted(Q,Qmin)
2789    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2790    Ib[:] = data['FltBack'][0]
2791    Ic[:] = 0
2792    Scale = data['Scale'][0]
2793    if data['Layer Seq'] == []:
2794        return
2795    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2796    Nlayers = len(laySeq)
2797    depth = np.zeros(Nlayers)
2798    rho = np.zeros(Nlayers)
2799    irho = np.zeros(Nlayers)
2800    sigma = np.zeros(Nlayers)
2801    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
2802        layer = data['Layers'][lay]
2803        name = layer['Name']
2804        if 'Thick' in layer:    #skips first & last layers
2805            depth[ilay] = layer['Thick'][0]
2806        if 'Rough' in layer:    #skips first layer
2807            sigma[ilay] = layer['Rough'][0]
2808        if 'unit scatter' == name:
2809            rho[ilay] = layer['DenMul'][0]
2810            irho[ilay] = layer['iDenMul'][0]
2811        else:
2812            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2813            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
2814        if 'Mag SLD' in layer:
2815            rho[ilay] += layer['Mag SLD'][0]
2816    if data['dQ type'] == 'None':
2817        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
2818    elif 'const' in data['dQ type']:
2819        res = data['Resolution'][0]/(100.*sateln2)
2820        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
2821    else:       #dQ/Q in data
2822        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
2823    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
2824
2825def abeles(kz, depth, rho, irho=0, sigma=0):
2826    """
2827    Optical matrix form of the reflectivity calculation.
2828    O.S. Heavens, Optical Properties of Thin Solid Films
2829   
2830    Reflectometry as a function of kz for a set of slabs.
2831
2832    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
2833        This is :math:`\\tfrac12 Q_z`.       
2834    :param depth:  float[m] (Ang).
2835        thickness of each layer.  The thickness of the incident medium
2836        and substrate are ignored.
2837    :param rho:  float[n,k] (1e-6/Ang^2)
2838        Real scattering length density for each layer for each kz
2839    :param irho:  float[n,k] (1e-6/Ang^2)
2840        Imaginary scattering length density for each layer for each kz
2841        Note: absorption cross section mu = 2 irho/lambda for neutrons
2842    :param sigma: float[m-1] (Ang)
2843        interfacial roughness.  This is the roughness between a layer
2844        and the previous layer. The sigma array should have m-1 entries.
2845
2846    Slabs are ordered with the surface SLD at index 0 and substrate at
2847    index -1, or reversed if kz < 0.
2848    """
2849    def calc(kz, depth, rho, irho, sigma):
2850        if len(kz) == 0: return kz
2851   
2852        # Complex index of refraction is relative to the incident medium.
2853        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
2854        # in place of kz^2, and ignoring rho_o
2855        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
2856        k = kz
2857   
2858        # According to Heavens, the initial matrix should be [ 1 F; F 1],
2859        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
2860        # multiply versus some coding convenience.
2861        B11 = 1
2862        B22 = 1
2863        B21 = 0
2864        B12 = 0
2865        for i in range(0, len(depth)-1):
2866            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
2867            F = (k - k_next) / (k + k_next)
2868            F *= np.exp(-2*k*k_next*sigma[i]**2)
2869            #print "==== layer",i
2870            #print "kz:", kz
2871            #print "k:", k
2872            #print "k_next:",k_next
2873            #print "F:",F
2874            #print "rho:",rho[:,i+1]
2875            #print "irho:",irho[:,i+1]
2876            #print "d:",depth[i],"sigma:",sigma[i]
2877            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
2878            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
2879            M21 = F*M11
2880            M12 = F*M22
2881            C1 = B11*M11 + B21*M12
2882            C2 = B11*M21 + B21*M22
2883            B11 = C1
2884            B21 = C2
2885            C1 = B12*M11 + B22*M12
2886            C2 = B12*M21 + B22*M22
2887            B12 = C1
2888            B22 = C2
2889            k = k_next
2890   
2891        r = B12/B11
2892        return np.absolute(r)**2
2893
2894    if np.isscalar(kz): kz = np.asarray([kz], 'd')
2895
2896    m = len(depth)
2897
2898    # Make everything into arrays
2899    depth = np.asarray(depth,'d')
2900    rho = np.asarray(rho,'d')
2901    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
2902    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
2903
2904    # Repeat rho,irho columns as needed
2905    if len(rho.shape) == 1:
2906        rho = rho[None,:]
2907        irho = irho[None,:]
2908
2909    return calc(kz, depth, rho, irho, sigma)
2910   
2911def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
2912    y = abeles(kz, depth, rho, irho, sigma)
2913    s = dq/2.
2914    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
2915    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
2916    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
2917    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
2918    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
2919    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
2920    y *= 0.137023
2921    return y
2922       
2923def makeRefdFFT(Limits,Profile):
2924    G2fil.G2Print ('make fft')
2925    Q,Io = Profile[:2]
2926    Qmin = Limits[1][0]
2927    Qmax = Limits[1][1]
2928    iBeg = np.searchsorted(Q,Qmin)
2929    iFin = np.searchsorted(Q,Qmax)+1    #include last point
2930    Qf = np.linspace(0.,Q[iFin-1],5000)
2931    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
2932    If = QI(Qf)*Qf**4
2933    R = np.linspace(0.,5000.,5000)
2934    F = fft.rfft(If)
2935    return R,F
2936
2937   
2938################################################################################
2939#### Stacking fault simulation codes
2940################################################################################
2941
2942def GetStackParms(Layers):
2943   
2944    Parms = []
2945#cell parms
2946    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
2947        Parms.append('cellA')
2948        Parms.append('cellC')
2949    else:
2950        Parms.append('cellA')
2951        Parms.append('cellB')
2952        Parms.append('cellC')
2953        if Layers['Laue'] != 'mmm':
2954            Parms.append('cellG')
2955#Transition parms
2956    for iY in range(len(Layers['Layers'])):
2957        for iX in range(len(Layers['Layers'])):
2958            Parms.append('TransP;%d;%d'%(iY,iX))
2959            Parms.append('TransX;%d;%d'%(iY,iX))
2960            Parms.append('TransY;%d;%d'%(iY,iX))
2961            Parms.append('TransZ;%d;%d'%(iY,iX))
2962    return Parms
2963
2964def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
2965    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
2966   
2967    :param dict Layers: dict with following items
2968
2969      ::
2970
2971       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
2972       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
2973       'Layers':[],'Stacking':[],'Transitions':[]}
2974       
2975    :param str ctrls: controls string to be written on DIFFaX controls.dif file
2976    :param float scale: scale factor
2977    :param dict background: background parameters
2978    :param list limits: min/max 2-theta to be calculated
2979    :param dict inst: instrument parameters dictionary
2980    :param list profile: powder pattern data
2981   
2982    Note that parameters all updated in place   
2983    '''
2984    import atmdata
2985    path = sys.path
2986    for name in path:
2987        if 'bin' in name:
2988            DIFFaX = name+'/DIFFaX.exe'
2989            G2fil.G2Print (' Execute '+DIFFaX)
2990            break
2991    # make form factor file that DIFFaX wants - atom types are GSASII style
2992    sf = open('data.sfc','w')
2993    sf.write('GSASII special form factor file for DIFFaX\n\n')
2994    atTypes = list(Layers['AtInfo'].keys())
2995    if 'H' not in atTypes:
2996        atTypes.insert(0,'H')
2997    for atType in atTypes:
2998        if atType == 'H': 
2999            blen = -.3741
3000        else:
3001            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
3002        Adat = atmdata.XrayFF[atType]
3003        text = '%4s'%(atType.ljust(4))
3004        for i in range(4):
3005            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
3006        text += '%11.6f%11.6f'%(Adat['fc'],blen)
3007        text += '%3d\n'%(Adat['Z'])
3008        sf.write(text)
3009    sf.close()
3010    #make DIFFaX control.dif file - future use GUI to set some of these flags
3011    cf = open('control.dif','w')
3012    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
3013        x0 = profile[0]
3014        iBeg = np.searchsorted(x0,limits[0])
3015        iFin = np.searchsorted(x0,limits[1])+1
3016        if iFin-iBeg > 20000:
3017            iFin = iBeg+20000
3018        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3019        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3020        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
3021    else:
3022        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3023        inst = {'Type':['XSC','XSC',]}
3024    cf.close()
3025    #make DIFFaX data file
3026    df = open('GSASII-DIFFaX.dat','w')
3027    df.write('INSTRUMENTAL\n')
3028    if 'X' in inst['Type'][0]:
3029        df.write('X-RAY\n')
3030    elif 'N' in inst['Type'][0]:
3031        df.write('NEUTRON\n')
3032    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
3033        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
3034        U = ateln2*inst['U'][1]/10000.
3035        V = ateln2*inst['V'][1]/10000.
3036        W = ateln2*inst['W'][1]/10000.
3037        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
3038        HW = np.sqrt(np.mean(HWHM))
3039    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
3040        if 'Mean' in Layers['selInst']:
3041            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
3042        elif 'Gaussian' in Layers['selInst']:
3043            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
3044        else:
3045            df.write('None\n')
3046    else:
3047        df.write('0.10\nNone\n')
3048    df.write('STRUCTURAL\n')
3049    a,b,c = Layers['Cell'][1:4]
3050    gam = Layers['Cell'][6]
3051    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
3052    laue = Layers['Laue']
3053    if laue == '2/m(ab)':
3054        laue = '2/m(1)'
3055    elif laue == '2/m(c)':
3056        laue = '2/m(2)'
3057    if 'unknown' in Layers['Laue']:
3058        df.write('%s %.3f\n'%(laue,Layers['Toler']))
3059    else:   
3060        df.write('%s\n'%(laue))
3061    df.write('%d\n'%(len(Layers['Layers'])))
3062    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
3063        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
3064    layerNames = []
3065    for layer in Layers['Layers']:
3066        layerNames.append(layer['Name'])
3067    for il,layer in enumerate(Layers['Layers']):
3068        if layer['SameAs']:
3069            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
3070            continue
3071        df.write('LAYER %d\n'%(il+1))
3072        if '-1' in layer['Symm']:
3073            df.write('CENTROSYMMETRIC\n')
3074        else:
3075            df.write('NONE\n')
3076        for ia,atom in enumerate(layer['Atoms']):
3077            [name,atype,x,y,z,frac,Uiso] = atom
3078            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
3079                frac /= 2.
3080            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
3081    df.write('STACKING\n')
3082    df.write('%s\n'%(Layers['Stacking'][0]))
3083    if 'recursive' in Layers['Stacking'][0]:
3084        df.write('%s\n'%Layers['Stacking'][1])
3085    else:
3086        if 'list' in Layers['Stacking'][1]:
3087            Slen = len(Layers['Stacking'][2])
3088            iB = 0
3089            iF = 0
3090            while True:
3091                iF += 68
3092                if iF >= Slen:
3093                    break
3094                iF = min(iF,Slen)
3095                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
3096                iB = iF
3097        else:
3098            df.write('%s\n'%Layers['Stacking'][1])   
3099    df.write('TRANSITIONS\n')
3100    for iY in range(len(Layers['Layers'])):
3101        sumPx = 0.
3102        for iX in range(len(Layers['Layers'])):
3103            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
3104            p = round(p,3)
3105            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
3106            sumPx += p
3107        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
3108            G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
3109            df.close()
3110            os.remove('data.sfc')
3111            os.remove('control.dif')
3112            os.remove('GSASII-DIFFaX.dat')
3113            return       
3114    df.close()   
3115    time0 = time.time()
3116    try:
3117        subp.call(DIFFaX)
3118    except OSError:
3119        G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
3120    G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
3121    if os.path.exists('GSASII-DIFFaX.spc'):
3122        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
3123        iFin = iBeg+Xpat.shape[1]
3124        bakType,backDict,backVary = SetBackgroundParms(background)
3125        backDict['Lam1'] = G2mth.getWave(inst)
3126        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
3127        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
3128        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
3129            rv = st.poisson(profile[3][iBeg:iFin])
3130            profile[1][iBeg:iFin] = rv.rvs()
3131            Z = np.ones_like(profile[3][iBeg:iFin])
3132            Z[1::2] *= -1
3133            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3134            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3135        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3136    #cleanup files..
3137        os.remove('GSASII-DIFFaX.spc')
3138    elif os.path.exists('GSASII-DIFFaX.sadp'):
3139        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
3140        Sadp = np.reshape(Sadp,(256,-1))
3141        Layers['Sadp']['Img'] = Sadp
3142        os.remove('GSASII-DIFFaX.sadp')
3143    os.remove('data.sfc')
3144    os.remove('control.dif')
3145    os.remove('GSASII-DIFFaX.dat')
3146   
3147def SetPWDRscan(inst,limits,profile):
3148   
3149    wave = G2mth.getMeanWave(inst)
3150    x0 = profile[0]
3151    iBeg = np.searchsorted(x0,limits[0])
3152    iFin = np.searchsorted(x0,limits[1])
3153    if iFin-iBeg > 20000:
3154        iFin = iBeg+20000
3155    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3156    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
3157    return iFin-iBeg
3158       
3159def SetStackingSF(Layers,debug):
3160# Load scattering factors into DIFFaX arrays
3161    import atmdata
3162    atTypes = Layers['AtInfo'].keys()
3163    aTypes = []
3164    for atype in atTypes:
3165        aTypes.append('%4s'%(atype.ljust(4)))
3166    SFdat = []
3167    for atType in atTypes:
3168        Adat = atmdata.XrayFF[atType]
3169        SF = np.zeros(9)
3170        SF[:8:2] = Adat['fa']
3171        SF[1:8:2] = Adat['fb']
3172        SF[8] = Adat['fc']
3173        SFdat.append(SF)
3174    SFdat = np.array(SFdat)
3175    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
3176   
3177def SetStackingClay(Layers,Type):
3178# Controls
3179    rand.seed()
3180    ranSeed = rand.randint(1,2**16-1)
3181    try:   
3182        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
3183            '6/m','6/mmm'].index(Layers['Laue'])+1
3184    except ValueError:  #for 'unknown'
3185        laueId = -1
3186    if 'SADP' in Type:
3187        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
3188        lmax = int(Layers['Sadp']['Lmax'])
3189    else:
3190        planeId = 0
3191        lmax = 0
3192# Sequences
3193    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
3194    try:
3195        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
3196    except ValueError:
3197        StkParm = -1
3198    if StkParm == 2:    #list
3199        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
3200        Nstk = len(StkSeq)
3201    else:
3202        Nstk = 1
3203        StkSeq = [0,]
3204    if StkParm == -1:
3205        StkParm = int(Layers['Stacking'][1])
3206    Wdth = Layers['Width'][0]
3207    mult = 1
3208    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
3209    LaueSym = Layers['Laue'].ljust(12)
3210    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
3211    return laueId,controls
3212   
3213def SetCellAtoms(Layers):
3214    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
3215# atoms in layers
3216    atTypes = list(Layers['AtInfo'].keys())
3217    AtomXOU = []
3218    AtomTp = []
3219    LayerSymm = []
3220    LayerNum = []
3221    layerNames = []
3222    Natm = 0
3223    Nuniq = 0
3224    for layer in Layers['Layers']:
3225        layerNames.append(layer['Name'])
3226    for il,layer in enumerate(Layers['Layers']):
3227        if layer['SameAs']:
3228            LayerNum.append(layerNames.index(layer['SameAs'])+1)
3229            continue
3230        else:
3231            LayerNum.append(il+1)
3232            Nuniq += 1
3233        if '-1' in layer['Symm']:
3234            LayerSymm.append(1)
3235        else:
3236            LayerSymm.append(0)
3237        for ia,atom in enumerate(layer['Atoms']):
3238            [name,atype,x,y,z,frac,Uiso] = atom
3239            Natm += 1
3240            AtomTp.append('%4s'%(atype.ljust(4)))
3241            Ta = atTypes.index(atype)+1
3242            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
3243    AtomXOU = np.array(AtomXOU)
3244    Nlayers = len(layerNames)
3245    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
3246    return Nlayers
3247   
3248def SetStackingTrans(Layers,Nlayers):
3249# Transitions
3250    TransX = []
3251    TransP = []
3252    for Ytrans in Layers['Transitions']:
3253        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
3254        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
3255    TransP = np.array(TransP,dtype='float').T
3256    TransX = np.array(TransX,dtype='float')
3257#    GSASIIpath.IPyBreak()
3258    pyx.pygettrans(Nlayers,TransP,TransX)
3259   
3260def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
3261# Scattering factors
3262    SetStackingSF(Layers,debug)
3263# Controls & sequences
3264    laueId,controls = SetStackingClay(Layers,'PWDR')
3265# cell & atoms
3266    Nlayers = SetCellAtoms(Layers)
3267    Volume = Layers['Cell'][7]   
3268# Transitions
3269    SetStackingTrans(Layers,Nlayers)
3270# PWDR scan
3271    Nsteps = SetPWDRscan(inst,limits,profile)
3272# result as Spec
3273    x0 = profile[0]
3274    profile[3] = np.zeros(len(profile[0]))
3275    profile[4] = np.zeros(len(profile[0]))
3276    profile[5] = np.zeros(len(profile[0]))
3277    iBeg = np.searchsorted(x0,limits[0])
3278    iFin = np.searchsorted(x0,limits[1])+1
3279    if iFin-iBeg > 20000:
3280        iFin = iBeg+20000
3281    Nspec = 20001       
3282    spec = np.zeros(Nspec,dtype='double')   
3283    time0 = time.time()
3284    pyx.pygetspc(controls,Nspec,spec)
3285    G2fil.G2Print (' GETSPC time = %.2fs'%(time.time()-time0))
3286    time0 = time.time()
3287    U = ateln2*inst['U'][1]/10000.
3288    V = ateln2*inst['V'][1]/10000.
3289    W = ateln2*inst['W'][1]/10000.
3290    HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
3291    HW = np.sqrt(np.mean(HWHM))
3292    BrdSpec = np.zeros(Nsteps)
3293    if 'Mean' in Layers['selInst']:
3294        pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
3295    elif 'Gaussian' in Layers['selInst']:
3296        pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
3297    else:
3298        BrdSpec = spec[:Nsteps]
3299    BrdSpec /= Volume
3300    iFin = iBeg+Nsteps
3301    bakType,backDict,backVary = SetBackgroundParms(background)
3302    backDict['Lam1'] = G2mth.getWave(inst)
3303    profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
3304    profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
3305    if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
3306        try:
3307            rv = st.poisson(profile[3][iBeg:iFin])
3308            profile[1][iBeg:iFin] = rv.rvs()
3309        except ValueError:
3310            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
3311        Z = np.ones_like(profile[3][iBeg:iFin])
3312        Z[1::2] *= -1
3313        profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3314        profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3315    profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3316    G2fil.G2Print (' Broadening time = %.2fs'%(time.time()-time0))
3317   
3318def CalcStackingSADP(Layers,debug):
3319   
3320# Scattering factors
3321    SetStackingSF(Layers,debug)
3322# Controls & sequences
3323    laueId,controls = SetStackingClay(Layers,'SADP')
3324# cell & atoms
3325    Nlayers = SetCellAtoms(Layers)   
3326# Transitions
3327    SetStackingTrans(Layers,Nlayers)
3328# result as Sadp
3329    Nspec = 20001       
3330    spec = np.zeros(Nspec,dtype='double')   
3331    time0 = time.time()
3332    hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
3333    Sapd = np.zeros((256,256))
3334    iB = 0
3335    for i in range(hkLim):
3336        iF = iB+Nblk
3337        p1 = 127+int(i*Incr)
3338        p2 = 128-int(i*Incr)
3339        if Nblk == 128:
3340            if i:
3341                Sapd[128:,p1] = spec[iB:iF]
3342                Sapd[:128,p1] = spec[iF:iB:-1]
3343            Sapd[128:,p2] = spec[iB:iF]
3344            Sapd[:128,p2] = spec[iF:iB:-1]
3345        else:
3346            if i:
3347                Sapd[:,p1] = spec[iB:iF]
3348            Sapd[:,p2] = spec[iB:iF]
3349        iB += Nblk
3350    Layers['Sadp']['Img'] = Sapd
3351    G2fil.G2Print (' GETSAD time = %.2fs'%(time.time()-time0))
3352   
3353###############################################################################
3354#### Maximum Entropy Method - Dysnomia
3355###############################################################################
3356   
3357def makePRFfile(data,MEMtype):
3358    ''' makes Dysnomia .prf control file from Dysnomia GUI controls
3359   
3360    :param dict data: GSAS-II phase data
3361    :param int MEMtype: 1 for neutron data with negative scattering lengths
3362                        0 otherwise
3363    :returns str: name of Dysnomia control file
3364    '''
3365
3366    generalData = data['General']
3367    pName = generalData['Name'].replace(' ','_')
3368    DysData = data['Dysnomia']
3369    prfName = pName+'.prf'
3370    prf = open(prfName,'w')
3371    prf.write('$PREFERENCES\n')
3372    prf.write(pName+'.mem\n') #or .fos?
3373    prf.write(pName+'.out\n')
3374    prf.write(pName+'.pgrid\n')
3375    prf.write(pName+'.fba\n')
3376    prf.write(pName+'_eps.raw\n')
3377    prf.write('%d\n'%MEMtype)
3378    if DysData['DenStart'] == 'uniform':
3379        prf.write('0\n')
3380    else:
3381        prf.write('1\n')
3382    if DysData['Optimize'] == 'ZSPA':
3383        prf.write('0\n')
3384    else:
3385        prf.write('1\n')
3386    prf.write('1\n')
3387    if DysData['Lagrange'][0] == 'user':
3388        prf.write('0\n')
3389    else:
3390        prf.write('1\n')
3391    prf.write('%.4f %d\n'%(DysData['Lagrange'][1],DysData['wt pwr']))
3392    prf.write('%.3f\n'%DysData['Lagrange'][2])
3393    prf.write('%.2f\n'%DysData['E_factor'])
3394    prf.write('1\n')
3395    prf.write('0\n')
3396    prf.write('%d\n'%DysData['Ncyc'])
3397    prf.write('1\n')
3398    prf.write('1 0 0 0 0 0 0 0\n')
3399    if DysData['prior'] == 'uniform':
3400        prf.write('0\n')
3401    else:
3402        prf.write('1\n')
3403    prf.close()
3404    return prfName
3405
3406def makeMEMfile(data,reflData,MEMtype,DYSNOMIA):
3407    ''' make Dysnomia .mem file of reflection data, etc.
3408
3409    :param dict data: GSAS-II phase data
3410    :param list reflData: GSAS-II reflection data
3411    :param int MEMtype: 1 for neutron data with negative scattering lengths
3412                        0 otherwise
3413    :param str DYSNOMIA: path to dysnomia.exe
3414    '''
3415   
3416    DysData = data['Dysnomia']
3417    generalData = data['General']
3418    cell = generalData['Cell'][1:7]
3419    A = G2lat.cell2A(cell)
3420    SGData = generalData['SGData']
3421    pName = generalData['Name'].replace(' ','_')
3422    memName = pName+'.mem'
3423    Map = generalData['Map']
3424    Type = Map['Type']
3425    UseList = Map['RefList']
3426    mem = open(memName,'w')
3427    mem.write('%s\n'%(generalData['Name']+' from '+UseList[0]))
3428    a,b,c,alp,bet,gam = cell
3429    mem.write('%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n'%(a,b,c,alp,bet,gam))
3430    mem.write('      0.0000000      0.0000000     -1    0    0    0     P\n')   #dummy PO stuff
3431    SGSym = generalData['SGData']['SpGrp']
3432    try:
3433        SGId = G2spc.spgbyNum.index(SGSym)
3434    except ValueError:
3435        return False
3436    org = 1
3437    if SGSym in G2spc.spg2origins:
3438        org = 2
3439    mapsize = Map['rho'].shape
3440    sumZ = 0.
3441    sumpos = 0.
3442    sumneg = 0.
3443    mem.write('%5d%5d%5d%5d%5d\n'%(SGId,org,mapsize[0],mapsize[1],mapsize[2]))
3444    for atm in generalData['NoAtoms']:
3445        Nat = generalData['NoAtoms'][atm]
3446        AtInfo = G2elem.GetAtomInfo(atm)
3447        sumZ += Nat*AtInfo['Z']
3448        isotope = generalData['Isotope'][atm]
3449        blen = generalData['Isotopes'][atm][isotope]['SL'][0]
3450        if blen < 0.:
3451            sumneg += blen*Nat
3452        else:
3453            sumpos += blen*Nat
3454    if 'X' in Type:
3455        mem.write('%10.2f  0.001\n'%sumZ)
3456    elif 'N' in Type and MEMtype:
3457        mem.write('%10.3f%10.3f 0.001\n'%(sumpos,sumneg))
3458    else:
3459        mem.write('%10.3f 0.001\n'%sumpos)
3460       
3461    dmin = DysData['MEMdmin']
3462    TOFlam = 2.0*dmin*npsind(80.0)
3463    refSet = G2lat.GenHLaue(dmin,SGData,A)      #list of h,k,l,d
3464    refDict = {'%d %d %d'%(ref[0],ref[1],ref[2]):ref for ref in refSet}
3465       
3466    refs = []
3467    prevpos = 0.
3468    for ref in reflData:
3469        if ref[3] < 0:
3470            continue
3471        if 'T' in Type:
3472            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,x,x,x,prfo = ref[:16]
3473            s = np.sqrt(max(sig,0.0001))   #var -> sig in deg
3474            FWHM = getgamFW(gam,s)
3475            if dsp < dmin:
3476                continue
3477            theta = npasind(TOFlam/(2.*dsp))
3478            FWHM *= nptand(theta)/pos
3479            pos = 2.*theta
3480        else:
3481            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,prfo = ref[:13]
3482            g = gam/100.    #centideg -> deg
3483            s = np.sqrt(max(sig,0.0001))/100.   #var -> sig in deg
3484            FWHM = getgamFW(g,s)
3485        delt = pos-prevpos
3486        refs.append([h,k,l,mult,pos,FWHM,Fobs,phase,delt])
3487        prevpos = pos
3488           
3489    ovlp = DysData['overlap']
3490    refs1 = []
3491    refs2 = []
3492    nref2 = 0
3493    iref = 0
3494    Nref = len(refs)
3495    start = False
3496    while iref < Nref-1:
3497        if refs[iref+1][-1] < ovlp*refs[iref][5]:
3498            if refs[iref][-1] > ovlp*refs[iref][5]:
3499                refs2.append([])
3500                start = True
3501            if nref2 == len(refs2):
3502                refs2.append([])
3503            refs2[nref2].append(refs[iref])
3504        else:
3505            if start:
3506                refs2[nref2].append(refs[iref])
3507                start = False
3508                nref2 += 1
3509            else:
3510                refs1.append(refs[iref])
3511        iref += 1
3512    if start:
3513        refs2[nref2].append(refs[iref])
3514    else:
3515        refs1.append(refs[iref])
3516   
3517    mem.write('%5d\n'%len(refs1))
3518    for ref in refs1:
3519        h,k,l = ref[:3]
3520        hkl = '%d %d %d'%(h,k,l)
3521        if hkl in refDict:
3522            del refDict[hkl]
3523        Fobs = np.sqrt(ref[6])
3524        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)))
3525    while True and nref2:
3526        if not len(refs2[-1]):
3527            del refs2[-1]
3528        else:
3529            break
3530    mem.write('%5d\n'%len(refs2))
3531    for iref2,ref2 in enumerate(refs2):
3532        mem.write('#%5d\n'%iref2)
3533        mem.write('%5d\n'%len(ref2))
3534        Gsum = 0.
3535        Msum = 0
3536        for ref in ref2:
3537            Gsum += ref[6]*ref[3]
3538            Msum += ref[3]
3539        G = np.sqrt(Gsum/Msum)
3540        h,k,l = ref2[0][:3]
3541        hkl = '%d %d %d'%(h,k,l)
3542        if hkl in refDict:
3543            del refDict[hkl]
3544        mem.write('%5d%5d%5d%10.3f%10.3f%5d\n'%(h,k,l,G,max(0.01*G,0.1),ref2[0][3]))
3545        for ref in ref2[1:]:
3546            h,k,l,m = ref[:4]
3547            mem.write('%5d%5d%5d%5d\n'%(h,k,l,m))
3548            hkl = '%d %d %d'%(h,k,l)
3549            if hkl in refDict:
3550                del refDict[hkl]
3551    if len(refDict):
3552        mem.write('%d\n'%len(refDict))
3553        for hkl in list(refDict.keys()):
3554            h,k,l = refDict[hkl][:3]
3555            mem.write('%5d%5d%5d\n'%(h,k,l))
3556    else:
3557        mem.write('0\n')
3558    mem.close()
3559    return True
3560
3561def MEMupdateReflData(prfName,data,reflData):
3562    ''' Update reflection data with new Fosq, phase result from Dysnomia
3563
3564    :param str prfName: phase.mem file name
3565    :param list reflData: GSAS-II reflection data
3566    '''
3567   
3568    generalData = data['General']
3569    Map = generalData['Map']
3570    Type = Map['Type']
3571    cell = generalData['Cell'][1:7]
3572    A = G2lat.cell2A(cell)
3573    reflDict = {}
3574    newRefs = []
3575    for iref,ref in enumerate(reflData):
3576        if ref[3] > 0:
3577            newRefs.append(ref)
3578            reflDict[hash('%5d%5d%5d'%(ref[0],ref[1],ref[2]))] = iref
3579    fbaName = os.path.splitext(prfName)[0]+'.fba'
3580    try: # patch for FileNotFoundError not in Python 2.7
3581        FileNotFoundError
3582    except NameError:
3583        FileNotFoundError = Exception
3584    try:
3585        fba = open(fbaName,'r')
3586    except FileNotFoundError:
3587        return False
3588    fba.readline()
3589    Nref = int(fba.readline()[:-1])
3590    fbalines = fba.readlines()
3591    for line in fbalines[:Nref]:
3592        info = line.split()
3593        h = int(info[0])
3594        k = int(info[1])
3595        l = int(info[2])
3596        FoR = float(info[3])
3597        FoI = float(info[4])
3598        Fosq = FoR**2+FoI**2
3599        phase = npatan2d(FoI,FoR)
3600        try:
3601            refId = reflDict[hash('%5d%5d%5d'%(h,k,l))]
3602        except KeyError:    #added reflections at end skipped
3603            d = float(1/np.sqrt(G2lat.calc_rDsq([h,k,l],A)))
3604            if 'T' in Type:
3605                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])
3606            else:
3607                newRefs.append([h,k,l,-1,d,0.,0.01,1.0,Fosq,Fosq,phase,1.0,1.0,1.0,1.0])
3608            continue
3609        newRefs[refId][8] = Fosq
3610        newRefs[refId][10] = phase
3611    newRefs = np.array(newRefs)
3612    return True,newRefs
3613   
3614#### testing data
3615NeedTestData = True
3616def TestData():
3617    'needs a doc string'
3618#    global NeedTestData
3619    global bakType
3620    bakType = 'chebyschev'
3621    global xdata
3622    xdata = np.linspace(4.0,40.0,36000)
3623    global parmDict0
3624    parmDict0 = {
3625        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
3626        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
3627        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
3628        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
3629        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'Z':0.0,'SH/L':0.002,
3630        'Back0':5.384,'Back1':-0.015,'Back2':.004,
3631        }
3632    global parmDict1
3633    parmDict1 = {
3634        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
3635        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
3636        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
3637        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
3638        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
3639        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
3640        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'Z':0.0,'SH/L':0.002,
3641        'Back0':36.897,'Back1':-0.508,'Back2':.006,
3642        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
3643        }
3644    global parmDict2
3645    parmDict2 = {
3646        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
3647        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'Z':0.0,'SH/L':0.02,
3648        'Back0':5.,'Back1':-0.02,'Back2':.004,
3649#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
3650        }
3651    global varyList
3652    varyList = []
3653
3654def test0():
3655    if NeedTestData: TestData()
3656    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
3657    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
3658    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
3659    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
3660    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
3661    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
3662   
3663def test1():
3664    if NeedTestData: TestData()
3665    time0 = time.time()
3666    for i in range(100):
3667        getPeakProfile(parmDict1,xdata,varyList,bakType)
3668    G2fil.G2Print ('100+6*Ka1-2 peaks=1200 peaks %.2f'%time.time()-time0)
3669   
3670def test2(name,delt):
3671    if NeedTestData: TestData()
3672    varyList = [name,]
3673    xdata = np.linspace(5.6,5.8,400)
3674    hplot = plotter.add('derivatives test for '+name).gca()
3675    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
3676    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
3677    parmDict2[name] += delt
3678    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
3679    hplot.plot(xdata,(y1-y0)/delt,'r+')
3680   
3681def test3(name,delt):
3682    if NeedTestData: TestData()
3683    names = ['pos','sig','gam','shl']
3684    idx = names.index(name)
3685    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
3686    xdata = np.linspace(5.6,5.8,800)
3687    dx = xdata[1]-xdata[0]
3688    hplot = plotter.add('derivatives test for '+name).gca()
3689    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
3690    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
3691    myDict[name] += delt
3692    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
3693    hplot.plot(xdata,(y1-y0)/delt,'r+')
3694
3695if __name__ == '__main__':
3696    import GSASIItestplot as plot
3697    global plotter
3698    plotter = plot.PlotNotebook()
3699#    test0()
3700#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)']:
3701    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
3702        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['Z',0.0001],['SH/L',0.00005]]:
3703        test2(name,shft)
3704    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
3705        test3(name,shft)
3706    G2fil.G2Print ("OK")
3707    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.