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

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since d6e5b31 was d6e5b31, checked in by Piotr Rozyczko <rozyczko@…>, 6 years ago

Distance objects for resolution should really contain two parameters

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