source: trunk/GSASIIpwd.py @ 4671

Last change on this file since 4671 was 4671, checked in by vondreele, 10 months ago

add refinable multiplier for a fixed background; multiplier is now normally > 0
usable for peak fitting and Rietved refinement
add the fixed background entry to all background defaults
fix bug in GetDetectorXY fo when cursor outside image - returns [0,0] not None; changes elsewhere to use this
GetTthAzmDsp? now returns explicit list not assumed tuple
put the abs in the nl.qr test for singularities in HessianLSQ
Add 'BF mult' to name list in G2obj
put a try - except TypeError? around setting plot style stuff in PlotPatterns?
Remove the setting of Pattern[0]BackFile? - this was redundant for PWDR
remove picker/pickradius from linescan plot - failed
remove the alternate fixed background definition ('_fixedVary', etc.)
clear the PhaseReaderClass?Drawing? dictionary upon import of phase from a gpx file

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