source: trunk/NIST_profile/atan_windowed_FP_profile.py @ 4413

Last change on this file since 4413 was 4413, checked in by toby, 20 months ago

add texture docs and fix build problem

  • Property svn:eol-style set to native
File size: 12.6 KB
Line 
1# -*- coding: utf-8 -*-
2"""
3Code to convert FP_Profile to use a windowed spectrum,
4and to generate real-space and Fourier IPF
5MHM October 2019
6"""
7
8from __future__ import print_function
9
10import numpy
11import math
12import sys
13
14from . import profile_functions_class
15
16class FP_atan_windowed_convolver:
17
18    ## @brief compute the emission spectrum with a window function
19    #  @param self self
20    #  @return the convolver for the emission
21    # requires an array of emiss_wavelengths, emiss_lor_widths, emiss_intensities,
22    # and a function emiss_passband_func which takes an argument in wavelength which
23    # will be offset by emiss_window_offset before being applied.
24    # @note this tramples the regular conv_emission, so this class must be to the left
25    #    of the main class in the final class creation e.g.:
26    #    class my_fp(FP_atan_windowed_convolver, profile_functions_class.FP_profile):
27    #        pass
28
29    def conv_emission(self):
30        """handle Lorentzian-sum emission spectrum with window function"""
31        me=self.get_function_name() #the name of this convolver,as a string
32        kwargs={}
33        kwargs.update(self.param_dicts[me]) #get all of our parameters
34        kwargs.update(self.param_dicts["conv_global"])
35 
36        #convert arrays to lists for key checking
37        key={}
38        key.update(kwargs)
39        for k,v in key.items():
40            if hasattr(v,'tolist') and hasattr(v,'__len__') and len(v) > 1:
41                key[k]=tuple(v.tolist())
42   
43        flag, emiss = self.get_conv(me, key, numpy.complex)
44        if flag: return emiss #already up to date
45       
46        xx=type("data", (), kwargs) #make it dot-notation accessible
47   
48        theta=xx.twotheta0/2
49        mono_width = (xx.b_mono/xx.diffractometer_radius) * math.tan(xx.two_theta_mono*math.pi/360)
50        dispersion=(2*math.tan(theta)+mono_width)
51        relative_passwidth = xx.ibm_source_width/(xx.a_mono*math.tan(xx.two_theta_mono*math.pi/360)) #FW delta-lambda/lambda eq. 6 in paper
52        passwidth=dispersion*relative_passwidth/2 #to match Topas macro notation
53        #and compute wavelength-scaled passband parameters from this
54        center_wavelength=xx.dominant_wavelength #fiducial wavelength from instrument global for mistune, etc.
55        center_angle=math.asin(center_wavelength/(2.0 * xx.d)) #nominal bragg angle
56
57        self.abs_window_pm=relative_passwidth*center_wavelength*1e12
58
59        X=self.twothetasamples-2*center_angle #to match notation in Topas macro
60        real_emiss=numpy.zeros_like(X) #emission in real space
61
62        for wid, lam0, intens in zip(xx.emiss_lor_widths, xx.emiss_wavelengths, xx.emiss_intensities):
63            #wid comes in as an FWHM, convert to HWHM
64            xvals=X-dispersion*(lam0-center_wavelength)/center_wavelength
65            hw=dispersion*(wid/center_wavelength)/2
66            real_emiss += (hw*intens/math.pi)/(xvals*xvals+hw*hw)
67
68        left_passband  = numpy.arctan((X+passwidth-xx.passband_mistune*passwidth)/(xx.passband_shoulder*passwidth))
69        right_passband = numpy.arctan((X-passwidth-xx.passband_mistune*passwidth)/(xx.passband_shoulder*passwidth))
70
71        self.raw_emission_shape=numpy.array(real_emiss)
72       
73        self.passband_window=(left_passband-right_passband)
74        real_emiss *= (left_passband-right_passband)
75        real_emiss /= real_emiss.sum() #normalize result
76       
77        self.full_emission_shape=real_emiss
78       
79        emiss[:]=profile_functions_class.best_rfft(real_emiss)
80       
81        return emiss
82
83    ## @brief compute the emission spectrum and (for convenience) the particle size widths
84    #  @param self self
85    #  @return the convolver for the emission and particle sizes
86    #  @note the particle size and strain stuff here is just to be consistent with \a Topas
87    #  and to be vaguely efficient about the computation, since all of these
88    #  have the same general shape.
89    def conv_crystallite_size(self):
90        """handle crystallite size broadening separately from windowed emission"""
91        me=self.get_function_name() #the name of this convolver,as a string
92        kwargs={}
93        kwargs.update(self.param_dicts["conv_emission"]) #steal our parameters from conv_emission
94        kwargs.update(self.param_dicts["conv_global"])
95        #if the crystallite size and strain parameters are not set, set them to
96        #values that make their corrections disappear
97        kwargs.setdefault("crystallite_size_lor",1e10)
98        kwargs.setdefault("crystallite_size_gauss",1e10)
99        kwargs.setdefault("strain_lor",0)
100        kwargs.setdefault("strain_gauss",0)
101
102        #convert arrays to lists for key checking
103        key={}
104        key.update(kwargs)
105        for k,v in key.items():
106            if hasattr(v,'tolist') and hasattr(v,'__len__') and len(v) > 1:
107                key[k]=tuple(v.tolist()) #tuples are better keys
108   
109        flag, emiss = self.get_conv(me, key, numpy.complex)
110        if flag: return emiss #already up to date
111       
112        xx=type("data", (), kwargs) #make it dot-notation accessible
113       
114        theta=xx.twotheta0/2
115        ## Emission profile FWHM + crystallite broadening (scale factors are Topas choice!) (Lorentzian)
116        #note:  the strain broadenings in Topas are expressed in degrees 2theta, must convert to radians(theta) with pi/360
117        widths = (
118            xx.strain_lor*(math.pi/360) * math.tan(theta) +
119            (xx.emiss_wavelengths / (2*xx.crystallite_size_lor * math.cos(theta)))
120        )
121        #save weighted average width for future reference in periodicity fixer
122        self.lor_widths[me]=sum(widths*xx.emiss_intensities)/sum(xx.emiss_intensities)
123        #gaussian bits add in quadrature
124        gfwhm2s=(
125            (xx.strain_gauss*(math.pi/360) * math.tan(theta))**2 +
126            (xx.emiss_wavelengths / (xx.crystallite_size_gauss * math.cos(theta)))**2
127        )
128       
129        # note that the Fourier transform of a lorentzian with FWHM 2a
130        # is exp(-abs(a omega))
131        omega_vals=self.omega_vals
132        for wid, gfwhm2, intens in zip(widths, gfwhm2s, xx.emiss_intensities):
133            xvals=numpy.clip(omega_vals*(-wid),-100,0)
134            sig2=gfwhm2/(8*math.log(2.0)) #convert fwhm**2 to sigma**2
135            gxv=numpy.clip((sig2/-2.0)*omega_vals*omega_vals,-100,0)
136            emiss += numpy.exp(xvals+gxv)*intens
137        return emiss
138
139
140#cu_ka_spectdata=numpy.array(( #each cluster is wavelength/m, intensity, fwhm/m, from Cu kalpha paper
141#        (0.15405925, 3.91, 0.0436e-3), #ka11
142#        (0.15410769, 0.474, 0.0558e-3), #ka12
143#        (0.15443873, 1.53, 0.0487e-3), #ka21
144#        (0.15446782, 0.754, 0.0630e-3), #ka22
145#    ))*(1e-9,1,1e-9)
146# replaced due to Sphinx problem with scaled values:
147cu_ka_spectdata=numpy.array(( #each cluster is wavelength/m, intensity, fwhm/m, from Cu kalpha paper
148        (0.15405925e-9, 3.91, 0.0436e-12), #ka11
149        (0.15410769e-9, 0.474, 0.0558e-12), #ka12
150        (0.15443873e-9, 1.53, 0.0487e-12), #ka21
151        (0.15446782e-9, 0.754, 0.0630e-12), #ka22
152    ))
153
154if __name__ == "__main__":
155    ## fixed parameters
156
157    class FP_windowed(FP_atan_windowed_convolver, profile_functions_class.FP_profile):
158        pass
159
160    def fixphase(xfrm, second_try_bin=None):
161        """adjust phases to a transform is close to centered. Warning: done in place!"""
162        phase=xfrm[1].conjugate()/abs(xfrm[1])
163        phase0=phase
164        for i in range(1,len(xfrm)):
165            xfrm[i]*=phase
166            phase *= phase0
167        if second_try_bin:
168            zz=xfrm[second_try_bin]
169            theta=-math.atan2(zz.imag, zz.real)/second_try_bin #better average phase shift
170            phase=complex(math.cos(theta), math.sin(theta))
171            phase0=phase
172            for i in range(1,len(xfrm)):
173                xfrm[i]*=phase
174                phase *= phase0
175   
176        return xfrm
177   
178    from diffraction_constants import omega_length_scale_deg
179   
180    from matplotlib import pyplot as plt
181    plt.style.use('posters_and_pubs')
182   
183    emiss_wavelengths, emiss_intensities, emiss_lor_widths = cu_ka_spectdata.transpose()
184   
185    p=FP_windowed(anglemode="twotheta",
186            output_gaussian_smoother_bins_sigma=1.0,
187            oversampling=2
188        )
189
190    #parameters from dbd_ibm_setup_atan_20181116.inc
191    #    prm !pixedge  0.08562_0.00399 'in mm
192    #    hat  = Rad pixedge/Rs; :  0.0168`_0.00016 num_hats 1 'rectangular pixels
193    #
194    #    #ifndef REFINE_SOURCE_WIDTH
195    #    prm !ibm_source_width  0.04828_0.00032 min 0.0002 max 1
196    #    prm !passband_shoulder  0.25772_0.00982 min 0.01 max 1.5  'width ratio of Lorentzian component to flattop
197    #    #endif
198    #
199    #    prm !two_theta_mono 27.27 'degrees twotheta, of the Ge crystal for a peak at CuKa1
200    #    prm !b_mono 217.5 'happens to be the same as the DBD radius
201    #
202   
203    #set parameters for each convolver
204    p.set_parameters(convolver="emission",
205        emiss_wavelengths=emiss_wavelengths,
206        emiss_intensities=emiss_intensities,
207        emiss_lor_widths=emiss_lor_widths,
208        emiss_gauss_widths=[0,0,0,0],
209        ibm_source_width=0.05102e-3, #51 um
210        passband_shoulder=0.08711,
211        two_theta_mono=27.27,
212        a_mono=119.e-3,
213        b_mono=217.5e-3,
214        passband_mistune=-0.14496,
215    )
216   
217    p.set_parameters(convolver="axial",
218        axDiv="full", slit_length_source=8e-3,
219        slit_length_target=12e-3,
220        length_sample=15e-3,
221        n_integral_points=20,
222        angI_deg=2.73, angD_deg=2.73,
223    )
224   
225    p.set_parameters(convolver="absorption",
226        absorption_coefficient=50.0*100*10000, #like SRM 1979, in m^(-1)
227    )
228
229    p.set_parameters(convolver="receiver_slit",
230     slit_width=0.11562e-3, #convert mm to m
231    )
232
233    import time
234    t_start=time.time()
235   
236    profiles=[]
237    for twotheta_x in (31.769, 66.37, 135,):
238#    for dstar in (4., 7., 11.):
239#        twotheta_x=360/math.pi*math.asin(0.15409*dstar/2)
240        #put the compute window in the right place and clear all histories
241        if twotheta_x > 90: window=2.
242        else: window=1.   
243        p.set_window(twotheta_output_points=2000,
244            twotheta_window_center_deg=twotheta_x,
245            twotheta_window_fullwidth_deg=window,
246        )
247        #set parameters which are shared by many things
248        p.set_parameters(twotheta0_deg=twotheta_x,
249            dominant_wavelength=emiss_wavelengths[0],
250            diffractometer_radius=217.5e-3)
251   
252        p.set_parameters(equatorial_divergence_deg=1.05)
253       
254        result=p.compute_line_profile(return_convolver=True)
255        print(p,file=sys.stderr)
256        abs_window_pm=p.abs_window_pm
257       
258        plt.figure("real {0:.0f}".format(twotheta_x))
259        pk=result.peak/result.peak.max()
260        emiss=p.full_emission_shape[::2]/p.full_emission_shape.max()
261        raw_emiss=p.raw_emission_shape[::2]/p.raw_emission_shape.max()
262        passband=p.passband_window[::2]/p.passband_window.max()
263        plt.plot(result.twotheta_deg, pk)
264        plt.plot(result.twotheta_deg, emiss)
265        plt.plot(result.twotheta_deg, raw_emiss)
266        plt.plot(result.twotheta_deg, passband)
267        plt.xlabel(r"$2\theta$ / degree")
268        plt.savefig("realspace_window_{abs_window_pm:.2f}_tth_{twotheta_x:.0f}.pdf".format(**locals()),
269                    bbox_inches='tight', transparent=True, frameon=False)
270        #plt.legend()
271        #plt.show()
272        dstarplot=numpy.array((result.twotheta_deg, 2*numpy.sin(result.twotheta/2)/0.15409, pk)).transpose()
273        numpy.savetxt("dstar_{abs_window_pm:.2f}_tth_{twotheta_x:.0f}.dat".format(**locals()),
274            dstarplot)
275       
276        omega=result.omega_inv_deg*omega_length_scale_deg(0.15409, twotheta_x)
277        omega_max=numpy.searchsorted(omega, 300)
278        fft=numpy.array(result.convolver)
279        fft /= fft[0]
280        fft[1::2]*=-1 #center it
281        fixphase(fft, second_try_bin=5) #adjust to best centering
282        numpy.savetxt("fourier_{abs_window_pm:.2f}_tth_{twotheta_x:.0f}.dat".format(**locals()),
283            numpy.transpose((omega[:omega_max], numpy.abs(fft[:omega_max])))
284        )
285        print(fft[:10], file=sys.stderr)
286        plt.figure("fourier log {0:.0f}".format(twotheta_x))
287        plt.semilogy(omega[:omega_max], numpy.abs(fft[:omega_max]))     
288        plt.xlabel(r"$\omega$ / nm")
289        plt.savefig("ft_log_window_{abs_window_pm:.2f}_tth_{twotheta_x:.0f}.pdf".format(**locals()),
290                    bbox_inches='tight', transparent=True, frameon=False)
291        plt.figure("fourier lin {0:.0f}".format(twotheta_x))
292        plt.plot(omega[:omega_max], numpy.abs(fft[:omega_max]))       
293        plt.xlabel(r"$\omega$ / nm")
294        plt.savefig("ft_lin_window_{abs_window_pm:.2f}_tth_{twotheta_x:.0f}.pdf".format(**locals()),
295                    bbox_inches='tight', transparent=True, frameon=False)
296        #plt.legend()
297        #plt.show()
Note: See TracBrowser for help on using the repository browser.