"""
This object is a small tool to allow user to quickly
determine the variance in q from the
instrumental parameters.
"""
import sys
from math import pi, sqrt
import math
import logging
import numpy as np
from .instrument import Sample
from .instrument import Detector
from .instrument import TOF as Neutron
from .instrument import Aperture
logger = logging.getLogger(__name__)
#Plank's constant in cgs unit
_PLANK_H = 6.62606896E-27
#Gravitational acc. in cgs unit
_GRAVITY = 981.0
[docs]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_lamd = 0
# x-dir (no lamda)
self.sigma_1 = 0
#y-dir (no lamda)
self.sigma_2 = 0
# 1D total
self.sigma_1d = 0
self.gravity_phi = None
# 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
# 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 = []
self.get_all_instrument_params()
# max q range for all lambdas
self.qxrange = []
self.qyrange = []
[docs] 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 = []
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)
# cehck if tof
if num_lamda > 1:
tof = True
else:
tof = False
# compute 2d resolution
_, _, sigma_1, sigma_2, sigma_r, sigma1d = \
self.compute(lam, dlam, qx_value, qy_value, coord, tof)
# 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)
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
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)
[docs] 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
[docs] def compute(self, wavelength, wavelength_spread, qx_value, qy_value,
coord='cartesian', tof=False):
"""
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
# the shape of wavelength distribution
if tof:
# rectangular
tof_factor = 2
else:
# triangular
tof_factor = 1
# 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 for 1d
sigma_1grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb,
lamb_spread, phi, comp1, 'on') / tof_factor
# for wavelength spread
# reserve for 1d calculation
A_value = self._cal_A_value(lamb, l_ssa, l_sad)
sigma_wave_1, sigma_wave_1_1d = self.get_variance_wave(A_value,
radius, l_two, lamb_spread,
phi, 'radial', 'on')
sigma_wave_1 /= tof_factor
sigma_wave_1_1d /= tof_factor
# for 1d
variance_1d_1 = (sigma_1 + sigma_1grav1d) / 2 + sigma_wave_1_1d
# 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 / (tof_factor *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 for 1d
sigma_2grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb,
lamb_spread, phi, comp2, 'on') / tof_factor
# for wavelength spread
# reserve for 1d calculation
sigma_wave_2, sigma_wave_2_1d = self.get_variance_wave(A_value,
radius, l_two, lamb_spread,
phi, 'phi', 'on')
sigma_wave_2 /= tof_factor
sigma_wave_2_1d /= tof_factor
# for 1d
variance_1d_2 = (sigma_2 + sigma_2grav1d) / 2 + sigma_wave_2_1d
# 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
[docs] 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...
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
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 = np.arange(self.qx_min, self.qx_max, dx_size)
y_val = np.arange(self.qy_max, self.qy_min, -dy_size)
q_1, q_2 = np.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
[docs] 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
[docs] def reset_image(self):
"""
Reset image to default (=[])
"""
self.image = []
[docs] 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
[docs] def get_variance_wave(self, A_value, 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 for 2d, sigma^2 for 1d [tuple]
"""
if switch.lower() == 'off':
return 0, 0
# check the singular point
if distance == 0 or comp == 'phi':
return 0, 0
else:
# calculate sigma^2 for 1d
sigma1d = 2 * math.pow(radius/distance*spread, 2)
if comp == 'x':
sigma1d *= (math.cos(phi)*math.cos(phi))
elif comp == 'y':
sigma1d *= (math.sin(phi)*math.sin(phi))
else:
sigma1d *= 1
# sigma^2 for 2d
# shift the coordinate due to the gravitational shift
rad_x = radius * math.cos(phi)
rad_y = A_value - radius * math.sin(phi)
radius = math.sqrt(rad_x * rad_x + rad_y * rad_y)
# new phi
phi = math.atan2(-rad_y, rad_x)
self.gravity_phi = phi
# 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, sigma1d
[docs] 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:
a_value = self._cal_A_value(None, s_distance, d_distance)
# calculate sigma^2
sigma = math.pow(a_value / d_distance, 2)
sigma *= math.pow(wavelength, 4)
sigma *= math.pow(spread, 2)
sigma *= 8
return sigma
def _cal_A_value(self, lamda, s_distance, d_distance):
"""
Calculate A value for gravity
: s_distance: source to sample distance
: d_distance: sample to detector distance
"""
# 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
# if lamda is give (broad meanning of A) return 2* lamda^2 * A
if lamda is not None:
a_value *= (4 * lamda * lamda)
return a_value
[docs] def get_intensity(self):
"""
Get intensity
"""
return self.wave.intensity
[docs] def get_wavelength(self):
"""
Get wavelength
"""
return self.wave.wavelength
[docs] def get_default_spectrum(self):
"""
Get default_spectrum
"""
return self.wave.get_default_spectrum()
[docs] def get_spectrum(self):
"""
Get _spectrum
"""
return self.wave.get_spectrum()
[docs] def get_wavelength_spread(self):
"""
Get wavelength spread
"""
return self.wave.wavelength_spread
[docs] def get_neutron_mass(self):
"""
Get Neutron mass
"""
return self.wave.mass
[docs] def get_source_aperture_size(self):
"""
Get source aperture size
"""
return self.aperture.source_size
[docs] def get_sample_aperture_size(self):
"""
Get sample aperture size
"""
return self.aperture.sample_size
[docs] def get_detector_pix_size(self):
"""
Get detector pixel size
"""
return self.detector.pix_size
[docs] def get_detector_size(self):
"""
Get detector size
"""
return self.detector.size
[docs] def get_source2sample_distance(self):
"""
Get detector source2sample_distance
"""
return self.aperture.sample_distance
[docs] def get_sample2sample_distance(self):
"""
Get detector sampleslitsample_distance
"""
return self.sample.distance
[docs] def get_sample2detector_distance(self):
"""
Get detector sample2detector_distance
"""
return self.detector.distance
[docs] def set_intensity(self, intensity):
"""
Set intensity
"""
self.wave.set_intensity(intensity)
[docs] 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 TypeError("invalid wavlength---should be list or float")
[docs] 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])
else:
raise TypeError("invalid wavelength spread---should be list or float")
[docs] def set_wavelength(self, wavelength):
"""
Set wavelength
"""
self.wavelength = wavelength
self.wave.set_wavelength(wavelength)
[docs] def set_spectrum(self, spectrum):
"""
Set spectrum
"""
self.spectrum = spectrum
self.wave.set_spectrum(spectrum)
[docs] def set_wavelength_spread(self, wavelength_spread):
"""
Set wavelength spread
"""
self.wavelength_spread = wavelength_spread
self.wave.set_wavelength_spread(wavelength_spread)
[docs] def set_wave_list(self, wavelength_list, wavelengthspread_list):
"""
Set wavelength and its spread list
"""
self.wave.set_wave_list(wavelength_list)
self.wave.set_wave_spread_list(wavelengthspread_list)
[docs] def get_wave_list(self):
"""
Set wavelength spread
"""
return self.wave.get_wave_list()
[docs] def get_intensity_list(self):
"""
Set wavelength spread
"""
return self.wave.get_intensity_list()
[docs] 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)
[docs] def set_neutron_mass(self, mass):
"""
Set Neutron mass
"""
self.wave.set_mass(mass)
self.mass = mass
[docs] 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)
[docs] def set_detector_pix_size(self, size):
"""
Set detector pixel size
"""
self.detector.set_pix_size(size)
[docs] 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)
[docs] 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)
[docs] 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)
[docs] 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)
[docs] def get_all_instrument_params(self):
"""
Get all instrumental parameters
"""
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()
[docs] 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 = np.arctan2(y_val, x_val)
# phi correction due to the gravity shift (in phi)
phi_0 = math.atan2(y0_val, x0_val)
phi_i = phi_i - phi_0 + self.gravity_phi
sin_phi = np.sin(self.gravity_phi)
cos_phi = np.cos(self.gravity_phi)
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 = np.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 *= np.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
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 Exception as exc:
logger.error(exc)
# 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 = np.arange(detector_pix_nums_x)
detector_ind_y = np.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 = np.zeros(len(detector_ind_x))
qy_value = np.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 = np.min(qx_value)
self.qx_max = np.max(qx_value)
self.qy_min = np.min(qy_value)
self.qy_max = np.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 sas.sascalc.dataloader.data_info import Data2D
output = Data2D()
inten = np.zeros_like(qx_value)
output.data = inten
output.qx_data = qx_value
output.qy_data = qy_value
except Exception as exc:
logger.error(exc)
return output
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 = np.arctan(plane_dist / det_dist)
qx_value = (2.0 * pi / wavelength) * np.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
sampletodetector = self.sample2detector_distance[0] - \
self.sample2sample_distance[0]
delta_y *= sampletodetector
delta_y *= (self.source2sample_distance[0] + self.sample2detector_distance[0])
delta_y /= (velocity * velocity)
return delta_y