source: trunk/NIST_profile/profile_functions_class.py

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

Add FPA peak fitting

  • Property svn:eol-style set to native
File size: 66.7 KB
Line 
1## @file
2#  @brief The file containing the core definitions for the XRD Fundamental Parameneters Model (FPA) computation in python
3#  @mainpage
4#  @author Marcus H. Mendenhall (marcus.mendenhall@nist.gov)
5#  @date March, 2015, updated August 2018
6#  @copyright
7#  The "Fundamental Parameters Python Code" ("software") is provided by the National Institute of Standards and Technology (NIST),
8#  an agency of the United States Department of Commerce, as a public service.
9#  This software is to be used for non-commercial research purposes only and is expressly provided "AS IS."
10#  Use of this software is subject to your acceptance of these terms and conditions.
11#  \n \n
12#  NIST MAKES NO WARRANTY OF ANY KIND, EXPRESS, IMPLIED, IN FACT OR ARISING BY OPERATION OF LAW,
13#  INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE,
14#  NON-INFRINGEMENT AND DATA ACCURACY.
15#  NIST NEITHER REPRESENTS NOR WARRANTS THAT THE OPERATION OF THE SOFTWARE WILL BE UNINTERRUPTED
16#  OR ERROR-FREE, OR THAT ANY DEFECTS WILL BE CORRECTED.
17#  NIST DOES NOT WARRANT OR MAKE ANY REPRESENTATIONS REGARDING THE USE OF THE SOFTWARE OR THE RESULTS THEREOF,
18#  INCLUDING BUT NOT LIMITED TO THE CORRECTNESS, ACCURACY, RELIABILITY, OR USEFULNESS OF THE SOFTWARE.
19#  \n \n
20#  NIST SHALL NOT BE LIABLE AND YOU HEREBY RELEASE NIST FROM LIABILITY FOR ANY INDIRECT, CONSEQUENTIAL,
21#  SPECIAL, OR INCIDENTAL DAMAGES (INCLUDING DAMAGES FOR LOSS OF BUSINESS PROFITS, BUSINESS INTERRUPTION,
22#  LOSS OF BUSINESS INFORMATION, AND THE LIKE), WHETHER ARISING IN TORT, CONTRACT, OR OTHERWISE,
23#  ARISING FROM OR RELATING TO THE SOFTWARE (OR THE USE OF OR INABILITY TO USE THIS SOFTWARE),
24#  EVEN IF NIST HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
25#  \n \n
26#  Software developed by NIST employees is not subject to copyright protection within the United States.
27#  By using this software, creating derivative works or by incorporating this software into another product,
28#  you agree that you will use the software only for non-commercial research purposes
29#  and will indemnify and hold harmless the United States Government
30#  for any and all damages or liabilities that arise out of any use by you.
31#
32#
33#  @file
34#  This file contains the main computational class FP_profile, which stores cached information to allow it
35#  to efficiently recompute profiles when parameters have been modified.  It also includes some
36#  helper classes, and a function @ref fourier_line_profile(), which has too many parameters,
37#  which computes a profile in a single call.
38#
39#  @section intro_sec Introduction
40#  Compute diffractometer line profile functions by methods
41#  from Cheary & Coelho 1998 and Mullen & Cline paper and 'R' package.
42#  Accumulate all convolutions in Fourier space, for efficiency,
43#  except for axial divergence, which needs to be weighted in real space for I3 integral.
44#
45# changelog:
46#
47# 15 August, 2018 -- MHM
48# fixed apparent error in flip of eps0 across twotheta=90 around line 549
49#
50# @file
51
52from __future__ import print_function
53
54import os, sys, math
55import numpy
56## @brief figure out which FFT package we have, and import it
57try:
58    from pyfftw.interfaces import numpy_fft, cache
59    ## recorded variant of real fft that we will use
60    best_rfft=numpy_fft.rfft
61    ## recorded variant of inverse real fft that we will use
62    best_irfft=numpy_fft.irfft
63    ## @brief create a table of nice factorizations for the FFT package
64    #
65    #  this is built once, and shared by all instances
66    #  fftw can handle a variety of transform factorizations
67    ft_factors=[
68        2*2**i*3**j*5**k for i in xrange(20) for j in range(10) for k in range(8)
69        if 2*2**i*3**j*5**k <= 1000000
70    ]
71    cache.enable()
72    cache.set_keepalive_time(1.0)
73except ImportError:
74    best_rfft=numpy.fft.rfft
75    best_irfft=numpy.fft.irfft
76    ## @brief create a table of nice factorizations for the FFT package
77    #
78    # this is built once, and shared by all instances
79    # numpy fft is not too efficient for things other than a power of two,
80    # although my own measurements says it really does fine.  For now, leave all factors available
81    ft_factors=[
82        2*2**i*3**j*5**k for i in range(20) for j in range(10) for k in range(8)
83        if 2*2**i*3**j*5**k <= 1000000
84    ]
85
86ft_factors.sort()
87ft_factors=numpy.array(ft_factors, numpy.int)
88
89## @brief used for debugging moments from FP_profile.axial_helper().
90moment_list=[]
91## @brief if this is \a True, compute and save moment errors
92collect_moment_errors=False
93
94## @brief a skeleton class which makes a combined dict and namespace interface for easy pickling and data passing
95class profile_data(object):
96    ## @brief initialize the class
97    #  @param self self
98    #  @param kwargs keyword=value list to pre-populate the class
99    def __init__(self, **kwargs):
100        mydict={}
101        mydict.update(kwargs)
102        for k,v in mydict.items():
103            setattr(self,k,v)
104        ## a dictionary which shadows our attributes.
105        self.dictionary=mydict
106
107    ## @brief add new symbols to both the attributes and dictionary for the class
108    #  @param self self
109    #  @param kwargs keyword=value pairs
110    def add_symbol(self, **kwargs):
111        self.dictionary.update(kwargs)
112        for k,v in kwargs.items():
113            setattr(self,k,v)
114
115## @brief the main class, which handles a single reflection
116#
117#  This class is designed to be highly extensible by inheriting new convolvers.  When it is initialized,
118#  it scans its namespace for specially formatted names, which can come from mixin classes.
119#  If it finds a function name of the form conv_xxx, it will call this funtion to create a convolver.
120#  If it finds a name of the form info_xxx it will associate the dictionary with that convolver, which
121#  can be used in UI generation, for example.  The class, as it stands, does nothing significant with it.
122#  If it finds str_xxx, it will use that function to format a printout of the current state of the
123#  convolver conv_xxx, to allow improved report generation for convolvers.
124#
125#  When it is asked to generate a profile, it calls all known convolvers.  Each convolver returns the
126#  Fourier transform of its convolvution.  The transforms are multiplied together, inverse transformed,
127#  and after fixing the periodicity issue, subsampled, smoothed and returned.
128#
129#  If a convolver returns \a None, it is not multipled into the product.
130class FP_profile:
131
132    ## @brief the number of histories to cache; can be overridden if memory is an issue.
133    max_history_length=5 #max number of histories for each convolver
134   
135    ## @brief length_scale_m sets scaling for nice printing of parameters.
136    #   if the units are in mm everywhere, set it to 0.001, e.g.
137    #   convolvers which implement their own str_xxx method may use this to format their results,
138    #   especially if 'natural' units are not meters.  Typical is wavelengths and lattices in nm or angstroms,
139    #   for example.
140    length_scale_m = 1.0
141   
142    ## @brief initialize the instance
143    #
144    #  @param anglemode = 'd' if setup will be in terms of a d-spacing,
145    #  otherwise 'twotheta' if setup will be at a fixed 2theta value.
146    #  @param output_gaussian_smoother_bins_sigma the number of bins for post-smoothing of data. 1.0 is good.
147    #  @param oversampling the number of bins internally which will get computed for each bin the the final result.
148    def __init__(self, anglemode,
149            output_gaussian_smoother_bins_sigma=1.0,
150            oversampling=10
151        ):
152        if anglemode not in ("d","twotheta"):
153            raise Exception("invalid angle mode %s, must be 'd' or 'twotheta'"%anglemode)
154        ## set to either 'd' for d-spacing based position, or 'twotheta' for angle-based position
155        self.anglemode=anglemode
156        ## sigma, in units of bins, for the final smoother. \a None means no final smoothing step.
157        self.output_gaussian_smoother_bins_sigma=output_gaussian_smoother_bins_sigma
158        ## the number of internal bins computed for each bin in the final output. 5-10 is usually plenty.
159        self.oversampling=oversampling
160        ## List of our convolvers, found by introspection of names beginning with 'conv_'
161        self.convolvers=convolvers=[x for x in dir(self) if x.startswith("conv_")]
162        ## A dictionary which will store all the parameters local to each convolution
163        self.param_dicts=dict([(c,{}) for c in convolvers])
164        #add global parameters, associated with no specific convolver
165        ## A dictionary of bound functions to call to compute convolutions
166        self.convolver_funcs=dict([(x, getattr(self,x)) for x in convolvers])
167        ## If \a True, print cache hit information
168        self.debug_cache=False
169
170    ## @brief return the name of the function that called this.
171    #  @param self self
172    #  @return name of calling function
173    #
174    #  Useful for convolvers to identify themselves
175    def get_function_name(self):
176        return sys._getframe(1).f_code.co_name
177
178    ## @brief add a numpy array to the list of objects that can be thrown away on pickling.
179    #  @param self self
180    #  @param b the buffer to add to the list
181    #  @return the same buffer, to make nesting easy.
182    def add_buffer(self, b):
183        self._clean_on_pickle.add(id(b))
184        return b
185   
186    ## @brief move the computation window to a new position, and clear history.
187    #  @param self self
188    #  @param twotheta_window_center_deg the center position of the middle bin of the window, in degrees
189    #  @param twotheta_window_fullwidth_deg the full width of the window, in degrees
190    #  @param twotheta_output_points the number of bins in the final output
191    #
192    #  sets up many arrays
193    def set_window(self, twotheta_window_center_deg, twotheta_window_fullwidth_deg,
194            twotheta_output_points
195        ):
196        """move the compute window to a new location and compute grids,
197            without resetting all parameters.  Clears convolution history.
198        """
199       
200        ## the saved width of the window, in degrees
201        self.twotheta_window_fullwidth_deg=twotheta_window_fullwidth_deg
202        ## the saved center of the window, in degrees
203        self.twotheta_window_center_deg=twotheta_window_center_deg
204        ## the width of the window, in radians
205        self.window_fullwidth=window_fullwidth=twotheta_window_fullwidth_deg*math.pi/180
206        ## the center of the window, in radians
207        self.twotheta_window_center=twotheta=twotheta_window_center_deg*math.pi/180.0
208        ## the number of points to compute in the final results
209        self.twotheta_output_points=twotheta_output_points
210        ## the number of points in Fourier space to compute
211        self.n_omega_points=nn=self.oversampling*twotheta_output_points//2+1
212        #build all the arrays
213        x=self._clean_on_pickle=set() #keep a record of things we don't keep when pickled
214        b=self.add_buffer #shortcut
215       
216        ## a real-format scratch buffer
217        self._rb1=b(numpy.zeros(nn,numpy.float))
218        ## a real-format scratch buffer
219        self._rb2=b(numpy.zeros(nn,numpy.float))
220        ## a real-format scratch buffer
221        self._rb3=b(numpy.zeros(nn,numpy.float))
222        ## a complex-format scratch buffer
223        self._cb1=b(numpy.zeros(nn,numpy.complex))
224        ## a scratch buffer used by the axial helper
225        self._f0buf=b(numpy.zeros(self.oversampling*twotheta_output_points, numpy.float))
226        ## a scratch buffer used for axial divergence
227        self._epsb2=b(numpy.zeros(self.oversampling*twotheta_output_points, numpy.float))
228        ## the I2+ buffer
229        self._I2p=b(numpy.zeros(self.oversampling*twotheta_output_points, numpy.float))
230        ## the I2- buffer
231        self._I2m=b(numpy.zeros(self.oversampling*twotheta_output_points, numpy.float))
232        ## another buffer used for axial divergence
233        self._axial=b(numpy.zeros(self.oversampling*twotheta_output_points, numpy.float))
234        ## the largest frequency in Fourier space
235        omega_max=self.n_omega_points*2*math.pi/window_fullwidth
236        #build the x grid and the complex array that is the convolver
237        #omega is in inverse radians in twotheta space (!)
238        #i.e. if a transform is written
239        #I(ds) = integral(A(L) exp(2 pi i L ds) dL where L is a real-space length,
240        #and s=2 sin(twotheta/2)/lambda
241        #then ds=2*pi*omega*cos(twotheta/2)/lambda (double check this!)
242        ## The grid in Fourier space, in inverse radians
243        self.omega_vals=b(numpy.linspace(0, omega_max, self.n_omega_points, endpoint=True))
244        ## The grid in Fourier space, in inverse degrees
245        self.omega_inv_deg=b(self.omega_vals*(math.pi/180))
246        ## The grid in real space, in radians, with full oversampling
247        self.twothetasamples=b(numpy.linspace(
248            twotheta-window_fullwidth/2.0,twotheta+window_fullwidth/2.0,
249            self.twotheta_output_points*self.oversampling, endpoint=False))
250        ## The grid in real space, in degrees, with full oversampling
251        self.twothetasamples_deg=b(numpy.linspace(
252            twotheta_window_center_deg-twotheta_window_fullwidth_deg/2.0,
253            twotheta_window_center_deg+twotheta_window_fullwidth_deg/2.0,
254            self.twotheta_output_points*self.oversampling, endpoint=False))
255        ## Offsets around the center of the window, in radians
256        self.epsilon=b(self.twothetasamples-twotheta)
257       
258        ## A dictionary in which we collect recent state for each convolution.
259        ## whenever the window gets reset, all of these get cleared
260        self.convolution_history=dict([(x, []) for x in self.convolvers])
261
262        ## A dictionary of Lorentz widths, used for de-periodizing the final result.
263        self.lor_widths={}
264
265    ## @brief find a bin count close to what we need, which works well for Fourier transforms.
266    #  @param self self
267    #  @param count a number of bins.
268    #  @return a bin count somewhat larger than \a count which is efficient for FFT
269    def get_good_bin_count(self, count):
270        return ft_factors[ft_factors.searchsorted(count)]
271
272    ## @brief given a window center, and approximate width, and an exact bin spacing,
273    #  set up with a really good, slightly wider window for efficient FFT.
274    #  @param self self
275    #  @param twotheta_window_center_deg exact position of center bin, in degrees
276    #  @param twotheta_approx_window_fullwidth_deg approximate desired width
277    #  @param twotheta_exact_bin_spacing_deg the exact bin spacing to use
278    def set_optimized_window(self, twotheta_window_center_deg,
279        twotheta_approx_window_fullwidth_deg,
280        twotheta_exact_bin_spacing_deg):
281        """pick a bin count which factors cleanly for FFT, and adjust the window width
282        to preserve the exact center and bin spacing"""
283        bins=self.get_good_bin_count(int(1+
284            twotheta_approx_window_fullwidth_deg/twotheta_exact_bin_spacing_deg))
285        window_actwidth=twotheta_exact_bin_spacing_deg*bins
286        self.set_window(twotheta_window_center_deg=twotheta_window_center_deg,
287            twotheta_window_fullwidth_deg=window_actwidth,
288            twotheta_output_points=bins)
289   
290    ## @brief update the dictionary of parameters associated with the given convolver
291    #  @param self self
292    #  @param convolver the name of the convolver.  name 'global', e.g., attaches to function 'conv_global'
293    #  @param kwargs keyword-value pairs to update the convolvers dictionary.
294    def set_parameters(self, convolver="global", **kwargs):
295        """update the arguments for a specific convolver, by name.  no convolver -> global parameters"""
296        self.param_dicts["conv_"+convolver].update(kwargs)
297   
298    ## @brief get a cached, pre-computed convolver associated with the given parameters,
299    #  or a newly zeroed convolver if the cache doesn't contain it. Recycles old cache entries.
300    #  @param name the name of the convolver to seek
301    #  @param key any hashable object which identifies the parameters for the computation
302    #  @param format the type of the array to create, if one is not found.
303    #  @return \a flag, which is \a True if valid data were found, or \a False if the returned array is zero,
304    #  and \a array, which must be computed by the convolver if \a flag was \a False.
305    #
306    #  This takes advantage of the mutability of arrays.
307    #  When the contents of the array are changed by the convolver, the
308    #  cached copy is implicitly updated, so thsat the next time this is called with the same
309    #  parameters, it will return the previous array.
310    def get_conv(self, name, key, format=numpy.float):
311        #either get an old buffer associated with key, with valid convolution kernel,
312        #or create a new, empty buffer. If we already have too many buffers, swap one out.
313        history=self.convolution_history[name] #previous computed values as a list
314        for idx, (k,b) in enumerate(history):
315            if k==key:
316                history.insert(0,history.pop(idx)) #move to front to mark recently used
317                if self.debug_cache: print >> sys.stderr, name, True
318                return True, b #True says we got a buffer with valid data
319        if len(history)==self.max_history_length:
320            buf=history.pop(-1)[1] #re-use oldest buffer
321            buf[:]=0
322        else:
323            buf=numpy.zeros(self.n_omega_points, format)
324        history.insert(0,(key,buf))
325        if self.debug_cache: print >> sys.stderr, name, False
326        return False, buf #False says buffer is empty, need to recompute
327
328    ## @brief scan for functions named conv_xxx, and associated info_xxx entries.
329    #  @param self self
330    #  @return list of (convolver_xxx, info_xxx) pairs
331    def get_convolver_information(self):
332        """return a list of convolvers, and what we know about them"""
333        info_list=[]
334        for k, f in self.convolver_funcs.items():
335            info=getattr(self, "info_"+k[5:], {})
336            info["docstring"]=f.__doc__
337            info_list.append((k,info))
338
339        return info_list
340
341    ## A dictionary of default parameters for the global namespace,
342    #  used to seed a GUI which can harvest this for names, descriptions, and initial values
343    info_global=dict(
344        group_name="Global parameters",
345        help="this should be help information",
346        param_info=dict(
347            twotheta0_deg=("Bragg center of peak (degrees)", 30.0),
348            d=("d spacing (m)", 4.00e-10),
349            dominant_wavelength=("wavelength of most intense line (m)", 1.5e-10)
350        )
351    )
352   
353    ## @brief return a nicely formatted report descibing the current state of this class
354    #  @brief self self
355    #  @return string of formatted information
356    #
357    #  this looks for an str_xxx function associated with each conv_xxx name.  If it is found, that
358    #  function if called to get the state of conv_xxx.  Otherwise,
359    #  this simply formats the dictionary of parameters for the convolver, and uses that.
360    def __str__(self):
361        #scan for all convolvers, and find if they have str_xxx methods matching
362        #if they do, call it for custom output formatting.  Otherwise,
363        #just print the dict for each convolver.
364        keys=list(self.convolver_funcs.keys())
365        keys.sort() #always return info in the same order... maybe later some priority list
366        keys.insert(0, keys.pop(keys.index('conv_global'))) #global is always first, anyways!
367        strings=["", "***convolver id 0x%08x:"%id(self)]
368        for k in keys:
369            strfn="str_"+k[5:]
370            if hasattr(self, strfn):
371                strings.append(getattr(self, strfn)())
372            else:
373                dd=self.param_dicts["conv_"+k[5:]]
374                if dd:
375                    strings.append(k[5:]+": "+str(dd))
376        return '\n'.join(strings)
377       
378    ## @brief returns a string representation for the global context.
379    #  @param self self
380    #  @return report on global parameters.
381    def str_global(self):
382        self.param_dicts["conv_global"].setdefault("d",0) #in case it's not initialized
383        return "global: peak center=%(twotheta0_deg).4f, d=%(d).8g, eq. div=%(equatorial_divergence_deg).3f" % self.param_dicts["conv_global"]
384   
385    ## @brief the global context isn't really a convolver, returning \a None means ignore result
386    #  @param self self
387    #  @return \a None, always
388    def conv_global(self):
389        """a dummy convolver to hold global variables and information"""
390        return None
391   
392    ## @brief the function F0 from the paper.
393    #
394    #  compute k/sqrt(peakpos-x)+y0 nonzero between outer & inner (inner is closer to peak)
395    #  or k/sqrt(x-peakpos)+y0 if reversed (i.e. if outer > peak)
396    #  fully evaluated on a specified eps grid, and stuff into destination
397    #  @param self self
398    #  @param outerbound the edge of the function farthest from the singularity, referenced to epsvals
399    #  @param innerbound the edge closest to the singularity, referenced to epsvals
400    #  @param epsvals the array of two-theta values or offsets
401    #  @param[in,out] destination an array into which final results are summed.  modified in place!
402    #  @param peakpos the position of the singularity, referenced to epsvals.
403    #  @param y0 the constant offset
404    #  @param k the scale factor
405    #  @return (\a lower_index, \a upper_index ) python style bounds
406    #    for region of \a destination which has been modified.
407    def axial_helper(self, outerbound, innerbound, epsvals, destination, peakpos=0, y0=0, k=0):
408        if k==0:
409            return len(epsvals)//2,len(epsvals)//2+1 #nothing to do, point at the middle
410       
411        dx=epsvals[1]-epsvals[0] #bin width for normalizer so sum(result*dx)=exact integral
412        flip=outerbound > peakpos #flag for whether tail is to the left or right.
413
414        delta1=abs(innerbound-peakpos)
415        delta2=abs(outerbound-peakpos)
416        #this is the analytic area the function must have, integral(1/sqrt(eps0-eps)) from lower to upper
417        exactintegral=2*k*(math.sqrt(delta2)-math.sqrt(delta1))
418        exactintegral += y0*(delta2-delta1)
419        #exactintegral=max(0,exactintegral) #can never be < 0, beta out of range.
420        exactintegral *= 1/dx #normalize so sum is right
421
422        #compute the exact centroid we need for this
423        if abs(delta2-delta1) < 1e-12:
424            exact_moment1=0
425        else:
426            exact_moment1=( #simplified from Mathematica FortranForm
427                (4*k*(delta2**1.5-delta1**1.5) + 3*y0*(delta2**2-delta1**2))/
428                (6.*(2*k*(math.sqrt(delta2)-math.sqrt(delta1)) + y0*(delta2-delta1) ))
429            )
430            if not flip: exact_moment1=-exact_moment1
431        exact_moment1+=peakpos
432
433        # note: because of the way the search is done, this produces a result
434        # with a bias of 1/2 channel to the left of where it should be.
435        # this is fixed by shifting all the parameters up 1/2 channel
436        outerbound+=dx/2; innerbound+=dx/2; peakpos+=dx/2 #fix 1/2 channel average bias from search
437        #note: searchsorted(side="left") always returns the position of the bin to the right of the match, or exact bin
438        idx0, idx1=epsvals.searchsorted((outerbound, innerbound), side='left')
439       
440        if abs(outerbound-innerbound)< (2*dx) or abs(idx1-idx0) < 2: #peak has been squeezed out, nothing to do
441            #preserve the exact centroid: requires summing into two channels
442            #for a peak this narrow, no attempt to preserve the width.
443            #note that x1 (1-f1) + (x1+dx) f1 = mu has solution (mu - x1) / dx  = f1
444            #thus, we want to sum into a channel that has x1<mu by less than dx, and the one to its right
445            idx0 = min(idx0,idx1)-1 #pick left edge and make sure we are past it
446            while exact_moment1-epsvals[idx0] > dx:
447                #normally only one step max, but make it a loop in case of corner case
448                idx0 += 1
449            f1=(exact_moment1-epsvals[idx0])/dx
450            res=(exactintegral*(1-f1),exactintegral*f1)
451            destination[idx0:idx0+2]+=res
452            if collect_moment_errors:
453                centroid2=(res*epsvals[idx0:idx0+2]).sum()/sum(res)
454                moment_list.append((centroid2-exact_moment1)/dx)
455            return [idx0, idx0+2] #return collapsed bounds
456
457        if not flip:
458            if epsvals[idx0] != outerbound: idx0=max(idx0-1,0)
459            idx1=min(idx1+1, len(epsvals))
460            sign=1
461            deps=self._f0buf[idx0:idx1]
462            deps[:]=peakpos
463            deps-=epsvals[idx0:idx1]
464            deps[-1]=peakpos-min(innerbound, peakpos)
465            deps[0]=peakpos-outerbound
466        else:
467            idx0, idx1 = idx1, idx0
468            if epsvals[idx0] != innerbound: idx0=max(idx0-1,0)
469            idx1=min(idx1+1, len(epsvals))
470            sign=-1
471            deps=self._f0buf[idx0:idx1]
472            deps[:]=epsvals[idx0:idx1]
473            deps-=peakpos
474            deps[0]=max(innerbound, peakpos)-peakpos
475            deps[-1]=outerbound-peakpos
476
477        dx0=abs(deps[1]-deps[0])
478        dx1=abs(deps[-1]-deps[-2])
479       
480        #make the numerics accurate: compute average on each bin, which is
481        #integral of 1/sqrt = 2*sqrt, then difference integral
482        intg=numpy.sqrt(deps,deps) #do it in place, return value is actually deps too
483        intg *= 2*k*sign
484
485        intg[:-1]-=intg[1:] #do difference in place, running forward to avoid self-trampling
486        intg[1:-2] += y0*dx #add constant
487        #handle narrowed bins on ends carefully
488        intg[0] += y0*dx0
489        intg[-2] += y0*dx1
490       
491        if min(intg[:-1]) < -1e-10*max(intg[:-1]): #intensities are never less than zero!
492            print( "bad parameters:", (5*"%10.4e ") %(peakpos, innerbound, outerbound, k, y0))
493            print( len(intg), intg[:-1])
494            raise ValueError("Bad axial helper parameters")
495               
496        #now, make sure the underlying area is the exactly correct
497        #integral, without bumps due to discretizing the grid.
498        intg *= (exactintegral/(intg[:-1].sum()))
499
500        destination[idx0:idx1-1]+=intg[:-1]
501       
502        ## This is purely for debugging.  If collect_moment_errors is \a True,
503        #  compute exact vs. approximate moments.
504        if collect_moment_errors:
505            centroid2=(intg[:-1]*epsvals[idx0:idx1-1]).sum()/intg[:-1].sum()
506            moment_list.append((centroid2-exact_moment1)/dx)
507
508        return [idx0, idx1-1] #useful info for peak position
509
510    ## @brief return the \a I2 function
511    #  @param self self
512    #  @param Lx length of the xray filament
513    #  @param Ls length of the sample
514    #  @param Lr length of the receiver slit
515    #  @param R diffractometer length, assumed symmetrical
516    #  @param twotheta angle, in radians, of the center of the computation
517    #  @param beta offset angle
518    #  @param epsvals array of offsets from center of computation, in radians
519    #  @return ( \a epsvals, \a idxmin, \a idxmax, \a I2p, \a I2m ).
520    #   \a idxmin and \a idxmax are the full python-style bounds of the non-zero region of \a I2p and \a I2m.
521    #   \a I2p and \a I2m are I2+ and I2- from the paper, the contributions to the intensity.
522    def full_axdiv_I2(self, Lx=None, Ls=None, Lr=None, R=None, twotheta=None, beta=None,
523            epsvals=None
524        ):
525       
526        from math import sin, cos, tan, pi
527        beta1=(Ls-Lx)/(2*R) #Ch&Co after eq. 15abcd
528        beta2=(Ls+Lx)/(2*R) #Ch&Co after eq. 15abcd, corrected by KM
529       
530        eps0=beta*beta*math.tan(twotheta)/2 #after eq. 26 in Ch&Co
531       
532        if -beta2 <= beta < beta1:
533            z0p=Lx/2+beta*R*(1+1/cos(twotheta))
534        elif beta1 <= beta <= beta2:
535            z0p=Ls/2+beta*R/cos(twotheta)
536       
537        if -beta2 <= beta <= -beta1:
538            z0m=-Ls/2+beta*R/cos(twotheta)
539        elif -beta1 < beta <= beta2:
540            z0m=-Lx/2+beta*R*(1+1/cos(twotheta))
541
542        epsscale=tan(pi/2-twotheta)/(2*R*R) #=cotan(twotheta)...
543
544        eps1p=(eps0-epsscale*((Lr/2)-z0p)**2) #Ch&Co 18a&18b, KM sign correction
545        eps2p=(eps0-epsscale*((Lr/2)-z0m)**2)
546        eps2m=(eps0-epsscale*((Lr/2)+z0p)**2) #reversed eps2m and eps1m per KM R
547        eps1m=(eps0-epsscale*((Lr/2)+z0m)**2) #flip all epsilons
548
549        if twotheta>pi/2: #this set of inversions from KM 'R' code, simplified here
550            eps1p, eps2p, eps1m, eps2m = eps1m, eps2m, eps1p, eps2p
551
552        #identify ranges per Ch&Co 4.2.2 and table 1 and select parameters
553        #note table 1 is full of typos, but the minimized
554        #tests from 4.2.2 with redundancies removed seem fine.
555        if Lr > z0p-z0m:
556            if  z0p <= Lr/2 and z0m > -Lr/2: #beam entirely within slit
557                rng=1; ea=eps1p; eb=eps2p; ec=eps1m; ed=eps2m
558            elif (z0p > Lr/2 and z0m < Lr/2) or (z0m < -Lr/2 and z0p > -Lr/2):
559                rng=2; ea=eps2p; eb=eps1p; ec=eps1m; ed=eps2m
560            else:
561                rng=3; ea=eps2p; eb=eps1p; ec=eps1m; ed=eps2m
562        else:
563            if z0m < -Lr/2 and z0p > Lr/2: #beam hanging off both ends of slit, peak centered
564                rng=1; ea=eps1m; eb=eps2p; ec=eps1p; ed=eps2m
565            elif (-Lr/2 < z0m < Lr/2 and z0p > Lr/2) or (-Lr/2 < z0p < Lr/2 and z0m < -Lr/2): #one edge of beam within slit
566                rng=2; ea=eps2p; eb=eps1m; ec=eps1p; ed=eps2m
567            else:
568                rng=3; ea=eps2p; eb=eps1m; ec=eps1p; ed=eps2m
569
570        #now, evaluate function on bounds in table 1 based on ranges
571        #note: because of a sign convention in epsilon, the bounds all get switched
572       
573        def F1(dst, lower, upper, eea, eeb): #define them in our namespace so they inherit ea, eb, ec, ed, etc.
574            return self.axial_helper(destination=dst,
575                    innerbound=upper, outerbound=lower,
576                    epsvals=epsvals, peakpos=eps0,
577                    k=math.sqrt(abs(eps0-eeb))-math.sqrt(abs(eps0-eea)), y0=0
578                )
579        def F2(dst, lower, upper, eea):
580            return self.axial_helper(destination=dst,
581                    innerbound=upper, outerbound=lower,
582                    epsvals=epsvals, peakpos=eps0,
583                    k=math.sqrt(abs(eps0-eea)), y0=-1
584                )
585        def F3(dst, lower, upper, eea):
586            return self.axial_helper(destination=dst,
587                    innerbound=upper, outerbound=lower,
588                    epsvals=epsvals, peakpos=eps0,
589                    k=math.sqrt(abs(eps0-eea)), y0=+1
590                )
591        def F4(dst, lower, upper, eea):
592            #just like F2 but k and y0 negated
593            return self.axial_helper(destination=dst,
594                    innerbound=upper, outerbound=lower,
595                    epsvals=epsvals, peakpos=eps0,
596                    k=-math.sqrt(abs(eps0-eea)), y0=+1
597                )
598
599        I2p=self._I2p
600        I2p[:]=0
601        I2m=self._I2m
602        I2m[:]=0
603       
604        indices=[]
605        if rng==1:
606            indices+=F1(dst=I2p, lower=ea, upper=eps0, eea=ea, eeb=eb)
607            indices+=F2(dst=I2p, lower=eb, upper=ea,   eea=eb)
608            indices+=F1(dst=I2m, lower=ec, upper=eps0, eea=ec, eeb=ed)
609            indices+=F2(dst=I2m, lower=ed, upper=ec,   eea=ed)
610        elif rng==2:
611            indices+=F2(dst=I2p, lower=ea, upper=eps0, eea=ea)
612            indices+=F3(dst=I2m, lower=eb, upper=eps0, eea=ea)
613            indices+=F1(dst=I2m, lower=ec, upper=eb  , eea=ec, eeb=ed)
614            indices+=F2(dst=I2m, lower=ed, upper=ec  , eea=ed)
615        elif rng==3:
616            indices+=F4(dst=I2m, lower=eb, upper=ea  , eea=ea)
617            indices+=F1(dst=I2m, lower=ec, upper=eb  , eea=ec, eeb=ed)
618            indices+=F2(dst=I2m, lower=ed, upper=ec  , eea=ed)
619
620        idxmin=min(indices)
621        idxmax=max(indices)
622
623        return epsvals, idxmin, idxmax, I2p,I2m
624
625    ## @brief carry out the integral of \a I2 over \a beta and the Soller slits.
626    #  @param self self
627    #  @param Lx length of the xray filament
628    #  @param Ls length of the sample
629    #  @param Lr length of the receiver slit
630    #  @param R the (assumed symmetrical) diffractometer radius
631    #  @param twotheta angle, in radians, of the center of the computation
632    #  @param epsvals array of offsets from center of computation, in radians
633    #  @param sollerIdeg the full-width (both sides) cutoff angle of the incident Soller slit
634    #  @param sollerDdeg the full-width (both sides) cutoff angle of the detector Soller slit
635    #  @param nsteps the number of subdivisions for the integral
636    #  @param axDiv not used
637    #  @return the accumulated integral, a copy of a persistent buffer \a _axial
638    def full_axdiv_I3(self, Lx=None, Ls=None, Lr=None, R=None,
639        twotheta=None,
640        epsvals=None, sollerIdeg=None, sollerDdeg=None, nsteps=10, axDiv=""):
641        from math import sin, cos, tan, pi
642       
643        beta2=(Ls+Lx)/(2*R) #Ch&Co after eq. 15abcd, corrected by KM
644
645        if sollerIdeg is not None:
646            solIrad=sollerIdeg*pi/180/2
647            #solIfunc=numpy.frompyfunc(lambda x: max(0, 1.0-abs(x/solIrad)),1,1)
648            #note: solIfunc is only called with a scalar, no need for numpy-ized version, really
649            def solIfunc(x):
650                return numpy.clip(1.0-abs(x/solIrad),0,1)
651            beta2=min(beta2, solIrad) #no point going beyond Soller
652        else:
653            def solIfunc(x):
654                return numpy.ones_like(x)
655        if sollerDdeg is not None:
656            solDrad=sollerDdeg*pi/180/2
657            #solDfunc=numpy.frompyfunc(lambda x: max(0, 1.0-abs(x/solDrad)),1,1)
658            def solDfunc(x):
659                return numpy.clip(1.0-abs(x/solDrad),0,1)
660        else:
661            def solDfunc(x):
662                return numpy.ones_like(x)
663
664        accum=self._axial
665        accum[:]=0
666       
667        if twotheta > pi/2:
668            tth1=pi-twotheta
669        else:
670            tth1=twotheta
671
672        for iidx in range(nsteps):
673            beta=beta2*iidx/float(nsteps)
674           
675            eps, idxmin, idxmax, I2p, I2m=self.full_axdiv_I2(
676                Lx=Lx,
677                Lr=Lr,
678                Ls=Ls,
679                beta=beta, R=R,
680                twotheta=twotheta, epsvals=epsvals
681            )
682           
683            eps0=beta*beta*math.tan(twotheta)/2 #after eq. 26 in Ch&Co
684           
685            gamma0=beta/cos(tth1)
686            deps=self._f0buf[idxmin:idxmax]
687            deps[:]=eps0
688            deps-=epsvals[idxmin:idxmax]
689            deps *= 2*math.tan(twotheta)
690            #check two channels on each end for negative argument.
691            deps[-1]=max(deps[-1],0)
692            deps[0]=max(deps[0],0)
693            if len(deps) >= 2:
694                deps[-2]=max(deps[-2],0)
695                deps[1]=max(deps[1],0)
696           
697            gamarg=numpy.sqrt(deps, deps) #do sqrt in place for speed
698            #still need to convert these to in-place
699            gamp=gamma0+gamarg
700            gamm=gamma0-gamarg
701           
702            if iidx==0 or iidx==nsteps-1: weight=1.0 #trapezoidal rule weighting
703            else: weight=2.0
704
705            #sum into the accumulator only channels which can be non-zero
706            #do scaling in-place to save a  lot of slow array copying
707            I2p[idxmin:idxmax] *= solDfunc(gamp)
708            I2p[idxmin:idxmax] *=(weight*solIfunc(beta))
709            accum[idxmin:idxmax]+=I2p[idxmin:idxmax]
710            I2m[idxmin:idxmax] *= solDfunc(gamm)
711            I2m[idxmin:idxmax] *=(weight*solIfunc(beta))
712            accum[idxmin:idxmax]+=I2m[idxmin:idxmax]
713
714        #keep this normalized
715        K = 2 * R*R * abs(tan(twotheta))
716        accum *= K
717       
718        return accum
719   
720    ## @brief compute the Fourier transform of the axial divergence component
721    #  @param self self
722    #  @return the transform buffer, or \a None if this is being ignored
723    def conv_axial(self):
724        me=self.get_function_name() #the name of this convolver,as a string
725        if self.param_dicts[me].get("axDiv",None) is None:
726            return None
727        kwargs={}
728        kwargs.update(self.param_dicts[me]) #get all of our parameters
729        kwargs.update(self.param_dicts["conv_global"])
730        if "equatorial_divergence_deg" in kwargs:
731            del kwargs["equatorial_divergence_deg"] #not used
732       
733        flag, axfn = self.get_conv(me, kwargs, numpy.complex)
734        if flag: return axfn #already up to date if first return is True
735       
736        xx=type("data",(), kwargs)
737        if xx.axDiv!="full" or xx.twotheta0_deg==90.0: #no axial divergence, transform of delta fn
738            axfn[:]=1
739            return axfn
740        else:
741            axbuf=self.full_axdiv_I3(
742                nsteps=xx.n_integral_points,
743                epsvals=self.epsilon,
744                Lx=xx.slit_length_source,
745                Lr=xx.slit_length_target,
746                Ls=xx.length_sample,
747                sollerIdeg=xx.angI_deg,
748                sollerDdeg=xx.angD_deg,
749                R=xx.diffractometer_radius,
750                twotheta=xx.twotheta0
751            )
752        axfn[:]=best_rfft(axbuf)
753       
754        return axfn
755
756    ## @brief compute the Fourier transform of the rectangular tube tails function
757    #  @param self self
758    #  @return the transform buffer, or \a None if this is being ignored
759    def conv_tube_tails(self):
760        me=self.get_function_name() #the name of this convolver,as a string
761        kwargs={}
762        kwargs.update(self.param_dicts[me]) #get all of our parameters
763        if not kwargs:
764            return None #no convolver
765        #we also need the diffractometer radius from the global space
766        kwargs["diffractometer_radius"]=self.param_dicts["conv_global"]["diffractometer_radius"]
767        flag, tailfn = self.get_conv(me, kwargs, numpy.complex)
768        if flag: return tailfn #already up to date
769       
770        #tube_tails is (main width, left width, right width, intensity),
771        #so peak is raw peak + tophat centered at (left width+ right width)
772        #with area intensity*(right width-left width)/main_width
773        #check this normalization!
774        #note: this widths are as defined by Topas... really I think it should be
775        #x/(2*diffractometer_radius) since the detector is 2R from the source,
776        #but since this is just a fit parameter, we'll defin it as does Topas
777        xx=type("data",(), kwargs) #allow dotted notation
778       
779        tail_eps=(xx.tail_right-xx.tail_left)/xx.diffractometer_radius
780        main_eps=xx.main_width/xx.diffractometer_radius
781        tail_center=(xx.tail_right+xx.tail_left)/xx.diffractometer_radius/2.0
782        tail_area=xx.tail_intens*(xx.tail_right-xx.tail_left)/xx.main_width
783
784        cb1=self._cb1
785        rb1=self._rb1
786       
787        cb1.real=0
788        cb1.imag = self.omega_vals
789        cb1.imag *= tail_center #sign is consistent with Topas definition
790        numpy.exp(cb1, tailfn) #shifted center, computed into tailfn
791       
792        rb1[:]=self.omega_vals
793        rb1 *= (tail_eps/2/math.pi)
794        rb1=numpy.sinc(rb1)
795        tailfn *= rb1
796        tailfn *= tail_area #normalize
797       
798        rb1[:]=self.omega_vals
799        rb1 *= (main_eps/2/math.pi)
800        rb1=numpy.sinc(rb1)
801        tailfn += rb1 #add central peak
802        return tailfn
803   
804    ## @brief a utility to compute a transformed tophat function and save it in a convolver buffer
805    #  @param self self
806    #  @param name the name of the convolver cache buffer to update
807    #  @param width the width in 2-theta space of the tophat
808    #  @return the updated convolver buffer, or \a None if the width was \a None
809    def general_tophat(self, name="", width=None):
810        """handle all centered top-hats"""
811        if width is None:
812            return #no convolver
813        flag, conv = self.get_conv(name, width, numpy.float)
814        if flag: return conv #already up to date
815        rb1=self._rb1
816        rb1[:]=self.omega_vals
817        rb1 *= (width/2/math.pi)
818        conv[:]=numpy.sinc(rb1)
819        return conv
820
821    ## A dictionary of default parameters for conv_emissions,
822    #  used to seed a GUI which can harvest this for names, descriptions, and initial values
823    info_emission=dict(
824        group_name="Incident beam and crystal size",
825        help="this should be help information",
826        param_info=dict(
827            emiss_wavelengths=("wavelengths (m)", (1.58e-10,)),
828            emiss_intensities=("relative intensities", (1.00,)),
829            emiss_lor_widths=("Lorenztian emission fwhm (m)",(1e-13,)),
830            emiss_gauss_widths=("Gaussian emissions fwhm (m)", (1e-13,)),
831            crystallite_size_gauss=("Gaussian crystallite size fwhm (m)", 1e-6),
832            crystallite_size_lor=("Lorentzian crystallite size fwhm (m)", 1e-6),
833        )
834    )
835   
836    ## @brief format the emission spectrum and crystal size information
837    #  @param self self
838    #  @return the formatted information
839    def str_emission(self):
840        dd=self.param_dicts["conv_emission"]
841        if not dd: return "No emission spectrum"
842        dd.setdefault("crystallite_size_lor",1e10)
843        dd.setdefault("crystallite_size_gauss",1e10)
844        dd.setdefault("strain_lor",0)
845        dd.setdefault("strain_gauss",0)
846        xx=type("data",(), dd)
847        spect=numpy.array((
848            xx.emiss_wavelengths, xx.emiss_intensities,
849            xx.emiss_lor_widths, xx.emiss_gauss_widths))
850        spect[0]*=1e10*self.length_scale_m #convert to Angstroms, like Topas
851        spect[2]*=1e13*self.length_scale_m #milli-Angstroms
852        spect[3]*=1e13*self.length_scale_m #milli-Angstroms
853        nm=1e9*self.length_scale_m
854        items=["emission and broadening:"]
855        items.append("spectrum=\n"+str(spect.transpose()))
856        items.append("crystallite_size_lor (nm): %.5g" % (xx.crystallite_size_lor*nm))
857        items.append("crystallite_size_gauss (nm): %.5g" % (xx.crystallite_size_gauss*nm))
858        items.append("strain_lor: %.5g" % xx.strain_lor)
859        items.append("strain_gauss: %.5g" % xx.strain_gauss)
860        return '\n'.join(items)
861   
862    ## @brief compute the emission spectrum and (for convenience) the particle size widths
863    #  @param self self
864    #  @return the convolver for the emission and particle sizes
865    #  @note the particle size and strain stuff here is just to be consistent with \a Topas
866    #  and to be vaguely efficient about the computation, since all of these
867    #  have the same general shape.
868    def conv_emission(self):
869        """handle emission spectrum and crystal size together, since it makes sense"""
870        me=self.get_function_name() #the name of this convolver,as a string
871        kwargs={}
872        kwargs.update(self.param_dicts[me]) #get all of our parameters
873        kwargs.update(self.param_dicts["conv_global"])
874        #if the crystallite size and strain parameters are not set, set them to
875        #values that make their corrections disappear
876        kwargs.setdefault("crystallite_size_lor",1e10)
877        kwargs.setdefault("crystallite_size_gauss",1e10)
878        kwargs.setdefault("strain_lor",0)
879        kwargs.setdefault("strain_gauss",0)
880
881        #convert arrays to lists for key checking
882        key={}
883        key.update(kwargs)
884        for k,v in key.items():
885            if hasattr(v,'tolist'):
886                key[k]=v.tolist()
887   
888        flag, emiss = self.get_conv(me, key, numpy.complex)
889        if flag: return emiss #already up to date
890       
891        xx=type("data", (), kwargs) #make it dot-notation accessible
892       
893        epsilon0s = 2 * numpy.arcsin(xx.emiss_wavelengths / (2.0 * xx.d)) - xx.twotheta0
894        theta=xx.twotheta0/2
895        ## Emission profile FWHM + crystallite broadening (scale factors are Topas choice!) (Lorentzian)
896        #note:  the strain broadenings in Topas are expressed in degrees 2theta, must convert to radians(theta) with pi/360
897        widths = (
898            (xx.emiss_lor_widths/xx.emiss_wavelengths) * math.tan(theta) +
899            xx.strain_lor*(math.pi/360) * math.tan(theta) +
900            (xx.emiss_wavelengths / (2*xx.crystallite_size_lor * math.cos(theta)))
901        )
902        #save weighted average width for future reference in periodicity fixer
903        self.lor_widths[me]=sum(widths*xx.emiss_intensities)/sum(xx.emiss_intensities)
904        #gaussian bits add in quadrature
905        gfwhm2s=(
906            ((2*xx.emiss_gauss_widths/xx.emiss_wavelengths) * math.tan(theta))**2 +
907            (xx.strain_gauss*(math.pi/360) * math.tan(theta))**2 +
908            (xx.emiss_wavelengths / (xx.crystallite_size_gauss * math.cos(theta)))**2
909        )
910       
911        # note that the Fourier transform of a lorentzian with FWHM 2a
912        # is exp(-abs(a omega))
913        #now, the line profiles in Fourier space have to have phases
914        #carefully handled to put the lines in the right places.
915        #note that the transform of f(x+dx)=exp(i omega dx) f~(x)
916        omega_vals=self.omega_vals
917        for wid, gfwhm2, eps, intens in zip(widths, gfwhm2s, epsilon0s, xx.emiss_intensities):
918            xvals=numpy.clip(omega_vals*(-wid),-100,0)
919            sig2=gfwhm2/(8*math.log(2.0)) #convert fwhm**2 to sigma**2
920            gxv=numpy.clip((sig2/-2.0)*omega_vals*omega_vals,-100,0)
921            emiss += numpy.exp(xvals+gxv+complex(0,-eps)*omega_vals)*intens
922        return emiss
923   
924    ## @brief compute the convolver for the flat-specimen correction
925    #  @param self self
926    #  @return the convolver
927    def conv_flat_specimen(self):
928        """handle flat specimen"""
929        me=self.get_function_name() #the name of this convolver,as a string
930        equatorial_divergence_deg=self.param_dicts["conv_global"].get("equatorial_divergence_deg",None)
931        if not equatorial_divergence_deg: return None
932        twotheta0=self.param_dicts["conv_global"]["twotheta0"]
933        key=(twotheta0, equatorial_divergence_deg)
934        flag, conv = self.get_conv(me, key, numpy.complex)
935        if flag: return conv #already up to date
936       
937        #Flat-specimen error, from Cheary, Coelho & Cline 2004 NIST eq. 9 & 10
938        #compute epsm in radians from eq. divergence in degrees
939        #to make it easy to use the axial_helper to compute the function
940        epsm = ((equatorial_divergence_deg*math.pi/180)**2)/math.tan(twotheta0/ 2.0)/2.0
941        eqdiv=self._epsb2
942        eqdiv[:]=0
943        dtwoth=(self.twothetasamples[1]-self.twothetasamples[0])
944        idx0, idx1=self.axial_helper(destination=eqdiv,
945            outerbound=-epsm,
946            innerbound=0,
947            epsvals=self.epsilon,
948            peakpos=0, k=dtwoth/(2.0*math.sqrt(epsm)))
949
950        conv[:]=best_rfft(eqdiv)
951        conv[1::2] *= -1 #flip center
952        return conv
953   
954    ## @brief compute the sample transparency correction, including the finite-thickness version
955    #  @param self self
956    #  @return the convolver
957    def conv_absorption(self):
958        """handle transparency"""
959        me=self.get_function_name() #the name of this convolver,as a string
960        kwargs={}
961        kwargs.update(self.param_dicts[me]) #get all of our parameters
962        if not kwargs: return None
963        kwargs["twotheta0"]=self.param_dicts["conv_global"]["twotheta0"]
964        kwargs["diffractometer_radius"]=self.param_dicts["conv_global"]["diffractometer_radius"]
965
966        flag, conv = self.get_conv(me, kwargs, numpy.complex)
967        if flag: return conv #already up to date
968        xx=type("data", (), kwargs) #make it dot-notation accessible
969       
970        #absorption, from Cheary, Coelho & Cline 2004 NIST eq. 12,
971        #EXCEPT delta = 1/(2*mu*R) instead of 2/(mu*R)
972        #from Mathematica, unnormalized transform is
973        #(1-exp(epsmin*(i w + 1/delta)))/(i w + 1/delta)
974        delta = math.sin(xx.twotheta0)/(2*xx.absorption_coefficient*xx.diffractometer_radius)
975        #arg=(1/delta)+complex(0,-1)*omega_vals
976        cb=self._cb1
977        cb.imag=self.omega_vals
978        cb.imag*=-1
979        cb.real=1/delta
980        numpy.reciprocal(cb,conv) #limit for thick samples=1/(delta*arg)
981        conv *= 1.0/delta #normalize
982        if kwargs.get("sample_thickness", None) is not None: #rest of transform of function with cutoff
983            epsmin = -2.0*xx.sample_thickness*math.cos(xx.twotheta0/2.0)/xx.diffractometer_radius
984            cb*=epsmin
985            numpy.expm1(cb,cb)
986            cb*=-1
987            conv *= cb
988        return conv
989   
990    ## @brief compute the peak shift due to sample displacement and the \a 2theta zero
991    #  @param self self
992    #  @return the convolver
993    def conv_displacement(self):
994        """handle displacements and peak offset from center"""
995        me=self.get_function_name() #the name of this convolver,as a string
996        kwargs=self.param_dicts[me]
997        twotheta0=self.param_dicts["conv_global"]["twotheta0"]
998        diffractometer_radius=self.param_dicts["conv_global"]["diffractometer_radius"]
999        specimen_displacement=kwargs.get("specimen_displacement", 0.0)
1000        if specimen_displacement is None: specimen_displacement=0.0
1001        zero_error_deg=kwargs.get("zero_error_deg", 0.0)
1002        if zero_error_deg is None: zero_error_deg=0.0
1003               
1004        flag, conv = self.get_conv(me,
1005            (twotheta0, diffractometer_radius,
1006            specimen_displacement, zero_error_deg),
1007            numpy.complex)
1008        if flag: return conv#already up to date
1009
1010        delta=-2*math.cos(twotheta0/2.0)*specimen_displacement/diffractometer_radius
1011        conv.real=0
1012        conv.imag=self.omega_vals
1013        #convolver *= numpy.exp(complex(0,-delta-zero_error_deg*pi/180.0)*omega_vals)
1014        conv.imag *= (-delta-zero_error_deg*math.pi/180.0-(twotheta0-self.twotheta_window_center))
1015        numpy.exp(conv,conv)
1016        return conv
1017
1018    ## @brief compute the rectangular convolution for the receiver slit or SiPSD pixel size
1019    #  @param self self
1020    #  @return the convolver
1021    def conv_receiver_slit(self):
1022        """ receiver slit width or si psd pixel width"""
1023        me=self.get_function_name() #the name of this convolver,as a string
1024        #The receiver slit convolution is a top-hat of angular half-width
1025        #a=(slit_width/2)/diffractometer_radius
1026        #which has Fourier transform of sin(a omega)/(a omega)
1027        #NOTE! numpy's sinc(x) is sin(pi x)/(pi x), not sin(x)/x
1028        if self.param_dicts[me].get("slit_width",None) is None:
1029            return None
1030       
1031        epsr=(self.param_dicts["conv_receiver_slit"]["slit_width"]/
1032            self.param_dicts["conv_global"]["diffractometer_radius"])
1033        return self.general_tophat(me, epsr)
1034
1035    ## @brief compute the convolver for the integral of defocusing of the face of an Si PSD
1036    #  @param self self
1037    #  @return the convolver
1038    def conv_si_psd(self):
1039        #omega offset defocussing from Cheary, Coelho & Cline 2004 eq. 15
1040        #expressed in terms of a Si PSD looking at channels with vertical offset
1041        #from the center between psd_window_lower_offset and psd_window_upper_offset
1042        #do this last, because we may ultimately take a list of bounds,
1043        #and return a list of convolutions, for efficiency
1044        me=self.get_function_name() #the name of this convolver,as a string
1045        kwargs={}
1046        kwargs.update(self.param_dicts[me]) #get all of our parameters
1047        if not kwargs: return None
1048        kwargs.update(self.param_dicts["conv_global"])
1049
1050        flag, conv = self.get_conv(me, kwargs, numpy.float)
1051        if flag: return conv #already up to date
1052       
1053        xx=type("data",(), kwargs)
1054       
1055        if not xx.equatorial_divergence_deg or not xx.si_psd_window_bounds:
1056            #if either of these is zero or None, convolution is trivial
1057            conv[:]=1
1058            return conv
1059       
1060        psd_lower_window_pos, psd_upper_window_pos = xx.si_psd_window_bounds
1061        dthl=psd_lower_window_pos/xx.diffractometer_radius
1062        dthu=psd_upper_window_pos/xx.diffractometer_radius
1063        alpha=xx.equatorial_divergence_deg*math.pi/180
1064        argscale=alpha/(2.0*math.tan(xx.twotheta0/2))
1065        from scipy.special import sici #get the sine and cosine integral
1066        #WARNING si(x)=integral(sin(x)/x), not integral(sin(pi x)/(pi x))
1067        #i.e. they sinc function is not consistent with the si function
1068        #whence the missing pi in the denominator of argscale
1069        rb1=self._rb1
1070        rb2=self._rb2
1071        rb3=self._rb3
1072        rb1[:]=self.omega_vals
1073        rb1 *= argscale*dthu
1074        sici(rb1,conv,rb3) #gets both sine and cosine integral, si in conv
1075        if dthl:  #no need to do this if the lower bound is 0
1076            rb1[:]=self.omega_vals
1077            rb1 *= argscale*dthl
1078            sici(rb1,rb2,rb3) #gets both sine and cosine integral, si in rb2
1079            conv-=rb2
1080        conv[1:] /= self.omega_vals[1:]
1081        conv *= 1/argscale
1082        conv[0] = dthu-dthl #fix 0/0 with proper area
1083        return conv
1084   
1085    ## @brief compute the convolver to smooth the final result with a Gaussian before downsampling.
1086    #  @param self self
1087    #  @return the convolver
1088    def conv_smoother(self):
1089        #create a smoother for output result, independent of real physics, if wanted
1090        me=self.get_function_name() #the name of this convolver,as a string
1091        if not self.output_gaussian_smoother_bins_sigma: return # no smoothing
1092        flag, buf=self.get_conv(me, self.output_gaussian_smoother_bins_sigma,
1093            format=numpy.float)
1094        if flag: return buf #already computed
1095        buf[:]=self.omega_vals
1096        buf*=(self.output_gaussian_smoother_bins_sigma*(
1097            self.twothetasamples[1]-self.twothetasamples[0]))
1098        buf *= buf
1099        buf *=-0.5
1100        numpy.exp(buf, buf)
1101        return buf
1102   
1103    ## @brief execute all the convolutions; if convolver_names is None, use everything
1104    #    we have, otherwise, use named convolutions.
1105    #  @param self self
1106    #  @param convolver_names a list of convolvers to select. If \a None, use all found convolvers.
1107    #  @param compute_derivative if \a True, also return d/dx(function) for peak position fitting
1108    #  @return a profile_data object with much information about the peak
1109    def compute_line_profile(self, convolver_names=None, compute_derivative=False, return_convolver=False):
1110       
1111        #create a function which is the Fourier transform of the
1112        #combined convolutions of all the factors
1113        from numpy import sin, cos, tan, arcsin as asin, arccos as acos
1114        from math import pi
1115
1116        anglemode=self.anglemode #="d" if we are using 'd' space, "twotheta" if using twotheta
1117       
1118        ## the rough center of the spectrum, used for things which need it. Copied from global convolver.
1119        self.dominant_wavelength=dominant_wavelength=self.param_dicts["conv_global"].get("dominant_wavelength",None)
1120       
1121        if anglemode=="twotheta":
1122            twotheta0_deg=self.param_dicts["conv_global"]["twotheta0_deg"]
1123            twotheta0=math.radians(twotheta0_deg)
1124            d=dominant_wavelength/(2*sin(twotheta0/2.0))
1125        else:
1126            d=self.param_dicts["conv_global"]["d"]
1127            twotheta0 = 2 * asin(dominant_wavelength/ (2.0 * d))
1128            twotheta0_deg=math.degrees(twotheta0)
1129       
1130        self.set_parameters(d=d,twotheta0=twotheta0,twotheta0_deg=twotheta0_deg) #set these in global namespace
1131
1132        if convolver_names is None:
1133            convolver_names=self.convolver_funcs.keys() #get all names
1134       
1135        #run through the name list, and call the convolver to harvest its result
1136        conv_list=[self.convolver_funcs[x]() for x in convolver_names]
1137
1138        #now, multiply everything together
1139        convolver=self._cb1 #get a complex scratch buffer
1140        convolver[:]=1 #initialize
1141        for c in conv_list: #accumulate product
1142            if c is not None: convolver *= c
1143
1144        if convolver[1].real > 0: #recenter peak!
1145            convolver[1::2]*=-1
1146       
1147        peak=best_irfft(convolver)
1148
1149        #now, use the trick from Mendenhall JQSRT Voigt paper to remove periodic function correction
1150        #JQSRT 105 number 3 July 2007 p. 519 eq. 7
1151        emiss_intensities=self.param_dicts["conv_emission"]["emiss_intensities"]
1152        correction_width=2*sum(self.lor_widths.values()) #total lor widths, created by the various colvolvers
1153       
1154        d2p=2.0*pi/self.window_fullwidth
1155        alpha=correction_width/2.0 #be consistent with convolver
1156        mu=(peak*self.twothetasamples).sum()/peak.sum() #centroid
1157        dx=self.twothetasamples-mu
1158        eps_corr1=(math.sinh(d2p*alpha)/self.window_fullwidth)/(math.cosh(d2p*alpha)-numpy.cos(d2p*dx))
1159        eps_corr2=(alpha/pi)/(dx*dx+alpha*alpha)
1160        corr=(convolver[0].real/numpy.sum(eps_corr2))*(eps_corr1-eps_corr2)
1161        peak -= corr
1162       
1163        peak*=self.window_fullwidth/(self.twotheta_output_points/self.oversampling) #scale to area
1164       
1165        if compute_derivative:
1166            #this is useful
1167            convolver *= self.omega_vals
1168            convolver *= complex(0,1)
1169            deriv = best_irfft(convolver)
1170            deriv*=self.window_fullwidth/(self.twotheta_output_points/self.oversampling)
1171            deriv=deriv[::self.oversampling]
1172        else:
1173            deriv=None
1174       
1175        result=profile_data(twotheta0_deg=twotheta0*180/math.pi,
1176            twotheta=self.twothetasamples[::self.oversampling],
1177            omega_inv_deg=self.omega_inv_deg[:self.twotheta_output_points//2+1],
1178            twotheta_deg=self.twothetasamples_deg[::self.oversampling],
1179            peak=peak[::self.oversampling],
1180            derivative=deriv
1181        )
1182
1183        if return_convolver:
1184            result.add_symbol(convolver=convolver[:self.twotheta_output_points//2+1])
1185           
1186        return result
1187   
1188    ## @brief do some cleanup to make us more compact;
1189    #  Instance can no longer be used after doing this, but can be pickled.
1190    #  @param self self
1191    def self_clean(self):
1192        clean=self._clean_on_pickle
1193        pd=dict()
1194        pd.update(self.__dict__)       
1195        for thing in pd.keys():
1196            x=getattr(self, thing)
1197            if id(x) in clean:
1198                delattr(self,thing)
1199        #delete extra attributes cautiously, in case we have already been cleaned
1200        for k in ('convolver_funcs','convolvers','factors','convolution_history'):
1201            if pd.pop(k,None) is not None:
1202                delattr(self,k)
1203
1204    ## @brief return information for pickling
1205    #  @param self self
1206    #  @return dictionary of sufficient information to reanimate instance.
1207    #
1208    #  Removes transient data from cache of shadow copy so resulting object is fairly compact.
1209    #  This does not affect the state of the actual instance.
1210    def __getstate__(self):
1211        #do some cleanup on state before we get pickled
1212        #(even if main class not cleaned)
1213        clean=self._clean_on_pickle
1214        pd=dict()
1215        pd.update(self.__dict__)
1216        for thing in pd.keys():
1217            x=getattr(self, thing)
1218            if id(x) in clean:
1219                del pd[thing]
1220        #delete extra attributes cautiously, in case we have alread been cleaned
1221        for k in ('convolver_funcs','convolvers','factors','convolution_history'):
1222            pd.pop(k,None)
1223        return pd
1224
1225    ## @brief reconstruct class from pickled information
1226    #  @param self an empty class instance
1227    #  @param setdict dictionary from FP_profile.__getstate__()
1228    #
1229    #  This rebuilds the class instance so it is ready to use on unpickling.
1230    def __setstate__(self, setdict):
1231        self.__init__(anglemode=setdict["anglemode"],
1232            output_gaussian_smoother_bins_sigma=setdict["output_gaussian_smoother_bins_sigma"],
1233            oversampling=setdict["oversampling"]
1234        )
1235        for k,v in setdict.items():
1236            setattr(self,k,v)
1237        self.set_window(
1238            twotheta_window_center_deg=self.twotheta_window_center_deg,
1239            twotheta_window_fullwidth_deg=self.twotheta_window_fullwidth_deg,
1240            twotheta_output_points=self.twotheta_output_points
1241        )
1242        self.lor_widths=setdict["lor_widths"] #override clearing of this by set_window
1243
1244##
1245#  @brief A single function call interface to an instance of the class.
1246#  This is only intended for debugging and simple calculation.  Since
1247#  it re-creates the class each time, there is no cacheing of results.  It is
1248#  fairly inefficient to do stuff this way.
1249#  @param d d-spacing, if present.  Either a d-spacing and wavelength or an angle must be provided
1250#  @param twotheta_deg the peak center, if a d-spacing is not given.
1251#  @return an object with much information about a line profile.
1252def fourier_line_profile(
1253        d=None,  twotheta_deg=None,
1254        crystallite_size_lor=None,
1255        crystallite_size_gauss=None,
1256        strain_lor=0.0, strain_gauss=0.0,
1257        wavelength=None,
1258        slit_width=None,
1259        slit_length_target=None,
1260        slit_length_source=None,
1261        length_sample=None,
1262        angI_deg = None,
1263        angD_deg = None,
1264        axDiv = "simple",
1265        diffractometer_radius=None,
1266        si_psd_window_bounds=None,
1267        equatorial_divergence_deg=None, absorption_coefficient=None,
1268        sample_thickness=None,
1269        specimen_displacement=None, zero_error_deg=None, 
1270        target_width=None, specimen_tilt=None, defocus_delta_omega_deg=None,
1271        mat_wavelengths=None, mat_lor_widths=None, mat_gauss_widths=None, mat_intensities=None,
1272        tube_tails=None,
1273        window_fullwidth_deg=1.,
1274        twotheta_output_points=1000,
1275        n_integral_points=10,
1276        output_gaussian_smoother_bins_sigma=1.0,
1277        oversampling=10
1278    ):
1279    if twotheta_deg is not None: #give priority to angle
1280        anglemode="twotheta"
1281        d=0.0
1282    else:
1283        anglemode="d"
1284        twotheta_deg = 2 * math.degrees(math.asin(wavelength/ (2.0 * d)))
1285   
1286    p=FP_profile(anglemode=anglemode,
1287        output_gaussian_smoother_bins_sigma=output_gaussian_smoother_bins_sigma,
1288        oversampling=oversampling
1289    )
1290    #put the compute window in the right place, using old convention
1291    #centering window on line
1292    p.set_window(
1293        twotheta_output_points=twotheta_output_points,
1294        twotheta_window_center_deg=twotheta_deg,
1295        twotheta_window_fullwidth_deg=window_fullwidth_deg,
1296    )
1297
1298    #set parameters which are shared by many things
1299    p.set_parameters(d=d,twotheta0_deg=twotheta_deg, dominant_wavelength=wavelength,
1300        equatorial_divergence_deg=equatorial_divergence_deg,
1301        diffractometer_radius=diffractometer_radius)
1302    #set parameters for each convolver
1303    p.set_parameters(convolver="emission",
1304        emiss_wavelengths=mat_wavelengths,
1305        emiss_intensities=mat_intensities,
1306        emiss_gauss_widths=mat_gauss_widths,
1307        emiss_lor_widths=mat_lor_widths,
1308        crystallite_size_gauss=crystallite_size_gauss,
1309        crystallite_size_lor=crystallite_size_lor,
1310        strain_lor=strain_lor, strain_gauss=strain_gauss
1311    )
1312    p.set_parameters(convolver="displacement",
1313        zero_error_deg=zero_error_deg, specimen_displacement=specimen_displacement
1314    )
1315    if axDiv is not None:
1316        p.set_parameters(convolver="axial",
1317            angI_deg=angI_deg, angD_deg=angD_deg,
1318            axDiv=axDiv, slit_length_source=slit_length_source,
1319            slit_length_target=slit_length_target,
1320            length_sample=length_sample,
1321            n_integral_points=n_integral_points,
1322        )
1323    if absorption_coefficient is not None:
1324        p.set_parameters(convolver="absorption",
1325            absorption_coefficient=absorption_coefficient,
1326            sample_thickness=sample_thickness,
1327        )
1328    if si_psd_window_bounds is not None:
1329        p.set_parameters(convolver="si_psd",
1330         si_psd_window_bounds=si_psd_window_bounds
1331        )
1332    p.set_parameters(convolver="receiver_slit",
1333     slit_width=slit_width
1334    )
1335    if tube_tails is not None:
1336        main_width, tail_left, tail_right, tail_intens=tube_tails
1337        p.set_parameters(convolver="tube_tails",
1338            main_width=main_width, tail_left=tail_left,
1339            tail_right=tail_right, tail_intens=tail_intens
1340        )
1341
1342    res=p.compute_line_profile()
1343    res.add_symbol(profile_class=p)
1344
1345    return res
1346
1347if __name__=='__main__' and sys.argv.pop(-1)=='plot':
1348   
1349    ## fixed parameters
1350    mat_wavelengths = numpy.array((1.540596e-10, 1.540596e-10*1.001))
1351    mat_lor_widths = numpy.array((1.5e-13/5, 1.5e-13/5))
1352    mat_gauss_widths = numpy.array((1.5e-13/5, 1.5e-13/5))
1353    mat_intensities = numpy.array((1.0,0.0))
1354    collect_moment_errors=True
1355   
1356    try:
1357        from scardi_leoni_polydisperse import scardi_leoni_lognormal_polydisperse
1358        ## @brief subclass with scardi & leoni mixin for testing here;
1359        ## this shows how to trivially extend the base class with new modules
1360        class FP_with_poly(FP_profile, scardi_leoni_lognormal_polydisperse):
1361            pass
1362        useclass=FP_with_poly
1363    except:
1364        useclass=FP_profile
1365
1366    p=useclass(anglemode="twotheta",
1367            output_gaussian_smoother_bins_sigma=1.0,
1368            oversampling=2
1369        )
1370    p.debug_cache=False
1371    #set parameters for each convolver
1372    p.set_parameters(convolver="emission",
1373        emiss_wavelengths=mat_wavelengths,
1374        emiss_intensities=mat_intensities,
1375        emiss_gauss_widths=mat_gauss_widths,
1376        emiss_lor_widths=mat_lor_widths,
1377        crystallite_size_gauss=6000e-9,
1378        crystallite_size_lor=6000e-9
1379    )
1380    p.set_parameters(convolver="axial",
1381        axDiv="full", slit_length_source=15e-3,
1382        slit_length_target=10e-3,
1383        length_sample=12e-3,
1384        n_integral_points=10
1385    )
1386    p.set_parameters(convolver="absorption",
1387        absorption_coefficient=800.0*100, #like LaB6, in m^(-1)
1388        sample_thickness=1e-3,
1389    )
1390    p.set_parameters(convolver="si_psd",
1391     si_psd_window_bounds=(0.,2.5e-3)
1392    )
1393    p.set_parameters(convolver="receiver_slit",
1394     slit_width=75e-6,
1395    )
1396
1397    p.set_parameters(convolver="tube_tails",
1398        main_width=200e-6, tail_left=-1e-3,tail_right=1e-3, tail_intens=0.001
1399    )
1400
1401    if useclass != FP_profile: #we have scardi & leoni
1402        from scardi_leoni_polydisperse import H_params_hex
1403        p.set_parameters(convolver="scardi_leoni_lognormal_polydisperse",
1404            hkl=(1,0,0), shape_function=H_params_hex,
1405            crystallite_size=200.e-9, sigma=0.5, flat_is_100=True, 
1406            crystallite_height_over_edge=2., lattice_c_over_a=1.,
1407        )
1408
1409    import time
1410    t_start=time.time()
1411   
1412    profiles=[]
1413    for ww in (20,):
1414        for twotheta_x in range(20,161,20):
1415            #put the compute window in the right place and clear all histories
1416            p.set_window(twotheta_output_points=5000,
1417                twotheta_window_center_deg=twotheta_x,
1418                twotheta_window_fullwidth_deg=ww,
1419            )
1420            #set parameters which are shared by many things
1421            p.set_parameters(twotheta0_deg=twotheta_x,
1422                dominant_wavelength=mat_wavelengths[0],
1423                diffractometer_radius=217.5e-3)
1424       
1425            for anglescale in (1.,2.,3.):
1426                p.set_parameters(equatorial_divergence_deg=0.5*anglescale)
1427                p.set_parameters(convolver="axial",
1428                    angI_deg=2.5*anglescale, angD_deg=2.5*anglescale,
1429                )
1430                profiles.append((twotheta_x, 0.5*anglescale, 2.5*anglescale, p.compute_line_profile()))
1431
1432            print(p,file=sys.stderr)
1433    print("execution time=", time.time()-t_start, file=sys.stderr)
1434
1435    from matplotlib import pyplot as plt
1436   
1437    datasets=[]
1438    for idx, (twotheta, eqdiv, soller, result) in enumerate(profiles):
1439        #normalize both to integral=sum/dx
1440        dx1=(result.twotheta_deg[1]-result.twotheta_deg[0])
1441        peaknorm=40*result.peak/dx1
1442        plt.plot(result.twotheta_deg, peaknorm, label="%.1f : %.1f" % (eqdiv, soller) )
1443   
1444    plt.legend()
1445    plt.show()
1446   
1447    hh=numpy.histogram(moment_list, bins=20)
1448    print("error histogram=", hh)
1449    xl=(hh[1][1:]+hh[1][:-1])*0.5 #bin centers
1450    print("error mean, rms=", (xl*hh[0]).sum()/hh[0].sum(), math.sqrt((xl**2*hh[0]).sum()/hh[0].sum()))
Note: See TracBrowser for help on using the repository browser.