resolution_calculator.py @ bc202cb

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since bc202cb was b6627d9, checked in by Doucet, Mathieu <doucetm@…>, 10 years ago
pylint fixes
Property mode set to `100644`
File size: 38.4 KB

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

Note: See TracBrowser for help on using the repository browser.

SasView

source: sasview/src/sas/calculator/resolution_calculator.py @ bc202cb

Download in other formats: