source: sasview/src/sas/sascalc/calculator/resolution_calculator.py @ 574adc7

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalcmagnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 574adc7 was 574adc7, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

convert sascalc to python 2/3 syntax

  • Property mode set to 100644
File size: 38.5 KB
RevLine 
[51f14603]1"""
2This object is a small tool to allow user to quickly
3determine the variance in q  from the
4instrumental parameters.
5"""
[b6627d9]6import sys
[574adc7]7from math import pi, sqrt
8import math
[b6627d9]9import logging
[51f14603]10
[574adc7]11import numpy as np
12
13from .instrument import Sample
14from .instrument import Detector
15from .instrument import TOF as Neutron
16from .instrument import Aperture
17
[463e7ffc]18logger = logging.getLogger(__name__)
[c155a16]19
[51f14603]20#Plank's constant in cgs unit
21_PLANK_H = 6.62606896E-27
22#Gravitational acc. in cgs unit
23_GRAVITY = 981.0
24
25
26class ResolutionCalculator(object):
27    """
28    compute resolution in 2D
29    """
30    def __init__(self):
[6bd3a8d1]31
[51f14603]32        # wavelength
33        self.wave = Neutron()
34        # sample
35        self.sample = Sample()
36        # aperture
37        self.aperture = Aperture()
38        # detector
39        self.detector = Detector()
40        # 2d image of the resolution
41        self.image = []
42        self.image_lam = []
43        # resolutions
44        # lamda in r-direction
45        self.sigma_lamd = 0
46        # x-dir (no lamda)
47        self.sigma_1 = 0
48        #y-dir (no lamda)
49        self.sigma_2 = 0
50        # 1D total
51        self.sigma_1d = 0
52        self.gravity_phi = None
53        # q min and max
54        self.qx_min = -0.3
55        self.qx_max = 0.3
56        self.qy_min = -0.3
57        self.qy_max = 0.3
58        # q min and max of the detector
59        self.detector_qx_min = -0.3
60        self.detector_qx_max = 0.3
61        self.detector_qy_min = -0.3
62        self.detector_qy_max = 0.3
63        # possible max qrange
64        self.qxmin_limit = 0
65        self.qxmax_limit = 0
66        self.qymin_limit = 0
67        self.qymax_limit = 0
68
69        # plots
70        self.plot = None
71        # instrumental params defaults
72        self.mass = 0
73        self.intensity = 0
74        self.wavelength = 0
75        self.wavelength_spread = 0
76        self.source_aperture_size = []
77        self.source2sample_distance = []
78        self.sample2sample_distance = []
79        self.sample_aperture_size = []
80        self.sample2detector_distance = []
81        self.detector_pix_size = []
82        self.detector_size = []
83        self.get_all_instrument_params()
84        # max q range for all lambdas
85        self.qxrange = []
86        self.qyrange = []
[6bd3a8d1]87
[51f14603]88    def compute_and_plot(self, qx_value, qy_value, qx_min, qx_max,
[b6627d9]89                         qy_min, qy_max, coord='cartesian'):
[51f14603]90        """
91        Compute the resolution
92        : qx_value: x component of q
93        : qy_value: y component of q
94        """
95        # make sure to update all the variables need.
96        # except lambda, dlambda, and intensity
97        self.get_all_instrument_params()
98        # wavelength etc.
99        lamda_list, dlamb_list = self.get_wave_list()
100        intens_list = []
101        sig1_list = []
102        sig2_list = []
103        sigr_list = []
104        sigma1d_list = []
105        num_lamda = len(lamda_list)
106        for num in range(num_lamda):
107            lam = lamda_list[num]
108            # wavelength spread
109            dlam = dlamb_list[num]
110            intens = self.setup_tof(lam, dlam)
111            intens_list.append(intens)
112            # cehck if tof
113            if num_lamda > 1:
114                tof = True
115            else:
116                tof = False
117            # compute 2d resolution
118            _, _, sigma_1, sigma_2, sigma_r, sigma1d = \
119                        self.compute(lam, dlam, qx_value, qy_value, coord, tof)
120            # make image
121            image = self.get_image(qx_value, qy_value, sigma_1, sigma_2,
[b6627d9]122                                   sigma_r, qx_min, qx_max, qy_min, qy_max,
123                                   coord, False)
[51f14603]124            if qx_min > self.qx_min:
125                qx_min = self.qx_min
126            if qx_max < self.qx_max:
127                qx_max = self.qx_max
128            if qy_min > self.qy_min:
129                qy_min = self.qy_min
130            if qy_max < self.qy_max:
131                qy_max = self.qy_max
[6bd3a8d1]132
[51f14603]133            # set max qranges
134            self.qxrange = [qx_min, qx_max]
135            self.qyrange = [qy_min, qy_max]
136            sig1_list.append(sigma_1)
137            sig2_list.append(sigma_2)
138            sigr_list.append(sigma_r)
139            sigma1d_list.append(sigma1d)
140        # redraw image in global 2d q-space.
141        self.image_lam = []
142        total_intensity = 0
143        sigma_1 = 0
144        sigma_r = 0
145        sigma_2 = 0
146        sigma1d = 0
147        for ind in range(num_lamda):
148            lam = lamda_list[ind]
149            dlam = dlamb_list[ind]
150            intens = self.setup_tof(lam, dlam)
151            out = self.get_image(qx_value, qy_value, sig1_list[ind],
[b6627d9]152                                 sig2_list[ind], sigr_list[ind],
153                                 qx_min, qx_max, qy_min, qy_max, coord)
[51f14603]154            # this is the case of q being outside the detector
155            #if numpy.all(out==0.0):
156            #    continue
157            image = out
158            # set variance as sigmas
159            sigma_1 += sig1_list[ind] * sig1_list[ind] * self.intensity
160            sigma_r += sigr_list[ind] * sigr_list[ind] * self.intensity
161            sigma_2 += sig2_list[ind] * sig2_list[ind] * self.intensity
162            sigma1d += sigma1d_list[ind] * sigma1d_list[ind] * self.intensity
163            total_intensity += self.intensity
[6bd3a8d1]164
[51f14603]165        if total_intensity != 0:
166            # average variance
167            image_out = image / total_intensity
168            sigma_1 = sigma_1 / total_intensity
169            sigma_r = sigma_r / total_intensity
170            sigma_2 = sigma_2 / total_intensity
171            sigma1d = sigma1d / total_intensity
172            # set sigmas
173            self.sigma_1 = sqrt(sigma_1)
174            self.sigma_lamd = sqrt(sigma_r)
175            self.sigma_2 = sqrt(sigma_2)
176            self.sigma_1d = sqrt(sigma1d)
177            # rescale
178            max_im_val = 1
179            if max_im_val > 0:
180                image_out /= max_im_val
181        else:
182            image_out = image * 0.0
183            # Don't calculate sigmas nor set self.sigmas!
184            sigma_1 = 0
185            sigma_r = 0
186            sigma_2 = 0
187            sigma1d = 0
188        if len(self.image) > 0:
189            self.image += image_out
190        else:
191            self.image = image_out
[6bd3a8d1]192
[51f14603]193        # plot image
194        return self.plot_image(self.image)
[6bd3a8d1]195
[51f14603]196    def setup_tof(self, wavelength, wavelength_spread):
197        """
198        Setup all parameters in instrument
[6bd3a8d1]199
[51f14603]200        : param ind: index of lambda, etc
201        """
202
203        # set wave.wavelength
204        self.set_wavelength(wavelength)
205        self.set_wavelength_spread(wavelength_spread)
206        self.intensity = self.wave.get_intensity()
[6bd3a8d1]207
[51f14603]208        if wavelength == 0:
209            msg = "Can't compute the resolution: the wavelength is zero..."
[574adc7]210            raise RuntimeError(msg)
[51f14603]211        return self.intensity
[6bd3a8d1]212
[51f14603]213    def compute(self, wavelength, wavelength_spread, qx_value, qy_value,
214                coord='cartesian', tof=False):
215        """
216        Compute the Q resoltuion in || and + direction of 2D
217        : qx_value: x component of q
218        : qy_value: y component of q
219        """
220        coord = 'cartesian'
221        lamb = wavelength
222        lamb_spread = wavelength_spread
223        # the shape of wavelength distribution
224
225        if tof:
226            # rectangular
227            tof_factor = 2
228        else:
229            # triangular
230            tof_factor = 1
231        # Find polar values
232        qr_value, phi = self._get_polar_value(qx_value, qy_value)
233        # vacuum wave transfer
234        knot = 2*pi/lamb
235        # scattering angle theta; always true for plane detector
236        # aligned vertically to the ko direction
237        if qr_value > knot:
238            theta = pi/2
239        else:
240            theta = math.asin(qr_value/knot)
241        # source aperture size
242        rone = self.source_aperture_size
243        # sample aperture size
244        rtwo = self.sample_aperture_size
245        # detector pixel size
246        rthree = self.detector_pix_size
247        # source to sample(aperture) distance
248        l_ssa = self.source2sample_distance[0]
249        # sample(aperture) to detector distance
250        l_sad = self.sample2detector_distance[0]
251        # sample (aperture) to sample distance
252        l_sas = self.sample2sample_distance[0]
253        # source to sample distance
254        l_one = l_ssa + l_sas
255        # sample to detector distance
256        l_two = l_sad - l_sas
[6bd3a8d1]257
[51f14603]258        # Sample offset correction for l_one and Lp on variance calculation
259        l1_cor = (l_ssa * l_two) / (l_sas + l_two)
260        lp_cor = (l_ssa * l_two) / (l_one + l_two)
261        # the radial distance to the pixel from the center of the detector
262        radius = math.tan(theta) * l_two
263        #Lp = l_one*l_two/(l_one+l_two)
264        # default polar coordinate
265        comp1 = 'radial'
266        comp2 = 'phi'
267        # in the case of the cartesian coordinate
268        if coord == 'cartesian':
269            comp1 = 'x'
270            comp2 = 'y'
271
272        # sigma in the radial/x direction
273        # for source aperture
[b6627d9]274        sigma_1 = self.get_variance(rone, l1_cor, phi, comp1)
[51f14603]275        # for sample apperture
276        sigma_1 += self.get_variance(rtwo, lp_cor, phi, comp1)
277        # for detector pix
278        sigma_1 += self.get_variance(rthree, l_two, phi, comp1)
279        # for gravity term for 1d
280        sigma_1grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb,
[b6627d9]281                                                  lamb_spread, phi, comp1, 'on') / tof_factor
[51f14603]282        # for wavelength spread
283        # reserve for 1d calculation
284        A_value = self._cal_A_value(lamb, l_ssa, l_sad)
285        sigma_wave_1, sigma_wave_1_1d = self.get_variance_wave(A_value,
[b6627d9]286                                                               radius, l_two, lamb_spread,
287                                                               phi, 'radial', 'on')
[51f14603]288        sigma_wave_1 /= tof_factor
289        sigma_wave_1_1d /= tof_factor
290        # for 1d
291        variance_1d_1 = (sigma_1 + sigma_1grav1d) / 2 + sigma_wave_1_1d
292        # normalize
293        variance_1d_1 = knot * knot * variance_1d_1 / 12
[6bd3a8d1]294
[51f14603]295        # for 2d
296        #sigma_1 += sigma_wave_1
297        # normalize
298        sigma_1 = knot * sqrt(sigma_1 / 12)
299        sigma_r = knot * sqrt(sigma_wave_1 / (tof_factor *12))
300        # sigma in the phi/y direction
301        # for source apperture
[b6627d9]302        sigma_2 = self.get_variance(rone, l1_cor, phi, comp2)
[51f14603]303
304        # for sample apperture
305        sigma_2 += self.get_variance(rtwo, lp_cor, phi, comp2)
306
307        # for detector pix
308        sigma_2 += self.get_variance(rthree, l_two, phi, comp2)
309
310        # for gravity term for 1d
311        sigma_2grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb,
[b6627d9]312                                                  lamb_spread, phi, comp2, 'on') / tof_factor
[51f14603]313
314        # for wavelength spread
315        # reserve for 1d calculation
316        sigma_wave_2, sigma_wave_2_1d = self.get_variance_wave(A_value,
[b6627d9]317                                                               radius, l_two, lamb_spread,
318                                                               phi, 'phi', 'on')
[51f14603]319        sigma_wave_2 /= tof_factor
320        sigma_wave_2_1d /= tof_factor
321        # for 1d
322        variance_1d_2 = (sigma_2 + sigma_2grav1d) / 2 + sigma_wave_2_1d
323        # normalize
324        variance_1d_2 = knot * knot * variance_1d_2 / 12
[6bd3a8d1]325
[51f14603]326        # for 2d
327        #sigma_2 =  knot*sqrt(sigma_2/12)
328        #sigma_2 += sigma_wave_2
329        # normalize
330        sigma_2 = knot * sqrt(sigma_2 / 12)
331        sigma1d = sqrt(variance_1d_1 + variance_1d_2)
332        # set sigmas
333        self.sigma_1 = sigma_1
334        self.sigma_lamd = sigma_r
335        self.sigma_2 = sigma_2
336        self.sigma_1d = sigma1d
337        return qr_value, phi, sigma_1, sigma_2, sigma_r, sigma1d
[6bd3a8d1]338
[51f14603]339    def _within_detector_range(self, qx_value, qy_value):
340        """
341        check if qvalues are within detector range
342        """
343        # detector range
344        detector_qx_min = self.detector_qx_min
345        detector_qx_max = self.detector_qx_max
346        detector_qy_min = self.detector_qy_min
347        detector_qy_max = self.detector_qy_max
348        if self.qxmin_limit > detector_qx_min:
349            self.qxmin_limit = detector_qx_min
350        if self.qxmax_limit < detector_qx_max:
351            self.qxmax_limit = detector_qx_max
352        if self.qymin_limit > detector_qy_min:
353            self.qymin_limit = detector_qy_min
354        if self.qymax_limit < detector_qy_max:
355            self.qymax_limit = detector_qy_max
356        if qx_value < detector_qx_min or qx_value > detector_qx_max:
357            return False
358        if qy_value < detector_qy_min or qy_value > detector_qy_max:
359            return False
360        return True
[6bd3a8d1]361
[51f14603]362    def get_image(self, qx_value, qy_value, sigma_1, sigma_2, sigma_r,
363                  qx_min, qx_max, qy_min, qy_max,
364                  coord='cartesian', full_cal=True):
365        """
366        Get the resolution in polar coordinate ready to plot
367        : qx_value: qx_value value
368        : qy_value: qy_value value
369        : sigma_1: variance in r direction
370        : sigma_2: variance in phi direction
371        : coord: coordinate system of image, 'polar' or 'cartesian'
372        """
373        # Get  qx_max and qy_max...
374        self._get_detector_qxqy_pixels()
[6bd3a8d1]375
[51f14603]376        qr_value, phi = self._get_polar_value(qx_value, qy_value)
377
378        # Check whether the q value is within the detector range
379        if qx_min < self.qx_min:
380            self.qx_min = qx_min
[574adc7]381            #raise ValueError(msg)
[51f14603]382        if qx_max > self.qx_max:
383            self.qx_max = qx_max
[574adc7]384            #raise ValueError(msg)
[51f14603]385        if qy_min < self.qy_min:
386            self.qy_min = qy_min
[574adc7]387            #raise ValueError(msg)
[51f14603]388        if qy_max > self.qy_max:
389            self.qy_max = qy_max
[574adc7]390            #raise ValueError(msg)
[51f14603]391        if not full_cal:
392            return None
[6bd3a8d1]393
[51f14603]394        # Make an empty graph in the detector scale
395        dx_size = (self.qx_max - self.qx_min) / (1000 - 1)
396        dy_size = (self.qy_max - self.qy_min) / (1000 - 1)
[9a5097c]397        x_val = np.arange(self.qx_min, self.qx_max, dx_size)
398        y_val = np.arange(self.qy_max, self.qy_min, -dy_size)
399        q_1, q_2 = np.meshgrid(x_val, y_val)
[51f14603]400        #q_phi = numpy.arctan(q_1,q_2)
401        # check whether polar or cartesian
402        if coord == 'polar':
403            # Find polar values
404            qr_value, phi = self._get_polar_value(qx_value, qy_value)
405            q_1, q_2 = self._rotate_z(q_1, q_2, phi)
406            qc_1 = qr_value
407            qc_2 = 0.0
408            # Calculate the 2D Gaussian distribution image
409            image = self._gaussian2d_polar(q_1, q_2, qc_1, qc_2,
[b6627d9]410                                           sigma_1, sigma_2, sigma_r)
[51f14603]411        else:
412            # catesian coordinate
413            # qx_center
414            qc_1 = qx_value
415            # qy_center
416            qc_2 = qy_value
[6bd3a8d1]417
[51f14603]418            # Calculate the 2D Gaussian distribution image
419            image = self._gaussian2d(q_1, q_2, qc_1, qc_2,
420                                     sigma_1, sigma_2, sigma_r)
421        # out side of detector
422        if not self._within_detector_range(qx_value, qy_value):
423            image *= 0.0
424            self.intensity = 0.0
425            #return self.image
426
427        # Add it if there are more than one inputs.
428        if len(self.image_lam) > 0:
429            self.image_lam += image * self.intensity
430        else:
431            self.image_lam = image * self.intensity
[6bd3a8d1]432
[51f14603]433        return self.image_lam
[6bd3a8d1]434
[51f14603]435    def plot_image(self, image):
436        """
437        Plot image using pyplot
438        : image: 2d resolution image
[6bd3a8d1]439
[51f14603]440        : return plt: pylab object
441        """
442        import matplotlib.pyplot as plt
443
444        self.plot = plt
445        plt.xlabel('$\\rm{Q}_{x} [A^{-1}]$')
446        plt.ylabel('$\\rm{Q}_{y} [A^{-1}]$')
447        # Max value of the image
448        # max = numpy.max(image)
449        qx_min, qx_max, qy_min, qy_max = self.get_detector_qrange()
450
451        # Image
452        im = plt.imshow(image,
[b6627d9]453                        extent=[qx_min, qx_max, qy_min, qy_max])
[51f14603]454
455        # bilinear interpolation to make it smoother
456        im.set_interpolation('bilinear')
457
458        return plt
[6bd3a8d1]459
[51f14603]460    def reset_image(self):
461        """
462        Reset image to default (=[])
463        """
464        self.image = []
[6bd3a8d1]465
[51f14603]466    def get_variance(self, size=[], distance=0, phi=0, comp='radial'):
467        """
468        Get the variance when the slit/pinhole size is given
469        : size: list that can be one(diameter for circular) or two components(lengths for rectangular)
470        : distance: [z, x] where z along the incident beam, x // qx_value
471        : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial'
[6bd3a8d1]472
[51f14603]473        : return variance: sigma^2
474        """
475        # check the length of size (list)
476        len_size = len(size)
[6bd3a8d1]477
[51f14603]478        # define sigma component direction
479        if comp == 'radial':
480            phi_x = math.cos(phi)
481            phi_y = math.sin(phi)
482        elif comp == 'phi':
483            phi_x = math.sin(phi)
484            phi_y = math.cos(phi)
485        elif comp == 'x':
486            phi_x = 1
487            phi_y = 0
488        elif comp == 'y':
489            phi_x = 0
490            phi_y = 1
491        else:
492            phi_x = 0
493            phi_y = 0
494        # calculate each component
495        # for pinhole w/ radius = size[0]/2
496        if len_size == 1:
497            x_comp = (0.5 * size[0]) * sqrt(3)
498            y_comp = 0
499        # for rectangular slit
500        elif len_size == 2:
501            x_comp = size[0] * phi_x
502            y_comp = size[1] * phi_y
503        # otherwise
504        else:
[574adc7]505            raise ValueError(" Improper input...")
[51f14603]506        # get them squared
[b6627d9]507        sigma = x_comp * x_comp
[51f14603]508        sigma += y_comp * y_comp
509        # normalize by distance
510        sigma /= (distance * distance)
511
512        return sigma
513
514    def get_variance_wave(self, A_value, radius, distance, spread, phi,
515                          comp='radial', switch='on'):
516        """
517        Get the variance when the wavelength spread is given
[6bd3a8d1]518
[51f14603]519        : radius: the radial distance from the beam center to the pix of q
520        : distance: sample to detector distance
521        : spread: wavelength spread (ratio)
522        : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial'
[6bd3a8d1]523
[51f14603]524        : return variance: sigma^2 for 2d, sigma^2 for 1d [tuple]
525        """
526        if switch.lower() == 'off':
527            return 0, 0
528        # check the singular point
529        if distance == 0 or comp == 'phi':
530            return 0, 0
531        else:
532            # calculate sigma^2 for 1d
533            sigma1d = 2 * math.pow(radius/distance*spread, 2)
534            if comp == 'x':
535                sigma1d *= (math.cos(phi)*math.cos(phi))
536            elif comp == 'y':
537                sigma1d *= (math.sin(phi)*math.sin(phi))
538            else:
539                sigma1d *= 1
540            # sigma^2 for 2d
541            # shift the coordinate due to the gravitational shift
542            rad_x = radius * math.cos(phi)
543            rad_y = A_value - radius * math.sin(phi)
544            radius = math.sqrt(rad_x * rad_x + rad_y * rad_y)
545            # new phi
546            phi = math.atan2(-rad_y, rad_x)
547            self.gravity_phi = phi
548            # calculate sigma^2
549            sigma = 2 * math.pow(radius/distance*spread, 2)
550            if comp == 'x':
551                sigma *= (math.cos(phi)*math.cos(phi))
552            elif comp == 'y':
553                sigma *= (math.sin(phi)*math.sin(phi))
554            else:
555                sigma *= 1
[6bd3a8d1]556
[51f14603]557            return sigma, sigma1d
558
559    def get_variance_gravity(self, s_distance, d_distance, wavelength, spread,
560                             phi, comp='radial', switch='on'):
561        """
562        Get the variance from gravity when the wavelength spread is given
[6bd3a8d1]563
[51f14603]564        : s_distance: source to sample distance
565        : d_distance: sample to detector distance
566        : wavelength: wavelength
567        : spread: wavelength spread (ratio)
568        : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial'
[6bd3a8d1]569
[51f14603]570        : return variance: sigma^2
571        """
572        if switch.lower() == 'off':
573            return 0
574        if self.mass == 0.0:
575            return 0
576        # check the singular point
577        if d_distance == 0 or comp == 'x':
578            return 0
579        else:
580            a_value = self._cal_A_value(None, s_distance, d_distance)
581            # calculate sigma^2
582            sigma = math.pow(a_value / d_distance, 2)
583            sigma *= math.pow(wavelength, 4)
584            sigma *= math.pow(spread, 2)
585            sigma *= 8
586            return sigma
[6bd3a8d1]587
[51f14603]588    def _cal_A_value(self, lamda, s_distance, d_distance):
589        """
590        Calculate A value for gravity
[6bd3a8d1]591
[51f14603]592        : s_distance: source to sample distance
593        : d_distance: sample to detector distance
594        """
595        # neutron mass in cgs unit
596        self.mass = self.get_neutron_mass()
597        # plank constant in cgs unit
598        h_constant = _PLANK_H
599        # gravity in cgs unit
600        gravy = _GRAVITY
601        # m/h
602        m_over_h = self.mass / h_constant
603        # A value
604        a_value = d_distance * (s_distance + d_distance)
605        a_value *= math.pow(m_over_h / 2, 2)
606        a_value *= gravy
607        # unit correction (1/cm to 1/A) for A and d_distance below
608        a_value *= 1.0E-16
609        # if lamda is give (broad meanning of A)  return 2* lamda^2 * A
[7432acb]610        if lamda is not None:
[51f14603]611            a_value *= (4 * lamda * lamda)
612        return a_value
[6bd3a8d1]613
[51f14603]614    def get_intensity(self):
615        """
616        Get intensity
617        """
618        return self.wave.intensity
619
620    def get_wavelength(self):
621        """
622        Get wavelength
623        """
624        return self.wave.wavelength
[6bd3a8d1]625
[51f14603]626    def get_default_spectrum(self):
627        """
628        Get default_spectrum
629        """
630        return self.wave.get_default_spectrum()
[6bd3a8d1]631
[51f14603]632    def get_spectrum(self):
633        """
634        Get _spectrum
635        """
636        return self.wave.get_spectrum()
[6bd3a8d1]637
[51f14603]638    def get_wavelength_spread(self):
639        """
640        Get wavelength spread
641        """
642        return self.wave.wavelength_spread
[6bd3a8d1]643
[51f14603]644    def get_neutron_mass(self):
645        """
646        Get Neutron mass
647        """
648        return self.wave.mass
[6bd3a8d1]649
[51f14603]650    def get_source_aperture_size(self):
651        """
652        Get source aperture size
653        """
654        return self.aperture.source_size
[6bd3a8d1]655
[51f14603]656    def get_sample_aperture_size(self):
657        """
658        Get sample aperture size
659        """
660        return self.aperture.sample_size
[6bd3a8d1]661
[51f14603]662    def get_detector_pix_size(self):
663        """
664        Get detector pixel size
665        """
666        return self.detector.pix_size
[6bd3a8d1]667
[51f14603]668    def get_detector_size(self):
669        """
670        Get detector size
671        """
672        return self.detector.size
[6bd3a8d1]673
[51f14603]674    def get_source2sample_distance(self):
675        """
676        Get detector source2sample_distance
677        """
678        return self.aperture.sample_distance
[6bd3a8d1]679
[51f14603]680    def get_sample2sample_distance(self):
681        """
682        Get detector sampleslitsample_distance
683        """
684        return self.sample.distance
[6bd3a8d1]685
[51f14603]686    def get_sample2detector_distance(self):
687        """
688        Get detector sample2detector_distance
689        """
690        return self.detector.distance
[6bd3a8d1]691
[51f14603]692    def set_intensity(self, intensity):
693        """
694        Set intensity
695        """
696        self.wave.set_intensity(intensity)
[6bd3a8d1]697
[51f14603]698    def set_wave(self, wavelength):
699        """
700        Set wavelength list or wavelength
701        """
702        if wavelength.__class__.__name__ == 'list':
703            self.wave.set_wave_list(wavelength)
704        elif wavelength.__class__.__name__ == 'float':
705            self.wave.set_wave_list([wavelength])
706            #self.set_wavelength(wavelength)
707        else:
[574adc7]708            raise TypeError("invalid wavlength---should be list or float")
[6bd3a8d1]709
[51f14603]710    def set_wave_spread(self, wavelength_spread):
711        """
712        Set wavelength spread  or wavelength spread
713        """
714        if wavelength_spread.__class__.__name__ == 'list':
715            self.wave.set_wave_spread_list(wavelength_spread)
716        elif wavelength_spread.__class__.__name__ == 'float':
717            self.wave.set_wave_spread_list([wavelength_spread])
718        else:
[574adc7]719            raise TypeError("invalid wavelength spread---should be list or float")
[6bd3a8d1]720
[51f14603]721    def set_wavelength(self, wavelength):
722        """
723        Set wavelength
724        """
725        self.wavelength = wavelength
726        self.wave.set_wavelength(wavelength)
[6bd3a8d1]727
[51f14603]728    def set_spectrum(self, spectrum):
729        """
730        Set spectrum
731        """
732        self.spectrum = spectrum
733        self.wave.set_spectrum(spectrum)
[6bd3a8d1]734
[51f14603]735    def set_wavelength_spread(self, wavelength_spread):
736        """
737        Set wavelength spread
738        """
739        self.wavelength_spread = wavelength_spread
740        self.wave.set_wavelength_spread(wavelength_spread)
[6bd3a8d1]741
[51f14603]742    def set_wave_list(self, wavelength_list, wavelengthspread_list):
743        """
744        Set wavelength and its spread list
745        """
746        self.wave.set_wave_list(wavelength_list)
747        self.wave.set_wave_spread_list(wavelengthspread_list)
[6bd3a8d1]748
[51f14603]749    def get_wave_list(self):
750        """
751        Set wavelength spread
752        """
753        return self.wave.get_wave_list()
[6bd3a8d1]754
[51f14603]755    def get_intensity_list(self):
756        """
757        Set wavelength spread
758        """
759        return self.wave.get_intensity_list()
[6bd3a8d1]760
[51f14603]761    def set_source_aperture_size(self, size):
762        """
763        Set source aperture size
[6bd3a8d1]764
[51f14603]765        : param size: [dia_value] or [x_value, y_value]
766        """
767        if len(size) < 1 or len(size) > 2:
[574adc7]768            raise RuntimeError("The length of the size must be one or two.")
[51f14603]769        self.aperture.set_source_size(size)
[6bd3a8d1]770
[51f14603]771    def set_neutron_mass(self, mass):
772        """
773        Set Neutron mass
774        """
775        self.wave.set_mass(mass)
776        self.mass = mass
[6bd3a8d1]777
[51f14603]778    def set_sample_aperture_size(self, size):
779        """
780        Set sample aperture size
[6bd3a8d1]781
[51f14603]782        : param size: [dia_value] or [xheight_value, yheight_value]
783        """
784        if len(size) < 1 or len(size) > 2:
[574adc7]785            raise RuntimeError("The length of the size must be one or two.")
[51f14603]786        self.aperture.set_sample_size(size)
[6bd3a8d1]787
[51f14603]788    def set_detector_pix_size(self, size):
789        """
790        Set detector pixel size
791        """
792        self.detector.set_pix_size(size)
[6bd3a8d1]793
[51f14603]794    def set_detector_size(self, size):
795        """
796        Set detector size in number of pixels
797        : param size: [pixel_nums] or [x_pix_num, yx_pix_num]
798        """
799        self.detector.set_size(size)
[6bd3a8d1]800
[51f14603]801    def set_source2sample_distance(self, distance):
802        """
803        Set detector source2sample_distance
[6bd3a8d1]804
[51f14603]805        : param distance: [distance, x_offset]
806        """
807        if len(distance) < 1 or len(distance) > 2:
[574adc7]808            raise RuntimeError("The length of the size must be one or two.")
[51f14603]809        self.aperture.set_sample_distance(distance)
810
811    def set_sample2sample_distance(self, distance):
812        """
813        Set detector sample_slit2sample_distance
[6bd3a8d1]814
[51f14603]815        : param distance: [distance, x_offset]
816        """
817        if len(distance) < 1 or len(distance) > 2:
[574adc7]818            raise RuntimeError("The length of the size must be one or two.")
[51f14603]819        self.sample.set_distance(distance)
[6bd3a8d1]820
[51f14603]821    def set_sample2detector_distance(self, distance):
822        """
823        Set detector sample2detector_distance
[6bd3a8d1]824
[51f14603]825        : param distance: [distance, x_offset]
826        """
827        if len(distance) < 1 or len(distance) > 2:
[574adc7]828            raise RuntimeError("The length of the size must be one or two.")
[51f14603]829        self.detector.set_distance(distance)
[6bd3a8d1]830
[51f14603]831    def get_all_instrument_params(self):
832        """
833        Get all instrumental parameters
834        """
835        self.mass = self.get_neutron_mass()
836        self.spectrum = self.get_spectrum()
837        self.source_aperture_size = self.get_source_aperture_size()
838        self.sample_aperture_size = self.get_sample_aperture_size()
839        self.detector_pix_size = self.get_detector_pix_size()
840        self.detector_size = self.get_detector_size()
841        self.source2sample_distance = self.get_source2sample_distance()
842        self.sample2sample_distance = self.get_sample2sample_distance()
843        self.sample2detector_distance = self.get_sample2detector_distance()
[6bd3a8d1]844
[51f14603]845    def get_detector_qrange(self):
846        """
847        get max detector q ranges
[6bd3a8d1]848
[51f14603]849        : return: qx_min, qx_max, qy_min, qy_max tuple
850        """
851        if len(self.qxrange) != 2 or len(self.qyrange) != 2:
852            return None
853        qx_min = self.qxrange[0]
854        qx_max = self.qxrange[1]
855        qy_min = self.qyrange[0]
856        qy_max = self.qyrange[1]
[6bd3a8d1]857
[51f14603]858        return qx_min, qx_max, qy_min, qy_max
[6bd3a8d1]859
[51f14603]860    def _rotate_z(self, x_value, y_value, theta=0.0):
861        """
862        Rotate x-y cordinate around z-axis by theta
863        : x_value: numpy array of x values
864        : y_value: numpy array of y values
865        : theta: angle to rotate by in rad
[6bd3a8d1]866
[51f14603]867        :return: x_prime, y-prime
[6bd3a8d1]868        """
[51f14603]869        # rotate by theta
870        x_prime = x_value * math.cos(theta) + y_value * math.sin(theta)
871        y_prime = -x_value * math.sin(theta) + y_value * math.cos(theta)
[6bd3a8d1]872
[51f14603]873        return x_prime, y_prime
[6bd3a8d1]874
[51f14603]875    def _gaussian2d(self, x_val, y_val, x0_val, y0_val,
876                    sigma_x, sigma_y, sigma_r):
877        """
878        Calculate 2D Gaussian distribution
879        : x_val: x value
880        : y_val: y value
881        : x0_val: mean value in x-axis
882        : y0_val: mean value in y-axis
883        : sigma_x: variance in x-direction
884        : sigma_y: variance in y-direction
[6bd3a8d1]885
[51f14603]886        : return: gaussian (value)
887        """
888        # phi values at each points (not at the center)
889        x_value = x_val - x0_val
890        y_value = y_val - y0_val
[9a5097c]891        phi_i = np.arctan2(y_val, x_val)
[6bd3a8d1]892
[51f14603]893        # phi correction due to the gravity shift (in phi)
894        phi_0 = math.atan2(y0_val, x0_val)
895        phi_i = phi_i - phi_0 + self.gravity_phi
896
[9a5097c]897        sin_phi = np.sin(self.gravity_phi)
898        cos_phi = np.cos(self.gravity_phi)
[6bd3a8d1]899
[51f14603]900        x_p = x_value * cos_phi + y_value * sin_phi
901        y_p = -x_value * sin_phi + y_value * cos_phi
[6bd3a8d1]902
[51f14603]903        new_sig_x = sqrt(sigma_r * sigma_r / (sigma_x * sigma_x) + 1)
904        new_sig_y = sqrt(sigma_r * sigma_r / (sigma_y * sigma_y) + 1)
905        new_x = x_p * cos_phi / new_sig_x - y_p * sin_phi
906        new_x /= sigma_x
907        new_y = x_p * sin_phi / new_sig_y + y_p * cos_phi
908        new_y /= sigma_y
909
910        nu_value = -0.5 * (new_x * new_x + new_y * new_y)
911
[9a5097c]912        gaussian = np.exp(nu_value)
[51f14603]913        # normalizing factor correction
914        gaussian /= gaussian.sum()
915
916        return gaussian
917
918    def _gaussian2d_polar(self, x_val, y_val, x0_val, y0_val,
[b6627d9]919                          sigma_x, sigma_y, sigma_r):
[51f14603]920        """
921        Calculate 2D Gaussian distribution for polar coodinate
922        : x_val: x value
923        : y_val: y value
924        : x0_val: mean value in x-axis
925        : y0_val: mean value in y-axis
926        : sigma_x: variance in r-direction
927        : sigma_y: variance in phi-direction
928        : sigma_r: wavelength variance in r-direction
[6bd3a8d1]929
[51f14603]930        : return: gaussian (value)
931        """
932        sigma_x = sqrt(sigma_x * sigma_x + sigma_r * sigma_r)
933        # call gaussian1d
[b6627d9]934        gaussian = self._gaussian1d(x_val, x0_val, sigma_x)
[51f14603]935        gaussian *= self._gaussian1d(y_val, y0_val, sigma_y)
[6bd3a8d1]936
[51f14603]937        # normalizing factor correction
938        if sigma_x != 0 and sigma_y != 0:
939            gaussian *= sqrt(2 * pi)
940        return gaussian
[6bd3a8d1]941
[51f14603]942    def _gaussian1d(self, value, mean, sigma):
943        """
944        Calculate 1D Gaussian distribution
945        : value: value
946        : mean: mean value
947        : sigma: variance
[6bd3a8d1]948
[51f14603]949        : return: gaussian (value)
950        """
951        # default
952        gaussian = 1.0
953        if sigma != 0:
954            # get exponent
955            nu_value = (value - mean) / sigma
956            nu_value *= nu_value
957            nu_value *= -0.5
[9a5097c]958            gaussian *= np.exp(nu_value)
[51f14603]959            gaussian /= sigma
960            # normalize
961            gaussian /= sqrt(2 * pi)
[6bd3a8d1]962
[51f14603]963        return gaussian
[6bd3a8d1]964
[51f14603]965    def _atan_phi(self, qy_value, qx_value):
966        """
967        Find the angle phi of q on the detector plane for qx_value, qy_value given
968        : qx_value: x component of q
969        : qy_value: y component of q
[6bd3a8d1]970
[51f14603]971        : return phi: the azimuthal angle of q on x-y plane
972        """
973        phi = math.atan2(qy_value, qx_value)
974        return phi
975
976    def _get_detector_qxqy_pixels(self):
977        """
978        Get the pixel positions of the detector in the qx_value-qy_value space
979        """
[6bd3a8d1]980
[51f14603]981        # update all param values
982        self.get_all_instrument_params()
[6bd3a8d1]983
[51f14603]984        # wavelength
985        wavelength = self.wave.wavelength
986        # Gavity correction
987        delta_y = self._get_beamcenter_drop()  # in cm
[6bd3a8d1]988
[51f14603]989        # detector_pix size
990        detector_pix_size = self.detector_pix_size
991        # Square or circular pixel
992        if len(detector_pix_size) == 1:
993            pix_x_size = detector_pix_size[0]
994            pix_y_size = detector_pix_size[0]
995        # rectangular pixel pixel
996        elif len(detector_pix_size) == 2:
997            pix_x_size = detector_pix_size[0]
998            pix_y_size = detector_pix_size[1]
999        else:
[574adc7]1000            raise ValueError(" Input value format error...")
[51f14603]1001        # Sample to detector distance = sample slit to detector
1002        # minus sample offset
1003        sample2detector_distance = self.sample2detector_distance[0] - \
1004                                    self.sample2sample_distance[0]
1005        # detector offset in x-direction
1006        detector_offset = 0
1007        try:
1008            detector_offset = self.sample2detector_distance[1]
1009        except:
[c155a16]1010            logger.error(sys.exc_value)
[6bd3a8d1]1011
[51f14603]1012        # detector size in [no of pix_x,no of pix_y]
1013        detector_pix_nums_x = self.detector_size[0]
[6bd3a8d1]1014
[51f14603]1015        # get pix_y if it exists, otherwse take it from [0]
1016        try:
1017            detector_pix_nums_y = self.detector_size[1]
1018        except:
1019            detector_pix_nums_y = self.detector_size[0]
[6bd3a8d1]1020
[51f14603]1021        # detector offset in pix number
1022        offset_x = detector_offset / pix_x_size
1023        offset_y = delta_y / pix_y_size
[6bd3a8d1]1024
[51f14603]1025        # beam center position in pix number (start from 0)
1026        center_x, center_y = self._get_beamcenter_position(detector_pix_nums_x,
[b6627d9]1027                                                           detector_pix_nums_y,
1028                                                           offset_x, offset_y)
[51f14603]1029        # distance [cm] from the beam center on detector plane
[9a5097c]1030        detector_ind_x = np.arange(detector_pix_nums_x)
1031        detector_ind_y = np.arange(detector_pix_nums_y)
[51f14603]1032
1033        # shif 0.5 pixel so that pix position is at the center of the pixel
1034        detector_ind_x = detector_ind_x + 0.5
1035        detector_ind_y = detector_ind_y + 0.5
1036
1037        # the relative postion from the beam center
1038        detector_ind_x = detector_ind_x - center_x
1039        detector_ind_y = detector_ind_y - center_y
[6bd3a8d1]1040
[51f14603]1041        # unit correction in cm
1042        detector_ind_x = detector_ind_x * pix_x_size
1043        detector_ind_y = detector_ind_y * pix_y_size
[6bd3a8d1]1044
[9a5097c]1045        qx_value = np.zeros(len(detector_ind_x))
1046        qy_value = np.zeros(len(detector_ind_y))
[51f14603]1047        i = 0
1048
1049        for indx in detector_ind_x:
1050            qx_value[i] = self._get_qx(indx, sample2detector_distance, wavelength)
1051            i += 1
1052        i = 0
1053        for indy in detector_ind_y:
1054            qy_value[i] = self._get_qx(indy, sample2detector_distance, wavelength)
1055            i += 1
[6bd3a8d1]1056
[51f14603]1057        # qx_value and qy_value values in array
1058        qx_value = qx_value.repeat(detector_pix_nums_y)
1059        qx_value = qx_value.reshape(detector_pix_nums_x, detector_pix_nums_y)
1060        qy_value = qy_value.repeat(detector_pix_nums_x)
1061        qy_value = qy_value.reshape(detector_pix_nums_y, detector_pix_nums_x)
1062        qy_value = qy_value.transpose()
1063
1064        # p min and max values among the center of pixels
[9a5097c]1065        self.qx_min = np.min(qx_value)
1066        self.qx_max = np.max(qx_value)
1067        self.qy_min = np.min(qy_value)
1068        self.qy_max = np.max(qy_value)
[6bd3a8d1]1069
[51f14603]1070        # Appr. min and max values of the detector display limits
1071        # i.e., edges of the last pixels.
1072        self.qy_min += self._get_qx(-0.5 * pix_y_size,
[b6627d9]1073                                    sample2detector_distance, wavelength)
[51f14603]1074        self.qy_max += self._get_qx(0.5 * pix_y_size,
[b6627d9]1075                                    sample2detector_distance, wavelength)
[51f14603]1076        #if self.qx_min == self.qx_max:
1077        self.qx_min += self._get_qx(-0.5 * pix_x_size,
[b6627d9]1078                                    sample2detector_distance, wavelength)
[51f14603]1079        self.qx_max += self._get_qx(0.5 * pix_x_size,
1080                                    sample2detector_distance, wavelength)
[6bd3a8d1]1081
[51f14603]1082        # min and max values of detecter
1083        self.detector_qx_min = self.qx_min
1084        self.detector_qx_max = self.qx_max
1085        self.detector_qy_min = self.qy_min
1086        self.detector_qy_max = self.qy_max
[6bd3a8d1]1087
[51f14603]1088        # try to set it as a Data2D otherwise pass (not required for now)
1089        try:
[b699768]1090            from sas.sascalc.dataloader.data_info import Data2D
[51f14603]1091            output = Data2D()
[9a5097c]1092            inten = np.zeros_like(qx_value)
[b6627d9]1093            output.data = inten
[51f14603]1094            output.qx_data = qx_value
1095            output.qy_data = qy_value
1096        except:
[c155a16]1097            logger.error(sys.exc_value)
[6bd3a8d1]1098
[51f14603]1099        return output
[6bd3a8d1]1100
[51f14603]1101    def _get_qx(self, dx_size, det_dist, wavelength):
1102        """
1103        :param dx_size: x-distance from beam center [cm]
1104        :param det_dist: sample to detector distance [cm]
[6bd3a8d1]1105
[51f14603]1106        :return: q-value at the given position
1107        """
1108        # Distance from beam center in the plane of detector
1109        plane_dist = dx_size
1110        # full scattering angle on the x-axis
[9a5097c]1111        theta = np.arctan(plane_dist / det_dist)
1112        qx_value = (2.0 * pi / wavelength) * np.sin(theta)
[51f14603]1113        return qx_value
[6bd3a8d1]1114
[51f14603]1115    def _get_polar_value(self, qx_value, qy_value):
1116        """
1117        Find qr_value and phi from qx_value and qy_value values
[6bd3a8d1]1118
[51f14603]1119        : return qr_value, phi
1120        """
1121        # find |q| on detector plane
1122        qr_value = sqrt(qx_value*qx_value + qy_value*qy_value)
1123        # find angle phi
1124        phi = self._atan_phi(qy_value, qx_value)
[6bd3a8d1]1125
[51f14603]1126        return qr_value, phi
[6bd3a8d1]1127
[51f14603]1128    def _get_beamcenter_position(self, num_x, num_y, offset_x, offset_y):
1129        """
1130        :param num_x: number of pixel in x-direction
1131        :param num_y: number of pixel in y-direction
1132        :param offset: detector offset in x-direction in pix number
[6bd3a8d1]1133
[51f14603]1134        :return: pix number; pos_x, pos_y in pix index
1135        """
1136        # beam center position
1137        pos_x = num_x / 2
1138        pos_y = num_y / 2
1139
1140        # correction for offset
1141        pos_x += offset_x
1142        # correction for gravity that is always negative
1143        pos_y -= offset_y
1144
1145        return pos_x, pos_y
1146
1147    def _get_beamcenter_drop(self):
1148        """
1149        Get the beam center drop (delta y) in y diection due to gravity
[6bd3a8d1]1150
[51f14603]1151        :return delta y: the beam center drop in cm
1152        """
1153        # Check if mass == 0 (X-ray).
1154        if self.mass == 0:
1155            return 0
1156        # Covert unit from A to cm
1157        unit_cm = 1e-08
1158        # Velocity of neutron in horizontal direction (~ actual velocity)
1159        velocity = _PLANK_H / (self.mass * self.wave.wavelength * unit_cm)
1160        # Compute delta y
1161        delta_y = 0.5
1162        delta_y *= _GRAVITY
1163        sampletodetector = self.sample2detector_distance[0] - \
1164                                    self.sample2sample_distance[0]
1165        delta_y *= sampletodetector
1166        delta_y *= (self.source2sample_distance[0] + self.sample2detector_distance[0])
1167        delta_y /= (velocity * velocity)
1168
1169        return delta_y
Note: See TracBrowser for help on using the repository browser.