Changeset be6e99a in sasview for calculator/resolution_calculator.py

Timestamp:

Apr 8, 2011 1:19:10 PM (13 years ago)

Author:

Jae Cho <jhjcho@…>

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

ef73618

Parents:

cdb3a9a

Message:

fixed problems in the estimator

File:

• calculator/resolution_calculator.py (modified) (25 diffs)

calculator/resolution_calculator.py

@@ -11 +11 @@
 from math import pi
 from math import sqrt
-from math import cos
-from math import sin
-from math import atan
-from math import atan2
-from math import pow
-from math import asin
-from math import tan
-
+import math 
+import scipy
 import numpy
 
@@ -43 +37 @@
         self.image = []
         # resolutions
+        # lamda in r-direction
+        self.sigma_lamda = 0
+        # x-dir (no lamda)
         self.sigma_1 = 0
+        #y-dir (no lamda)
         self.sigma_2 = 0
+        # 1D total
         self.sigma_1d = 0
         # q min and max
@@ -74 +73 @@
         
     def compute_and_plot(self, qx_value, qy_value, qx_min, qx_max, 
-                          qy_min, qy_max, coord = 'polar'):
+                          qy_min, qy_max, coord = 'cartesian'):
         """
         Compute the resolution
@@ -81 +80 @@
         """
         # compute 2d resolution
-        _, _, sigma_1, sigma_2 = self.compute(qx_value, qy_value, coord)
+        _, _, sigma_1, sigma_2, sigma_r = \
+                            self.compute(qx_value, qy_value, coord)
         # make image
-        image = self.get_image(qx_value, qy_value, sigma_1, sigma_2, 
+        image = self.get_image(qx_value, qy_value, sigma_1, sigma_2, sigma_r, 
                                qx_min, qx_max, qy_min, qy_max, coord)
         # plot image
         return self.plot_image(image) 
 
-    def compute(self, qx_value, qy_value, coord = 'polar'):
+    def compute(self, qx_value, qy_value, coord = 'cartesian'):
         """
         Compute the Q resoltuion in || and + direction of 2D
@@ -94 +94 @@
         : qy_value: y component of q
         """
+        coord = 'cartesian'
         # make sure to update all the variables need.
         self.get_all_instrument_params()
@@ -113 +114 @@
             theta = pi/2
         else:
-            theta = asin(qr_value/knot)
+            theta = math.asin(qr_value/knot)
         # source aperture size
         rone = self.source_aperture_size
@@ -135 +136 @@
         lp_cor = (l_ssa * l_two) / (l_one + l_two)
         # the radial distance to the pixel from the center of the detector
-        radius = tan(theta)*l_two
+        radius = math.tan(theta)*l_two
         #Lp = l_one*l_two/(l_one+l_two)
         # default polar coordinate
@@ -158 +159 @@
         # reserve for 1d calculation
         sigma_wave_1 = self.get_variance_wave(radius, l_two, lamb_spread, 
-                                          phi, comp1, 'on')
+                                          phi, 'radial', 'on')
         # for 1d
         variance_1d_1 = sigma_1/2 +sigma_wave_1
@@ -165 +166 @@
         
         # for 2d
-        sigma_1 += sigma_wave_1
+        #sigma_1 += sigma_wave_1
         # normalize
         sigma_1 = knot*sqrt(sigma_1/12)
-
+        sigma_r = knot*sqrt(sigma_wave_1/12)
         # sigma in the phi/y direction
         # for source apperture
         sigma_2  = self.get_variance(rone, l1_cor, phi, comp2)
+
         # for sample apperture
         sigma_2 += self.get_variance(rtwo, lp_cor, phi, comp2)
+
         # for detector pix
         sigma_2 += self.get_variance(rthree, l_two, phi, comp2)
+
         # for gravity term
         sigma_2 +=  self.get_variance_gravity(l_ssa, l_sad, lamb, lamb_spread, 
                              phi, comp2, 'on')
+
+        
         # for wavelength spread
         # reserve for 1d calculation
         sigma_wave_2 = self.get_variance_wave(radius, l_two, lamb_spread, 
-                                          phi, comp2, 'on') 
+                                          phi, 'phi', 'on') 
         # for 1d
         variance_1d_2 = sigma_2/2 +sigma_wave_2
@@ -189 +195 @@
         
         # for 2d
-        sigma_2 += sigma_wave_2
+        #sigma_2 =  knot*sqrt(sigma_2/12)
+        #sigma_2 += sigma_wave_2
         # normalize
         sigma_2 =  knot*sqrt(sigma_2/12)
@@ -195 +202 @@
         # set sigmas
         self.sigma_1 = sigma_1
+        self.sigma_lamd = sigma_r
         self.sigma_2 = sigma_2
         
         self.sigma_1d = sqrt(variance_1d_1 + variance_1d_2)
-        return qr_value, phi, sigma_1, sigma_2
-    
-    def get_image(self, qx_value, qy_value, sigma_1, sigma_2,
-                  qx_min, qx_max, qy_min, qy_max, coord = 'polar'): 
+        return qr_value, phi, sigma_1, sigma_2, sigma_r
+    
+    def get_image(self, qx_value, qy_value, sigma_1, sigma_2, sigma_r,
+                  qx_min, qx_max, qy_min, qy_max, coord = 'cartesian'): 
         """
         Get the resolution in polar coordinate ready to plot
@@ -217 +225 @@
         #qy_min = self.qy_min
         #qy_max = self.qy_max
-
-        # Find polar values
         qr_value, phi = self._get_polar_value(qx_value, qy_value)
-        
+
         # Check whether the q value is within the detector range
         msg = "Invalid input: Q value out of the detector range..."
@@ -242 +248 @@
         y_val = numpy.arange(self.qy_max, self.qy_min, -dy_size)
         q_1, q_2 = numpy.meshgrid(x_val, y_val)
-
+        #q_phi = numpy.arctan(q_1,q_2)
         # check whether polar or cartesian
         if coord == 'polar':
+            # Find polar values
+            qr_value, phi = self._get_polar_value(qx_value, qy_value)
             q_1, q_2 = self._rotate_z(q_1, q_2, phi)
             qc_1 = qr_value
             qc_2 = 0.0
+            # Calculate the 2D Gaussian distribution image
+            image = self._gaussian2d_polar(q_1, q_2, qc_1, qc_2, 
+                                 sigma_1, sigma_2, sigma_r)
         else:
             # catesian coordinate
@@ -255 +266 @@
             qc_2 = qy_value
             
-        # Calculate the 2D Gaussian distribution image
-        image = self._gaussian2d(q_1, q_2, qc_1, qc_2, sigma_1, sigma_2)
+            # Calculate the 2D Gaussian distribution image
+            image = self._gaussian2d(q_1, q_2, qc_1, qc_2, 
+                                     sigma_1, sigma_2, sigma_r)
         # Add it if there are more than one inputs.
         if len(self.image) > 0:
@@ -309 +321 @@
         # define sigma component direction
         if comp == 'radial':
-            phi_x = cos(phi)
-            phi_y = sin(phi)
+            phi_x = math.cos(phi)
+            phi_y = math.sin(phi)
         elif comp == 'phi':
-            phi_x = sin(phi)
-            phi_y = cos(phi)
+            phi_x = math.sin(phi)
+            phi_y = math.cos(phi)
         elif comp == 'x':
             phi_x = 1
@@ -362 +374 @@
         else:
             # calculate sigma^2
-            sigma = 2 * pow(radius/distance*spread, 2)
+            sigma = 2 * math.pow(radius/distance*spread, 2)
             if comp == 'x':
-                sigma *= (cos(phi)*cos(phi))
+                sigma *= (math.cos(phi)*math.cos(phi))
             elif comp == 'y':
-                sigma *= (sin(phi)*sin(phi))
+                sigma *= (math.sin(phi)*math.sin(phi))
             else:
                 sigma *= 1          
@@ -403 +415 @@
             # A value
             a_value = d_distance * (s_distance + d_distance)
-            a_value *= pow(m_over_h / 2, 2)
+            a_value *= math.pow(m_over_h / 2, 2)
             a_value *= gravy
             # unit correction (1/cm to 1/A) for A and d_distance below
@@ -409 +421 @@
             
             # calculate sigma^2
-            sigma = pow(a_value / d_distance, 2)
-            sigma *= pow(wavelength, 4)
-            sigma *= pow(spread, 2)
+            sigma = math.pow(a_value / d_distance, 2)
+            sigma *= math.pow(wavelength, 4)
+            sigma *= math.pow(spread, 2)
             sigma *= 8
             
             # only for the polar coordinate
-            if comp == 'radial':
-                sigma *= (sin(phi) * sin(phi))
-            elif comp == 'phi':
-                sigma *= (cos(phi) * cos(phi))
+            #if comp == 'radial':
+            #    sigma *= (math.sin(phi) * math.sin(phi))
+            #elif comp == 'phi':
+            #    sigma *= (math.cos(phi) * math.cos(phi))
             
             return sigma
@@ -521 +533 @@
         """
         self.wave.set_mass(mass)
+        self.mass = mass
         
     def set_sample_aperture_size(self, size):
@@ -603 +616 @@
         """        
         # rotate by theta
-        x_prime = x_value * cos(theta) + y_value * sin(theta)
-        y_prime =  -x_value * sin(theta) + y_value * cos(theta)
+        x_prime = x_value * math.cos(theta) + y_value * math.sin(theta)
+        y_prime =  -x_value * math.sin(theta) + y_value * math.cos(theta)
     
         return x_prime, y_prime
     
-    def _gaussian2d(self, x_val, y_val, x0_val, y0_val, sigma_x, sigma_y):
+    def _gaussian2d(self, x_val, y_val, x0_val, y0_val, 
+                    sigma_x, sigma_y, sigma_r):
         """
         Calculate 2D Gaussian distribution
@@ -620 +634 @@
         : return: gaussian (value)
         """
+        # phi values at each points (not at the center)
+        x_value = x_val - x0_val
+        y_value = y_val - y0_val
+        phi_i = numpy.arctan2(y_val, x_val)
+
+        sin_phi = numpy.sin(phi_i)
+        cos_phi = numpy.cos(phi_i)
+        
+        x_p = x_value * cos_phi + y_value * sin_phi
+        y_p = -x_value * sin_phi + y_value * cos_phi
+        
+        new_x = x_p * cos_phi / (sigma_r / sigma_x + 1) - y_p * sin_phi
+        new_x /= sigma_x
+        new_y = x_p * sin_phi / (sigma_r / sigma_y + 1) + y_p * cos_phi
+        new_y /= sigma_y
+
+        nu_value = -0.5 *(new_x * new_x + new_y * new_y)
+
+        gaussian = numpy.exp(nu_value)
+        # normalizing factor correction
+        gaussian /= gaussian.sum()
+
+        return gaussian
+
+    def _gaussian2d_polar(self, x_val, y_val, x0_val, y0_val, 
+                        sigma_x, sigma_y, sigma_r):
+        """
+        Calculate 2D Gaussian distribution for polar coodinate 
+        : x_val: x value
+        : y_val: y value
+        : x0_val: mean value in x-axis
+        : y0_val: mean value in y-axis
+        : sigma_x: variance in r-direction
+        : sigma_y: variance in phi-direction
+        : sigma_r: wavelength variance in r-direction
+        
+        : return: gaussian (value)
+        """
+        sigma_x = sqrt(sigma_x * sigma_x + sigma_r * sigma_r)
         # call gaussian1d 
         gaussian  = self._gaussian1d(x_val, x0_val, sigma_x)
@@ -628 +681 @@
             gaussian *= sqrt(2 * pi)
         return gaussian
-
+    
     def _gaussian1d(self, value, mean, sigma):
         """
@@ -660 +713 @@
         : return phi: the azimuthal angle of q on x-y plane
         """
+        phi = math.atan2(qy_value, qx_value)
+        return phi
         # default
         phi = 0
@@ -674 +729 @@
         else:
             # the angle
-            phi = atan2(qy_value, qx_value)
+            phi = math.atan2(qy_value, qx_value)
 
         return phi
@@ -811 +866 @@
         plane_dist = dx_size
         # full scattering angle on the x-axis
-        theta  = atan(plane_dist / det_dist)
-        qx_value     = (2.0 * pi / wavelength) * sin(theta)
+        theta  = math.atan(plane_dist / det_dist)
+        qx_value     = (2.0 * pi / wavelength) * math.sin(theta)
         return qx_value