source: sasview/src/sas/calculator/resolution_calculator.py @ b6627d9

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since b6627d9 was b6627d9, checked in by Doucet, Mathieu <doucetm@…>, 9 years ago

pylint fixes

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