source: trunk/GSASIIpwd.py @ 5018

Last change on this file since 5018 was 5018, checked in by vondreele, 4 years ago

do np.nan_to_num on xye to avoid nans in plotting arrays.
remove extra reference to data in OnRelease? for reflection lists
move an import to top, fix a reference to 'B' in MakeRDF

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