source: sasview/sanscalculator/src/sans/calculator/resolution_calculator.py @ 440f020

"""
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 TOF as Neutron
from instrument import Aperture
# import math stuffs
from math import pi
from math import sqrt
import math 
import scipy
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 = []
        self.image_lam = []
        # 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
        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
        # possible max qrange
        self.qxmin_limit = 0
        self.qxmax_limit = 0
        self.qymin_limit = 0
        self.qymax_limit = 0
        # set sigmas
        self.sigma_1 = 0
        self.sigma_lamd = 0
        self.sigma_2 = 0
        self.sigma_1d = 0
        # 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.intensity = self.get_intensity()
        #self.wavelength = self.get_wavelength()
        #self.wavelength_spread = self.get_wavelength_spread()
        self.get_all_instrument_params()
        # max q range for all lambdas
        self.qxrange = []
        self.qyrange = []
        
    def compute_and_plot(self, qx_value, qy_value, qx_min, qx_max, 
                          qy_min, qy_max, coord = 'cartesian'):
        """
        Compute the resolution
        : qx_value: x component of q
        : qy_value: y component of q
        """
        # make sure to update all the variables need.
        # except lambda, dlambda, and intensity
        self.get_all_instrument_params()
        # wavelength etc.
        lamda_list, dlamb_list = self.get_wave_list()
        intens_list = []#self.get_intensity_list()

        sig1_list = []
        sig2_list = []
        sigr_list = []
        sigma1d_list = []
        num_lamda = len(lamda_list)
        for num in range(num_lamda):
            lam = lamda_list[num]
            # wavelength spread
            dlam = dlamb_list[num]
            intens = self.setup_tof(lam, dlam)
            intens_list.append(intens)
            # save q min max
            #qx_min = 0
            #qx_max = 0
            #qy_min = 0
            #qy_max = 0
            # compute 2d resolution
            _, _, sigma_1, sigma_2, sigma_r, sigma1d = \
                            self.compute(lam, dlam, qx_value, qy_value, coord)
            # make image
            image = self.get_image(qx_value, qy_value, sigma_1, sigma_2, 
                            sigma_r, qx_min, qx_max, qy_min, qy_max,
                            coord, False)
            
            # Non tof mode to be speed up 
            #if num_lamda < 2:
            #    return self.plot_image(image)
            
            if qx_min > self.qx_min:
                qx_min = self.qx_min
            if qx_max < self.qx_max:
                qx_max = self.qx_max
            if qy_min > self.qy_min:
                qy_min = self.qy_min
            if qy_max < self.qy_max:
                qy_max = self.qy_max
                
            # set max qranges
            self.qxrange = [qx_min, qx_max]
            self.qyrange = [qy_min, qy_max]
            
            sig1_list.append(sigma_1)
            sig2_list.append(sigma_2)
            sigr_list.append(sigma_r)
            sigma1d_list.append(sigma1d)
        # redraw image in global 2d q-space.   
        self.image_lam = []
        total_intensity = 0
        sigma_1 = 0
        sigma_r = 0
        sigma_2 = 0
        sigma1d = 0
        for ind in range(num_lamda):
            lam = lamda_list[ind]
            dlam = dlamb_list[ind]
            intens = self.setup_tof(lam, dlam)
            out = self.get_image(qx_value, qy_value, sig1_list[ind], 
                                   sig2_list[ind], sigr_list[ind], 
                                   qx_min, qx_max, qy_min, qy_max, coord)
            # this is the case of q being outside the detector
            #if numpy.all(out==0.0):
            #    continue
            image = out
            # set variance as sigmas
            sigma_1 += sig1_list[ind] *  sig1_list[ind] * self.intensity
            sigma_r += sigr_list[ind] * sigr_list[ind] * self.intensity
            sigma_2 += sig2_list[ind] * sig2_list[ind] * self.intensity
            sigma1d += sigma1d_list[ind] * sigma1d_list[ind] * self.intensity
            total_intensity += self.intensity
        
        if total_intensity != 0:
            # average variance
            image_out = image / total_intensity
            sigma_1 = sigma_1 / total_intensity
            sigma_r = sigma_r / total_intensity
            sigma_2 = sigma_2 / total_intensity
            sigma1d = sigma1d / total_intensity
            # set sigmas
            self.sigma_1 = sqrt(sigma_1)
            self.sigma_lamd = sqrt(sigma_r)
            self.sigma_2 = sqrt(sigma_2)
            self.sigma_1d = sqrt(sigma1d)
            # rescale 
            max_im_val = 1 #image_out.max()
            if max_im_val > 0:
                image_out /= max_im_val
        else:
            image_out = image * 0.0
            # Don't calculate sigmas nor set self.sigmas!
            sigma_1 = 0
            sigma_r = 0
            sigma_2 = 0
            sigma1d = 0
        if len(self.image) > 0:
            self.image += image_out
        else:
            self.image = image_out
        
        # plot image
        return self.plot_image(self.image) 
    
    def setup_tof(self, wavelength, wavelength_spread):
        """
        Setup all parameters in instrument
        
        : param ind: index of lambda, etc
        """

        # set wave.wavelength
        self.set_wavelength(wavelength)
        self.set_wavelength_spread(wavelength_spread)
        self.intensity = self.wave.get_intensity()
        
        if wavelength == 0:
            msg = "Can't compute the resolution: the wavelength is zero..." 
            raise RuntimeError, msg
        return self.intensity
        
    def compute(self, wavelength, wavelength_spread, qx_value, qy_value, 
                coord = 'cartesian'):
        """
        Compute the Q resoltuion in || and + direction of 2D
        : qx_value: x component of q
        : qy_value: y component of q
        """
        coord = 'cartesian'
        lamb = wavelength
        lamb_spread = 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 = math.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 = math.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, 'radial', '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_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, 'phi', '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 =  knot*sqrt(sigma_2/12)
        #sigma_2 += sigma_wave_2
        # normalize
        sigma_2 =  knot * sqrt(sigma_2 / 12)
        sigma1d = sqrt(variance_1d_1 + variance_1d_2)
        # set sigmas
        self.sigma_1 = sigma_1
        self.sigma_lamd = sigma_r
        self.sigma_2 = sigma_2
        self.sigma_1d = sigma1d
        return qr_value, phi, sigma_1, sigma_2, sigma_r, sigma1d
    
    def _within_detector_range(self,qx_value, qy_value):
        """
        check if qvalues are within detector range
        """
        # detector range
        detector_qx_min = self.detector_qx_min
        detector_qx_max = self.detector_qx_max
        detector_qy_min = self.detector_qy_min
        detector_qy_max = self.detector_qy_max
        if self.qxmin_limit > detector_qx_min:
            self.qxmin_limit = detector_qx_min
        if self.qxmax_limit < detector_qx_max:
            self.qxmax_limit = detector_qx_max
        if self.qymin_limit > detector_qy_min:
            self.qymin_limit = detector_qy_min
        if self.qymax_limit < detector_qy_max:
            self.qymax_limit = detector_qy_max
        if qx_value < detector_qx_min or qx_value > detector_qx_max:
            return False
        if qy_value < detector_qy_min or qy_value > detector_qy_max:
            return False
        return True
    
    def get_image(self, qx_value, qy_value, sigma_1, sigma_2, sigma_r,
                  qx_min, qx_max, qy_min, qy_max, 
                  coord = 'cartesian', full_cal=True): 
        """
        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()
       
        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
        if not full_cal:
            return None
 
        # 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)
        #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
            # 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, sigma_r)
        # out side of detector
        if not self._within_detector_range(qx_value, qy_value):
            image *= 0.0
            self.intensity = 0.0
            #return self.image

        # Add it if there are more than one inputs.
        if len(self.image_lam) > 0:
            self.image_lam += image * self.intensity
        else:
            self.image_lam = image * self.intensity
        
        return self.image_lam
    
    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)
        qx_min, qx_max, qy_min, qy_max = self.get_detector_qrange()

        # Image
        im = plt.imshow(image, 
                extent = [qx_min, qx_max, qy_min, 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 = math.cos(phi)
            phi_y = math.sin(phi)
        elif comp == 'phi':
            phi_x = math.sin(phi)
            phi_y = math.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 * math.pow(radius/distance*spread, 2)
            if comp == 'x':
                sigma *= (math.cos(phi)*math.cos(phi))
            elif comp == 'y':
                sigma *= (math.sin(phi)*math.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 *= math.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 = 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 *= (math.sin(phi) * math.sin(phi))
            #elif comp == 'phi':
            #    sigma *= (math.cos(phi) * math.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_spectrum(self):
        """
        Get spectrum
        """
        return self.wave.spectrum   
    
    def get_default_spectrum(self):
        """
        Get default_spectrum
        """
        return self.wave.get_default_spectrum()    
    
    def get_spectrum(self):
        """
        Get _spectrum
        """
        return self.wave.get_spectrum() 
         
    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_wave(self, wavelength):
        """
        Set wavelength list or wavelength
        """
        if wavelength.__class__.__name__ == 'list':
            self.wave.set_wave_list(wavelength)
        elif wavelength.__class__.__name__ == 'float':
            self.wave.set_wave_list([wavelength])
            #self.set_wavelength(wavelength)
        else:
            raise
    
    def set_wave_spread(self, wavelength_spread):
        """
        Set wavelength spread  or wavelength spread
        """
        if wavelength_spread.__class__.__name__ == 'list':
            self.wave.set_wave_spread_list(wavelength_spread)
        elif wavelength_spread.__class__.__name__ == 'float':
            self.wave.set_wave_spread_list([wavelength_spread])
            #self.set_wavelength_spread(wavelength_spread)
        else:
            raise
        
    def set_wavelength(self, wavelength):
        """
        Set wavelength
        """
        self.wavelength = wavelength
        self.wave.set_wavelength(wavelength)
        
    def set_spectrum(self, spectrum):
        """
        Set spectrum
        """
        self.spectrum = spectrum
        self.wave.set_spectrum(spectrum)  
          
    def set_wavelength_spread(self, wavelength_spread):
        """
        Set wavelength spread
        """
        self.wavelength_spread = wavelength_spread
        self.wave.set_wavelength_spread(wavelength_spread)
        
    def set_wave_list(self, wavelength_list, wavelengthspread_list):
        """
        Set wavelength and its spread list
        """
        self.wave.set_wave_list(wavelength_list, wavelengthspread_list)
    
    def get_wave_list(self):
        """
        Set wavelength spread
        """
        return self.wave.get_wave_list()
    
    def get_intensity_list(self):
        """
        Set wavelength spread
        """
        return self.wave.get_intensity_list()
           
    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)
        self.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.spectrum = self.get_spectrum()
        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 get_detector_qrange(self):
        """
        get max detector q ranges
        
        : return: qx_min, qx_max, qy_min, qy_max tuple
        """
        if len(self.qxrange) != 2 or len(self.qyrange) != 2:
            return None
        qx_min = self.qxrange[0]
        qx_max = self.qxrange[1]
        qy_min = self.qyrange[0]
        qy_max = self.qyrange[1]
        
        return qx_min, qx_max, qy_min, qy_max
    
    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 * 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, sigma_r):
        """
        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)
        """
        # 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_sig_x = sqrt(sigma_r * sigma_r / (sigma_x * sigma_x ) + 1)
        new_sig_y = sqrt(sigma_r * sigma_r / (sigma_y * sigma_y ) + 1)
        new_x = x_p * cos_phi / new_sig_x - y_p * sin_phi
        new_x /= sigma_x
        new_y = x_p * sin_phi / new_sig_y + 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)
        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
        """
        phi = math.atan2(qy_value, qx_value)
        return phi
        # 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 = math.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.wave.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 sans.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  = math.atan(plane_dist / det_dist)
        qx_value     = (2.0 * pi / wavelength) * math.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.wave.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