source: sasview/calculator/resolution_calculator.py @ 4b5bd73

"""
This object is a small tool to allow user to quickly
determine the variance in q  from the 
instrumental parameters.
"""
from instrument import Sample
from instrument import Detector
from instrument import Neutron
from instrument import Aperture
# import math stuffs
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 numpy

#Plank's constant in cgs unit
_PLANK_H = 6.62606896E-27
#Gravitational acc. in cgs unit
_GRAVITY = 981.0

class ResolutionCalculator(object):
    """
    compute resolution in 2D
    """
    def __init__(self):
        
        # wavelength
        self.wave = Neutron()
        # sample
        self.sample = Sample()
        # aperture
        self.aperture = Aperture()
        # detector
        self.detector = Detector()
        # 2d image of the resolution
        self.image = []
        # resolutions
        self.sigma_1 = 0
        self.sigma_2 = 0
        self.sigma_1d = 0
        # q min and max
        self.qx_min = -0.3
        self.qx_max = 0.3
        self.qy_min = -0.3
        self.qy_max = 0.3
        # q min and max of the detector
        self.detector_qx_min = -0.3
        self.detector_qx_max = 0.3
        self.detector_qy_min = -0.3
        self.detector_qy_max = 0.3
        # plots
        self.plot = None
        # instrumental params defaults
        self.mass = 0
        self.intensity = 0
        self.wavelength = 0
        self.wavelength_spread = 0
        self.source_aperture_size = []
        self.source2sample_distance = []
        self.sample2sample_distance = []
        self.sample_aperture_size = []
        self.sample2detector_distance = []
        self.detector_pix_size = []
        self.detector_size = []
        # get all the values of the instrumental parameters
        self.get_all_instrument_params()
        
    def compute_and_plot(self, qx_value, qy_value, qx_min, qx_max, 
                          qy_min, qy_max, coord = 'polar'):
        """
        Compute the resolution
        : qx_value: x component of q
        : qy_value: y component of q
        """
        # compute 2d resolution
        _, _, sigma_1, sigma_2 = self.compute(qx_value, qy_value, coord)
        # make image
        image = self.get_image(qx_value, qy_value, sigma_1, sigma_2, 
                               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'):
        """
        Compute the Q resoltuion in || and + direction of 2D
        : qx_value: x component of q
        : qy_value: y component of q
        """
        # make sure to update all the variables need.
        self.get_all_instrument_params()
        # wavelength
        lamb = self.wavelength
   
        if lamb == 0:
            msg = "Can't compute the resolution: the wavelength is zero..." 
            raise RuntimeError, msg
        # wavelength spread
        lamb_spread = self.wavelength_spread
        # Find polar values
        qr_value, phi = self._get_polar_value(qx_value, qy_value)
        # vacuum wave transfer
        knot = 2*pi/lamb
        # scattering angle theta; always true for plane detector
        # aligned vertically to the ko direction
        if qr_value > knot:
            theta = pi/2
        else:
            theta = asin(qr_value/knot)
        # source aperture size
        rone = self.source_aperture_size
        # sample aperture size
        rtwo = self.sample_aperture_size
        # detector pixel size
        rthree = self.detector_pix_size
        # source to sample(aperture) distance
        l_ssa = self.source2sample_distance[0]
        # sample(aperture) to detector distance
        l_sad = self.sample2detector_distance[0]
        # sample (aperture) to sample distance
        l_sas = self.sample2sample_distance[0]
        # source to sample distance
        l_one = l_ssa + l_sas
        # sample to detector distance
        l_two = l_sad - l_sas
        
        # Sample offset correction for l_one and Lp on variance calculation
        l1_cor = (l_ssa * l_two) / (l_sas + l_two)
        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
        #Lp = l_one*l_two/(l_one+l_two)
        # default polar coordinate
        comp1 = 'radial'
        comp2 = 'phi'
        # in the case of the cartesian coordinate
        if coord == 'cartesian':
            comp1 = 'x'
            comp2 = 'y'

        # sigma in the radial/x direction
        # for source aperture
        sigma_1  = self.get_variance(rone, l1_cor, phi, comp1)
        # for sample apperture
        sigma_1 += self.get_variance(rtwo, lp_cor, phi, comp1)
        # for detector pix
        sigma_1 += self.get_variance(rthree, l_two, phi, comp1)
        # for gravity term
        sigma_1 +=  self.get_variance_gravity(l_ssa, l_sad, lamb, lamb_spread, 
                             phi, comp1, 'on')  
        # for wavelength spread
        # reserve for 1d calculation
        sigma_wave_1 = self.get_variance_wave(radius, l_two, lamb_spread, 
                                          phi, comp1, 'on')
        # for 1d
        variance_1d_1 = sigma_1/2 +sigma_wave_1
        # normalize
        variance_1d_1 = knot*knot*variance_1d_1/12
        
        # for 2d
        sigma_1 += sigma_wave_1
        # normalize
        sigma_1 = knot*sqrt(sigma_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') 
        # for 1d
        variance_1d_2 = sigma_2/2 +sigma_wave_2
        # normalize
        variance_1d_2 = knot*knot*variance_1d_2/12
        
        # for 2d
        sigma_2 += sigma_wave_2
        # normalize
        sigma_2 =  knot*sqrt(sigma_2/12)

        # set sigmas
        self.sigma_1 = sigma_1
        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'): 
        """
        Get the resolution in polar coordinate ready to plot
        : qx_value: qx_value value
        : qy_value: qy_value value
        : sigma_1: variance in r direction
        : sigma_2: variance in phi direction
        : coord: coordinate system of image, 'polar' or 'cartesian'
        """
        # Get  qx_max and qy_max...
        output = self._get_detector_qxqy_pixels()
        # Set qx_value/qy_value min/max
        #qx_min = self.qx_min
        #qx_max = self.qx_max
        #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..."
        if qx_min < self.qx_min:
            self.qx_min = qx_min
            #raise ValueError, msg
        if qx_max > self.qx_max:
            self.qx_max = qx_max
            #raise ValueError, msg
        if qy_min < self.qy_min:
            self.qy_min = qy_min
            #raise ValueError, msg
        if qy_max > self.qy_max:
            self.qy_max = qy_max
            #raise ValueError, msg

        # Make an empty graph in the detector scale
        dx_size = (self.qx_max - self.qx_min) / (1000 - 1)
        dy_size = (self.qy_max - self.qy_min) / (1000 - 1)
        x_val = numpy.arange(self.qx_min, self.qx_max, dx_size)
        y_val = numpy.arange(self.qy_max, self.qy_min, -dy_size)
        q_1, q_2 = numpy.meshgrid(x_val, y_val)

        # check whether polar or cartesian
        if coord == 'polar':
            q_1, q_2 = self._rotate_z(q_1, q_2, phi)
            qc_1 = qr_value
            qc_2 = 0.0
        else:
            # catesian coordinate
            # qx_center
            qc_1 = qx_value
            # qy_center
            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)
        # Add it if there are more than one inputs.
        if len(self.image) > 0:
            self.image += image
        else:
            self.image = image

        return self.image
    
    def plot_image(self, image):
        """
        Plot image using pyplot
        : image: 2d resolution image 
        
        : return plt: pylab object
        """
        import matplotlib.pyplot as plt

        self.plot = plt
        plt.xlabel('$\\rm{Q}_{x} [A^{-1}]$')
        plt.ylabel('$\\rm{Q}_{y} [A^{-1}]$')
        # Max value of the image
        max = numpy.max(image)
        # Image
        im = plt.imshow(image, 
                extent = [self.qx_min, self.qx_max, self.qy_min, self.qy_max])

        # bilinear interpolation to make it smoother
        im.set_interpolation('bilinear')

        return plt
    
    def reset_image(self):
        """
        Reset image to default (=[])
        """
        self.image = []
        
    def get_variance(self, size = [], distance = 0, phi = 0, comp = 'radial'):
        """
        Get the variance when the slit/pinhole size is given
        : size: list that can be one(diameter for circular) 
                or two components(lengths for rectangular)
        : distance: [z, x] where z along the incident beam, x // qx_value
        : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial'
        
        : return variance: sigma^2
        """   
        # check the length of size (list)
        len_size = len(size)
        
        # define sigma component direction
        if comp == 'radial':
            phi_x = cos(phi)
            phi_y = sin(phi)
        elif comp == 'phi':
            phi_x = sin(phi)
            phi_y = cos(phi)
        elif comp == 'x':
            phi_x = 1
            phi_y = 0
        elif comp == 'y':
            phi_x = 0
            phi_y = 1
        else:
            phi_x = 0
            phi_y = 0
        # calculate each component
        # for pinhole w/ radius = size[0]/2
        if len_size == 1:
            x_comp = (0.5 * size[0]) * sqrt(3)
            y_comp = 0
        # for rectangular slit
        elif len_size == 2:
            x_comp = size[0] * phi_x 
            y_comp = size[1] * phi_y
        # otherwise
        else:
            raise ValueError, " Improper input..."
        # get them squared  
        sigma  = x_comp * x_comp 
        sigma += y_comp * y_comp
        # normalize by distance
        sigma /= (distance * distance)

        return sigma

    def get_variance_wave(self, radius, distance, spread, phi, 
                          comp = 'radial', switch = 'on'):
        """
        Get the variance when the wavelength spread is given
        
        : radius: the radial distance from the beam center to the pix of q
        : distance: sample to detector distance
        : spread: wavelength spread (ratio)
        : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial'
        
        : return variance: sigma^2
        """
        if switch.lower() == 'off':
            return 0
        # check the singular point
        if distance == 0 or comp == 'phi':
            return 0
        else:
            # calculate sigma^2
            sigma = 2 * pow(radius/distance*spread, 2)
            if comp == 'x':
                sigma *= (cos(phi)*cos(phi))
            elif comp == 'y':
                sigma *= (sin(phi)*sin(phi))
            else:
                sigma *= 1          
                
            return sigma

    def get_variance_gravity(self, s_distance, d_distance, wavelength, spread, 
                             phi, comp = 'radial', switch = 'on'):
        """
        Get the variance from gravity when the wavelength spread is given
        
        : s_distance: source to sample distance
        : d_distance: sample to detector distance
        : wavelength: wavelength
        : spread: wavelength spread (ratio)
        : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial'
        
        : return variance: sigma^2
        """
        if switch.lower() == 'off':
            return 0
        if self.mass == 0.0:
            return 0
        # check the singular point
        if d_distance == 0 or comp == 'x':
            return 0
        else:
            # neutron mass in cgs unit
            self.mass = self.get_neutron_mass()
            # plank constant in cgs unit
            h_constant = _PLANK_H
            # gravity in cgs unit
            gravy = _GRAVITY
            # m/h
            m_over_h = self.mass /h_constant
            # A value
            a_value = d_distance * (s_distance + d_distance)
            a_value *= pow(m_over_h / 2, 2)
            a_value *= gravy
            # unit correction (1/cm to 1/A) for A and d_distance below
            a_value *= 1.0E-16
            
            # calculate sigma^2
            sigma = pow(a_value / d_distance, 2)
            sigma *= pow(wavelength, 4)
            sigma *= 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))
            
            return sigma
        
    def get_intensity(self):
        """
        Get intensity
        """
        return self.wave.intensity

    def get_wavelength(self):
        """
        Get wavelength
        """
        return self.wave.wavelength
    
    def get_wavelength_spread(self):
        """
        Get wavelength spread
        """
        return self.wave.wavelength_spread
    
    def get_neutron_mass(self):
        """
        Get Neutron mass
        """
        return self.wave.mass
    
    def get_source_aperture_size(self):
        """
        Get source aperture size
        """
        return self.aperture.source_size
    
    def get_sample_aperture_size(self):
        """
        Get sample aperture size
        """
        return self.aperture.sample_size
    
    def get_detector_pix_size(self):
        """
        Get detector pixel size
        """
        return self.detector.pix_size
    
    def get_detector_size(self):
        """
        Get detector size
        """
        return self.detector.size   
     
    def get_source2sample_distance(self):
        """
        Get detector source2sample_distance
        """
        return self.aperture.sample_distance
    
    def get_sample2sample_distance(self):
        """
        Get detector sampleslitsample_distance
        """
        return self.sample.distance
    
    def get_sample2detector_distance(self):
        """
        Get detector sample2detector_distance
        """
        return self.detector.distance
    
    def set_intensity(self, intensity):
        """
        Set intensity
        """
        self.wave.set_intensity(intensity)
   
    def set_wavelength(self, wavelength):
        """
        Set wavelength
        """
        self.wave.set_wavelength(wavelength)
    
    def set_wavelength_spread(self, wavelength_spread):
        """
        Set wavelength spread
        """
        self.wave.set_wavelength_spread(wavelength_spread)
    
    def set_source_aperture_size(self, size):
        """
        Set source aperture size
        
        : param size: [dia_value] or [x_value, y_value]
        """        
        if len(size) < 1 or len(size) > 2:
            raise RuntimeError, "The length of the size must be one or two."
        self.aperture.set_source_size(size)
        
    def set_neutron_mass(self, mass):
        """
        Set Neutron mass
        """
        self.wave.set_mass(mass)
        
    def set_sample_aperture_size(self, size):
        """
        Set sample aperture size
        
        : param size: [dia_value] or [xheight_value, yheight_value]
        """
        if len(size) < 1 or len(size) > 2:
            raise RuntimeError, "The length of the size must be one or two."
        self.aperture.set_sample_size(size)
    
    def set_detector_pix_size(self, size):
        """
        Set detector pixel size
        """
        self.detector.set_pix_size(size)
        
    def set_detector_size(self, size):
        """
        Set detector size in number of pixels
        : param size: [pixel_nums] or [x_pix_num, yx_pix_num]
        """
        self.detector.set_size(size)
        
    def set_source2sample_distance(self, distance):
        """
        Set detector source2sample_distance
        
        : param distance: [distance, x_offset]
        """
        if len(distance) < 1 or len(distance) > 2:
            raise RuntimeError, "The length of the size must be one or two."
        self.aperture.set_sample_distance(distance)

    def set_sample2sample_distance(self, distance):
        """
        Set detector sample_slit2sample_distance
        
        : param distance: [distance, x_offset]
        """
        if len(distance) < 1 or len(distance) > 2:
            raise RuntimeError, "The length of the size must be one or two."
        self.sample.set_distance(distance)
   
    def set_sample2detector_distance(self, distance):
        """
        Set detector sample2detector_distance
        
        : param distance: [distance, x_offset]
        """
        if len(distance) < 1 or len(distance) > 2:
            raise RuntimeError, "The length of the size must be one or two."
        self.detector.set_distance(distance)
        
    def get_all_instrument_params(self):
        """
        Get all instrumental parameters
        """
        self.intensity = self.get_intensity()
        self.wavelength = self.get_wavelength()
        self.wavelength_spread = self.get_wavelength_spread()
        self.mass = self.get_neutron_mass()
        self.source_aperture_size = self.get_source_aperture_size()
        self.sample_aperture_size = self.get_sample_aperture_size()
        self.detector_pix_size = self.get_detector_pix_size()
        self.detector_size = self.get_detector_size()
        self.source2sample_distance = self.get_source2sample_distance()
        self.sample2sample_distance = self.get_sample2sample_distance()
        self.sample2detector_distance = self.get_sample2detector_distance()


          
    def _rotate_z(self, x_value, y_value, theta= 0.0):
        """
        Rotate x-y cordinate around z-axis by theta
        : x_value: numpy array of x values
        : y_value: numpy array of y values
        : theta: angle to rotate by in rad
        
        :return: x_prime, y-prime
        """        
        # rotate by theta
        x_prime = x_value * cos(theta) + y_value * sin(theta)
        y_prime =  -x_value * sin(theta) + y_value * cos(theta)
    
        return x_prime, y_prime
    
    def _gaussian2d(self, x_val, y_val, x0_val, y0_val, sigma_x, sigma_y):
        """
        Calculate 2D Gaussian distribution
        : 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 x-direction
        : sigma_y: variance in y-direction
        
        : return: gaussian (value)
        """
        # call gaussian1d 
        gaussian  = self._gaussian1d(x_val, x0_val, sigma_x)
        gaussian *= self._gaussian1d(y_val, y0_val, sigma_y)
 
        # normalizing factor correction
        if sigma_x != 0 and sigma_y != 0:
            gaussian *= sqrt(2 * pi)
        return gaussian

    def _gaussian1d(self, value, mean, sigma):
        """
        Calculate 1D Gaussian distribution
        : value: value
        : mean: mean value
        : sigma: variance 
        
        : return: gaussian (value)
        """
        # default
        gaussian = 1.0
        if sigma != 0:
            # get exponent
            nu_value = (value - mean) / sigma
            nu_value *= nu_value
            nu_value *= -0.5
            gaussian *= numpy.exp(nu_value)
            gaussian /= sigma
            # normalize
            gaussian /= sqrt(2 * pi)
            
        return gaussian
    
    def _atan_phi(self, qy_value, qx_value):
        """
        Find the angle phi of q on the detector plane for qx_value, qy_value given
        : qx_value: x component of q
        : qy_value: y component of q
        
        : return phi: the azimuthal angle of q on x-y plane
        """
        # default
        phi = 0
        # ToDo: This is misterious - sign???
        #qy_value = -qy_value
        # Take care of the singular point
        if qx_value == 0:
            if qy_value > 0:
                phi = pi / 2
            elif qy_value < 0:
                phi = -pi / 2
            else:
                phi = 0
        else:
            # the angle
            phi = atan2(qy_value, qx_value)

        return phi

    def _get_detector_qxqy_pixels(self):
        """
        Get the pixel positions of the detector in the qx_value-qy_value space
        """
        
        # update all param values
        self.get_all_instrument_params()
        
        # wavelength
        wavelength = self.wavelength
        # Gavity correction
        delta_y = self._get_beamcenter_drop() # in cm
        
        # detector_pix size
        detector_pix_size = self.detector_pix_size
        # Square or circular pixel
        if len(detector_pix_size) == 1:
            pix_x_size = detector_pix_size[0]
            pix_y_size = detector_pix_size[0]
        # rectangular pixel pixel
        elif len(detector_pix_size) == 2:
            pix_x_size = detector_pix_size[0]
            pix_y_size = detector_pix_size[1]
        else:
            raise ValueError, " Input value format error..."
        # Sample to detector distance = sample slit to detector 
        # minus sample offset
        sample2detector_distance = self.sample2detector_distance[0] - \
                                    self.sample2sample_distance[0]
        # detector offset in x-direction
        detector_offset = 0
        try:
            detector_offset = self.sample2detector_distance[1]
        except:
            pass
        
        # detector size in [no of pix_x,no of pix_y]
        detector_pix_nums_x = self.detector_size[0]
        
        # get pix_y if it exists, otherwse take it from [0]
        try:
            detector_pix_nums_y = self.detector_size[1]
        except:
            detector_pix_nums_y = self.detector_size[0]
        
        # detector offset in pix number
        offset_x = detector_offset / pix_x_size
        offset_y = delta_y / pix_y_size
        
        # beam center position in pix number (start from 0)
        center_x, center_y = self._get_beamcenter_position(detector_pix_nums_x,
                                    detector_pix_nums_y, offset_x, offset_y)
        # distance [cm] from the beam center on detector plane
        detector_ind_x = numpy.arange(detector_pix_nums_x)
        detector_ind_y = numpy.arange(detector_pix_nums_y)

        # shif 0.5 pixel so that pix position is at the center of the pixel
        detector_ind_x = detector_ind_x + 0.5
        detector_ind_y = detector_ind_y + 0.5

        # the relative postion from the beam center
        detector_ind_x = detector_ind_x - center_x
        detector_ind_y = detector_ind_y - center_y
        
        # unit correction in cm
        detector_ind_x = detector_ind_x * pix_x_size
        detector_ind_y = detector_ind_y * pix_y_size
        
        qx_value = numpy.zeros(len(detector_ind_x))
        qy_value = numpy.zeros(len(detector_ind_y))
        i = 0

        for indx in detector_ind_x:
            qx_value[i] = self._get_qx(indx, sample2detector_distance, wavelength)
            i += 1
        i = 0
        for indy in detector_ind_y:
            qy_value[i] = self._get_qx(indy, sample2detector_distance, wavelength)
            i += 1
            
        # qx_value and qy_value values in array
        qx_value = qx_value.repeat(detector_pix_nums_y) 
        qx_value = qx_value.reshape(detector_pix_nums_x, detector_pix_nums_y)
        qy_value = qy_value.repeat(detector_pix_nums_x)  
        qy_value = qy_value.reshape(detector_pix_nums_y, detector_pix_nums_x)
        qy_value = qy_value.transpose()

        # p min and max values among the center of pixels
        self.qx_min = numpy.min(qx_value)
        self.qx_max = numpy.max(qx_value)
        self.qy_min = numpy.min(qy_value)
        self.qy_max = numpy.max(qy_value)
                
        # Appr. min and max values of the detector display limits 
        # i.e., edges of the last pixels.
        self.qy_min += self._get_qx(-0.5 * pix_y_size, 
                                sample2detector_distance, wavelength)
        self.qy_max += self._get_qx(0.5 * pix_y_size, 
                                sample2detector_distance, wavelength)
        #if self.qx_min == self.qx_max:
        self.qx_min += self._get_qx(-0.5 * pix_x_size, 
                                sample2detector_distance, wavelength)
        self.qx_max += self._get_qx(0.5 * pix_x_size, 
                                    sample2detector_distance, wavelength)
        # min and max values of detecter
        self.detector_qx_min = self.qx_min
        self.detector_qx_max = self.qx_max
        self.detector_qy_min = self.qy_min
        self.detector_qy_max = self.qy_max

        # try to set it as a Data2D otherwise pass (not required for now)
        try:
            from DataLoader.data_info import Data2D
            output = Data2D()
            inten = numpy.zeros_like(qx_value)
            output.data     = inten
            output.qx_data  = qx_value
            output.qy_data  = qy_value            
        except:
            pass
        
        return output#qx_value,qy_value
        
    def _get_qx(self, dx_size, det_dist, wavelength):
        """
        :param dx_size: x-distance from beam center [cm]
        :param det_dist: sample to detector distance [cm]
        
        :return: q-value at the given position
        """
        # Distance from beam center in the plane of detector
        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)
        return qx_value  
    
    def _get_polar_value(self, qx_value, qy_value):
        """
        Find qr_value and phi from qx_value and qy_value values
        
        : return qr_value, phi
        """
        # find |q| on detector plane
        qr_value = sqrt(qx_value*qx_value + qy_value*qy_value)
        # find angle phi
        phi = self._atan_phi(qy_value, qx_value)
        
        return qr_value, phi
    
    def _get_beamcenter_position(self, num_x, num_y, offset_x, offset_y):
        """
        :param num_x: number of pixel in x-direction
        :param num_y: number of pixel in y-direction
        :param offset: detector offset in x-direction in pix number
        
        :return: pix number; pos_x, pos_y in pix index
        """
        # beam center position
        pos_x = num_x / 2
        pos_y = num_y / 2

        # correction for offset
        pos_x += offset_x
        # correction for gravity that is always negative
        pos_y -= offset_y

        return pos_x, pos_y     

    def _get_beamcenter_drop(self):
        """
        Get the beam center drop (delta y) in y diection due to gravity
        
        :return delta y: the beam center drop in cm
        """
        # Check if mass == 0 (X-ray). 
        if self.mass == 0:
            return 0
        # Covert unit from A to cm
        unit_cm = 1e-08
        # Velocity of neutron in horizontal direction (~ actual velocity)
        velocity = _PLANK_H / (self.mass * self.wavelength * unit_cm)
        # Compute delta y
        delta_y = 0.5
        delta_y *= _GRAVITY
        delta_y *= self.sample2detector_distance[0] 
        delta_y *= (self.source2sample_distance[0] + self.sample2detector_distance[0])
        delta_y /= (velocity * velocity)

        return delta_y