source: trunk/GSASIIpwd.py @ 4880

Last change on this file since 4880 was 4880, checked in by toby, 6 months ago

implement constraints in scriptable, part 3; misc docs fixes

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