source: trunk/GSASIIpwd.py @ 2189

Last change on this file since 2189 was 2189, checked in by vondreele, 6 years ago

comment out import pydiffax since g77 version woun't compile (argh)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 82.1 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASII powder calculation module*
5==================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2016-03-31 21:12:34 +0000 (Thu, 31 Mar 2016) $
10# $Author: vondreele $
11# $Revision: 2189 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 2189 2016-03-31 21:12:34Z vondreele $
14########### SVN repository information ###################
15import sys
16import math
17import time
18import os
19import subprocess as subp
20
21import numpy as np
22import scipy as sp
23import numpy.linalg as nl
24import random as rand
25from numpy.fft import ifft, fft, fftshift
26import scipy.interpolate as si
27import scipy.stats as st
28import scipy.optimize as so
29
30import GSASIIpath
31GSASIIpath.SetVersionNumber("$Revision: 2189 $")
32import GSASIIlattice as G2lat
33import GSASIIspc as G2spc
34import GSASIIElem as G2elem
35import GSASIIgrid as G2gd
36import GSASIIIO as G2IO
37import GSASIImath as G2mth
38import pypowder as pyd
39#import pydiffax as pyx
40
41# trig functions in degrees
42sind = lambda x: math.sin(x*math.pi/180.)
43asind = lambda x: 180.*math.asin(x)/math.pi
44tand = lambda x: math.tan(x*math.pi/180.)
45atand = lambda x: 180.*math.atan(x)/math.pi
46atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
47cosd = lambda x: math.cos(x*math.pi/180.)
48acosd = lambda x: 180.*math.acos(x)/math.pi
49rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
50#numpy versions
51npsind = lambda x: np.sin(x*np.pi/180.)
52npasind = lambda x: 180.*np.arcsin(x)/math.pi
53npcosd = lambda x: np.cos(x*math.pi/180.)
54npacosd = lambda x: 180.*np.arccos(x)/math.pi
55nptand = lambda x: np.tan(x*math.pi/180.)
56npatand = lambda x: 180.*np.arctan(x)/np.pi
57npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
58npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
59npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
60forln2 = 4.0*math.log(2.0)
61   
62#GSASII pdf calculation routines
63       
64def Transmission(Geometry,Abs,Diam):
65    '''
66    Calculate sample transmission
67
68    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
69    :param float Abs: absorption coeff in cm-1
70    :param float Diam: sample thickness/diameter in mm
71    '''
72    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
73        MuR = Abs*Diam/20.0
74        if MuR <= 3.0:
75            T0 = 16/(3.*math.pi)
76            T1 = -0.045780
77            T2 = -0.02489
78            T3 = 0.003045
79            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
80            if T < -20.:
81                return 2.06e-9
82            else:
83                return math.exp(T)
84        else:
85            T1 = 1.433902
86            T2 = 0.013869+0.337894
87            T3 = 1.933433+1.163198
88            T4 = 0.044365-0.04259
89            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
90            return T/100.
91    elif 'plate' in Geometry:
92        MuR = Abs*Diam/10.
93        return math.exp(-MuR)
94    elif 'Bragg' in Geometry:
95        return 0.0
96       
97def SurfaceRough(SRA,SRB,Tth):
98    ''' Suortti (J. Appl. Cryst, 5,325-331, 1972) surface roughness correction
99    :param float SRA: Suortti surface roughness parameter
100    :param float SRB: Suortti surface roughness parameter
101    :param float Tth: 2-theta(deg) - can be numpy array
102   
103    '''
104    sth = npsind(Tth/2.)
105    T1 = np.exp(-SRB/sth)
106    T2 = SRA+(1.-SRA)*np.exp(-SRB)
107    return (SRA+(1.-SRA)*T1)/T2
108   
109def SurfaceRoughDerv(SRA,SRB,Tth):
110    ''' Suortti surface roughness correction derivatives
111    :param float SRA: Suortti surface roughness parameter (dimensionless)
112    :param float SRB: Suortti surface roughness parameter (dimensionless)
113    :param float Tth: 2-theta(deg) - can be numpy array
114    :return list: [dydSRA,dydSRB] derivatives to be used for intensity derivative
115    '''
116    sth = npsind(Tth/2.)
117    T1 = np.exp(-SRB/sth)
118    T2 = SRA+(1.-SRA)*np.exp(-SRB)
119    Trans = (SRA+(1.-SRA)*T1)/T2
120    dydSRA = ((1.-T1)*T2-(1.-np.exp(-SRB))*Trans)/T2**2
121    dydSRB = ((SRA-1.)*T1*T2/sth-Trans*(SRA-T2))/T2**2
122    return [dydSRA,dydSRB]
123
124def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
125    '''Calculate sample absorption
126    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
127    :param float MuR: absorption coeff * sample thickness/2 or radius
128    :param Tth: 2-theta scattering angle - can be numpy array
129    :param float Phi: flat plate tilt angle - future
130    :param float Psi: flat plate tilt axis - future
131    '''
132   
133    def muRunder3(Sth2):
134        T0 = 16.0/(3.*np.pi)
135        T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
136            0.109561*np.sqrt(Sth2)-26.04556
137        T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
138            1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
139        T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
140        Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
141        return np.exp(Trns)
142   
143    def muRover3(Sth2):
144        T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
145            10.02088*Sth2**3-3.36778*Sth2**4
146        T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
147            (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
148        T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
149            0.13576*np.sqrt(Sth2)+1.163198
150        T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
151        Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
152        return Trns/100.
153       
154    Sth2 = npsind(Tth/2.0)**2
155    Cth2 = 1.-Sth2
156    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
157        if MuR <= 3.0:
158            return muRunder3(Sth2)
159        else:
160            return muRover3(Sth2)
161    elif 'Bragg' in Geometry:
162        return 1.0
163    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
164        # and only defined for 2theta < 90
165        MuT = 2.*MuR
166        T1 = np.exp(-MuT)
167        T2 = np.exp(-MuT/npcosd(Tth))
168        Tb = MuT-MuT/npcosd(Tth)
169        return (T2-T1)/Tb
170    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
171        MuT = 2.*MuR
172        cth = npcosd(Tth/2.0)
173        return np.exp(-MuT/cth)/cth
174       
175def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
176    'needs a doc string'
177    dA = 0.001
178    AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
179    if MuR:
180        AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
181        return (AbsP-AbsM)/(2.0*dA)
182    else:
183        return (AbsP-1.)/dA
184       
185def Polarization(Pola,Tth,Azm=0.0):
186    """   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
187
188    :param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
189    :param Azm: azimuthal angle e.g. 0.0 in plane of polarization
190    :param Tth: 2-theta scattering angle - can be numpy array
191      which (if either) of these is "right"?
192    :return: (pola, dpdPola)
193      * pola = ((1-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \
194        (1-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
195      * dpdPola: derivative needed for least squares
196
197    """
198    pola = ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
199        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
200    dpdPola = -npsind(Tth)**2*(npsind(Azm)**2-npcosd(Azm)**2)
201    return pola,dpdPola
202   
203def Oblique(ObCoeff,Tth):
204    'currently assumes detector is normal to beam'
205    if ObCoeff:
206        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
207    else:
208        return 1.0
209               
210def Ruland(RulCoff,wave,Q,Compton):
211    'needs a doc string'
212    C = 2.9978e8
213    D = 1.5e-3
214    hmc = 0.024262734687
215    sinth2 = (Q*wave/(4.0*np.pi))**2
216    dlam = (wave**2)*Compton*Q/C
217    dlam_c = 2.0*hmc*sinth2-D*wave**2
218    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
219   
220def LorchWeight(Q):
221    'needs a doc string'
222    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
223           
224def GetAsfMean(ElList,Sthl2):
225    '''Calculate various scattering factor terms for PDF calcs
226
227    :param dict ElList: element dictionary contains scattering factor coefficients, etc.
228    :param np.array Sthl2: numpy array of sin theta/lambda squared values
229    :returns: mean(f^2), mean(f)^2, mean(compton)
230    '''
231    sumNoAtoms = 0.0
232    FF = np.zeros_like(Sthl2)
233    FF2 = np.zeros_like(Sthl2)
234    CF = np.zeros_like(Sthl2)
235    for El in ElList:
236        sumNoAtoms += ElList[El]['FormulaNo']
237    for El in ElList:
238        el = ElList[El]
239        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
240        cf = G2elem.ComptonFac(el,Sthl2)
241        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
242        FF2 += ff2*el['FormulaNo']/sumNoAtoms
243        CF += cf*el['FormulaNo']/sumNoAtoms
244    return FF2,FF**2,CF
245   
246def GetNumDensity(ElList,Vol):
247    'needs a doc string'
248    sumNoAtoms = 0.0
249    for El in ElList:
250        sumNoAtoms += ElList[El]['FormulaNo']
251    return sumNoAtoms/Vol
252           
253def CalcPDF(data,inst,xydata):
254    'needs a doc string'
255    auxPlot = []
256    import copy
257    import scipy.fftpack as ft
258    #subtract backgrounds - if any
259    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
260    if data['Sample Bkg.']['Name']:
261        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
262            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
263    if data['Container']['Name']:
264        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
265        if data['Container Bkg.']['Name']:
266            xycontainer += (xydata['Container Bkg.'][1][1]+
267                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
268        xydata['IofQ'][1][1] += xycontainer
269    #get element data & absorption coeff.
270    ElList = data['ElList']
271    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
272    #Apply angle dependent corrections
273    Tth = xydata['Sample'][1][0]
274    dt = (Tth[1]-Tth[0])
275    MuR = Abs*data['Diam']/20.0
276    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],MuR,Tth)
277    if 'X' in inst['Type'][0]:
278        xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
279    if data['DetType'] == 'Image plate':
280        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
281    XY = xydata['IofQ'][1]   
282    #convert to Q
283    hc = 12.397639
284    wave = G2mth.getWave(inst)
285    keV = hc/wave
286    minQ = npT2q(Tth[0],wave)
287    maxQ = npT2q(Tth[-1],wave)   
288    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
289    dq = Qpoints[1]-Qpoints[0]
290    XY[0] = npT2q(XY[0],wave)   
291#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
292    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
293    Qdata = T(Qpoints)
294   
295    qLimits = data['QScaleLim']
296    minQ = np.searchsorted(Qpoints,qLimits[0])
297    maxQ = np.searchsorted(Qpoints,qLimits[1])
298    newdata = []
299    xydata['IofQ'][1][0] = Qpoints
300    xydata['IofQ'][1][1] = Qdata
301    for item in xydata['IofQ'][1]:
302        newdata.append(item[:maxQ])
303    xydata['IofQ'][1] = newdata
304   
305
306    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
307    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
308    Q = xydata['SofQ'][1][0]
309    ruland = Ruland(data['Ruland'],wave,Q,CF)
310#    auxPlot.append([Q,ruland,'Ruland'])     
311    CF *= ruland
312#    auxPlot.append([Q,CF,'CF-Corr'])
313    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
314    xydata['SofQ'][1][1] *= scale
315    xydata['SofQ'][1][1] -= CF
316    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
317    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
318    xydata['SofQ'][1][1] *= scale
319   
320    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
321    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
322    if data['Lorch']:
323        xydata['FofQ'][1][1] *= LorchWeight(Q)
324   
325    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
326    nR = len(xydata['GofR'][1][1])
327    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
328    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
329   
330       
331    return auxPlot
332
333################################################################################       
334#GSASII peak fitting routines: Finger, Cox & Jephcoat model       
335################################################################################
336
337def factorize(num):
338    ''' Provide prime number factors for integer num
339    :returns: dictionary of prime factors (keys) & power for each (data)
340    '''
341    factors = {}
342    orig = num
343
344    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
345    i, sqi = 2, 4
346    while sqi <= num:
347        while not num%i:
348            num /= i
349            factors[i] = factors.get(i, 0) + 1
350
351        sqi += 2*i + 1
352        i += 1
353
354    if num != 1 and num != orig:
355        factors[num] = factors.get(num, 0) + 1
356
357    if factors:
358        return factors
359    else:
360        return {num:1}          #a prime number!
361           
362def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
363    ''' Provide list of optimal data sizes for FFT calculations
364
365    :param int nmin: minimum data size >= 1
366    :param int nmax: maximum data size > nmin
367    :param int thresh: maximum prime factor allowed
368    :Returns: list of data sizes where the maximum prime factor is < thresh
369    ''' 
370    plist = []
371    nmin = max(1,nmin)
372    nmax = max(nmin+1,nmax)
373    for p in range(nmin,nmax):
374        if max(factorize(p).keys()) < thresh:
375            plist.append(p)
376    return plist
377
378np.seterr(divide='ignore')
379
380# Normal distribution
381
382# loc = mu, scale = std
383_norm_pdf_C = 1./math.sqrt(2*math.pi)
384class norm_gen(st.rv_continuous):
385    'needs a doc string'
386     
387    def pdf(self,x,*args,**kwds):
388        loc,scale=kwds['loc'],kwds['scale']
389        x = (x-loc)/scale
390        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
391       
392norm = norm_gen(name='norm',longname='A normal',extradoc="""
393
394Normal distribution
395
396The location (loc) keyword specifies the mean.
397The scale (scale) keyword specifies the standard deviation.
398
399normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
400""")
401
402## Cauchy
403
404# median = loc
405
406class cauchy_gen(st.rv_continuous):
407    'needs a doc string'
408
409    def pdf(self,x,*args,**kwds):
410        loc,scale=kwds['loc'],kwds['scale']
411        x = (x-loc)/scale
412        return 1.0/np.pi/(1.0+x*x) / scale
413       
414cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
415
416Cauchy distribution
417
418cauchy.pdf(x) = 1/(pi*(1+x**2))
419
420This is the t distribution with one degree of freedom.
421""")
422   
423   
424#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
425
426
427class fcjde_gen(st.rv_continuous):
428    """
429    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
430    Ref: J. Appl. Cryst. (1994) 27, 892-900.
431
432    :param x: array -1 to 1
433    :param t: 2-theta position of peak
434    :param s: sum(S/L,H/L); S: sample height, H: detector opening,
435      L: sample to detector opening distance
436    :param dx: 2-theta step size in deg
437
438    :returns: for fcj.pdf
439
440     * T = x*dx+t
441     * s = S/L+H/L
442     * if x < 0::
443
444        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|
445
446     * if x >= 0: fcj.pdf = 0   
447    """
448    def _pdf(self,x,t,s,dx):
449        T = dx*x+t
450        ax2 = abs(npcosd(T))
451        ax = ax2**2
452        bx = npcosd(t)**2
453        bx = np.where(ax>bx,bx,ax)
454        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
455        fx = np.where(fx > 0.,fx,0.0)
456        return fx
457             
458    def pdf(self,x,*args,**kwds):
459        loc=kwds['loc']
460        return self._pdf(x-loc,*args)
461       
462fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
463               
464def getWidthsCW(pos,sig,gam,shl):
465    '''Compute the peak widths used for computing the range of a peak
466    for constant wavelength data. On low-angle side, 10 FWHM are used,
467    on high-angle side 15 are used, low angle side extended by axial divergence
468    (for peaks above 90 deg, these are reversed.)
469    '''
470    widths = [np.sqrt(sig)/100.,gam/100.]
471    fwhm = 2.355*widths[0]+widths[1]
472    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
473    fmax = 15.0*fwhm
474    if pos > 90:
475        fmin,fmax = [fmax,fmin]         
476    return widths,fmin,fmax
477   
478def getWidthsTOF(pos,alp,bet,sig,gam):
479    '''Compute the peak widths used for computing the range of a peak
480    for constant wavelength data. 10 FWHM are used on both sides each
481    extended by exponential coeff.
482    '''
483    lnf = 3.3      # =log(0.001/2)
484    widths = [np.sqrt(sig),gam]
485    fwhm = 2.355*widths[0]+2.*widths[1]
486    fmin = 10.*fwhm*(1.+1./alp)   
487    fmax = 10.*fwhm*(1.+1./bet)
488    return widths,fmin,fmax
489   
490def getFWHM(pos,Inst):
491    '''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
492    via getgamFW(g,s).
493   
494    :param pos: float peak position in deg 2-theta or tof in musec
495    :param Inst: dict instrument parameters
496   
497    :returns float: total FWHM of pseudoVoigt in deg or musec
498    ''' 
499   
500    sig = lambda Th,U,V,W: np.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W))
501    sigTOF = lambda dsp,S0,S1,S2,Sq:  S0+S1*dsp**2+S2*dsp**4+Sq/dsp**2
502    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))
503    gamTOF = lambda dsp,X,Y: X*dsp+Y*dsp**2
504    if 'C' in Inst['Type'][0]:
505        s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
506        g = gam(pos/2.,Inst['X'][1],Inst['Y'][1])
507        return getgamFW(g,s)/100.  #returns FWHM in deg
508    else:
509        dsp = pos/Inst['difC'][0]
510        s = sigTOF(dsp,Inst['sig-0'][1],Inst['sig-1'][1],Inst['sig-2'][1],Inst['sig-q'][1])
511        g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1])
512        return getgamFW(g,s)
513   
514def getgamFW(g,s):
515    '''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
516    lambda fxn needs FWHM for both Gaussian & Lorentzian components
517   
518    :param g: float Lorentzian gamma = FWHM(L)
519    :param s: float Gaussian sig
520   
521    :returns float: total FWHM of pseudoVoigt
522    ''' 
523    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.)
524    return gamFW(2.35482*s,g)   #sqrt(8ln2)*sig = FWHM(G)
525               
526def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
527    'needs a doc string'
528    DX = xdata[1]-xdata[0]
529    widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
530    x = np.linspace(pos-fmin,pos+fmin,256)
531    dx = x[1]-x[0]
532    Norm = norm.pdf(x,loc=pos,scale=widths[0])
533    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
534    arg = [pos,shl/57.2958,dx,]
535    FCJ = fcjde.pdf(x,*arg,loc=pos)
536    if len(np.nonzero(FCJ)[0])>5:
537        z = np.column_stack([Norm,Cauchy,FCJ]).T
538        Z = fft(z)
539        Df = ifft(Z.prod(axis=0)).real
540    else:
541        z = np.column_stack([Norm,Cauchy]).T
542        Z = fft(z)
543        Df = fftshift(ifft(Z.prod(axis=0))).real
544    Df /= np.sum(Df)
545    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
546    return intens*Df(xdata)*DX/dx
547
548def getBackground(pfx,parmDict,bakType,dataType,xdata):
549    'needs a doc string'
550    if 'T' in dataType:
551        q = 2.*np.pi*parmDict[pfx+'difC']/xdata
552    elif 'C' in dataType:
553        wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
554        q = 2.*np.pi*npsind(xdata/2.)/wave
555    yb = np.zeros_like(xdata)
556    nBak = 0
557    cw = np.diff(xdata)
558    cw = np.append(cw,cw[-1])
559    sumBk = [0.,0.,0]
560    while True:
561        key = pfx+'Back;'+str(nBak)
562        if key in parmDict:
563            nBak += 1
564        else:
565            break
566#empirical functions
567    if bakType in ['chebyschev','cosine']:
568        dt = xdata[-1]-xdata[0]   
569        for iBak in range(nBak):
570            key = pfx+'Back;'+str(iBak)
571            if bakType == 'chebyschev':
572                ybi = parmDict[key]*(2.*(xdata-xdata[0])/dt-1.)**iBak
573            elif bakType == 'cosine':
574                ybi = parmDict[key]*npcosd(xdata*iBak)
575            yb += ybi
576        sumBk[0] = np.sum(yb)
577    elif bakType in ['Q^2 power series','Q^-2 power series']:
578        QT = 1.
579        yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
580        for iBak in range(nBak-1):
581            key = pfx+'Back;'+str(iBak+1)
582            if '-2' in bakType:
583                QT *= (iBak+1)*q**-2
584            else:
585                QT *= q**2/(iBak+1)
586            yb += QT*parmDict[key]
587        sumBk[0] = np.sum(yb)
588    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
589        if nBak == 1:
590            yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
591        elif nBak == 2:
592            dX = xdata[-1]-xdata[0]
593            T2 = (xdata-xdata[0])/dX
594            T1 = 1.0-T2
595            yb = parmDict[pfx+'Back;0']*T1+parmDict[pfx+'Back;1']*T2
596        else:
597            if bakType == 'lin interpolate':
598                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
599            elif bakType == 'inv interpolate':
600                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
601            elif bakType == 'log interpolate':
602                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
603            bakPos[0] = xdata[0]
604            bakPos[-1] = xdata[-1]
605            bakVals = np.zeros(nBak)
606            for i in range(nBak):
607                bakVals[i] = parmDict[pfx+'Back;'+str(i)]
608            bakInt = si.interp1d(bakPos,bakVals,'linear')
609            yb = bakInt(xdata)
610        sumBk[0] = np.sum(yb)
611#Debye function       
612    if pfx+'difC' in parmDict:
613        ff = 1.
614    else:       
615        try:
616            wave = parmDict[pfx+'Lam']
617        except KeyError:
618            wave = parmDict[pfx+'Lam1']
619        SQ = (q/(4.*np.pi))**2
620        FF = G2elem.GetFormFactorCoeff('Si')[0]
621        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
622    iD = 0       
623    while True:
624        try:
625            dbA = parmDict[pfx+'DebyeA;'+str(iD)]
626            dbR = parmDict[pfx+'DebyeR;'+str(iD)]
627            dbU = parmDict[pfx+'DebyeU;'+str(iD)]
628            ybi = ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
629            yb += ybi
630            sumBk[1] += np.sum(ybi)
631            iD += 1       
632        except KeyError:
633            break
634#peaks
635    iD = 0
636    while True:
637        try:
638            pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
639            pkI = parmDict[pfx+'BkPkint;'+str(iD)]
640            pkS = parmDict[pfx+'BkPksig;'+str(iD)]
641            pkG = parmDict[pfx+'BkPkgam;'+str(iD)]
642            if 'C' in dataType:
643                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
644            else: #'T'OF
645                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
646            iBeg = np.searchsorted(xdata,pkP-fmin)
647            iFin = np.searchsorted(xdata,pkP+fmax)
648            lenX = len(xdata)
649            if not iBeg:
650                iFin = np.searchsorted(xdata,pkP+fmax)
651            elif iBeg == lenX:
652                iFin = iBeg
653            else:
654                iFin = np.searchsorted(xdata,pkP+fmax)
655            if 'C' in dataType:
656                ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])
657                yb[iBeg:iFin] += ybi
658            else:   #'T'OF
659                ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
660                yb[iBeg:iFin] += ybi
661            sumBk[2] += np.sum(ybi)
662            iD += 1       
663        except KeyError:
664            break
665        except ValueError:
666            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
667            break       
668    return yb,sumBk
669   
670def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata):
671    'needs a doc string'
672    if 'T' in dataType:
673        q = 2.*np.pi*parmDict[hfx+'difC']/xdata
674    elif 'C' in dataType:
675        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
676        q = 2.*np.pi*npsind(xdata/2.)/wave
677    nBak = 0
678    while True:
679        key = hfx+'Back;'+str(nBak)
680        if key in parmDict:
681            nBak += 1
682        else:
683            break
684    dydb = np.zeros(shape=(nBak,len(xdata)))
685    dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
686    dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
687    cw = np.diff(xdata)
688    cw = np.append(cw,cw[-1])
689
690    if bakType in ['chebyschev','cosine']:
691        dt = xdata[-1]-xdata[0]   
692        for iBak in range(nBak):   
693            if bakType == 'chebyschev':
694                dydb[iBak] = (2.*(xdata-xdata[0])/dt-1.)**iBak
695            elif bakType == 'cosine':
696                dydb[iBak] = npcosd(xdata*iBak)
697    elif bakType in ['Q^2 power series','Q^-2 power series']:
698        QT = 1.
699        dydb[0] = np.ones_like(xdata)
700        for iBak in range(nBak-1):
701            if '-2' in bakType:
702                QT *= (iBak+1)*q**-2
703            else:
704                QT *= q**2/(iBak+1)
705            dydb[iBak+1] = QT
706    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
707        if nBak == 1:
708            dydb[0] = np.ones_like(xdata)
709        elif nBak == 2:
710            dX = xdata[-1]-xdata[0]
711            T2 = (xdata-xdata[0])/dX
712            T1 = 1.0-T2
713            dydb = [T1,T2]
714        else:
715            if bakType == 'lin interpolate':
716                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
717            elif bakType == 'inv interpolate':
718                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
719            elif bakType == 'log interpolate':
720                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
721            bakPos[0] = xdata[0]
722            bakPos[-1] = xdata[-1]
723            for i,pos in enumerate(bakPos):
724                if i == 0:
725                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
726                elif i == len(bakPos)-1:
727                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
728                else:
729                    dydb[i] = np.where(xdata>bakPos[i],
730                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
731                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
732    if hfx+'difC' in parmDict:
733        ff = 1.
734    else:
735        try:
736            wave = parmDict[hfx+'Lam']
737        except KeyError:
738            wave = parmDict[hfx+'Lam1']
739        q = 4.0*np.pi*npsind(xdata/2.0)/wave
740        SQ = (q/(4*np.pi))**2
741        FF = G2elem.GetFormFactorCoeff('Si')[0]
742        ff = np.array(G2elem.ScatFac(FF,SQ)[0])
743    iD = 0       
744    while True:
745        try:
746            if hfx+'difC' in parmDict:
747                q = 2*np.pi*parmDict[hfx+'difC']/xdata
748            dbA = parmDict[hfx+'DebyeA;'+str(iD)]
749            dbR = parmDict[hfx+'DebyeR;'+str(iD)]
750            dbU = parmDict[hfx+'DebyeU;'+str(iD)]
751            sqr = np.sin(q*dbR)/(q*dbR)
752            cqr = np.cos(q*dbR)
753            temp = np.exp(-dbU*q**2)
754            dyddb[3*iD] = ff*sqr*temp/(np.pi*cw)
755            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(np.pi*dbR*cw)
756            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2/(np.pi*cw)
757            iD += 1       #ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
758        except KeyError:
759            break
760    iD = 0
761    while True:
762        try:
763            pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
764            pkI = parmDict[hfx+'BkPkint;'+str(iD)]
765            pkS = parmDict[hfx+'BkPksig;'+str(iD)]
766            pkG = parmDict[hfx+'BkPkgam;'+str(iD)]
767            if 'C' in dataType:
768                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
769            else: #'T'OF
770                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
771            iBeg = np.searchsorted(xdata,pkP-fmin)
772            iFin = np.searchsorted(xdata,pkP+fmax)
773            lenX = len(xdata)
774            if not iBeg:
775                iFin = np.searchsorted(xdata,pkP+fmax)
776            elif iBeg == lenX:
777                iFin = iBeg
778            else:
779                iFin = np.searchsorted(xdata,pkP+fmax)
780            if 'C' in dataType:
781                Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
782                dydpk[4*iD][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdp
783                dydpk[4*iD+1][iBeg:iFin] += 100.*cw[iBeg:iFin]*Df
784                dydpk[4*iD+2][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFds
785                dydpk[4*iD+3][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdg
786            else:   #'T'OF
787                Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
788                dydpk[4*iD][iBeg:iFin] += pkI*dFdp
789                dydpk[4*iD+1][iBeg:iFin] += Df
790                dydpk[4*iD+2][iBeg:iFin] += pkI*dFds
791                dydpk[4*iD+3][iBeg:iFin] += pkI*dFdg
792            iD += 1       
793        except KeyError:
794            break
795        except ValueError:
796            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
797            break       
798    return dydb,dyddb,dydpk
799
800#use old fortran routine
801def getFCJVoigt3(pos,sig,gam,shl,xdata):
802    'needs a doc string'
803   
804    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
805#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
806    Df /= np.sum(Df)
807    return Df
808
809def getdFCJVoigt3(pos,sig,gam,shl,xdata):
810    'needs a doc string'
811   
812    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
813#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
814    sumDf = np.sum(Df)
815    return Df,dFdp,dFds,dFdg,dFdsh
816
817def getPsVoigt(pos,sig,gam,xdata):
818    'needs a doc string'
819   
820    Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
821    Df /= np.sum(Df)
822    return Df
823
824def getdPsVoigt(pos,sig,gam,xdata):
825    'needs a doc string'
826   
827    Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
828    sumDf = np.sum(Df)
829    return Df,dFdp,dFds,dFdg
830
831def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
832    'needs a doc string'
833    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
834    Df /= np.sum(Df)
835    return Df 
836   
837def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
838    'needs a doc string'
839    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
840    sumDf = np.sum(Df)
841    return Df,dFdp,dFda,dFdb,dFds,dFdg   
842
843def ellipseSize(H,Sij,GB):
844    'needs a doc string'
845    HX = np.inner(H.T,GB)
846    lenHX = np.sqrt(np.sum(HX**2))
847    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
848    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
849    lenR = np.sqrt(np.sum(R**2))
850    return lenR
851
852def ellipseSizeDerv(H,Sij,GB):
853    'needs a doc string'
854    lenR = ellipseSize(H,Sij,GB)
855    delt = 0.001
856    dRdS = np.zeros(6)
857    for i in range(6):
858        Sij[i] -= delt
859        lenM = ellipseSize(H,Sij,GB)
860        Sij[i] += 2.*delt
861        lenP = ellipseSize(H,Sij,GB)
862        Sij[i] -= delt
863        dRdS[i] = (lenP-lenM)/(2.*delt)
864    return lenR,dRdS
865
866def getHKLpeak(dmin,SGData,A,Inst=None):
867    'needs a doc string'
868    HKL = G2lat.GenHLaue(dmin,SGData,A)       
869    HKLs = []
870    for h,k,l,d in HKL:
871        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
872        if not ext:
873            if Inst == None:
874                HKLs.append([h,k,l,d,0,-1])
875            else:
876                HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
877    return HKLs
878
879def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
880    'needs a doc string'
881    HKLs = []
882    vec = np.array(Vec)
883    vstar = np.sqrt(G2lat.calc_rDsq(vec,A))     #find extra needed for -n SS reflections
884    dvec = 1./(maxH*vstar+1./dmin)
885    HKL = G2lat.GenHLaue(dvec,SGData,A)       
886    SSdH = [vec*h for h in range(-maxH,maxH+1)]
887    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
888    for h,k,l,d in HKL:
889        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
890        if not ext and d >= dmin:
891            HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
892        for dH in SSdH:
893            if dH:
894                DH = SSdH[dH]
895                H = [h+DH[0],k+DH[1],l+DH[2]]
896                d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
897                if d >= dmin:
898                    HKLM = np.array([h,k,l,dH])
899                    if G2spc.checkSSextc(HKLM,SSGData):
900                        HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])   
901#    GSASIIpath.IPyBreak()
902    return G2lat.sortHKLd(HKLs,True,True,True)
903
904def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
905    'needs a doc string'
906   
907    yb = getBackground('',parmDict,bakType,dataType,xdata)[0]
908    yc = np.zeros_like(yb)
909    cw = np.diff(xdata)
910    cw = np.append(cw,cw[-1])
911    if 'C' in dataType:
912        shl = max(parmDict['SH/L'],0.002)
913        Ka2 = False
914        if 'Lam1' in parmDict.keys():
915            Ka2 = True
916            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
917            kRatio = parmDict['I(L2)/I(L1)']
918        iPeak = 0
919        while True:
920            try:
921                pos = parmDict['pos'+str(iPeak)]
922                tth = (pos-parmDict['Zero'])
923                intens = parmDict['int'+str(iPeak)]
924                sigName = 'sig'+str(iPeak)
925                if sigName in varyList:
926                    sig = parmDict[sigName]
927                else:
928                    sig = G2mth.getCWsig(parmDict,tth)
929                sig = max(sig,0.001)          #avoid neg sigma^2
930                gamName = 'gam'+str(iPeak)
931                if gamName in varyList:
932                    gam = parmDict[gamName]
933                else:
934                    gam = G2mth.getCWgam(parmDict,tth)
935                gam = max(gam,0.001)             #avoid neg gamma
936                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
937                iBeg = np.searchsorted(xdata,pos-fmin)
938                iFin = np.searchsorted(xdata,pos+fmin)
939                if not iBeg+iFin:       #peak below low limit
940                    iPeak += 1
941                    continue
942                elif not iBeg-iFin:     #peak above high limit
943                    return yb+yc
944                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
945                if Ka2:
946                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
947                    iBeg = np.searchsorted(xdata,pos2-fmin)
948                    iFin = np.searchsorted(xdata,pos2+fmin)
949                    if iBeg-iFin:
950                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
951                iPeak += 1
952            except KeyError:        #no more peaks to process
953                return yb+yc
954    else:
955        Pdabc = parmDict['Pdabc']
956        difC = parmDict['difC']
957        iPeak = 0
958        while True:
959            try:
960                pos = parmDict['pos'+str(iPeak)]               
961                tof = pos-parmDict['Zero']
962                dsp = tof/difC
963                intens = parmDict['int'+str(iPeak)]
964                alpName = 'alp'+str(iPeak)
965                if alpName in varyList:
966                    alp = parmDict[alpName]
967                else:
968                    if len(Pdabc):
969                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
970                    else:
971                        alp = G2mth.getTOFalpha(parmDict,dsp)
972                alp = max(0.0001,alp)
973                betName = 'bet'+str(iPeak)
974                if betName in varyList:
975                    bet = parmDict[betName]
976                else:
977                    if len(Pdabc):
978                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
979                    else:
980                        bet = G2mth.getTOFbeta(parmDict,dsp)
981                bet = max(0.0001,bet)
982                sigName = 'sig'+str(iPeak)
983                if sigName in varyList:
984                    sig = parmDict[sigName]
985                else:
986                    sig = G2mth.getTOFsig(parmDict,dsp)
987                gamName = 'gam'+str(iPeak)
988                if gamName in varyList:
989                    gam = parmDict[gamName]
990                else:
991                    gam = G2mth.getTOFgamma(parmDict,dsp)
992                gam = max(gam,0.001)             #avoid neg gamma
993                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
994                iBeg = np.searchsorted(xdata,pos-fmin)
995                iFin = np.searchsorted(xdata,pos+fmax)
996                lenX = len(xdata)
997                if not iBeg:
998                    iFin = np.searchsorted(xdata,pos+fmax)
999                elif iBeg == lenX:
1000                    iFin = iBeg
1001                else:
1002                    iFin = np.searchsorted(xdata,pos+fmax)
1003                if not iBeg+iFin:       #peak below low limit
1004                    iPeak += 1
1005                    continue
1006                elif not iBeg-iFin:     #peak above high limit
1007                    return yb+yc
1008                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1009                iPeak += 1
1010            except KeyError:        #no more peaks to process
1011                return yb+yc
1012           
1013def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
1014    'needs a doc string'
1015# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1016    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1017    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,dataType,xdata)
1018    if 'Back;0' in varyList:            #background derivs are in front if present
1019        dMdv[0:len(dMdb)] = dMdb
1020    names = ['DebyeA','DebyeR','DebyeU']
1021    for name in varyList:
1022        if 'Debye' in name:
1023            parm,id = name.split(';')
1024            ip = names.index(parm)
1025            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
1026    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1027    for name in varyList:
1028        if 'BkPk' in name:
1029            parm,id = name.split(';')
1030            ip = names.index(parm)
1031            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
1032    cw = np.diff(xdata)
1033    cw = np.append(cw,cw[-1])
1034    if 'C' in dataType:
1035        shl = max(parmDict['SH/L'],0.002)
1036        Ka2 = False
1037        if 'Lam1' in parmDict.keys():
1038            Ka2 = True
1039            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1040            kRatio = parmDict['I(L2)/I(L1)']
1041        iPeak = 0
1042        while True:
1043            try:
1044                pos = parmDict['pos'+str(iPeak)]
1045                tth = (pos-parmDict['Zero'])
1046                intens = parmDict['int'+str(iPeak)]
1047                sigName = 'sig'+str(iPeak)
1048                if sigName in varyList:
1049                    sig = parmDict[sigName]
1050                    dsdU = dsdV = dsdW = 0
1051                else:
1052                    sig = G2mth.getCWsig(parmDict,tth)
1053                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1054                sig = max(sig,0.001)          #avoid neg sigma
1055                gamName = 'gam'+str(iPeak)
1056                if gamName in varyList:
1057                    gam = parmDict[gamName]
1058                    dgdX = dgdY = 0
1059                else:
1060                    gam = G2mth.getCWgam(parmDict,tth)
1061                    dgdX,dgdY = G2mth.getCWgamDeriv(tth)
1062                gam = max(gam,0.001)             #avoid neg gamma
1063                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1064                iBeg = np.searchsorted(xdata,pos-fmin)
1065                iFin = np.searchsorted(xdata,pos+fmin)
1066                if not iBeg+iFin:       #peak below low limit
1067                    iPeak += 1
1068                    continue
1069                elif not iBeg-iFin:     #peak above high limit
1070                    break
1071                dMdpk = np.zeros(shape=(6,len(xdata)))
1072                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1073                for i in range(1,5):
1074                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
1075                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
1076                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1077                if Ka2:
1078                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1079                    iBeg = np.searchsorted(xdata,pos2-fmin)
1080                    iFin = np.searchsorted(xdata,pos2+fmin)
1081                    if iBeg-iFin:
1082                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1083                        for i in range(1,5):
1084                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
1085                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
1086                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
1087                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1088                for parmName in ['pos','int','sig','gam']:
1089                    try:
1090                        idx = varyList.index(parmName+str(iPeak))
1091                        dMdv[idx] = dervDict[parmName]
1092                    except ValueError:
1093                        pass
1094                if 'U' in varyList:
1095                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1096                if 'V' in varyList:
1097                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1098                if 'W' in varyList:
1099                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1100                if 'X' in varyList:
1101                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1102                if 'Y' in varyList:
1103                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1104                if 'SH/L' in varyList:
1105                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
1106                if 'I(L2)/I(L1)' in varyList:
1107                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1108                iPeak += 1
1109            except KeyError:        #no more peaks to process
1110                break
1111    else:
1112        Pdabc = parmDict['Pdabc']
1113        difC = parmDict['difC']
1114        iPeak = 0
1115        while True:
1116            try:
1117                pos = parmDict['pos'+str(iPeak)]               
1118                tof = pos-parmDict['Zero']
1119                dsp = tof/difC
1120                intens = parmDict['int'+str(iPeak)]
1121                alpName = 'alp'+str(iPeak)
1122                if alpName in varyList:
1123                    alp = parmDict[alpName]
1124                else:
1125                    if len(Pdabc):
1126                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1127                        dad0 = 0
1128                    else:
1129                        alp = G2mth.getTOFalpha(parmDict,dsp)
1130                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1131                betName = 'bet'+str(iPeak)
1132                if betName in varyList:
1133                    bet = parmDict[betName]
1134                else:
1135                    if len(Pdabc):
1136                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1137                        dbdb0 = dbdb1 = dbdb2 = 0
1138                    else:
1139                        bet = G2mth.getTOFbeta(parmDict,dsp)
1140                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1141                sigName = 'sig'+str(iPeak)
1142                if sigName in varyList:
1143                    sig = parmDict[sigName]
1144                    dsds0 = dsds1 = dsds2 = dsds3 = 0
1145                else:
1146                    sig = G2mth.getTOFsig(parmDict,dsp)
1147                    dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1148                gamName = 'gam'+str(iPeak)
1149                if gamName in varyList:
1150                    gam = parmDict[gamName]
1151                    dsdX = dsdY = 0
1152                else:
1153                    gam = G2mth.getTOFgamma(parmDict,dsp)
1154                    dsdX,dsdY = G2mth.getTOFgammaDeriv(dsp)
1155                gam = max(gam,0.001)             #avoid neg gamma
1156                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1157                iBeg = np.searchsorted(xdata,pos-fmin)
1158                lenX = len(xdata)
1159                if not iBeg:
1160                    iFin = np.searchsorted(xdata,pos+fmax)
1161                elif iBeg == lenX:
1162                    iFin = iBeg
1163                else:
1164                    iFin = np.searchsorted(xdata,pos+fmax)
1165                if not iBeg+iFin:       #peak below low limit
1166                    iPeak += 1
1167                    continue
1168                elif not iBeg-iFin:     #peak above high limit
1169                    break
1170                dMdpk = np.zeros(shape=(7,len(xdata)))
1171                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1172                for i in range(1,6):
1173                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1174                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1175                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1176                for parmName in ['pos','int','alp','bet','sig','gam']:
1177                    try:
1178                        idx = varyList.index(parmName+str(iPeak))
1179                        dMdv[idx] = dervDict[parmName]
1180                    except ValueError:
1181                        pass
1182                if 'alpha' in varyList:
1183                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1184                if 'beta-0' in varyList:
1185                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1186                if 'beta-1' in varyList:
1187                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1188                if 'beta-q' in varyList:
1189                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1190                if 'sig-0' in varyList:
1191                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1192                if 'sig-1' in varyList:
1193                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1194                if 'sig-2' in varyList:
1195                    dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1196                if 'sig-q' in varyList:
1197                    dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1198                if 'X' in varyList:
1199                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1200                if 'Y' in varyList:
1201                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
1202                iPeak += 1
1203            except KeyError:        #no more peaks to process
1204                break
1205    return dMdv
1206       
1207def Dict2Values(parmdict, varylist):
1208    '''Use before call to leastsq to setup list of values for the parameters
1209    in parmdict, as selected by key in varylist'''
1210    return [parmdict[key] for key in varylist] 
1211   
1212def Values2Dict(parmdict, varylist, values):
1213    ''' Use after call to leastsq to update the parameter dictionary with
1214    values corresponding to keys in varylist'''
1215    parmdict.update(zip(varylist,values))
1216   
1217def SetBackgroundParms(Background):
1218    'needs a doc string'
1219    if len(Background) == 1:            # fix up old backgrounds
1220        BackGround.append({'nDebye':0,'debyeTerms':[]})
1221    bakType,bakFlag = Background[0][:2]
1222    backVals = Background[0][3:]
1223    backNames = ['Back;'+str(i) for i in range(len(backVals))]
1224    Debye = Background[1]           #also has background peaks stuff
1225    backDict = dict(zip(backNames,backVals))
1226    backVary = []
1227    if bakFlag:
1228        backVary = backNames
1229
1230    backDict['nDebye'] = Debye['nDebye']
1231    debyeDict = {}
1232    debyeList = []
1233    for i in range(Debye['nDebye']):
1234        debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1235        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1236        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1237    debyeVary = []
1238    for item in debyeList:
1239        if item[1]:
1240            debyeVary.append(item[0])
1241    backDict.update(debyeDict)
1242    backVary += debyeVary
1243
1244    backDict['nPeaks'] = Debye['nPeaks']
1245    peaksDict = {}
1246    peaksList = []
1247    for i in range(Debye['nPeaks']):
1248        peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1249        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1250        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1251    peaksVary = []
1252    for item in peaksList:
1253        if item[1]:
1254            peaksVary.append(item[0])
1255    backDict.update(peaksDict)
1256    backVary += peaksVary   
1257    return bakType,backDict,backVary
1258   
1259def DoCalibInst(IndexPeaks,Inst):
1260   
1261    def SetInstParms():
1262        dataType = Inst['Type'][0]
1263        insVary = []
1264        insNames = []
1265        insVals = []
1266        for parm in Inst:
1267            insNames.append(parm)
1268            insVals.append(Inst[parm][1])
1269            if parm in ['Lam','difC','difA','difB','Zero',]:
1270                if Inst[parm][2]:
1271                    insVary.append(parm)
1272        instDict = dict(zip(insNames,insVals))
1273        return dataType,instDict,insVary
1274       
1275    def GetInstParms(parmDict,Inst,varyList):
1276        for name in Inst:
1277            Inst[name][1] = parmDict[name]
1278       
1279    def InstPrint(Inst,sigDict):
1280        print 'Instrument Parameters:'
1281        if 'C' in Inst['Type'][0]:
1282            ptfmt = "%12.6f"
1283        else:
1284            ptfmt = "%12.3f"
1285        ptlbls = 'names :'
1286        ptstr =  'values:'
1287        sigstr = 'esds  :'
1288        for parm in Inst:
1289            if parm in  ['Lam','difC','difA','difB','Zero',]:
1290                ptlbls += "%s" % (parm.center(12))
1291                ptstr += ptfmt % (Inst[parm][1])
1292                if parm in sigDict:
1293                    sigstr += ptfmt % (sigDict[parm])
1294                else:
1295                    sigstr += 12*' '
1296        print ptlbls
1297        print ptstr
1298        print sigstr
1299       
1300    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1301        parmDict.update(zip(varyList,values))
1302        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1303
1304    peakPos = []
1305    peakDsp = []
1306    peakWt = []
1307    for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1308        if peak[2] and peak[3] and sig > 0.:
1309            peakPos.append(peak[0])
1310            peakDsp.append(peak[-1])    #d-calc
1311            peakWt.append(peak[-1]**2/sig**2)   #weight by d**2
1312    peakPos = np.array(peakPos)
1313    peakDsp = np.array(peakDsp)
1314    peakWt = np.array(peakWt)
1315    dataType,insDict,insVary = SetInstParms()
1316    parmDict = {}
1317    parmDict.update(insDict)
1318    varyList = insVary
1319    if not len(varyList):
1320        print '**** ERROR - nothing to refine! ****'
1321        return False
1322    while True:
1323        begin = time.time()
1324        values =  np.array(Dict2Values(parmDict, varyList))
1325        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1326            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1327        ncyc = int(result[2]['nfev']/2)
1328        runtime = time.time()-begin   
1329        chisq = np.sum(result[2]['fvec']**2)
1330        Values2Dict(parmDict, varyList, result[0])
1331        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
1332        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(peakPos),' Number of parameters: ',len(varyList)
1333        print 'calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1334        print 'chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)
1335        try:
1336            sig = np.sqrt(np.diag(result[1])*GOF)
1337            if np.any(np.isnan(sig)):
1338                print '*** Least squares aborted - some invalid esds possible ***'
1339            break                   #refinement succeeded - finish up!
1340        except ValueError:          #result[1] is None on singular matrix
1341            print '**** Refinement failed - singular matrix ****'
1342       
1343    sigDict = dict(zip(varyList,sig))
1344    GetInstParms(parmDict,Inst,varyList)
1345    InstPrint(Inst,sigDict)
1346    return True
1347           
1348def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
1349    '''Called to perform a peak fit, refining the selected items in the peak
1350    table as well as selected items in the background.
1351
1352    :param str FitPgm: type of fit to perform. At present this is ignored.
1353    :param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1354      four values followed by a refine flag where the values are: position, intensity,
1355      sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1356      tree entry, dict item "peaks"
1357    :param list Background: describes the background. List with two items.
1358      Item 0 specifies a background model and coefficients. Item 1 is a dict.
1359      From the Histogram/Background tree entry.
1360    :param list Limits: min and max x-value to use
1361    :param dict Inst: Instrument parameters
1362    :param dict Inst2: more Instrument parameters
1363    :param numpy.array data: a 5xn array. data[0] is the x-values,
1364      data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1365      calc, background and difference intensities, respectively.
1366    :param list prevVaryList: Used in sequential refinements to override the
1367      variable list. Defaults as an empty list.
1368    :param bool oneCycle: True if only one cycle of fitting should be performed
1369    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1370      and derivType = controls['deriv type']. If None default values are used.
1371    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1372      Defaults to None, which means no updates are done.
1373    '''
1374    def GetBackgroundParms(parmList,Background):
1375        iBak = 0
1376        while True:
1377            try:
1378                bakName = 'Back;'+str(iBak)
1379                Background[0][iBak+3] = parmList[bakName]
1380                iBak += 1
1381            except KeyError:
1382                break
1383        iDb = 0
1384        while True:
1385            names = ['DebyeA;','DebyeR;','DebyeU;']
1386            try:
1387                for i,name in enumerate(names):
1388                    val = parmList[name+str(iDb)]
1389                    Background[1]['debyeTerms'][iDb][2*i] = val
1390                iDb += 1
1391            except KeyError:
1392                break
1393        iDb = 0
1394        while True:
1395            names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1396            try:
1397                for i,name in enumerate(names):
1398                    val = parmList[name+str(iDb)]
1399                    Background[1]['peaksList'][iDb][2*i] = val
1400                iDb += 1
1401            except KeyError:
1402                break
1403               
1404    def BackgroundPrint(Background,sigDict):
1405        print 'Background coefficients for',Background[0][0],'function'
1406        ptfmt = "%12.5f"
1407        ptstr =  'value: '
1408        sigstr = 'esd  : '
1409        for i,back in enumerate(Background[0][3:]):
1410            ptstr += ptfmt % (back)
1411            if Background[0][1]:
1412                sigstr += ptfmt % (sigDict['Back;'+str(i)])
1413            if len(ptstr) > 75:
1414                print ptstr
1415                if Background[0][1]: print sigstr
1416                ptstr =  'value: '
1417                sigstr = 'esd  : '
1418        if len(ptstr) > 8:
1419            print ptstr
1420            if Background[0][1]: print sigstr
1421
1422        if Background[1]['nDebye']:
1423            parms = ['DebyeA;','DebyeR;','DebyeU;']
1424            print 'Debye diffuse scattering coefficients'
1425            ptfmt = "%12.5f"
1426            print ' term       DebyeA       esd        DebyeR       esd        DebyeU        esd'
1427            for term in range(Background[1]['nDebye']):
1428                line = ' term %d'%(term)
1429                for ip,name in enumerate(parms):
1430                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1431                    if name+str(term) in sigDict:
1432                        line += ptfmt%(sigDict[name+str(term)])
1433                print line
1434        if Background[1]['nPeaks']:
1435            print 'Coefficients for Background Peaks'
1436            ptfmt = "%15.3f"
1437            for j,pl in enumerate(Background[1]['peaksList']):
1438                names =  'peak %3d:'%(j+1)
1439                ptstr =  'values  :'
1440                sigstr = 'esds    :'
1441                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1442                    val = pl[2*i]
1443                    prm = lbl+";"+str(j)
1444                    names += '%15s'%(prm)
1445                    ptstr += ptfmt%(val)
1446                    if prm in sigDict:
1447                        sigstr += ptfmt%(sigDict[prm])
1448                    else:
1449                        sigstr += " "*15
1450                print names
1451                print ptstr
1452                print sigstr
1453                           
1454    def SetInstParms(Inst):
1455        dataType = Inst['Type'][0]
1456        insVary = []
1457        insNames = []
1458        insVals = []
1459        for parm in Inst:
1460            insNames.append(parm)
1461            insVals.append(Inst[parm][1])
1462            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1463                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1464                    insVary.append(parm)
1465        instDict = dict(zip(insNames,insVals))
1466        instDict['X'] = max(instDict['X'],0.01)
1467        instDict['Y'] = max(instDict['Y'],0.01)
1468        if 'SH/L' in instDict:
1469            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1470        return dataType,instDict,insVary
1471       
1472    def GetInstParms(parmDict,Inst,varyList,Peaks):
1473        for name in Inst:
1474            Inst[name][1] = parmDict[name]
1475        iPeak = 0
1476        while True:
1477            try:
1478                sigName = 'sig'+str(iPeak)
1479                pos = parmDict['pos'+str(iPeak)]
1480                if sigName not in varyList:
1481                    if 'C' in Inst['Type'][0]:
1482                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1483                    else:
1484                        dsp = G2lat.Pos2dsp(Inst,pos)
1485                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1486                gamName = 'gam'+str(iPeak)
1487                if gamName not in varyList:
1488                    if 'C' in Inst['Type'][0]:
1489                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1490                    else:
1491                        dsp = G2lat.Pos2dsp(Inst,pos)
1492                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1493                iPeak += 1
1494            except KeyError:
1495                break
1496       
1497    def InstPrint(Inst,sigDict):
1498        print 'Instrument Parameters:'
1499        ptfmt = "%12.6f"
1500        ptlbls = 'names :'
1501        ptstr =  'values:'
1502        sigstr = 'esds  :'
1503        for parm in Inst:
1504            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1505                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1506                ptlbls += "%s" % (parm.center(12))
1507                ptstr += ptfmt % (Inst[parm][1])
1508                if parm in sigDict:
1509                    sigstr += ptfmt % (sigDict[parm])
1510                else:
1511                    sigstr += 12*' '
1512        print ptlbls
1513        print ptstr
1514        print sigstr
1515
1516    def SetPeaksParms(dataType,Peaks):
1517        peakNames = []
1518        peakVary = []
1519        peakVals = []
1520        if 'C' in dataType:
1521            names = ['pos','int','sig','gam']
1522        else:
1523            names = ['pos','int','alp','bet','sig','gam']
1524        for i,peak in enumerate(Peaks):
1525            for j,name in enumerate(names):
1526                peakVals.append(peak[2*j])
1527                parName = name+str(i)
1528                peakNames.append(parName)
1529                if peak[2*j+1]:
1530                    peakVary.append(parName)
1531        return dict(zip(peakNames,peakVals)),peakVary
1532               
1533    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1534        if 'C' in Inst['Type'][0]:
1535            names = ['pos','int','sig','gam']
1536        else:   #'T'
1537            names = ['pos','int','alp','bet','sig','gam']
1538        for i,peak in enumerate(Peaks):
1539            pos = parmDict['pos'+str(i)]
1540            if 'difC' in Inst:
1541                dsp = pos/Inst['difC'][1]
1542            for j in range(len(names)):
1543                parName = names[j]+str(i)
1544                if parName in varyList:
1545                    peak[2*j] = parmDict[parName]
1546                elif 'alpha' in parName:
1547                    peak[2*j] = parmDict['alpha']/dsp
1548                elif 'beta' in parName:
1549                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1550                elif 'sig' in parName:
1551                    if 'C' in Inst['Type'][0]:
1552                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1553                    else:
1554                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1555                elif 'gam' in parName:
1556                    if 'C' in Inst['Type'][0]:
1557                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1558                    else:
1559                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1560                       
1561    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1562        print 'Peak coefficients:'
1563        if 'C' in dataType:
1564            names = ['pos','int','sig','gam']
1565        else:   #'T'
1566            names = ['pos','int','alp','bet','sig','gam']           
1567        head = 13*' '
1568        for name in names:
1569            if name in ['alp','bet']:
1570                head += name.center(8)+'esd'.center(8)
1571            else:
1572                head += name.center(10)+'esd'.center(10)
1573        print head
1574        if 'C' in dataType:
1575            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1576        else:
1577            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1578        for i,peak in enumerate(Peaks):
1579            ptstr =  ':'
1580            for j in range(len(names)):
1581                name = names[j]
1582                parName = name+str(i)
1583                ptstr += ptfmt[name] % (parmDict[parName])
1584                if parName in varyList:
1585                    ptstr += ptfmt[name] % (sigDict[parName])
1586                else:
1587                    if name in ['alp','bet']:
1588                        ptstr += 8*' '
1589                    else:
1590                        ptstr += 10*' '
1591            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1592               
1593    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1594        parmdict.update(zip(varylist,values))
1595        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1596           
1597    def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):       
1598        parmdict.update(zip(varylist,values))
1599        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1600        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1601        if dlg:
1602            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1603            if not GoOn:
1604                return -M           #abort!!
1605        return M
1606       
1607    if controls:
1608        Ftol = controls['min dM/M']
1609        derivType = controls['deriv type']
1610    else:
1611        Ftol = 0.0001
1612        derivType = 'analytic'
1613    if oneCycle:
1614        Ftol = 1.0
1615    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1616    yc *= 0.                            #set calcd ones to zero
1617    yb *= 0.
1618    yd *= 0.
1619    xBeg = np.searchsorted(x,Limits[0])
1620    xFin = np.searchsorted(x,Limits[1])
1621    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1622    dataType,insDict,insVary = SetInstParms(Inst)
1623    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1624    parmDict = {}
1625    parmDict.update(bakDict)
1626    parmDict.update(insDict)
1627    parmDict.update(peakDict)
1628    parmDict['Pdabc'] = []      #dummy Pdabc
1629    parmDict.update(Inst2)      #put in real one if there
1630    if prevVaryList:
1631        varyList = prevVaryList[:]
1632    else:
1633        varyList = bakVary+insVary+peakVary
1634    fullvaryList = varyList[:]
1635    while True:
1636        begin = time.time()
1637        values =  np.array(Dict2Values(parmDict, varyList))
1638        Rvals = {}
1639        badVary = []
1640        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1641               args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1642        ncyc = int(result[2]['nfev']/2)
1643        runtime = time.time()-begin   
1644        chisq = np.sum(result[2]['fvec']**2)
1645        Values2Dict(parmDict, varyList, result[0])
1646        Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1647        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1648        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1649        if ncyc:
1650            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1651        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1652        sig = [0]*len(varyList)
1653        if len(varyList) == 0: break  # if nothing was refined
1654        try:
1655            sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1656            if np.any(np.isnan(sig)):
1657                print '*** Least squares aborted - some invalid esds possible ***'
1658            break                   #refinement succeeded - finish up!
1659        except ValueError:          #result[1] is None on singular matrix
1660            print '**** Refinement failed - singular matrix ****'
1661            Ipvt = result[2]['ipvt']
1662            for i,ipvt in enumerate(Ipvt):
1663                if not np.sum(result[2]['fjac'],axis=1)[i]:
1664                    print 'Removing parameter: ',varyList[ipvt-1]
1665                    badVary.append(varyList[ipvt-1])
1666                    del(varyList[ipvt-1])
1667                    break
1668            else: # nothing removed
1669                break
1670    if dlg: dlg.Destroy()
1671    sigDict = dict(zip(varyList,sig))
1672    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]
1673    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1674    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1675    GetBackgroundParms(parmDict,Background)
1676    if bakVary: BackgroundPrint(Background,sigDict)
1677    GetInstParms(parmDict,Inst,varyList,Peaks)
1678    if insVary: InstPrint(Inst,sigDict)
1679    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1680    if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList)
1681    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
1682
1683def calcIncident(Iparm,xdata):
1684    'needs a doc string'
1685
1686    def IfunAdv(Iparm,xdata):
1687        Itype = Iparm['Itype']
1688        Icoef = Iparm['Icoeff']
1689        DYI = np.ones((12,xdata.shape[0]))
1690        YI = np.ones_like(xdata)*Icoef[0]
1691       
1692        x = xdata/1000.                 #expressions are in ms
1693        if Itype == 'Exponential':
1694            for i in [1,3,5,7,9]:
1695                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1696                YI += Icoef[i]*Eterm
1697                DYI[i] *= Eterm
1698                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
1699        elif 'Maxwell'in Itype:
1700            Eterm = np.exp(-Icoef[2]/x**2)
1701            DYI[1] = Eterm/x**5
1702            DYI[2] = -Icoef[1]*DYI[1]/x**2
1703            YI += (Icoef[1]*Eterm/x**5)
1704            if 'Exponential' in Itype:
1705                for i in range(3,12,2):
1706                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1707                    YI += Icoef[i]*Eterm
1708                    DYI[i] *= Eterm
1709                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
1710            else:   #Chebyschev
1711                T = (2./x)-1.
1712                Ccof = np.ones((12,xdata.shape[0]))
1713                Ccof[1] = T
1714                for i in range(2,12):
1715                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1716                for i in range(1,10):
1717                    YI += Ccof[i]*Icoef[i+2]
1718                    DYI[i+2] =Ccof[i]
1719        return YI,DYI
1720       
1721    Iesd = np.array(Iparm['Iesd'])
1722    Icovar = Iparm['Icovar']
1723    YI,DYI = IfunAdv(Iparm,xdata)
1724    YI = np.where(YI>0,YI,1.)
1725    WYI = np.zeros_like(xdata)
1726    vcov = np.zeros((12,12))
1727    k = 0
1728    for i in range(12):
1729        for j in range(i,12):
1730            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1731            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1732            k += 1
1733    M = np.inner(vcov,DYI.T)
1734    WYI = np.sum(M*DYI,axis=0)
1735    WYI = np.where(WYI>0.,WYI,0.)
1736    return YI,WYI
1737   
1738################################################################################
1739# Stacking fault simulation codes
1740################################################################################
1741
1742def GetStackParms(Layers):
1743   
1744    Parms = []
1745#cell parms
1746    cell = Layers['Cell'] 
1747    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
1748        Parms.append('cellA')
1749        Parms.append('cellC')
1750    else:
1751        Parms.append('cellA')
1752        Parms.append('cellB')
1753        Parms.append('cellC')
1754        if Layers['Laue'] != 'mmm':
1755            Parms.append('cellG')
1756#Transition parms
1757    Trans = Layers['Transitions']
1758    for iY in range(len(Layers['Layers'])):
1759        for iX in range(len(Layers['Layers'])):
1760            Parms.append('TransP;%d;%d'%(iY,iX))
1761            Parms.append('TransX;%d;%d'%(iY,iX))
1762            Parms.append('TransY;%d;%d'%(iY,iX))
1763            Parms.append('TransZ;%d;%d'%(iY,iX))
1764    return Parms
1765
1766def StackSim(Layers,ctrls,HistName='',scale=0.,background={},limits=[],inst={},profile=[]):
1767    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
1768   
1769    param: Layers dict: 'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
1770                        'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
1771                        'Layers':[],'Stacking':[],'Transitions':[]}
1772    param: HistName str: histogram name to simulate 'PWDR...'
1773    param: scale float: scale factor
1774    param: background dict: background parameters
1775    param: limits list: min/max 2-theta to be calculated
1776    param: inst dict: instrument parameters dictionary
1777    param: profile list: powder pattern data
1778   
1779    all updated in place   
1780    '''
1781    import atmdata
1782    path = sys.path
1783    for name in path:
1784        if 'bin' in name:
1785            DIFFaX = name+'/DIFFaX.exe'
1786            print ' Execute ',DIFFaX
1787            break
1788    # make form factor file that DIFFaX wants - atom types are GSASII style
1789    sf = open('data.sfc','w')
1790    sf.write('GSASII special form factor file for DIFFaX\n\n')
1791    atTypes = Layers['AtInfo'].keys()
1792    if 'H' not in atTypes:
1793        atTypes.insert(0,'H')
1794    for atType in atTypes:
1795        if atType == 'H': 
1796            blen = -.3741
1797        else:
1798            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
1799        Adat = atmdata.XrayFF[atType]
1800        text = '%4s'%(atType.ljust(4))
1801        for i in range(4):
1802            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
1803        text += '%11.6f%11.6f'%(Adat['fc'],blen)
1804        text += '%3d\n'%(Adat['Z'])
1805        sf.write(text)
1806    sf.close()
1807    #make DIFFaX control.dif file - future use GUI to set some of these flags
1808    cf = open('control.dif','w')
1809    if ctrls == '0\n0\n3\n': 
1810        x0 = profile[0]
1811        iBeg = np.searchsorted(x0,limits[0])
1812        iFin = np.searchsorted(x0,limits[1])
1813        if iFin-iBeg > 20000:
1814            iFin = iBeg+20000
1815        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
1816        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
1817        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
1818    else:
1819        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
1820        inst = {'Type':['XSC','XSC',]}
1821    cf.close()
1822    #make DIFFaX data file
1823    df = open('GSASII-DIFFaX.dat','w')
1824    df.write('INSTRUMENTAL\n')
1825    if 'X' in inst['Type'][0]:
1826        df.write('X-RAY\n')
1827    elif 'N' in inst['Type'][0]:
1828        df.write('NEUTRON\n')
1829    if ctrls == '0\n0\n3\n': 
1830        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
1831        U = forln2*inst['U'][1]/10000.
1832        V = forln2*inst['V'][1]/10000.
1833        W = forln2*inst['W'][1]/10000.
1834        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
1835        HW = np.mean(HWHM)
1836    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
1837    #    df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
1838        df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
1839    else:
1840        df.write('0.10\nNone\n')
1841    df.write('STRUCTURAL\n')
1842    a,b,c = Layers['Cell'][1:4]
1843    gam = Layers['Cell'][6]
1844    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
1845    laue = Layers['Laue']
1846    if laue == '2/m(ab)':
1847        laue = '2/m(1)'
1848    elif laue == '2/m(c)':
1849        laue = '2/m(2)'
1850    if 'unknown' in Layers['Laue']:
1851        df.write('%s %.3f\n'%(laue,Layers['Toler']))
1852    else:   
1853        df.write('%s\n'%(laue))
1854    df.write('%d\n'%(len(Layers['Layers'])))
1855    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
1856        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
1857    layerNames = []
1858    for layer in Layers['Layers']:
1859        layerNames.append(layer['Name'])
1860    for il,layer in enumerate(Layers['Layers']):
1861        if layer['SameAs']:
1862            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
1863            continue
1864        df.write('LAYER %d\n'%(il+1))
1865        if '-1' in layer['Symm']:
1866            df.write('CENTROSYMMETRIC\n')
1867        else:
1868            df.write('NONE\n')
1869        for ia,atom in enumerate(layer['Atoms']):
1870            [name,atype,x,y,z,frac,Uiso] = atom
1871            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
1872                frac /= 2.
1873            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
1874    df.write('STACKING\n')
1875    df.write('%s\n'%(Layers['Stacking'][0]))
1876    if 'recursive' in Layers['Stacking'][0]:
1877        df.write('%s\n'%Layers['Stacking'][1])
1878    else:
1879        if 'list' in Layers['Stacking'][1]:
1880            Slen = len(Layers['Stacking'][2])
1881            iB = 0
1882            iF = 0
1883            while True:
1884                iF += 68
1885                if iF >= Slen:
1886                    break
1887                iF = min(iF,Slen)
1888                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
1889                iB = iF
1890        else:
1891            df.write('%s\n'%Layers['Stacking'][1])   
1892    df.write('TRANSITIONS\n')
1893    for iY in range(len(Layers['Layers'])):
1894        sumPx = 0.
1895        for iX in range(len(Layers['Layers'])):
1896            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
1897            p = round(p,3)
1898            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
1899            sumPx += p
1900        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
1901            print 'ERROR - Layer probabilities sum to ',sumPx,' DIFFaX will insist it = 1.0'
1902            df.close()
1903            os.remove('data.sfc')
1904            os.remove('control.dif')
1905            os.remove('GSASII-DIFFaX.dat')
1906            return       
1907    df.close()   
1908    time0 = time.time()
1909    subp.call(DIFFaX)
1910    print 'DIFFaX time = %.2fs'%(time.time()-time0)
1911    if os.path.exists('GSASII-DIFFaX.spc'):
1912        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
1913        iFin = iBeg+Xpat.shape[1]
1914        bakType,backDict,backVary = SetBackgroundParms(background)
1915        backDict['Lam1'] = G2mth.getWave(inst)
1916        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
1917        profile[3][iBeg:iFin] = Xpat[2]*scale+profile[4][iBeg:iFin]
1918        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
1919            rv = st.poisson(profile[3][iBeg:iFin])
1920            profile[1][iBeg:iFin] = rv.rvs()
1921            Z = np.ones_like(profile[3][iBeg:iFin])
1922            Z[1::2] *= -1
1923            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
1924            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
1925        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
1926    #cleanup files..
1927        os.remove('GSASII-DIFFaX.spc')
1928    elif os.path.exists('GSASII-DIFFaX.sadp'):
1929        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
1930        Sadp = np.reshape(Sadp,(256,-1))
1931        Layers['Sadp']['Img'] = Sadp
1932        os.remove('GSASII-DIFFaX.sadp')
1933    os.remove('data.sfc')
1934    os.remove('control.dif')
1935    os.remove('GSASII-DIFFaX.dat')
1936   
1937def CalcStackingSADP(Layers):
1938   
1939    rand.seed()
1940    ranSeed = rand.randint(1,2**16-1)
1941# Scattering factors
1942    import atmdata
1943    atTypes = Layers['AtInfo'].keys()
1944    aTypes = []
1945    for atype in atTypes:
1946        aTypes.append('%4s'%(atype.ljust(4)))
1947    SFdat = []
1948    for atType in atTypes:
1949        if atType == 'H': 
1950            blen = -.3741
1951        else:
1952            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
1953        Adat = atmdata.XrayFF[atType]
1954        SF = np.zeros(9)
1955        SF[:8:2] = Adat['fa']
1956        SF[1:8:2] = Adat['fb']
1957        SF[8] = Adat['fc']
1958        SFdat.append(SF)
1959    SFdat = np.array(SFdat)
1960    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T)
1961# Controls
1962    try:   
1963        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
1964            '6/m','6/mmm'].index(Layers['Laue'])+1
1965    except ValueError:
1966        laueId = -1
1967    planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
1968    lmax = int(Layers['Sadp']['Lmax'])
1969    mult = int(Layers['Sadp']['Mult'])
1970# Sequences
1971    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
1972    try:
1973        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
1974    except ValueError:
1975        StkParm = -1
1976    if StkParm == 2:    #list
1977        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
1978        Nstk = len(StkSeq)
1979    else:
1980        Nstk = 1
1981        StkSeq = [0,]
1982    if StkParm == -1:
1983        StkParm = int(Layers['Stacking'][1])
1984    Wdth = Layers['Width'][0]
1985    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
1986    LaueSym = Layers['Laue'].ljust(12)
1987    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
1988   
1989    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
1990# atoms in layers
1991    AtomXOU = []
1992    AtomTp = []
1993    LayerSymm = []
1994    LayerNum = []
1995    layerNames = []
1996    Natm = 0
1997    Nuniq = 0
1998    for layer in Layers['Layers']:
1999        layerNames.append(layer['Name'])
2000    for il,layer in enumerate(Layers['Layers']):
2001        if layer['SameAs']:
2002            LayerNum.append(layerNames.index(layer['SameAs'])+1)
2003            continue
2004        else:
2005            LayerNum.append(il+1)
2006            Nuniq += 1
2007        if '-1' in layer['Symm']:
2008            LayerSymm.append(1)
2009        else:
2010            LayerSymm.append(0)
2011        for ia,atom in enumerate(layer['Atoms']):
2012            [name,atype,x,y,z,frac,Uiso] = atom
2013            Natm += 1
2014            AtomTp.append('%4s'%(atype.ljust(4)))
2015            Ta = atTypes.index(atype)+1
2016            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
2017    AtomXOU = np.array(AtomXOU)
2018    Nlayers = len(layerNames)
2019    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
2020# Transitions
2021    TransX = []
2022    TransP = []
2023    for Ytrans in Layers['Transitions']:
2024        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
2025        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
2026    TransP = np.array(TransP,dtype='float')
2027    TransX = np.array(TransX,dtype='float')
2028    pyx.pygettrans(Nlayers,TransP,TransX)
2029# result as Sadp
2030    mirror = laueId in [2,3,4,7,8,9,10]
2031    Nspec = 20001       
2032    spec = np.zeros(Nspec,dtype='double')   
2033    time0 = time.time()
2034    hkLim,Incr = pyx.pygetsadp(controls,Nspec,spec)
2035#    GSASIIpath.IPyBreak()
2036    Sapd = np.zeros((256,256))
2037    maxInt = np.max(spec[1:])
2038    Scale = mult*32767./maxInt
2039    iB = 0
2040    for i in range(hkLim):
2041        iF = iB+128
2042        p1 = 128+int(i*Incr)
2043        Sapd[:128,p1] = spec[iB:iF]
2044        Sapd[128:,p1] = spec[iF:iB:-1]
2045        if mirror:
2046            p2 = 128-int(i*Incr)
2047            Sapd[:128,p2] = spec[iB:iF]
2048            Sapd[128:,p2] = spec[iF:iB:-1]
2049        iB += 128
2050    Sapd *= Scale
2051    Sapd = np.where(Sapd<32767.,Sapd,32767.)
2052    Layers['Sadp']['Img'] = Sapd
2053    print 'GETSAD time = %.2fs'%(time.time()-time0)
2054   
2055#testing data
2056NeedTestData = True
2057def TestData():
2058    'needs a doc string'
2059#    global NeedTestData
2060    NeedTestData = False
2061    global bakType
2062    bakType = 'chebyschev'
2063    global xdata
2064    xdata = np.linspace(4.0,40.0,36000)
2065    global parmDict0
2066    parmDict0 = {
2067        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
2068        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
2069        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
2070        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
2071        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
2072        'Back0':5.384,'Back1':-0.015,'Back2':.004,
2073        }
2074    global parmDict1
2075    parmDict1 = {
2076        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
2077        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
2078        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
2079        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
2080        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
2081        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
2082        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
2083        'Back0':36.897,'Back1':-0.508,'Back2':.006,
2084        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2085        }
2086    global parmDict2
2087    parmDict2 = {
2088        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
2089        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
2090        'Back0':5.,'Back1':-0.02,'Back2':.004,
2091#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2092        }
2093    global varyList
2094    varyList = []
2095
2096def test0():
2097    if NeedTestData: TestData()
2098    msg = 'test '
2099    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
2100    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
2101    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
2102    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
2103    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
2104    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
2105   
2106def test1():
2107    if NeedTestData: TestData()
2108    time0 = time.time()
2109    for i in range(100):
2110        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
2111    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
2112   
2113def test2(name,delt):
2114    if NeedTestData: TestData()
2115    varyList = [name,]
2116    xdata = np.linspace(5.6,5.8,400)
2117    hplot = plotter.add('derivatives test for '+name).gca()
2118    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
2119    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2120    parmDict2[name] += delt
2121    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2122    hplot.plot(xdata,(y1-y0)/delt,'r+')
2123   
2124def test3(name,delt):
2125    if NeedTestData: TestData()
2126    names = ['pos','sig','gam','shl']
2127    idx = names.index(name)
2128    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
2129    xdata = np.linspace(5.6,5.8,800)
2130    dx = xdata[1]-xdata[0]
2131    hplot = plotter.add('derivatives test for '+name).gca()
2132    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
2133    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2134    myDict[name] += delt
2135    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2136    hplot.plot(xdata,(y1-y0)/delt,'r+')
2137
2138if __name__ == '__main__':
2139    import GSASIItestplot as plot
2140    global plotter
2141    plotter = plot.PlotNotebook()
2142#    test0()
2143#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
2144    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
2145        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
2146        test2(name,shft)
2147    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
2148        test3(name,shft)
2149    print "OK"
2150    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.