source: sasview/calculator/resolution_calculator.py @ 7bc88bf

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 7bc88bf was 19637b1, checked in by Jae Cho <jhjcho@…>, 14 years ago

fixed a small bug and removed unnecessary raise
  • Property mode set to 100644
File size: 29.1 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 Neutron
9from instrument import Aperture
10# import math stuffs
11from math import pi
12from math import sqrt
13from math import cos
14from math import sin
15from math import atan
16from math import atan2
17from math import pow
18from math import asin
19from math import tan
20
21import numpy
22
23#Plank's constant in cgs unit
24_PLANK_H = 6.62606896E-27
25#Gravitational acc. in cgs unit
26_GRAVITY = 981.0
27
28class ResolutionCalculator(object):
29    """
30    compute resolution in 2D
31    """
32    def __init__(self):
33       
34        # wavelength
35        self.wave = Neutron()
36        # sample
37        self.sample = Sample()
38        # aperture
39        self.aperture = Aperture()
40        # detector
41        self.detector = Detector()
42        # 2d image of the resolution
43        self.image = []
44        # resolutions
45        self.sigma_1 = 0
46        self.sigma_2 = 0
47        self.sigma_1d = 0
48        # q min and max
49        self.qx_min = -0.3
50        self.qx_max = 0.3
51        self.qy_min = -0.3
52        self.qy_max = 0.3
53        # q min and max of the detector
54        self.detector_qx_min = -0.3
55        self.detector_qx_max = 0.3
56        self.detector_qy_min = -0.3
57        self.detector_qy_max = 0.3
58        # plots
59        self.plot = None
60        # instrumental params defaults
61        self.mass = 0
62        self.intensity = 0
63        self.wavelength = 0
64        self.wavelength_spread = 0
65        self.source_aperture_size = []
66        self.source2sample_distance = []
67        self.sample2sample_distance = []
68        self.sample_aperture_size = []
69        self.sample2detector_distance = []
70        self.detector_pix_size = []
71        self.detector_size = []
72        # get all the values of the instrumental parameters
73        self.get_all_instrument_params()
74       
75    def compute_and_plot(self, qx_value, qy_value, qx_min, qx_max, 
76                          qy_min, qy_max, coord = 'polar'):
77        """
78        Compute the resolution
79        : qx_value: x component of q
80        : qy_value: y component of q
81        """
82        # compute 2d resolution
83        _, _, sigma_1, sigma_2 = self.compute(qx_value, qy_value, coord)
84        # make image
85        image = self.get_image(qx_value, qy_value, sigma_1, sigma_2, 
86                               qx_min, qx_max, qy_min, qy_max, coord)
87        # plot image
88        return self.plot_image(image) 
89
90    def compute(self, qx_value, qy_value, coord = 'polar'):
91        """
92        Compute the Q resoltuion in || and + direction of 2D
93        : qx_value: x component of q
94        : qy_value: y component of q
95        """
96        # make sure to update all the variables need.
97        self.get_all_instrument_params()
98        # wavelength
99        lamb = self.wavelength
100   
101        if lamb == 0:
102            msg = "Can't compute the resolution: the wavelength is zero..." 
103            raise RuntimeError, msg
104        # wavelength spread
105        lamb_spread = self.wavelength_spread
106        # Find polar values
107        qr_value, phi = self._get_polar_value(qx_value, qy_value)
108        # vacuum wave transfer
109        knot = 2*pi/lamb
110        # scattering angle theta; always true for plane detector
111        # aligned vertically to the ko direction
112        if qr_value > knot:
113            theta = pi/2
114        else:
115            theta = asin(qr_value/knot)
116        # source aperture size
117        rone = self.source_aperture_size
118        # sample aperture size
119        rtwo = self.sample_aperture_size
120        # detector pixel size
121        rthree = self.detector_pix_size
122        # source to sample(aperture) distance
123        l_ssa = self.source2sample_distance[0]
124        # sample(aperture) to detector distance
125        l_sad = self.sample2detector_distance[0]
126        # sample (aperture) to sample distance
127        l_sas = self.sample2sample_distance[0]
128        # source to sample distance
129        l_one = l_ssa + l_sas
130        # sample to detector distance
131        l_two = l_sad - l_sas
132       
133        # Sample offset correction for l_one and Lp on variance calculation
134        l1_cor = (l_ssa * l_two) / (l_sas + l_two)
135        lp_cor = (l_ssa * l_two) / (l_one + l_two)
136        # the radial distance to the pixel from the center of the detector
137        radius = tan(theta)*l_two
138        #Lp = l_one*l_two/(l_one+l_two)
139        # default polar coordinate
140        comp1 = 'radial'
141        comp2 = 'phi'
142        # in the case of the cartesian coordinate
143        if coord == 'cartesian':
144            comp1 = 'x'
145            comp2 = 'y'
146
147        # sigma in the radial/x direction
148        # for source aperture
149        sigma_1  = self.get_variance(rone, l1_cor, phi, comp1)
150        # for sample apperture
151        sigma_1 += self.get_variance(rtwo, lp_cor, phi, comp1)
152        # for detector pix
153        sigma_1 += self.get_variance(rthree, l_two, phi, comp1)
154        # for gravity term
155        sigma_1 +=  self.get_variance_gravity(l_ssa, l_sad, lamb, lamb_spread, 
156                             phi, comp1, 'on') 
157        # for wavelength spread
158        # reserve for 1d calculation
159        sigma_wave_1 = self.get_variance_wave(radius, l_two, lamb_spread, 
160                                          phi, comp1, 'on')
161        # for 1d
162        variance_1d_1 = sigma_1/2 +sigma_wave_1
163        # normalize
164        variance_1d_1 = knot*knot*variance_1d_1/12
165       
166        # for 2d
167        sigma_1 += sigma_wave_1
168        # normalize
169        sigma_1 = knot*sqrt(sigma_1/12)
170
171        # sigma in the phi/y direction
172        # for source apperture
173        sigma_2  = self.get_variance(rone, l1_cor, phi, comp2)
174        # for sample apperture
175        sigma_2 += self.get_variance(rtwo, lp_cor, phi, comp2)
176        # for detector pix
177        sigma_2 += self.get_variance(rthree, l_two, phi, comp2)
178        # for gravity term
179        sigma_2 +=  self.get_variance_gravity(l_ssa, l_sad, lamb, lamb_spread, 
180                             phi, comp2, 'on')
181        # for wavelength spread
182        # reserve for 1d calculation
183        sigma_wave_2 = self.get_variance_wave(radius, l_two, lamb_spread, 
184                                          phi, comp2, 'on') 
185        # for 1d
186        variance_1d_2 = sigma_2/2 +sigma_wave_2
187        # normalize
188        variance_1d_2 = knot*knot*variance_1d_2/12
189       
190        # for 2d
191        sigma_2 += sigma_wave_2
192        # normalize
193        sigma_2 =  knot*sqrt(sigma_2/12)
194
195        # set sigmas
196        self.sigma_1 = sigma_1
197        self.sigma_2 = sigma_2
198       
199        self.sigma_1d = sqrt(variance_1d_1 + variance_1d_2)
200        return qr_value, phi, sigma_1, sigma_2
201   
202    def get_image(self, qx_value, qy_value, sigma_1, sigma_2,
203                  qx_min, qx_max, qy_min, qy_max, coord = 'polar'): 
204        """
205        Get the resolution in polar coordinate ready to plot
206        : qx_value: qx_value value
207        : qy_value: qy_value value
208        : sigma_1: variance in r direction
209        : sigma_2: variance in phi direction
210        : coord: coordinate system of image, 'polar' or 'cartesian'
211        """
212        # Get  qx_max and qy_max...
213        output = self._get_detector_qxqy_pixels()
214        # Set qx_value/qy_value min/max
215        #qx_min = self.qx_min
216        #qx_max = self.qx_max
217        #qy_min = self.qy_min
218        #qy_max = self.qy_max
219
220        # Find polar values
221        qr_value, phi = self._get_polar_value(qx_value, qy_value)
222       
223        # Check whether the q value is within the detector range
224        msg = "Invalid input: Q value out of the detector range..."
225        if qx_min < self.qx_min:
226            self.qx_min = qx_min
227            #raise ValueError, msg
228        if qx_max > self.qx_max:
229            self.qx_max = qx_max
230            #raise ValueError, msg
231        if qy_min < self.qy_min:
232            self.qy_min = qy_min
233            #raise ValueError, msg
234        if qy_max > self.qy_max:
235            self.qy_max = qy_max
236            #raise ValueError, msg
237
238        # Make an empty graph in the detector scale
239        dx_size = (self.qx_max - self.qx_min) / (1000 - 1)
240        dy_size = (self.qy_max - self.qy_min) / (1000 - 1)
241        x_val = numpy.arange(self.qx_min, self.qx_max, dx_size)
242        y_val = numpy.arange(self.qy_max, self.qy_min, -dy_size)
243        q_1, q_2 = numpy.meshgrid(x_val, y_val)
244
245        # check whether polar or cartesian
246        if coord == 'polar':
247            q_1, q_2 = self._rotate_z(q_1, q_2, phi)
248            qc_1 = qr_value
249            qc_2 = 0.0
250        else:
251            # catesian coordinate
252            # qx_center
253            qc_1 = qx_value
254            # qy_center
255            qc_2 = qy_value
256           
257        # Calculate the 2D Gaussian distribution image
258        image = self._gaussian2d(q_1, q_2, qc_1, qc_2, sigma_1, sigma_2)
259        # Add it if there are more than one inputs.
260        if len(self.image) > 0:
261            self.image += image
262        else:
263            self.image = image
264
265        return self.image
266   
267    def plot_image(self, image):
268        """
269        Plot image using pyplot
270        : image: 2d resolution image
271       
272        : return plt: pylab object
273        """
274        import matplotlib.pyplot as plt
275
276        self.plot = plt
277        plt.xlabel('$\\rm{Q}_{x} [A^{-1}]$')
278        plt.ylabel('$\\rm{Q}_{y} [A^{-1}]$')
279        # Max value of the image
280        max = numpy.max(image)
281        # Image
282        im = plt.imshow(image, 
283                extent = [self.qx_min, self.qx_max, self.qy_min, self.qy_max])
284
285        # bilinear interpolation to make it smoother
286        im.set_interpolation('bilinear')
287
288        return plt
289   
290    def reset_image(self):
291        """
292        Reset image to default (=[])
293        """
294        self.image = []
295       
296    def get_variance(self, size = [], distance = 0, phi = 0, comp = 'radial'):
297        """
298        Get the variance when the slit/pinhole size is given
299        : size: list that can be one(diameter for circular)
300                or two components(lengths for rectangular)
301        : distance: [z, x] where z along the incident beam, x // qx_value
302        : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial'
303       
304        : return variance: sigma^2
305        """   
306        # check the length of size (list)
307        len_size = len(size)
308       
309        # define sigma component direction
310        if comp == 'radial':
311            phi_x = cos(phi)
312            phi_y = sin(phi)
313        elif comp == 'phi':
314            phi_x = sin(phi)
315            phi_y = cos(phi)
316        elif comp == 'x':
317            phi_x = 1
318            phi_y = 0
319        elif comp == 'y':
320            phi_x = 0
321            phi_y = 1
322        else:
323            phi_x = 0
324            phi_y = 0
325        # calculate each component
326        # for pinhole w/ radius = size[0]/2
327        if len_size == 1:
328            x_comp = (0.5 * size[0]) * sqrt(3)
329            y_comp = 0
330        # for rectangular slit
331        elif len_size == 2:
332            x_comp = size[0] * phi_x
333            y_comp = size[1] * phi_y
334        # otherwise
335        else:
336            raise ValueError, " Improper input..."
337        # get them squared 
338        sigma  = x_comp * x_comp
339        sigma += y_comp * y_comp
340        # normalize by distance
341        sigma /= (distance * distance)
342
343        return sigma
344
345    def get_variance_wave(self, radius, distance, spread, phi, 
346                          comp = 'radial', switch = 'on'):
347        """
348        Get the variance when the wavelength spread is given
349       
350        : radius: the radial distance from the beam center to the pix of q
351        : distance: sample to detector distance
352        : spread: wavelength spread (ratio)
353        : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial'
354       
355        : return variance: sigma^2
356        """
357        if switch.lower() == 'off':
358            return 0
359        # check the singular point
360        if distance == 0 or comp == 'phi':
361            return 0
362        else:
363            # calculate sigma^2
364            sigma = 2 * pow(radius/distance*spread, 2)
365            if comp == 'x':
366                sigma *= (cos(phi)*cos(phi))
367            elif comp == 'y':
368                sigma *= (sin(phi)*sin(phi))
369            else:
370                sigma *= 1         
371               
372            return sigma
373
374    def get_variance_gravity(self, s_distance, d_distance, wavelength, spread, 
375                             phi, comp = 'radial', switch = 'on'):
376        """
377        Get the variance from gravity when the wavelength spread is given
378       
379        : s_distance: source to sample distance
380        : d_distance: sample to detector distance
381        : wavelength: wavelength
382        : spread: wavelength spread (ratio)
383        : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial'
384       
385        : return variance: sigma^2
386        """
387        if switch.lower() == 'off':
388            return 0
389        if self.mass == 0.0:
390            return 0
391        # check the singular point
392        if d_distance == 0 or comp == 'x':
393            return 0
394        else:
395            # neutron mass in cgs unit
396            self.mass = self.get_neutron_mass()
397            # plank constant in cgs unit
398            h_constant = _PLANK_H
399            # gravity in cgs unit
400            gravy = _GRAVITY
401            # m/h
402            m_over_h = self.mass /h_constant
403            # A value
404            a_value = d_distance * (s_distance + d_distance)
405            a_value *= pow(m_over_h / 2, 2)
406            a_value *= gravy
407            # unit correction (1/cm to 1/A) for A and d_distance below
408            a_value *= 1.0E-16
409           
410            # calculate sigma^2
411            sigma = pow(a_value / d_distance, 2)
412            sigma *= pow(wavelength, 4)
413            sigma *= pow(spread, 2)
414            sigma *= 8
415           
416            # only for the polar coordinate
417            if comp == 'radial':
418                sigma *= (sin(phi) * sin(phi))
419            elif comp == 'phi':
420                sigma *= (cos(phi) * cos(phi))
421           
422            return sigma
423       
424    def get_intensity(self):
425        """
426        Get intensity
427        """
428        return self.wave.intensity
429
430    def get_wavelength(self):
431        """
432        Get wavelength
433        """
434        return self.wave.wavelength
435   
436    def get_wavelength_spread(self):
437        """
438        Get wavelength spread
439        """
440        return self.wave.wavelength_spread
441   
442    def get_neutron_mass(self):
443        """
444        Get Neutron mass
445        """
446        return self.wave.mass
447   
448    def get_source_aperture_size(self):
449        """
450        Get source aperture size
451        """
452        return self.aperture.source_size
453   
454    def get_sample_aperture_size(self):
455        """
456        Get sample aperture size
457        """
458        return self.aperture.sample_size
459   
460    def get_detector_pix_size(self):
461        """
462        Get detector pixel size
463        """
464        return self.detector.pix_size
465   
466    def get_detector_size(self):
467        """
468        Get detector size
469        """
470        return self.detector.size   
471     
472    def get_source2sample_distance(self):
473        """
474        Get detector source2sample_distance
475        """
476        return self.aperture.sample_distance
477   
478    def get_sample2sample_distance(self):
479        """
480        Get detector sampleslitsample_distance
481        """
482        return self.sample.distance
483   
484    def get_sample2detector_distance(self):
485        """
486        Get detector sample2detector_distance
487        """
488        return self.detector.distance
489   
490    def set_intensity(self, intensity):
491        """
492        Set intensity
493        """
494        self.wave.set_intensity(intensity)
495   
496    def set_wavelength(self, wavelength):
497        """
498        Set wavelength
499        """
500        self.wave.set_wavelength(wavelength)
501   
502    def set_wavelength_spread(self, wavelength_spread):
503        """
504        Set wavelength spread
505        """
506        self.wave.set_wavelength_spread(wavelength_spread)
507   
508    def set_source_aperture_size(self, size):
509        """
510        Set source aperture size
511       
512        : param size: [dia_value] or [x_value, y_value]
513        """       
514        if len(size) < 1 or len(size) > 2:
515            raise RuntimeError, "The length of the size must be one or two."
516        self.aperture.set_source_size(size)
517       
518    def set_neutron_mass(self, mass):
519        """
520        Set Neutron mass
521        """
522        self.wave.set_mass(mass)
523       
524    def set_sample_aperture_size(self, size):
525        """
526        Set sample aperture size
527       
528        : param size: [dia_value] or [xheight_value, yheight_value]
529        """
530        if len(size) < 1 or len(size) > 2:
531            raise RuntimeError, "The length of the size must be one or two."
532        self.aperture.set_sample_size(size)
533   
534    def set_detector_pix_size(self, size):
535        """
536        Set detector pixel size
537        """
538        self.detector.set_pix_size(size)
539       
540    def set_detector_size(self, size):
541        """
542        Set detector size in number of pixels
543        : param size: [pixel_nums] or [x_pix_num, yx_pix_num]
544        """
545        self.detector.set_size(size)
546       
547    def set_source2sample_distance(self, distance):
548        """
549        Set detector source2sample_distance
550       
551        : param distance: [distance, x_offset]
552        """
553        if len(distance) < 1 or len(distance) > 2:
554            raise RuntimeError, "The length of the size must be one or two."
555        self.aperture.set_sample_distance(distance)
556
557    def set_sample2sample_distance(self, distance):
558        """
559        Set detector sample_slit2sample_distance
560       
561        : param distance: [distance, x_offset]
562        """
563        if len(distance) < 1 or len(distance) > 2:
564            raise RuntimeError, "The length of the size must be one or two."
565        self.sample.set_distance(distance)
566   
567    def set_sample2detector_distance(self, distance):
568        """
569        Set detector sample2detector_distance
570       
571        : param distance: [distance, x_offset]
572        """
573        if len(distance) < 1 or len(distance) > 2:
574            raise RuntimeError, "The length of the size must be one or two."
575        self.detector.set_distance(distance)
576       
577    def get_all_instrument_params(self):
578        """
579        Get all instrumental parameters
580        """
581        self.intensity = self.get_intensity()
582        self.wavelength = self.get_wavelength()
583        self.wavelength_spread = self.get_wavelength_spread()
584        self.mass = self.get_neutron_mass()
585        self.source_aperture_size = self.get_source_aperture_size()
586        self.sample_aperture_size = self.get_sample_aperture_size()
587        self.detector_pix_size = self.get_detector_pix_size()
588        self.detector_size = self.get_detector_size()
589        self.source2sample_distance = self.get_source2sample_distance()
590        self.sample2sample_distance = self.get_sample2sample_distance()
591        self.sample2detector_distance = self.get_sample2detector_distance()
592
593
594         
595    def _rotate_z(self, x_value, y_value, theta= 0.0):
596        """
597        Rotate x-y cordinate around z-axis by theta
598        : x_value: numpy array of x values
599        : y_value: numpy array of y values
600        : theta: angle to rotate by in rad
601       
602        :return: x_prime, y-prime
603        """       
604        # rotate by theta
605        x_prime = x_value * cos(theta) + y_value * sin(theta)
606        y_prime =  -x_value * sin(theta) + y_value * cos(theta)
607   
608        return x_prime, y_prime
609   
610    def _gaussian2d(self, x_val, y_val, x0_val, y0_val, sigma_x, sigma_y):
611        """
612        Calculate 2D Gaussian distribution
613        : x_val: x value
614        : y_val: y value
615        : x0_val: mean value in x-axis
616        : y0_val: mean value in y-axis
617        : sigma_x: variance in x-direction
618        : sigma_y: variance in y-direction
619       
620        : return: gaussian (value)
621        """
622        # call gaussian1d
623        gaussian  = self._gaussian1d(x_val, x0_val, sigma_x)
624        gaussian *= self._gaussian1d(y_val, y0_val, sigma_y)
625 
626        # normalizing factor correction
627        if sigma_x != 0 and sigma_y != 0:
628            gaussian *= sqrt(2 * pi)
629        return gaussian
630
631    def _gaussian1d(self, value, mean, sigma):
632        """
633        Calculate 1D Gaussian distribution
634        : value: value
635        : mean: mean value
636        : sigma: variance
637       
638        : return: gaussian (value)
639        """
640        # default
641        gaussian = 1.0
642        if sigma != 0:
643            # get exponent
644            nu_value = (value - mean) / sigma
645            nu_value *= nu_value
646            nu_value *= -0.5
647            gaussian *= numpy.exp(nu_value)
648            gaussian /= sigma
649            # normalize
650            gaussian /= sqrt(2 * pi)
651           
652        return gaussian
653   
654    def _atan_phi(self, qy_value, qx_value):
655        """
656        Find the angle phi of q on the detector plane for qx_value, qy_value given
657        : qx_value: x component of q
658        : qy_value: y component of q
659       
660        : return phi: the azimuthal angle of q on x-y plane
661        """
662        # default
663        phi = 0
664        # ToDo: This is misterious - sign???
665        #qy_value = -qy_value
666        # Take care of the singular point
667        if qx_value == 0:
668            if qy_value > 0:
669                phi = pi / 2
670            elif qy_value < 0:
671                phi = -pi / 2
672            else:
673                phi = 0
674        else:
675            # the angle
676            phi = atan2(qy_value, qx_value)
677
678        return phi
679
680    def _get_detector_qxqy_pixels(self):
681        """
682        Get the pixel positions of the detector in the qx_value-qy_value space
683        """
684       
685        # update all param values
686        self.get_all_instrument_params()
687       
688        # wavelength
689        wavelength = self.wavelength
690        # Gavity correction
691        delta_y = self._get_beamcenter_drop() # in cm
692       
693        # detector_pix size
694        detector_pix_size = self.detector_pix_size
695        # Square or circular pixel
696        if len(detector_pix_size) == 1:
697            pix_x_size = detector_pix_size[0]
698            pix_y_size = detector_pix_size[0]
699        # rectangular pixel pixel
700        elif len(detector_pix_size) == 2:
701            pix_x_size = detector_pix_size[0]
702            pix_y_size = detector_pix_size[1]
703        else:
704            raise ValueError, " Input value format error..."
705        # Sample to detector distance = sample slit to detector
706        # minus sample offset
707        sample2detector_distance = self.sample2detector_distance[0] - \
708                                    self.sample2sample_distance[0]
709        # detector offset in x-direction
710        detector_offset = 0
711        try:
712            detector_offset = self.sample2detector_distance[1]
713        except:
714            pass
715       
716        # detector size in [no of pix_x,no of pix_y]
717        detector_pix_nums_x = self.detector_size[0]
718       
719        # get pix_y if it exists, otherwse take it from [0]
720        try:
721            detector_pix_nums_y = self.detector_size[1]
722        except:
723            detector_pix_nums_y = self.detector_size[0]
724       
725        # detector offset in pix number
726        offset_x = detector_offset / pix_x_size
727        offset_y = delta_y / pix_y_size
728       
729        # beam center position in pix number (start from 0)
730        center_x, center_y = self._get_beamcenter_position(detector_pix_nums_x,
731                                    detector_pix_nums_y, offset_x, offset_y)
732        # distance [cm] from the beam center on detector plane
733        detector_ind_x = numpy.arange(detector_pix_nums_x)
734        detector_ind_y = numpy.arange(detector_pix_nums_y)
735
736        # shif 0.5 pixel so that pix position is at the center of the pixel
737        detector_ind_x = detector_ind_x + 0.5
738        detector_ind_y = detector_ind_y + 0.5
739
740        # the relative postion from the beam center
741        detector_ind_x = detector_ind_x - center_x
742        detector_ind_y = detector_ind_y - center_y
743       
744        # unit correction in cm
745        detector_ind_x = detector_ind_x * pix_x_size
746        detector_ind_y = detector_ind_y * pix_y_size
747       
748        qx_value = numpy.zeros(len(detector_ind_x))
749        qy_value = numpy.zeros(len(detector_ind_y))
750        i = 0
751
752        for indx in detector_ind_x:
753            qx_value[i] = self._get_qx(indx, sample2detector_distance, wavelength)
754            i += 1
755        i = 0
756        for indy in detector_ind_y:
757            qy_value[i] = self._get_qx(indy, sample2detector_distance, wavelength)
758            i += 1
759           
760        # qx_value and qy_value values in array
761        qx_value = qx_value.repeat(detector_pix_nums_y) 
762        qx_value = qx_value.reshape(detector_pix_nums_x, detector_pix_nums_y)
763        qy_value = qy_value.repeat(detector_pix_nums_x) 
764        qy_value = qy_value.reshape(detector_pix_nums_y, detector_pix_nums_x)
765        qy_value = qy_value.transpose()
766
767        # p min and max values among the center of pixels
768        self.qx_min = numpy.min(qx_value)
769        self.qx_max = numpy.max(qx_value)
770        self.qy_min = numpy.min(qy_value)
771        self.qy_max = numpy.max(qy_value)
772               
773        # Appr. min and max values of the detector display limits
774        # i.e., edges of the last pixels.
775        self.qy_min += self._get_qx(-0.5 * pix_y_size, 
776                                sample2detector_distance, wavelength)
777        self.qy_max += self._get_qx(0.5 * pix_y_size, 
778                                sample2detector_distance, wavelength)
779        #if self.qx_min == self.qx_max:
780        self.qx_min += self._get_qx(-0.5 * pix_x_size, 
781                                sample2detector_distance, wavelength)
782        self.qx_max += self._get_qx(0.5 * pix_x_size, 
783                                    sample2detector_distance, wavelength)
784        # min and max values of detecter
785        self.detector_qx_min = self.qx_min
786        self.detector_qx_max = self.qx_max
787        self.detector_qy_min = self.qy_min
788        self.detector_qy_max = self.qy_max
789
790        # try to set it as a Data2D otherwise pass (not required for now)
791        try:
792            from DataLoader.data_info import Data2D
793            output = Data2D()
794            inten = numpy.zeros_like(qx_value)
795            output.data     = inten
796            output.qx_data  = qx_value
797            output.qy_data  = qy_value           
798        except:
799            pass
800       
801        return output#qx_value,qy_value
802       
803    def _get_qx(self, dx_size, det_dist, wavelength):
804        """
805        :param dx_size: x-distance from beam center [cm]
806        :param det_dist: sample to detector distance [cm]
807       
808        :return: q-value at the given position
809        """
810        # Distance from beam center in the plane of detector
811        plane_dist = dx_size
812        # full scattering angle on the x-axis
813        theta  = atan(plane_dist / det_dist)
814        qx_value     = (2.0 * pi / wavelength) * sin(theta)
815        return qx_value 
816   
817    def _get_polar_value(self, qx_value, qy_value):
818        """
819        Find qr_value and phi from qx_value and qy_value values
820       
821        : return qr_value, phi
822        """
823        # find |q| on detector plane
824        qr_value = sqrt(qx_value*qx_value + qy_value*qy_value)
825        # find angle phi
826        phi = self._atan_phi(qy_value, qx_value)
827       
828        return qr_value, phi
829   
830    def _get_beamcenter_position(self, num_x, num_y, offset_x, offset_y):
831        """
832        :param num_x: number of pixel in x-direction
833        :param num_y: number of pixel in y-direction
834        :param offset: detector offset in x-direction in pix number
835       
836        :return: pix number; pos_x, pos_y in pix index
837        """
838        # beam center position
839        pos_x = num_x / 2
840        pos_y = num_y / 2
841
842        # correction for offset
843        pos_x += offset_x
844        # correction for gravity that is always negative
845        pos_y -= offset_y
846
847        return pos_x, pos_y     
848
849    def _get_beamcenter_drop(self):
850        """
851        Get the beam center drop (delta y) in y diection due to gravity
852       
853        :return delta y: the beam center drop in cm
854        """
855        # Check if mass == 0 (X-ray).
856        if self.mass == 0:
857            return 0
858        # Covert unit from A to cm
859        unit_cm = 1e-08
860        # Velocity of neutron in horizontal direction (~ actual velocity)
861        velocity = _PLANK_H / (self.mass * self.wavelength * unit_cm)
862        # Compute delta y
863        delta_y = 0.5
864        delta_y *= _GRAVITY
865        delta_y *= self.sample2detector_distance[0] 
866        delta_y *= (self.source2sample_distance[0] + self.sample2detector_distance[0])
867        delta_y /= (velocity * velocity)
868
869        return delta_y     
870
871   
Note: See TracBrowser for help on using the repository browser.