Changeset 8a621ac in sasview for sanscalculator
- Timestamp:
- Apr 27, 2012 9:22:40 AM (13 years ago)
- Branches:
- master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
- Children:
- 10bfeb3
- Parents:
- 34f3ad0
- Location:
- sanscalculator/src/sans/calculator
- Files:
-
- 4 edited
-
instrument.py (modified) (22 diffs)
-
kiessig_calculator.py (modified) (6 diffs)
-
resolution_calculator.py (modified) (53 diffs)
-
slit_length_calculator.py (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
sanscalculator/src/sans/calculator/instrument.py
rf93b6a1e r8a621ac 1 1 """ 2 This module is a small tool to allow user to 3 control instrumental parameters 2 This module is a small tool to allow user to 3 control instrumental parameters 4 4 """ 5 5 import numpy … … 7 7 # defaults in cgs unit 8 8 _SAMPLE_A_SIZE = [1.27] 9 _SOURCE_A_SIZE = [3.81] 10 _SAMPLE_DISTANCE = [1627, 0] 11 _SAMPLE_OFFSET = [0, 0] 12 _SAMPLE_SIZE = [2.54] 13 _SAMPLE_THICKNESS = 0.2 14 _D_DISTANCE = [1000, 0] 9 _SOURCE_A_SIZE = [3.81] 10 _SAMPLE_DISTANCE = [1627, 0] 11 _SAMPLE_OFFSET = [0, 0] 12 _SAMPLE_SIZE = [2.54] 13 _SAMPLE_THICKNESS = 0.2 14 _D_DISTANCE = [1000, 0] 15 15 _D_SIZE = [128, 128] 16 _D_PIX_SIZE = [0.5, 0.5] 16 _D_PIX_SIZE = [0.5, 0.5] 17 17 18 18 _MIN = 0.0 … … 21 21 _WAVE_LENGTH = 6.0 22 22 _WAVE_SPREAD = 0.125 23 _MASS = 1.67492729E-24 #[gr] 24 _LAMBDA_ARRAY = [[0, 1e+16],[_INTENSITY, _INTENSITY]] 23 _MASS = 1.67492729E-24 # [gr] 24 _LAMBDA_ARRAY = [[0, 1e+16], [_INTENSITY, _INTENSITY]] 25 25 26 26 27 class Aperture(object): … … 36 37 self.sample_distance = _SAMPLE_DISTANCE 37 38 38 def set_source_size(self, size 39 def set_source_size(self, size=[]): 39 40 """ 40 41 Set the source aperture size 41 42 """ 42 43 if len(size) == 0: 43 self.source_size = 0.0 44 self.source_size = 0.0 44 45 else: 45 46 self.source_size = size 46 47 validate(size[0]) 47 def set_sample_size(self, size =[]): 48 49 def set_sample_size(self, size=[]): 48 50 """ 49 51 Set the sample aperture size 50 52 """ 51 53 if len(size) == 0: 52 self.sample_size = 0.0 54 self.sample_size = 0.0 53 55 else: 54 56 self.sample_size = size 55 57 validate(size[0]) 56 58 57 def set_sample_distance(self, distance =[]):59 def set_sample_distance(self, distance=[]): 58 60 """ 59 61 Set the sample aperture distance 60 62 """ 61 63 if len(distance) == 0: 62 self.sample_distance = 0.0 64 self.sample_distance = 0.0 63 65 else: 64 66 self.sample_distance = distance … … 78 80 self.thickness = _SAMPLE_THICKNESS 79 81 80 81 def set_size(self, size =[]): 82 def set_size(self, size=[]): 82 83 """ 83 84 Set the sample size 84 85 """ 85 86 if len(size) == 0: 86 self.sample_size = 0.0 87 self.sample_size = 0.0 87 88 else: 88 89 self.sample_size = size 89 90 validate(size[0]) 90 91 91 def set_thickness(self, thickness =0.0):92 def set_thickness(self, thickness=0.0): 92 93 """ 93 94 Set the sample thickness … … 96 97 validate(thickness) 97 98 98 def set_distance(self, distance =[]):99 def set_distance(self, distance=[]): 99 100 """ 100 101 Set the sample distance 101 102 """ 102 103 if len(distance) == 0: 103 self.distance = 0.0 104 self.distance = 0.0 104 105 else: 105 106 self.distance = distance … … 122 123 123 124 124 def set_size(self, size 125 def set_size(self, size=[]): 125 126 """ 126 127 Set the detector size … … 132 133 validate(size[0]) 133 134 134 def set_pix_size(self, size =[]):135 def set_pix_size(self, size=[]): 135 136 """ 136 137 Set the detector pix_size 137 138 """ 138 139 if len(size) == 0: 139 self.pix_size = 0 140 self.pix_size = 0 140 141 else: 141 142 self.pix_size = size 142 143 validate(size[0]) 143 144 144 def set_distance(self, distance =[]):145 def set_distance(self, distance=[]): 145 146 """ 146 147 Set the detector distance … … 166 167 # wavelength spread (FWHM) 167 168 self.wavelength_spread = _WAVE_SPREAD 168 # wavelength spectrum 169 self.spectrum = self.get_default_spectrum() 169 # wavelength spectrum 170 self.spectrum = self.get_default_spectrum() 170 171 # intensity in counts/sec 171 self.intensity = numpy.interp(self.wavelength, 172 self.intensity = numpy.interp(self.wavelength, 172 173 self.spectrum[0], 173 174 self.spectrum[1], … … 188 189 """ 189 190 self.band = self.spectrum 190 def set_spectrum(self, spectrum): 191 192 def set_spectrum(self, spectrum): 191 193 """ 192 194 Set spectrum … … 194 196 :param spectrum: numpy array 195 197 """ 196 self.spectrum = spectrum 198 self.spectrum = spectrum 197 199 self.setup_spectrum() 198 200 199 201 def setup_spectrum(self): 200 202 """ 201 To set the wavelength spectrum, and intensity, assumes 203 To set the wavelength spectrum, and intensity, assumes 202 204 wavelength is already within the spectrum 203 205 """ 204 206 spectrum = self.spectrum 205 intensity = numpy.interp(self.wavelength, 207 intensity = numpy.interp(self.wavelength, 206 208 spectrum[0], 207 209 spectrum[1], … … 213 215 self.max = max(self.spectrum[0]) 214 216 # set default band 215 self.set_band([self.min, self.max])217 self.set_band([self.min, self.max]) 216 218 217 219 def set_band(self, band=[]): … … 224 226 if min(band) < self.min or\ 225 227 max(band) > self.max: 226 raise 228 raise 227 229 self.band = band 228 229 230 230 def set_intensity(self, intensity =368428):231 def set_intensity(self, intensity=368428): 231 232 """ 232 233 Sets the intensity in counts/sec 233 234 """ 234 self.intensity = intensity 235 validate(intensity) 236 237 def set_wavelength(self, wavelength =_WAVE_LENGTH):235 self.intensity = intensity 236 validate(intensity) 237 238 def set_wavelength(self, wavelength=_WAVE_LENGTH): 238 239 """ 239 240 Sets the wavelength … … 242 243 if wavelength < min(self.band) or\ 243 244 wavelength > max(self.band): 244 raise 245 raise 245 246 self.wavelength = wavelength 246 247 validate(wavelength) 247 self.intensity = numpy.interp(self.wavelength, 248 self.intensity = numpy.interp(self.wavelength, 248 249 self.spectrum[0], 249 250 self.spectrum[1], … … 251 252 0.0) 252 253 253 254 def set_mass(self, mass = _MASS): 254 def set_mass(self, mass=_MASS): 255 255 """ 256 256 Sets the wavelength … … 259 259 validate(mass) 260 260 261 def set_wavelength_spread(self, spread =_WAVE_SPREAD):261 def set_wavelength_spread(self, spread=_WAVE_SPREAD): 262 262 """ 263 263 Sets the wavelength spread … … 302 302 """ 303 303 return self.spectrum 304 304 305 def get_default_spectrum(self): 305 306 """ … … 307 308 """ 308 309 return numpy.array(_LAMBDA_ARRAY) 310 309 311 def get_band(self): 310 312 """ 311 313 To get the wavelength band 312 314 """ 313 return self.band 315 return self.band 314 316 315 317 def plot_spectrum(self): … … 320 322 try: 321 323 import matplotlib.pyplot as plt 322 plt.plot(self.spectrum[0], self.spectrum[1], linewidth = 2, color ='r')323 plt.legend(['Spectrum'], loc ='best')324 plt.plot(self.spectrum[0], self.spectrum[1], linewidth=2, color='r') 325 plt.legend(['Spectrum'], loc='best') 324 326 plt.show() 325 327 except: … … 345 347 get list of the intensity wrt wavelength_list 346 348 """ 347 out = numpy.interp(self.wavelength_list, 349 out = numpy.interp(self.wavelength_list, 348 350 self.spectrum[0], 349 351 self.spectrum[1], … … 375 377 376 378 377 def validate(value =None):379 def validate(value=None): 378 380 """ 379 381 Check if the value is folat > 0.0 … … 389 391 except: 390 392 val = False 391 #if not val:392 # raise ValueError, "Got improper value..." -
sanscalculator/src/sans/calculator/kiessig_calculator.py
rf93b6a1e r8a621ac 1 1 """ 2 2 This module is a small tool to allow user to quickly 3 determine the size value in real space from the 3 determine the size value in real space from the 4 4 fringe width in q space. 5 5 """ 6 6 from math import pi, fabs 7 7 _DQ_DEFAULT = 0.05 8 9 8 10 class KiessigThicknessCalculator(object): 9 11 """ … … 13 15 14 16 # dq value 15 self.deltaq = _DQ_DEFAULT 17 self.deltaq = _DQ_DEFAULT 16 18 # thickenss value 17 19 self.thickness = None … … 26 28 """ 27 29 # set dq 28 self.deltaq 30 self.deltaq = dq 29 31 30 32 def get_deltaq(self): … … 33 35 """ 34 36 # return dq 35 return self.deltaq 37 return self.deltaq 36 38 37 39 def compute_thickness(self): … … 51 53 else: 52 54 # calculate thickness 53 thickness = 2*pi/fabs(dq) 54 # return thickness value 55 thickness = 2*pi/fabs(dq) 56 # return thickness value 55 57 return thickness 56 58 57 def get_thickness_unit(self): 59 def get_thickness_unit(self): 58 60 """ 59 61 :return: the thickness unit. … … 61 63 # unit of thickness 62 64 return self.thickness_unit 63 64 -
sanscalculator/src/sans/calculator/resolution_calculator.py
r5f3164c r8a621ac 1 1 """ 2 2 This object is a small tool to allow user to quickly 3 determine the variance in q from the 3 determine the variance in q from the 4 4 instrumental parameters. 5 5 """ … … 11 11 from math import pi 12 12 from math import sqrt 13 import math 14 import scipy 13 import math 15 14 import numpy 16 15 … … 19 18 #Gravitational acc. in cgs unit 20 19 _GRAVITY = 981.0 20 21 21 22 22 class ResolutionCalculator(object): … … 86 86 self.qyrange = [] 87 87 88 def compute_and_plot(self, qx_value, qy_value, qx_min, qx_max, 89 qy_min, qy_max, coord ='cartesian'):88 def compute_and_plot(self, qx_value, qy_value, qx_min, qx_max, 89 qy_min, qy_max, coord='cartesian'): 90 90 """ 91 91 Compute the resolution … … 98 98 # wavelength etc. 99 99 lamda_list, dlamb_list = self.get_wave_list() 100 intens_list = []#self.get_intensity_list() 101 100 intens_list = [] 102 101 sig1_list = [] 103 102 sig2_list = [] … … 120 119 self.compute(lam, dlam, qx_value, qy_value, coord, tof) 121 120 # make image 122 image = self.get_image(qx_value, qy_value, sigma_1, sigma_2, 121 image = self.get_image(qx_value, qy_value, sigma_1, sigma_2, 123 122 sigma_r, qx_min, qx_max, qy_min, qy_max, 124 123 coord, False) 125 124 126 # Non tof mode to be speed up 125 # Non tof mode to be speed up 127 126 #if num_lamda < 2: 128 127 # return self.plot_image(image) … … 140 139 self.qxrange = [qx_min, qx_max] 141 140 self.qyrange = [qy_min, qy_max] 142 #print qy_max+qy_min,qy_max,qy_min143 141 sig1_list.append(sigma_1) 144 142 sig2_list.append(sigma_2) 145 143 sigr_list.append(sigma_r) 146 144 sigma1d_list.append(sigma1d) 147 # redraw image in global 2d q-space. 145 # redraw image in global 2d q-space. 148 146 self.image_lam = [] 149 147 total_intensity = 0 … … 156 154 dlam = dlamb_list[ind] 157 155 intens = self.setup_tof(lam, dlam) 158 out = self.get_image(qx_value, qy_value, sig1_list[ind], 159 sig2_list[ind], sigr_list[ind], 156 out = self.get_image(qx_value, qy_value, sig1_list[ind], 157 sig2_list[ind], sigr_list[ind], 160 158 qx_min, qx_max, qy_min, qy_max, coord) 161 159 # this is the case of q being outside the detector … … 164 162 image = out 165 163 # set variance as sigmas 166 sigma_1 += sig1_list[ind] * 164 sigma_1 += sig1_list[ind] * sig1_list[ind] * self.intensity 167 165 sigma_r += sigr_list[ind] * sigr_list[ind] * self.intensity 168 166 sigma_2 += sig2_list[ind] * sig2_list[ind] * self.intensity … … 182 180 self.sigma_2 = sqrt(sigma_2) 183 181 self.sigma_1d = sqrt(sigma1d) 184 # rescale 185 max_im_val = 1 #image_out.max()182 # rescale 183 max_im_val = 1 186 184 if max_im_val > 0: 187 185 image_out /= max_im_val … … 199 197 200 198 # plot image 201 return self.plot_image(self.image) 199 return self.plot_image(self.image) 202 200 203 201 def setup_tof(self, wavelength, wavelength_spread): … … 214 212 215 213 if wavelength == 0: 216 msg = "Can't compute the resolution: the wavelength is zero..." 214 msg = "Can't compute the resolution: the wavelength is zero..." 217 215 raise RuntimeError, msg 218 216 return self.intensity 219 217 220 def compute(self, wavelength, wavelength_spread, qx_value, qy_value, 221 coord ='cartesian', tof=False):218 def compute(self, wavelength, wavelength_spread, qx_value, qy_value, 219 coord='cartesian', tof=False): 222 220 """ 223 221 Compute the Q resoltuion in || and + direction of 2D … … 232 230 if tof: 233 231 # rectangular 234 tof_factor = 232 tof_factor = 2 235 233 else: 236 234 # triangular 237 tof_factor = 235 tof_factor = 1 238 236 # Find polar values 239 237 qr_value, phi = self._get_polar_value(qx_value, qy_value) … … 285 283 sigma_1 += self.get_variance(rthree, l_two, phi, comp1) 286 284 # for gravity term for 1d 287 sigma_1grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, lamb_spread,288 phi, comp1, 'on') / tof_factor285 sigma_1grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, 286 lamb_spread, phi, comp1, 'on') / tof_factor 289 287 # for wavelength spread 290 288 # reserve for 1d calculation 291 289 A_value = self._cal_A_value(lamb, l_ssa, l_sad) 292 sigma_wave_1, sigma_wave_1_1d = self.get_variance_wave(A_value, 293 radius, l_two, lamb_spread, 290 sigma_wave_1, sigma_wave_1_1d = self.get_variance_wave(A_value, 291 radius, l_two, lamb_spread, 294 292 phi, 'radial', 'on') 295 293 sigma_wave_1 /= tof_factor 296 sigma_wave_1_1d /= 294 sigma_wave_1_1d /= tof_factor 297 295 # for 1d 298 variance_1d_1 = (sigma_1 + sigma_1grav1d) / 2 + sigma_wave_1_1d296 variance_1d_1 = (sigma_1 + sigma_1grav1d) / 2 + sigma_wave_1_1d 299 297 # normalize 300 298 variance_1d_1 = knot * knot * variance_1d_1 / 12 … … 303 301 #sigma_1 += sigma_wave_1 304 302 # normalize 305 sigma_1 = knot *sqrt(sigma_1 / 12)306 sigma_r = knot *sqrt(sigma_wave_1 / (tof_factor *12))303 sigma_1 = knot * sqrt(sigma_1 / 12) 304 sigma_r = knot * sqrt(sigma_wave_1 / (tof_factor *12)) 307 305 # sigma in the phi/y direction 308 306 # for source apperture … … 316 314 317 315 # for gravity term for 1d 318 sigma_2grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, lamb_spread, 319 phi, comp2, 'on') / tof_factor 320 321 316 sigma_2grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, 317 lamb_spread, phi, comp2, 'on') / tof_factor 318 322 319 # for wavelength spread 323 320 # reserve for 1d calculation 324 sigma_wave_2, sigma_wave_2_1d = self.get_variance_wave(A_value, 325 radius, l_two, lamb_spread, 326 phi, 'phi', 'on') 327 sigma_wave_2 /= 328 sigma_wave_2_1d /= 321 sigma_wave_2, sigma_wave_2_1d = self.get_variance_wave(A_value, 322 radius, l_two, lamb_spread, 323 phi, 'phi', 'on') 324 sigma_wave_2 /= tof_factor 325 sigma_wave_2_1d /= tof_factor 329 326 # for 1d 330 327 variance_1d_2 = (sigma_2 + sigma_2grav1d) / 2 + sigma_wave_2_1d 331 328 # normalize 332 variance_1d_2 = knot *knot*variance_1d_2 / 12329 variance_1d_2 = knot * knot * variance_1d_2 / 12 333 330 334 331 # for 2d … … 336 333 #sigma_2 += sigma_wave_2 337 334 # normalize 338 sigma_2 = 335 sigma_2 = knot * sqrt(sigma_2 / 12) 339 336 sigma1d = sqrt(variance_1d_1 + variance_1d_2) 340 337 # set sigmas … … 345 342 return qr_value, phi, sigma_1, sigma_2, sigma_r, sigma1d 346 343 347 def _within_detector_range(self, qx_value, qy_value):344 def _within_detector_range(self, qx_value, qy_value): 348 345 """ 349 346 check if qvalues are within detector range … … 369 366 370 367 def get_image(self, qx_value, qy_value, sigma_1, sigma_2, sigma_r, 371 qx_min, qx_max, qy_min, qy_max, 372 coord = 'cartesian', full_cal=True):368 qx_min, qx_max, qy_min, qy_max, 369 coord='cartesian', full_cal=True): 373 370 """ 374 371 Get the resolution in polar coordinate ready to plot … … 380 377 """ 381 378 # Get qx_max and qy_max... 382 output =self._get_detector_qxqy_pixels()379 self._get_detector_qxqy_pixels() 383 380 384 381 qr_value, phi = self._get_polar_value(qx_value, qy_value) 385 382 386 383 # Check whether the q value is within the detector range 387 msg = "Invalid input: Q value out of the detector range..."384 #msg = "Invalid input: Q value out of the detector range..." 388 385 if qx_min < self.qx_min: 389 386 self.qx_min = qx_min … … 416 413 qc_2 = 0.0 417 414 # Calculate the 2D Gaussian distribution image 418 image = self._gaussian2d_polar(q_1, q_2, qc_1, qc_2, 415 image = self._gaussian2d_polar(q_1, q_2, qc_1, qc_2, 419 416 sigma_1, sigma_2, sigma_r) 420 417 else: … … 426 423 427 424 # Calculate the 2D Gaussian distribution image 428 image = self._gaussian2d(q_1, q_2, qc_1, qc_2, 425 image = self._gaussian2d(q_1, q_2, qc_1, qc_2, 429 426 sigma_1, sigma_2, sigma_r) 430 427 # out side of detector … … 445 442 """ 446 443 Plot image using pyplot 447 : image: 2d resolution image 444 : image: 2d resolution image 448 445 449 446 : return plt: pylab object … … 455 452 plt.ylabel('$\\rm{Q}_{y} [A^{-1}]$') 456 453 # Max value of the image 457 max = numpy.max(image)454 # max = numpy.max(image) 458 455 qx_min, qx_max, qy_min, qy_max = self.get_detector_qrange() 459 456 460 457 # Image 461 im = plt.imshow(image, 462 extent =[qx_min, qx_max, qy_min, qy_max])458 im = plt.imshow(image, 459 extent=[qx_min, qx_max, qy_min, qy_max]) 463 460 464 461 # bilinear interpolation to make it smoother … … 473 470 self.image = [] 474 471 475 def get_variance(self, size = [], distance = 0, phi = 0, comp ='radial'):472 def get_variance(self, size=[], distance=0, phi=0, comp='radial'): 476 473 """ 477 474 Get the variance when the slit/pinhole size is given 478 : size: list that can be one(diameter for circular) 475 : size: list that can be one(diameter for circular) 479 476 or two components(lengths for rectangular) 480 477 : distance: [z, x] where z along the incident beam, x // qx_value … … 482 479 483 480 : return variance: sigma^2 484 """ 481 """ 485 482 # check the length of size (list) 486 483 len_size = len(size) … … 509 506 # for rectangular slit 510 507 elif len_size == 2: 511 x_comp = size[0] * phi_x 508 x_comp = size[0] * phi_x 512 509 y_comp = size[1] * phi_y 513 510 # otherwise 514 511 else: 515 512 raise ValueError, " Improper input..." 516 # get them squared 517 sigma = x_comp * x_comp 513 # get them squared 514 sigma = x_comp * x_comp 518 515 sigma += y_comp * y_comp 519 516 # normalize by distance … … 522 519 return sigma 523 520 524 def get_variance_wave(self, A_value, radius, distance, spread, phi, 525 comp = 'radial', switch ='on'):521 def get_variance_wave(self, A_value, radius, distance, spread, phi, 522 comp='radial', switch='on'): 526 523 """ 527 524 Get the variance when the wavelength spread is given … … 547 544 sigma1d *= (math.sin(phi)*math.sin(phi)) 548 545 else: 549 sigma1d *= 1 550 # sigma^2 for 2d 546 sigma1d *= 1 547 # sigma^2 for 2d 551 548 # shift the coordinate due to the gravitational shift 552 549 rad_x = radius * math.cos(phi) 553 550 rad_y = A_value - radius * math.sin(phi) 554 551 radius = math.sqrt(rad_x * rad_x + rad_y * rad_y) 555 # new phi 552 # new phi 556 553 phi = math.atan2(-rad_y, rad_x) 557 self.gravity_phi = phi 554 self.gravity_phi = phi 558 555 # calculate sigma^2 559 556 sigma = 2 * math.pow(radius/distance*spread, 2) … … 563 560 sigma *= (math.sin(phi)*math.sin(phi)) 564 561 else: 565 sigma *= 1 562 sigma *= 1 566 563 567 564 return sigma, sigma1d 568 565 569 def get_variance_gravity(self, s_distance, d_distance, wavelength, spread, 570 phi, comp = 'radial', switch ='on'):566 def get_variance_gravity(self, s_distance, d_distance, wavelength, spread, 567 phi, comp='radial', switch='on'): 571 568 """ 572 569 Get the variance from gravity when the wavelength spread is given … … 603 600 return sigma 604 601 605 def _cal_A_value(self, lamda, s_distance, d_distance): 602 def _cal_A_value(self, lamda, s_distance, d_distance): 606 603 """ 607 604 Calculate A value for gravity … … 617 614 gravy = _GRAVITY 618 615 # m/h 619 m_over_h = self.mass / h_constant616 m_over_h = self.mass / h_constant 620 617 # A value 621 618 a_value = d_distance * (s_distance + d_distance) … … 641 638 return self.wave.wavelength 642 639 640 #TODO: why was this method duplicated? 641 #def get_spectrum(self): 642 # """ 643 # Get spectrum 644 # """ 645 # return self.wave.spectrum 646 647 def get_default_spectrum(self): 648 """ 649 Get default_spectrum 650 """ 651 return self.wave.get_default_spectrum() 652 643 653 def get_spectrum(self): 644 654 """ 645 Get spectrum646 """647 return self.wave.spectrum648 649 def get_default_spectrum(self):650 """651 Get default_spectrum652 """653 return self.wave.get_default_spectrum()654 655 def get_spectrum(self):656 """657 655 Get _spectrum 658 656 """ 659 return self.wave.get_spectrum() 657 return self.wave.get_spectrum() 660 658 661 659 def get_wavelength_spread(self): … … 693 691 Get detector size 694 692 """ 695 return self.detector.size 693 return self.detector.size 696 694 697 695 def get_source2sample_distance(self): … … 755 753 """ 756 754 self.spectrum = spectrum 757 self.wave.set_spectrum(spectrum) 755 self.wave.set_spectrum(spectrum) 758 756 759 757 def set_wavelength_spread(self, wavelength_spread): … … 787 785 788 786 : param size: [dia_value] or [x_value, y_value] 789 """ 787 """ 790 788 if len(size) < 1 or len(size) > 2: 791 789 raise RuntimeError, "The length of the size must be one or two." … … 884 882 return qx_min, qx_max, qy_min, qy_max 885 883 886 def _rotate_z(self, x_value, y_value, theta= 884 def _rotate_z(self, x_value, y_value, theta=0.0): 887 885 """ 888 886 Rotate x-y cordinate around z-axis by theta … … 895 893 # rotate by theta 896 894 x_prime = x_value * math.cos(theta) + y_value * math.sin(theta) 897 y_prime = 895 y_prime = -x_value * math.sin(theta) + y_value * math.cos(theta) 898 896 899 897 return x_prime, y_prime 900 898 901 def _gaussian2d(self, x_val, y_val, x0_val, y0_val, 899 def _gaussian2d(self, x_val, y_val, x0_val, y0_val, 902 900 sigma_x, sigma_y, sigma_r): 903 901 """ … … 927 925 y_p = -x_value * sin_phi + y_value * cos_phi 928 926 929 new_sig_x = sqrt(sigma_r * sigma_r / (sigma_x * sigma_x 930 new_sig_y = sqrt(sigma_r * sigma_r / (sigma_y * sigma_y 927 new_sig_x = sqrt(sigma_r * sigma_r / (sigma_x * sigma_x) + 1) 928 new_sig_y = sqrt(sigma_r * sigma_r / (sigma_y * sigma_y) + 1) 931 929 new_x = x_p * cos_phi / new_sig_x - y_p * sin_phi 932 930 new_x /= sigma_x … … 934 932 new_y /= sigma_y 935 933 936 nu_value = -0.5 * (new_x * new_x + new_y * new_y)934 nu_value = -0.5 * (new_x * new_x + new_y * new_y) 937 935 938 936 gaussian = numpy.exp(nu_value) … … 942 940 return gaussian 943 941 944 def _gaussian2d_polar(self, x_val, y_val, x0_val, y0_val, 942 def _gaussian2d_polar(self, x_val, y_val, x0_val, y0_val, 945 943 sigma_x, sigma_y, sigma_r): 946 944 """ 947 Calculate 2D Gaussian distribution for polar coodinate 945 Calculate 2D Gaussian distribution for polar coodinate 948 946 : x_val: x value 949 947 : y_val: y value … … 957 955 """ 958 956 sigma_x = sqrt(sigma_x * sigma_x + sigma_r * sigma_r) 959 # call gaussian1d 957 # call gaussian1d 960 958 gaussian = self._gaussian1d(x_val, x0_val, sigma_x) 961 959 gaussian *= self._gaussian1d(y_val, y0_val, sigma_y) … … 971 969 : value: value 972 970 : mean: mean value 973 : sigma: variance 971 : sigma: variance 974 972 975 973 : return: gaussian (value) … … 1028 1026 wavelength = self.wave.wavelength 1029 1027 # Gavity correction 1030 delta_y = self._get_beamcenter_drop() # in cm1028 delta_y = self._get_beamcenter_drop() # in cm 1031 1029 1032 1030 # detector_pix size … … 1042 1040 else: 1043 1041 raise ValueError, " Input value format error..." 1044 # Sample to detector distance = sample slit to detector 1042 # Sample to detector distance = sample slit to detector 1045 1043 # minus sample offset 1046 1044 sample2detector_distance = self.sample2detector_distance[0] - \ … … 1098 1096 1099 1097 # qx_value and qy_value values in array 1100 qx_value = qx_value.repeat(detector_pix_nums_y) 1098 qx_value = qx_value.repeat(detector_pix_nums_y) 1101 1099 qx_value = qx_value.reshape(detector_pix_nums_x, detector_pix_nums_y) 1102 qy_value = qy_value.repeat(detector_pix_nums_x) 1100 qy_value = qy_value.repeat(detector_pix_nums_x) 1103 1101 qy_value = qy_value.reshape(detector_pix_nums_y, detector_pix_nums_x) 1104 1102 qy_value = qy_value.transpose() … … 1110 1108 self.qy_max = numpy.max(qy_value) 1111 1109 1112 # Appr. min and max values of the detector display limits 1110 # Appr. min and max values of the detector display limits 1113 1111 # i.e., edges of the last pixels. 1114 self.qy_min += self._get_qx(-0.5 * pix_y_size, 1112 self.qy_min += self._get_qx(-0.5 * pix_y_size, 1115 1113 sample2detector_distance, wavelength) 1116 self.qy_max += self._get_qx(0.5 * pix_y_size, 1114 self.qy_max += self._get_qx(0.5 * pix_y_size, 1117 1115 sample2detector_distance, wavelength) 1118 1116 #if self.qx_min == self.qx_max: 1119 self.qx_min += self._get_qx(-0.5 * pix_x_size, 1117 self.qx_min += self._get_qx(-0.5 * pix_x_size, 1120 1118 sample2detector_distance, wavelength) 1121 self.qx_max += self._get_qx(0.5 * pix_x_size, 1119 self.qx_max += self._get_qx(0.5 * pix_x_size, 1122 1120 sample2detector_distance, wavelength) 1123 1121 … … 1128 1126 self.detector_qy_max = self.qy_max 1129 1127 1130 1131 1132 1128 # try to set it as a Data2D otherwise pass (not required for now) 1133 1129 try: … … 1135 1131 output = Data2D() 1136 1132 inten = numpy.zeros_like(qx_value) 1137 output.data 1138 output.qx_data 1139 output.qy_data = qy_value1133 output.data = inten 1134 output.qx_data = qx_value 1135 output.qy_data = qy_value 1140 1136 except: 1141 1137 pass 1142 1138 1143 return output #qx_value,qy_value1139 return output 1144 1140 1145 1141 def _get_qx(self, dx_size, det_dist, wavelength): … … 1153 1149 plane_dist = dx_size 1154 1150 # full scattering angle on the x-axis 1155 theta 1156 qx_value 1157 return qx_value 1151 theta = numpy.arctan(plane_dist / det_dist) 1152 qx_value = (2.0 * pi / wavelength) * numpy.sin(theta) 1153 return qx_value 1158 1154 1159 1155 def _get_polar_value(self, qx_value, qy_value): … … 1187 1183 pos_y -= offset_y 1188 1184 1189 return pos_x, pos_y 1185 return pos_x, pos_y 1190 1186 1191 1187 def _get_beamcenter_drop(self): … … 1195 1191 :return delta y: the beam center drop in cm 1196 1192 """ 1197 # Check if mass == 0 (X-ray). 1193 # Check if mass == 0 (X-ray). 1198 1194 if self.mass == 0: 1199 1195 return 0 … … 1211 1207 delta_y /= (velocity * velocity) 1212 1208 1213 return delta_y 1214 1215 1209 return delta_y -
sanscalculator/src/sans/calculator/slit_length_calculator.py
rf93b6a1e r8a621ac 3 3 determine the slit length value of data. 4 4 """ 5 5 6 6 7 class SlitlengthCalculator(object): … … 18 19 self.slit_length = 0.0 19 20 20 # The unit is unknown from SAXSess profile: 21 # It seems 1/nm but it could be not fixed, 21 # The unit is unknown from SAXSess profile: 22 # It seems 1/nm but it could be not fixed, 22 23 # so users should be notified to determine the unit by themselves. 23 24 self.slit_length_unit = "unknown" … … 25 26 def set_data(self, x=None, y=None): 26 27 """ 27 Receive two vector x, y and prepare the slit calculator for 28 Receive two vector x, y and prepare the slit calculator for 28 29 computation. 29 30 … … 50 51 max_y = y.max() 51 52 52 # initial values 53 # initial values 53 54 y_sum = 0.0 54 55 y_max = 0.0 55 56 ind = 0.0 56 57 57 # sum 10 or more y values until getting max_y, 58 # sum 10 or more y values until getting max_y, 58 59 while (True): 59 60 if ind >= 10 and y_max == max_y: 60 61 break 61 62 y_sum = y_sum + y[ind] 62 if y[ind] > y_max: 63 if y[ind] > y_max: 63 64 y_max = y[ind] 64 65 ind += 1 65 66 66 # find the average value/2 of the top values 67 # find the average value/2 of the top values 67 68 y_half = y_sum/(2.0*ind) 68 69 … … 73 74 while (True): 74 75 # no need to check when ind == 0 75 ind += 1 76 # y value and ind just after passed the spot of the half height 77 y_half_d = y[ind] 78 if y[ind] < y_half: 76 ind += 1 77 # y value and ind just after passed the spot of the half height 78 y_half_d = y[ind] 79 if y[ind] < y_half: 79 80 break 80 81 81 82 # y value and ind just before passed the spot of the half height 82 y_half_u = y[ind-1] 83 y_half_u = y[ind-1] 83 84 84 85 # get corresponding x values 85 86 x_half_d = x[ind] 86 x_half_u = x[ind-1] 87 x_half_u = x[ind-1] 87 88 88 89 # calculate x at y = y_half using linear interpolation … … 91 92 else: 92 93 x_half = (x_half_u * (y_half - y_half_d) \ 93 + x_half_d *(y_half_u-y_half)) \94 / (y_half_u - y_half_d)94 + x_half_d * (y_half_u - y_half)) \ 95 / (y_half_u - y_half_d) 95 96 96 97 # Our slit length is half width, so just give half beam value … … 101 102 return self.slit_length 102 103 103 def get_slit_length_unit(self): 104 def get_slit_length_unit(self): 104 105 """ 105 106 :return: the slit length unit.
Note: See TracChangeset
for help on using the changeset viewer.
