source: trunk/GSASIIpwd.py @ 4196

Last change on this file since 4196 was 4196, checked in by vondreele, 23 months ago

ad fullrmc dialog for set up
comment out & skip W plot update for d-plots & q-plots (they crash for TOF data)
put lower bounds on background single peak coefficients - avoids crash if refine negative
add MakePDB routine for fullrmc setup

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