source: trunk/GSASIIpwd.py @ 4333

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

work on spot masking; change some text, menu names, etc.
Add menu item for Auto search for selected images.

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