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

ESS_GUIESS_GUI_batch_fittingESS_GUI_iss879ESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since b8080e1 was b8080e1, checked in by Piotr Rozyczko <rozyczko@…>, 7 months ago

cherry picking sascalc changes from master SASVIEW-996
minor unit test fixes

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