source: branch/2frame/GSASIImpsubs.py @ 2944

Last change on this file since 2944 was 2944, checked in by toby, 4 years ago

Finally: multiprocessing for refinement
  • Property svn:eol-style set to native
File size: 21.3 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIImpsubs - routines used in multiprocessing*
4-------------------------------------------------
5
6The routines here are called either directly when GSAS-II is used without multiprocessing
7or in separate cores when multiprocessing is used.
8
9These routines are designed to be used in one of two ways:
10
11 * when multiprocessing is
12   enabled (see global variable useMP) the computational routines are called in
13   separate Python interpreter that is created and then deleted after use.
14
15 * when useMP is False, these routines are called directly from the main "thread".
16
17Note that :func:`GSASIImpsubs.InitMP` should be called before any of the other routines
18in this module are used.
19'''
20########### SVN repository information ###################
21# $Date: $
22# $Author: $
23# $Revision: $
24# $URL: $
25# $Id: $
26########### SVN repository information ###################
27import multiprocessing as mp
28import numpy as np
29import numpy.ma as ma
30import numpy.linalg as nl
31import GSASIIpath
32GSASIIpath.SetVersionNumber("$Revision: 2895 $")
33import GSASIIpwd as G2pwd
34import GSASIIstrMath as G2stMth
35
36sind = lambda x: np.sin(x*np.pi/180.)
37cosd = lambda x: np.cos(x*np.pi/180.)
38tand = lambda x: np.tan(x*np.pi/180.)
39#asind = lambda x: 180.*np.arcsin(x)/np.pi
40#acosd = lambda x: 180.*np.arccos(x)/np.pi
41#atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
42   
43ncores = None
44
45def InitMP(allowMP=True):
46    '''Called in routines to initialize use of Multiprocessing
47    '''
48    global useMP,ncores
49    if ncores is not None: return
50    useMP = False
51    if not allowMP:
52        print('Multiprocessing disabled')
53        ncores = 0
54        return
55    ncores = GSASIIpath.GetConfigValue('Multiprocessing_cores',-1)
56    if ncores < 0: ncores = mp.cpu_count()
57    if ncores > 1:
58        useMP = True
59    #if GSASIIpath.GetConfigValue('debug'):
60    if True:
61        print('Multiprocessing with {} cores enabled'.format(ncores))
62
63################################################################################
64# derivative computation
65################################################################################       
66def InitDerivGlobals(im1,calcControls1,SGMT1,hfx1,phfx1,pfx1,G1,GB1,g1,SGData1,
67                     parmDict1,wave1,shl1,x1,cw1,Ka21,A1,varylist1,dependentVars1,
68                     dFdvDict1,lamRatio1,kRatio1,doPawley1,pawleyLookup1):
69    '''Initialize for the computation of derivatives. Puts lots of junk into the global
70    namespace in this module, including the arrays for derivatives (when needed.)
71    '''
72    global im,calcControls,SGMT,hfx,phfx,pfx,G,GB,g,SGData,parmDict,wave,shl,x,cw,Ka2,A
73    global varylist,dependentVars,dFdvDict,lamRatio,kRatio,doPawley,pawleyLookup
74    im = im1
75    calcControls = calcControls1
76    SGMT = SGMT1
77    hfx = hfx1
78    phfx = phfx1
79    pfx = pfx1
80    G = G1
81    GB = GB1
82    g = g1
83    SGData = SGData1
84    parmDict = parmDict1
85    wave = wave1
86    shl = shl1
87    x = ma.getdata(x1)
88    cw = cw1
89    Ka2 = Ka21
90    A = A1
91    varylist = varylist1
92    dependentVars = dependentVars1
93    dFdvDict = dFdvDict1
94    lamRatio = lamRatio1
95    kRatio = kRatio1
96    doPawley = doPawley1
97    pawleyLookup = pawleyLookup1
98    # determine the parameters that will have derivatives computed only at end
99    global nonatomvarylist
100    nonatomvarylist = []
101    for name in varylist:
102        if '::RBV;' not in name:
103            try:
104                aname = name.split(pfx)[1][:2]
105                if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
106                    'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
107            except IndexError:
108                continue
109        nonatomvarylist.append(name)
110    global nonatomdependentVars
111    nonatomdependentVars = []
112    for name in dependentVars:
113        if '::RBV;' not in name:
114            try:
115                aname = name.split(pfx)[1][:2]
116                if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
117                    'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
118            except IndexError:
119                continue
120        nonatomdependentVars.append(name)
121    # create local copies of derivative arrays, if multiprocessing will be used
122    if useMP:
123        global dMdv
124        dMdv = np.zeros(shape=(len(varylist),len(x)))
125        global depDerivDict
126        depDerivDict = {j:np.zeros(shape=len(x)) for j in dependentVars}
127           
128
129def cellVaryDerv(pfx,SGData,dpdA): 
130    if SGData['SGLaue'] in ['-1',]:
131        return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
132            [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
133    elif SGData['SGLaue'] in ['2/m',]:
134        if SGData['SGUniq'] == 'a':
135            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
136        elif SGData['SGUniq'] == 'b':
137            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
138        else:
139            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
140    elif SGData['SGLaue'] in ['mmm',]:
141        return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
142    elif SGData['SGLaue'] in ['4/m','4/mmm']:
143        return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
144    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
145        return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
146    elif SGData['SGLaue'] in ['3R', '3mR']:
147        return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
148    elif SGData['SGLaue'] in ['m3m','m3']:
149        return [[pfx+'A0',dpdA[0]]]
150
151def ComputeDerivMPbatch(reflsList):
152    '''Computes a the derivatives for a batch of reflections and sums the them into
153    global arrays dMdv & depDerivDict. These arrays are returned once the computation
154    is completed.
155    '''
156    for refl,iref,fmin,fmax,iBeg,iFin in reflsList:
157        if ComputeDeriv(refl,iref,fmin,fmax,iBeg,iFin,dMdv,depDerivDict): break
158    return dMdv,depDerivDict
159
160def ComputeDeriv(refl,iref,fmin,fmax,iBeg,iFin,dMdv,depDerivDict):
161    '''Compute the parameter derivatives for a single reflection and add the results
162    into either array dMdv or depDerivDict
163    '''
164    global wave
165    if im:
166        h,k,l,m = refl[:4]
167    else:
168        h,k,l = refl[:3]
169    Uniq = np.inner(refl[:3],SGMT)
170    if 'T' in calcControls[hfx+'histType']:
171        wave = refl[14+im]
172    dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = G2stMth.GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
173    pos = refl[5+im]
174    calcKa2 = False
175    if 'C' in calcControls[hfx+'histType']:
176        tanth = tand(pos/2.0)
177        costh = cosd(pos/2.0)
178        if Ka2:
179            pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
180            iBeg2 = np.searchsorted(x,pos2-fmin)
181            iFin2 = np.searchsorted(x,pos2+fmax)
182            if iBeg2-iFin2:
183                calcKa2 = True
184                iFin = iFin2       
185        lenBF = iFin-iBeg
186        dMdpk = np.zeros(shape=(6,lenBF))
187        dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,x[iBeg:iFin])
188        for i in range(5):
189            dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
190        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
191        if calcKa2: 
192            dMdpk2 = np.zeros(shape=(6,lenBF))
193            dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,x[iBeg:iFin])
194            for i in range(5):
195                dMdpk2[i] = 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
196            dMdpk2[5] = 100.*cw[iBeg:iFin]*refl[11+im]*dMdipk2[0]
197            dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
198    else:   #'T'OF
199        lenBF = iFin-iBeg
200        if lenBF < 0: return True  #bad peak coeff
201        dMdpk = np.zeros(shape=(6,lenBF))
202        dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])
203        for i in range(6):
204            dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
205        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
206    if doPawley:
207        try:
208            if im:
209                pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
210            else:
211                pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
212            idx = varylist.index(pIdx)
213            dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9+im]
214            if Ka2: #not for TOF either
215                dMdv[idx][iBeg:iFin] += dervDict2['int']/refl[9+im]
216        except: # ValueError:
217            pass
218    if 'C' in calcControls[hfx+'histType']:
219        dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = G2stMth.GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
220        names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
221            hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
222            hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
223            hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
224            hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
225            hfx+'DisplaceY':[dpdY,'pos'],}
226        if 'Bragg' in calcControls[hfx+'instType']:
227            names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
228                hfx+'SurfRoughB':[dFdAb[1],'int'],})
229        else:
230            names.update({hfx+'Absorption':[dFdAb,'int'],})
231    else:   #'T'OF
232        dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = G2stMth.GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
233        names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
234            hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
235            hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
236            hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
237            hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
238            hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
239            hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
240    for name in names:
241        item = names[name]
242        if name in varylist:
243            dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
244            if calcKa2:
245                dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict2[item[1]]
246        elif name in dependentVars:
247            depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
248            if calcKa2:
249                depDerivDict[name][iBeg:iFin] += item[0]*dervDict2[item[1]]
250    for iPO in dIdPO:
251        if iPO in varylist:
252            dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
253            if calcKa2:
254                dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict2['int']
255        elif iPO in dependentVars:
256            depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
257            if calcKa2:
258                depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict2['int']
259    for i,name in enumerate(['omega','chi','phi']):
260        aname = pfx+'SH '+name
261        if aname in varylist:
262            dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
263            if calcKa2:
264                dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict2['int']
265        elif aname in dependentVars:
266            depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
267            if calcKa2:
268                depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict2['int']
269    for iSH in dFdODF:
270        if iSH in varylist:
271            dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
272            if calcKa2:
273                dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict2['int']
274        elif iSH in dependentVars:
275            depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
276            if calcKa2:
277                depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict2['int']
278    cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
279    for name,dpdA in cellDervNames:
280        if name in varylist:
281            dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
282            if calcKa2:
283                dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict2['pos']
284        elif name in dependentVars: #need to scale for mixed phase constraints?
285            depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
286            if calcKa2:
287                depDerivDict[name][iBeg:iFin] += dpdA*dervDict2['pos']
288    dDijDict = G2stMth.GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
289    for name in dDijDict:
290        if name in varylist:
291            dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
292            if calcKa2:
293                dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict2['pos']
294        elif name in dependentVars:
295            depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
296            if calcKa2:
297                depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict2['pos']
298    for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
299        if name in varylist:
300            dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
301            if calcKa2:
302                dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict2['pos']
303        elif name in dependentVars:
304            depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
305            if calcKa2:
306                depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict2['pos']
307    if 'C' in calcControls[hfx+'histType']:
308        sigDict,gamDict = G2stMth.GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
309    else:   #'T'OF
310        sigDict,gamDict = G2stMth.GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
311    for name in gamDict:
312        if name in varylist:
313            dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
314            if calcKa2:
315                dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict2['gam']
316        elif name in dependentVars:
317            depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
318            if calcKa2:
319                depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict2['gam']
320    for name in sigDict:
321        if name in varylist:
322            dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
323            if calcKa2:
324                dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict2['sig']
325        elif name in dependentVars:
326            depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
327            if calcKa2:
328                depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict2['sig']
329    for name in ['BabA','BabU']:
330        if refl[9+im]:
331            if phfx+name in varylist:
332                dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
333                if calcKa2:
334                    dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]
335            elif phfx+name in dependentVars:                   
336                depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
337                if calcKa2:
338                    depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]                 
339    if not doPawley and not parmDict[phfx+'LeBail']:
340        #do atom derivatives -  for RB,F,X & U so far - how do I scale mixed phase constraints?
341        corr = 0.
342        #corr2 = 0.
343        if refl[9+im]:             
344            corr = dervDict['int']/refl[9+im]
345            #if calcKa2:  # commented out in Bob's code. Why?
346            #    corr2 = dervDict2['int']/refl[9+im]
347        for name in nonatomvarylist:
348            dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
349            #if calcKa2:
350            #   dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr2 # unneeded w/o above
351        for name in nonatomdependentVars:
352           depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
353           #if calcKa2:
354           #    depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr2
355
356################################################################################
357# Fobs Squared computation
358################################################################################       
359#x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2
360def InitFobsSqGlobals(x1,ratio1,shl1,xB1,xF1,im1,lamRatio1,kRatio1,xMask1,Ka21):
361    '''Initialize for the computation of Fobs Squared for powder histograms.
362    Puts lots of junk into the global namespace in this module.
363    '''
364    global x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2
365    x = ma.getdata(x1)
366    ratio = ratio1
367    shl = shl1
368    xB = xB1
369    xF = xF1
370    im = im1
371    lamRatio = lamRatio1
372    kRatio = kRatio1
373    xMask = xMask1
374    Ka2 = Ka21
375
376def ComputeFobsSqCWbatch(profList):
377    sInt = 0
378    resList = []
379    for refl,iref in profList:
380        icod = ComputeFobsSqCW(refl,iref)
381        if type(icod) is tuple:
382            resList.append((icod[0],iref))
383            sInt += icod[1]
384        elif icod == -1:
385            res.append((None,iref))
386        elif icod == -2:
387            break
388    return sInt,resList
389
390def ComputeFobsSqTOFbatch(profList):
391    sInt = 0
392    resList = []
393    for refl,iref in profList:
394        icod = ComputeFobsSqTOF(refl,iref)
395        if type(icod) is tuple:
396            resList.append((icod[0],iref))
397            sInt += icod[1]
398        elif icod == -1:
399            res.append((None,iref))
400        elif icod == -2:
401            break
402    return sInt,resList
403       
404def ComputeFobsSqCW(refl,iref):
405    yp = np.zeros(len(x)) # not masked
406    sInt = 0
407    refl8im = 0
408    Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
409    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
410    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
411    iFin2 = iFin
412    if not iBeg+iFin:       #peak below low limit - skip peak
413        return 0
414    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
415        return -1
416    elif not iBeg-iFin:     #peak above high limit - done
417        return -2
418    elif iBeg < iFin:
419        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(
420            refl[5+im],refl[6+im],refl[7+im],shl,x[iBeg:iFin])
421        sInt = refl[11+im]*refl[9+im]
422        if Ka2:
423            pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
424            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
425            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
426            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
427            if iFin2 > iBeg2: 
428                yp[iBeg2:iFin2] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(
429                    pos2,refl[6+im],refl[7+im],shl,x[iBeg2:iFin2])
430                sInt *= 1.+kRatio
431    refl8im = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
432    return refl8im,sInt
433
434def ComputeFobsSqTOF(refl,iref):
435    yp = np.zeros(len(x)) # not masked
436    refl8im = 0
437    Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
438    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
439    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
440    if not iBeg+iFin:       #peak below low limit - skip peak
441        return 0
442    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
443        return -1
444    elif not iBeg-iFin:     #peak above high limit - done
445        return -2
446    if iBeg < iFin:
447        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(
448            refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])
449    refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
450    return refl8im,refl[11+im]*refl[9+im]
451################################################################################
452# Powder Profile computation
453################################################################################       
454def InitPwdrProfGlobals(im1,shl1,x1):
455    '''Initialize for the computation of Fobs Squared for powder histograms.
456    Puts lots of junk into the global namespace in this module.
457    '''
458    global im,shl,x
459    im = im1
460    shl = shl1
461    x = ma.getdata(x1)
462    global cw
463    cw = np.diff(x)
464    cw = np.append(cw,cw[-1])
465    # create local copies of ycalc array
466    if useMP:
467        global yc
468        yc = np.zeros_like(x1)
469
470
471def ComputePwdrProfCW(profList):
472    'Compute the peaks profile for a set of CW peaks and add into the yc array'
473    for pos,refl,iBeg,iFin,kRatio in profList:
474        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(
475            pos,refl[6+im],refl[7+im],shl,x[iBeg:iFin])
476    return yc
477
478def ComputePwdrProfTOF(profList):
479    'Compute the peaks profile for a set of TOF peaks and add into the yc array'
480    for pos,refl,iBeg,iFin in profList:
481        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(
482            pos,refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])/cw[iBeg:iFin]
483    return yc
484   
Note: See TracBrowser for help on using the repository browser.