source: trunk/GSASIImpsubs.py @ 4919

Last change on this file since 4919 was 4919, checked in by vondreele, 2 years ago

change powder peak shape functions to include effect of varying step size.
Impacts scale factors & peak intensities.
New notification system invoked to let users know of change.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 8.6 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: 2021-06-03 15:12:41 +0000 (Thu, 03 Jun 2021) $
22# $Author: vondreele $
23# $Revision: 4919 $
24# $URL: trunk/GSASIImpsubs.py $
25# $Id: GSASIImpsubs.py 4919 2021-06-03 15:12:41Z vondreele $
26########### SVN repository information ###################
27from __future__ import division, print_function
28import multiprocessing as mp
29import numpy as np
30import numpy.ma as ma
31import GSASIIpath
32GSASIIpath.SetVersionNumber("$Revision: 4919 $")
33import GSASIIpwd as G2pwd
34import GSASIIfiles as G2fil
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 ResetMP():
46    '''Call after changing Config var 'Multiprocessing_cores' to force a resetting
47    of the useMP from the parameter.
48    '''
49    global ncores
50    ncores = None
51   
52def InitMP(allowMP=True):
53    '''Called to initialize use of Multiprocessing
54    '''
55    global useMP,ncores
56    if ncores is not None: return useMP,ncores
57    useMP = False
58    if not allowMP:
59        G2fil.G2Print('Multiprocessing disabled')
60        ncores = 0
61        return useMP,ncores
62    ncores = GSASIIpath.GetConfigValue('Multiprocessing_cores',0)
63    if ncores < 0: ncores = mp.cpu_count()//2
64    if ncores > 1:
65        useMP = True
66    if useMP:
67        G2fil.G2Print('Multiprocessing with {} cores enabled'.format(ncores))
68    return useMP,ncores
69
70################################################################################
71# Fobs Squared computation
72################################################################################       
73def InitFobsSqGlobals(x1,ratio1,shl1,xB1,xF1,im1,lamRatio1,kRatio1,xMask1,Ka21):
74    '''Initialize for the computation of Fobs Squared for powder histograms.
75    Puts lots of junk into the global namespace in this module.
76    '''
77    global x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2,cw
78    x = ma.getdata(x1)
79    cw = np.diff(x)
80    cw = np.append(cw,cw[-1])
81    ratio = ratio1
82    shl = shl1
83    xB = xB1
84    xF = xF1
85    im = im1
86    lamRatio = lamRatio1
87    kRatio = kRatio1
88    xMask = xMask1
89    Ka2 = Ka21
90
91def ComputeFobsSqCWbatch(profList):
92    sInt = 0
93    resList = []
94    for refl,iref in profList:
95        icod = ComputeFobsSqCW(refl,iref)
96        if type(icod) is tuple:
97            resList.append((icod[0],iref))
98            sInt += icod[1]
99        elif icod == -1:
100            resList.append((None,iref))
101        elif icod == -2:
102            break
103    return sInt,resList
104
105def ComputeFobsSqTOFbatch(profList):
106    sInt = 0
107    resList = []
108    for refl,iref in profList:
109        icod = ComputeFobsSqTOF(refl,iref)
110        if type(icod) is tuple:
111            resList.append((icod[0],iref))
112            sInt += icod[1]
113        elif icod == -1:
114            resList.append((None,iref))
115        elif icod == -2:
116            break
117    return sInt,resList
118       
119def ComputeFobsSqPinkbatch(profList):
120    sInt = 0
121    resList = []
122    for refl,iref in profList:
123        icod = ComputeFobsSqPink(refl,iref)
124        if type(icod) is tuple:
125            resList.append((icod[0],iref))
126            sInt += icod[1]
127        elif icod == -1:
128            resList.append((None,iref))
129        elif icod == -2:
130            break
131    return sInt,resList
132
133def ComputeFobsSqCW(refl,iref):
134    yp = np.zeros(len(x)) # not masked
135    sInt = 0
136    refl8im = 0
137    Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
138    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
139    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
140    iFin2 = iFin
141    if not iBeg+iFin:       #peak below low limit - skip peak
142        return 0
143    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
144        return -1
145    elif not iBeg-iFin:     #peak above high limit - done
146        return -2
147    elif iBeg < iFin:
148        fp,sumfp = G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,x[iBeg:iFin])
149        yp[iBeg:iFin] = 100.*refl[11+im]*refl[9+im]*fp*cw[iBeg:iFin]/sumfp
150        sInt = refl[11+im]*refl[9+im]
151        if Ka2:
152            pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
153            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
154            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
155            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
156            if iFin2 > iBeg2:
157                fp2,sumfp2 = G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,x[iBeg2:iFin2])
158                yp[iBeg2:iFin2] += 100.*refl[11+im]*refl[9+im]*kRatio*fp2*cw[iBeg2:iFin2]/sumfp2
159                sInt *= 1.+kRatio
160    refl8im = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
161    return refl8im,sInt
162
163def ComputeFobsSqTOF(refl,iref):
164    yp = np.zeros(len(x)) # not masked
165    refl8im = 0
166    Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
167    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
168    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
169    if not iBeg+iFin:       #peak below low limit - skip peak
170        return 0
171    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
172        return -1
173    elif not iBeg-iFin:     #peak above high limit - done
174        return -2
175    if iBeg < iFin:
176        fp,sumfp = G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])
177        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*fp*cw[iBeg:iFin]/sumfp
178    refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
179    return refl8im,refl[11+im]*refl[9+im]
180
181def ComputeFobsSqPink(refl,iref):
182    yp = np.zeros(len(x)) # not masked
183    refl8im = 0
184    Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im]/1.e4,refl[7+im]/100.)
185    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
186    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
187    if not iBeg+iFin:       #peak below low limit - skip peak
188        return 0
189    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
190        return -1
191    elif not iBeg-iFin:     #peak above high limit - done
192        return -2
193    if iBeg < iFin:
194        fp,sumfp = G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im]/1.e4,refl[7+im]/100.,x[iBeg:iFin])
195        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*fp*cw[iBeg:iFin]/sumfp
196    refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
197    return refl8im,refl[11+im]*refl[9+im]
198
199################################################################################
200# Powder Profile computation
201################################################################################       
202def InitPwdrProfGlobals(im1,shl1,x1):
203    '''Initialize for the computation of Fobs Squared for powder histograms.
204    Puts lots of junk into the global namespace in this module.
205    '''
206    global im,shl,x,cw,yc
207    im = im1
208    shl = shl1
209    x = ma.getdata(x1)
210    cw = np.diff(x)
211    cw = np.append(cw,cw[-1])
212    # create local copies of ycalc array
213    yc = np.zeros_like(x)
214
215def ComputePwdrProfCW(profList):
216    'Compute the peaks profile for a set of CW peaks and add into the yc array'
217    for pos,refl,iBeg,iFin,kRatio in profList:
218        fp = G2pwd.getFCJVoigt3(pos,refl[6+im],refl[7+im],shl,x[iBeg:iFin])[0]
219        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*fp/cw[iBeg:iFin]
220    return yc
221
222def ComputePwdrProfTOF(profList):
223    'Compute the peaks profile for a set of TOF peaks and add into the yc array'
224    for pos,refl,iBeg,iFin in profList:
225        fp = G2pwd.getEpsVoigt(pos,refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])[0]
226        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*fp*cw[iBeg:iFin]
227   
228def ComputePwdrProfPink(profList):
229    'Compute the peaks profile for a set of TOF peaks and add into the yc array'
230    for pos,refl,iBeg,iFin in profList:
231        fp = G2pwd.getEpsVoigt(pos,refl[12+im],refl[13+im],refl[6+im]/1.e4,refl[7+im]/100.,x[iBeg:iFin])[0]
232        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*fp
233    return yc
Note: See TracBrowser for help on using the repository browser.