Commits
38 38 | self.antenna_block_m = 0.0 |
39 39 | self.taper = 10 |
40 40 | self.ref_freq = -1. |
41 41 | self.kernel_type = "" |
42 42 | self.kernel_param = {} |
43 43 | self.pa = "0.0deg" |
44 44 | self.sampling_arcsec = [] |
45 45 | self.cell_arcsec = [] |
46 46 | |
47 47 | def __to_arcsec_list(self, angle): |
48 - | """return a list of angles in arcsec (value only without unit)""" |
48 + | """Return a list of angles in arcsec (value only without unit).""" |
49 49 | if type(angle) not in [list, tuple, np.ndarray]: |
50 50 | angle = [angle] |
51 51 | return [self.__to_arcsec(val) for val in angle] |
52 52 | |
53 53 | def __to_arcsec(self, angle): |
54 - | """convert angle to arcsec and return the value without unit.""" |
54 + | """Convert angle to arcsec and return the value without unit.""" |
55 55 | my_qa = quanta() |
56 56 | if my_qa.isangle(angle): |
57 57 | return my_qa.getvalue(my_qa.convert(angle, "arcsec"))[0] |
58 58 | elif my_qa.getunit(angle) == '': |
59 59 | return float(angle) |
60 60 | else: |
61 61 | raise ValueError("Invalid angle: %s" % (str(angle))) |
62 62 | |
63 63 | def __parse_width(self, val, cell_size_arcsec): |
64 - | """ |
65 - | Convert value in unit of arcsec |
66 - | if val is angle, returns a float value in unit of arcsec. |
64 + | """Convert value in unit of arcsec. |
65 + | |
66 + | If val is angle, returns a float value in unit of arcsec. |
67 67 | else the unit is assumed to be pixel and multiplied by cell_size_arcsec |
68 68 | """ |
69 69 | my_qa = quanta() |
70 70 | if my_qa.isangle(val): |
71 71 | return self.__to_arcsec(val) |
72 72 | elif my_qa.getunit(val) in ('', 'pixel'): |
73 73 | return my_qa.getvalue(val) * cell_size_arcsec |
74 74 | else: |
75 75 | raise ValueError("Invalid width %s" % str(val)) |
76 76 | |
77 77 | def set_antenna(self, diam, blockage="0.0m", taper=10): |
78 - | """ |
79 - | set parameters to construct antenna beam |
78 + | """Set parameters to construct antenna beam. |
79 + | |
80 80 | antenna diameter and blockage |
81 81 | taper: the illumination taper in dB |
82 82 | """ |
83 83 | # try quantity |
84 84 | my_qa = quanta() |
85 85 | self.antenna_diam_m = my_qa.getvalue(my_qa.convert(diam, "m"))[0] |
86 86 | self.antenna_block_m = my_qa.getvalue(my_qa.convert(blockage, "m"))[0] |
87 87 | self.taper = taper |
88 88 | self.is_antenna_set = True |
89 89 | |
90 90 | def set_sampling(self, intervals, pa="0deg"): |
91 - | """ |
92 - | Aet sampling interval of observation |
91 + | """Set sampling interval of observation. |
92 + | |
93 93 | intervals: pointing inteval of observation (['10arcsec','10arcsec']) |
94 94 | pa: position angle (NOT USED) |
95 95 | """ |
96 96 | self.pa = pa |
97 97 | self.sampling_arcsec = [abs(a) for a in |
98 98 | self.__to_arcsec_list(intervals)] |
99 99 | self.is_sampling_set = True |
100 100 | |
101 101 | def set_image_param(self, cell, ref_freq, gridfunction, |
102 102 | convsupport, truncate, gwidth, jwidth, is_alma=False): |
103 - | """ |
104 - | Set imaging parameters |
103 + | """Set imaging parameters. |
104 + | |
105 105 | cell: image pixel size |
106 106 | ref_freq: reference frequency of image to calculate beam size |
107 107 | gridfunction, convsupport, truncate, gwidth, jwidth: |
108 108 | parameters passed to imager |
109 109 | is_alma: valid only for PB kernel to use 10.7m |
110 110 | """ |
111 111 | self.ref_freq = ref_freq |
112 112 | self.cell_arcsec = [abs(a) for a in |
113 113 | self.__to_arcsec_list(cell)] |
114 114 | if gridfunction.upper() == "SF": |
117 117 | self.__set_gjinc_kernel(truncate, gwidth, jwidth) |
118 118 | elif gridfunction.upper() == "GAUSS": |
119 119 | self.__set_gauss_kernel(truncate, gwidth) |
120 120 | elif gridfunction.upper() == "BOX": |
121 121 | self.__set_box_kernel(self.cell_arcsec[0]) |
122 122 | elif gridfunction.upper() == "PB": |
123 123 | self.__set_pb_kernel(is_alma) |
124 124 | self.is_kernel_set = True |
125 125 | |
126 126 | def __set_sf_kernel(self, convsupport): |
127 - | """Set SF kernel parameter to the class""" |
127 + | """Set SF kernel parameter to the class.""" |
128 128 | self.kernel_type = "SF" |
129 129 | self.kernel_param = dict(convsupport=convsupport) |
130 130 | |
131 131 | def __set_gjinc_kernel(self, truncate, gwidth, jwidth): |
132 - | """Set GJINC kernel parameter to the class""" |
132 + | """Set GJINC kernel parameter to the class.""" |
133 133 | self.kernel_type = "GJINC" |
134 134 | self.kernel_param = dict(truncate=truncate, gwidth=gwidth, jwidth=jwidth) |
135 135 | |
136 136 | def __set_gauss_kernel(self, truncate, gwidth): |
137 - | """Set GAUSS kernel parameter to the class""" |
137 + | """Set GAUSS kernel parameter to the class.""" |
138 138 | self.kernel_type = "GAUSS" |
139 139 | self.kernel_param = dict(truncate=truncate, gwidth=gwidth) |
140 140 | |
141 141 | def __set_box_kernel(self, width): |
142 - | """Set BOX kernel parameter to the class""" |
142 + | """Set BOX kernel parameter to the class.""" |
143 143 | self.kernel_type = "BOX" |
144 144 | self.kernel_param = dict(width=width) |
145 145 | |
146 146 | def __set_pb_kernel(self, alma=False): |
147 - | """Set PB kernel parameter to the class""" |
147 + | """Set PB kernel parameter to the class.""" |
148 148 | self.kernel_type = "PB" |
149 149 | self.kernel_param = dict(alma=alma) |
150 150 | |
151 151 | def get_beamsize_image(self): |
152 - | """ |
153 - | Returns FWHM of theoretical beam size in image. |
152 + | """Return FWHM of theoretical beam size in image. |
153 + | |
154 154 | The FWHM is derived by fitting gridded beam with Gaussian function. |
155 155 | """ |
156 156 | # assert necessary information is set |
157 157 | self.__assert_antenna() |
158 158 | self.__assert_kernel() |
159 159 | self.__assert_sampling() |
160 160 | casalog.post("Calculating theoretical beam size of the image") |
161 161 | # construct theoretic beam for image |
162 162 | axis, beam = self.get_antenna_beam() |
163 163 | casalog.post( |
210 210 | fwhm1) |
211 211 | del result |
212 212 | # clear-up axes |
213 213 | del axis, beam, kernel, sampling, gridded |
214 214 | |
215 215 | return dict(major="%farcsec" % (fwhm_geo_arcsec), |
216 216 | minor="%farcsec" % (fwhm_geo_arcsec), |
217 217 | pa=self.pa) |
218 218 | |
219 219 | def __assert_antenna(self): |
220 - | """Raise an error if antenna information is not set""" |
220 + | """Raise an error if antenna information is not set.""" |
221 221 | if not self.is_antenna_set: |
222 222 | raise RuntimeError("Antenna is not set") |
223 223 | |
224 224 | def __assert_kernel(self): |
225 - | """Raise an error if imaging parameters are not set""" |
225 + | """Raise an error if imaging parameters are not set.""" |
226 226 | if not self.is_kernel_set: |
227 227 | raise RuntimeError("Kernel is not set.") |
228 228 | |
229 229 | def __assert_sampling(self): |
230 - | """Raise an error if sampling information is not set""" |
230 + | """Raise an error if sampling information is not set.""" |
231 231 | if not self.is_sampling_set: |
232 232 | raise RuntimeError("Sampling information is not set") |
233 233 | |
234 234 | def get_antenna_beam(self): |
235 - | """ |
236 - | Returns arrays of antenna beam response and it's horizontal axis |
235 + | """Return arrays of antenna beam response and it's horizontal axis. |
236 + | |
237 237 | Picked from AnalysisUtils (revision 1.2204, 2015/02/18) |
238 238 | """ |
239 239 | self.__assert_antenna() |
240 240 | self.__assert_kernel() |
241 241 | |
242 242 | fwhm_arcsec = primaryBeamArcsec(frequency=self.ref_freq, |
243 243 | diameter=self.antenna_diam_m, |
244 244 | taper=self.taper, showEquation=True, |
245 245 | obscuration=self.antenna_block_m, |
246 246 | fwhmfactor=None) |
270 270 | diameter=self.antenna_diam_m) |
271 271 | else: |
272 272 | # Gaussian beam |
273 273 | myxaxis = np.arange(-3 * fwhm_arcsec, |
274 274 | 3 * fwhm_arcsec + 0.5 * convolutionPixelSize, |
275 275 | convolutionPixelSize) |
276 276 | myfunction = self.gauss(myxaxis, [fwhm_arcsec, truncate]) |
277 277 | return myxaxis, myfunction |
278 278 | |
279 279 | def get_kernel(self, axis): |
280 - | """Returns imaging kernel array""" |
280 + | """Return imaging kernel array.""" |
281 281 | self.__assert_kernel() |
282 282 | if self.kernel_type == "SF": |
283 283 | # TODO: what to do for cell[0]!=cell[1]??? |
284 284 | return self.get_sf_kernel(axis, self.kernel_param['convsupport'], |
285 285 | self.cell_arcsec[0]) |
286 286 | elif self.kernel_type == "GJINC": |
287 287 | return self.get_gjinc_kernel(axis, self.kernel_param['truncate'], |
288 288 | self.kernel_param['gwidth'], |
289 289 | self.kernel_param['jwidth'], |
290 290 | self.cell_arcsec[0]) |
301 301 | casalog.post( |
302 302 | f"Using effective antenna diameter {diam}m for " |
303 303 | f"{self.kernel_type} kernel of ALMA antennas") |
304 304 | epsilon = self.antenna_block_m / diam |
305 305 | return self.get_pb_kernel(axis, diam, self.ref_freq, epsilon=epsilon) |
306 306 | # return (self.rootAiryIntensity(axis, epsilon))**2 |
307 307 | else: |
308 308 | raise RuntimeError("Invalid kernel: %s" % self.kernel_type) |
309 309 | |
310 310 | def summary(self): |
311 - | """Print summary of parameters set to the class""" |
311 + | """Print summary of parameters set to the class.""" |
312 312 | casalog.post("=" * 40) |
313 313 | casalog.post("Summary of Image Beam Parameters") |
314 314 | casalog.post("=" * 40) |
315 315 | casalog.post("[Antenna]") |
316 316 | self.__antenna_summary() |
317 317 | |
318 318 | casalog.post("\n[Imaging Parameters]") |
319 319 | self.__kernel_summary() |
320 320 | |
321 321 | casalog.post("\n[Sampling]") |
322 322 | self.__sampling_summary() |
323 323 | |
324 324 | def __notset_message(self): |
325 325 | casalog.post("Not set.") |
326 326 | |
327 327 | def __antenna_summary(self): |
328 - | """Print summary of antenna setup""" |
328 + | """Print summary of antenna setup.""" |
329 329 | if not self.is_antenna_set: |
330 330 | self.__notset_message() |
331 331 | return |
332 332 | casalog.post("diameter: %f m" % (self.antenna_diam_m)) |
333 333 | casalog.post("blockage: %f m" % (self.antenna_block_m)) |
334 334 | |
335 335 | def __kernel_summary(self): |
336 - | """Print summary of imaging parameter setup""" |
336 + | """Print summary of imaging parameter setup.""" |
337 337 | if not self.is_kernel_set: |
338 338 | self.__notset_message() |
339 339 | return |
340 340 | casalog.post("reference frequency: %s" % str(self.ref_freq)) |
341 341 | casalog.post("cell size: %s arcsec" % str(self.cell_arcsec)) |
342 342 | casalog.post("kernel type: %s" % self.kernel_type) |
343 343 | for key, val in self.kernel_param.items(): |
344 344 | casalog.post("%s: %s" % (key, str(val))) |
345 345 | |
346 346 | def __sampling_summary(self): |
347 - | """Print summary of sampling setup""" |
347 + | """Print summary of sampling setup.""" |
348 348 | if not self.is_sampling_set: |
349 349 | self.__notset_message() |
350 350 | return |
351 351 | casalog.post("sampling interval: %s arcsec" % str(self.sampling_arcsec)) |
352 352 | casalog.post("position angle: %s" % (self.pa)) |
353 353 | |
354 354 | # |
355 355 | # Construct Kernel arrays |
356 356 | # |
357 357 | |
358 358 | # |
359 359 | # BOX |
360 360 | # |
361 361 | def get_box_kernel(self, axis, width): |
362 - | """ |
363 - | Returns a box kernel array with specified width. |
362 + | """Return a box kernel array with specified width. |
363 + | |
364 364 | axis: an array of xaxis values |
365 365 | width: kernel width |
366 366 | output array |
367 367 | out[i] = 1.0 (-width/2.0 <= x[i] <= width/2.0) |
368 368 | = 0.0 (else) |
369 369 | """ |
370 370 | data = np.zeros(len(axis)) |
371 371 | indices = np.where(abs(axis) <= width / 2.0) |
372 372 | data[indices] = 1.0 |
373 373 | return data |
374 374 | |
375 375 | # |
376 376 | # SF |
377 377 | # |
378 378 | def get_sf_kernel(self, axis, convsupport, cell_size): |
379 - | """ |
380 - | Returns spheroidal kernel array |
379 + | """Return spheroidal kernel array. |
380 + | |
381 381 | axis: an array of xaxis values |
382 382 | convsupport: the truncation radius of SF kernel in pixel unit. |
383 383 | cell_size: image pixel size |
384 384 | |
385 385 | Modified version of one in AnalysisUtils.sfBeamPredict (revision 1.2204, 2015/02/18) |
386 386 | """ |
387 387 | convsupport = 3 if convsupport == -1 else convsupport |
388 388 | supportwidth = (convsupport * 1.0 + 0.0) |
389 389 | # value obtained by matching Fred's grdsf.f output with scipy(m=0,n=0) |
390 390 | c = 5.356 * np.pi / 2.0 |
391 391 | sfaxis = axis / float(supportwidth * cell_size * 1.0) |
392 392 | indices = np.where(abs(sfaxis) < 1)[0] |
393 393 | centralRegion = sfaxis[indices] |
394 394 | centralRegionY = self.spheroidalWaveFunction(centralRegion, 0, 0, c, 1) |
395 395 | mysf = np.zeros(len(axis)) |
396 396 | mysf[indices] += centralRegionY / max(centralRegionY) |
397 397 | return mysf |
398 398 | |
399 399 | def spheroidalWaveFunction(self, x, m=0, n=0, c=0, alpha=0): |
400 - | """ |
401 - | Returns spheroidal wave function |
400 + | """Return spheroidal wave function. |
401 + | |
402 402 | Migrated from AnalysisUtils (revision 1.2204, 2015/02/18) |
403 403 | """ |
404 404 | if (type(x) != list and type(x) != np.ndarray): |
405 405 | returnScalar = True |
406 406 | x = [x] |
407 407 | else: |
408 408 | returnScalar = False |
409 409 | cv = scipy.special.pro_cv(m, n, c) # get the eigenvalue |
410 410 | result = scipy.special.pro_ang1_cv(m, n, c, cv, x)[0] |
411 411 | for i in range(len(x)): |
414 414 | # The peak of this function is about 10000 for m=0,n=0,c=6 |
415 415 | if (returnScalar): |
416 416 | return result[0] |
417 417 | else: |
418 418 | return result |
419 419 | |
420 420 | # |
421 421 | # GAUSS |
422 422 | # |
423 423 | def get_gauss_kernel(self, axis, truncate, gwidth, cell_arcsec): |
424 - | """ |
425 - | Returns a gaussian kernel array |
424 + | """Return a gaussian kernel array. |
425 + | |
426 426 | axis : an array of xaxis values |
427 427 | truncate : truncation radius |
428 428 | NOTE definition is different from truncate in gauss()! |
429 429 | gwidth : kernel gwidth |
430 430 | cell_arcsec : image pixel size in unit of arcsec |
431 431 | """ |
432 432 | if gwidth == -1: |
433 433 | gwidth = np.sqrt(np.log(2.0)) |
434 434 | gwidth_arcsec = self.__parse_width(gwidth, cell_arcsec) |
435 435 | # get gauss for full axis |
436 436 | result = self.gauss(axis, [gwidth_arcsec]) |
437 437 | # truncate kernel outside the truncation radius |
438 438 | if truncate == -1: |
439 439 | trunc_arcsec = gwidth_arcsec * 1.5 |
440 440 | elif truncate is not None: |
441 441 | trunc_arcsec = self.__parse_width(truncate, cell_arcsec) |
442 442 | idx = np.where(abs(axis) > trunc_arcsec) |
443 443 | result[idx] = 0. |
444 444 | return result |
445 445 | |
446 446 | def gauss(self, x, parameters): |
447 - | """ |
448 - | Computes the value of the Gaussian function at the specified |
447 + | """Compute the value of the Gaussian function. |
448 + | |
449 + | Compute the value of the Gaussian function at the specified |
449 450 | location(s) with respect to the peak (which is assumed to be at x=0). |
450 451 | truncate: if not None, then set result to zero if below this value. |
451 452 | -Todd Hunter |
452 453 | Migrated from AnalysisUtils (revision 1.2204, 2015/02/18) |
453 454 | """ |
454 455 | if (type(parameters) != np.ndarray and type(parameters) != list): |
455 456 | parameters = np.array([parameters]) |
456 457 | if (len(parameters) < 2): |
457 458 | parameters = np.array([parameters[0], 0]) |
458 459 | fwhm = parameters[0] |
459 460 | x = np.asarray(x, dtype=np.float64) |
460 461 | sigma = fwhm / 2.3548201 |
461 462 | result = np.exp(-(x**2 / (2.0 * sigma**2))) |
462 463 | idx = np.where(result < parameters[1])[0] |
463 464 | result[idx] = 0 |
464 465 | return result |
465 466 | |
466 467 | # |
467 468 | # GJINC |
468 469 | # |
469 470 | def get_gjinc_kernel(self, axis, truncate, gwidth, jwidth, cell_arcsec): |
470 - | """ |
471 - | Returns a GJinc kernel array |
471 + | """Return a GJinc kernel array. |
472 + | |
472 473 | axis : an array of xaxis values |
473 474 | truncate : truncation radius (NOT SUPPORTED YET) |
474 475 | gwidth jwidth : kernel gwidth |
475 476 | cell_arcsec : image pixel size in unit of arcsec |
476 477 | """ |
477 478 | if gwidth == -1: |
478 479 | gwidth = 2.52 * np.sqrt(np.log(2.0)) |
479 480 | if jwidth == -1: |
480 481 | jwidth = 1.55 |
481 482 | gwidth_arcsec = self.__parse_width(gwidth, cell_arcsec) |
482 483 | jwidth_arcsec = self.__parse_width(jwidth, cell_arcsec) |
483 484 | mygjinc = self.trunc(self.gjinc(axis, gwidth=gwidth_arcsec, |
484 485 | jwidth=jwidth_arcsec, |
485 486 | useCasaJinc=True, normalize=False)) |
486 487 | return mygjinc |
487 488 | |
488 489 | def gjinc(self, x, gwidth, jwidth, useCasaJinc=False, normalize=False): |
489 - | """Migrated from AnalysisUtils (revision 1.2204, 2015/02/18)""" |
490 + | """Migrated from AnalysisUtils (revision 1.2204, 2015/02/18).""" |
490 491 | if (useCasaJinc): |
491 492 | result = self.grdjinc1(x, jwidth, normalize) * self.gjincGauss(x, gwidth) |
492 493 | else: |
493 494 | result = self.jinc(x, jwidth) * self.gjincGauss(x, gwidth) |
494 495 | return result |
495 496 | |
496 497 | def grdjinc1(self, val, c, normalize=True): |
497 - | """Migrated from AnalysisUtils (revision 1.2204, 2015/02/18)""" |
498 + | """Migrated from AnalysisUtils (revision 1.2204, 2015/02/18).""" |
498 499 | # Casa's function |
499 500 | # // Calculate J_1(x) using approximate formula |
500 501 | xs = np.pi * val / c |
501 502 | result = [] |
502 503 | for x in xs: |
503 504 | x = abs(x) # I added this to make it symmetric |
504 505 | ax = abs(x) |
505 506 | if (ax < 8.0): |
506 507 | y = x * x |
507 508 | ans1 = x * (72362614232.0 + y * |
527 528 | if (x == 0.0): |
528 529 | out = 0.5 |
529 530 | else: |
530 531 | out = ans / x |
531 532 | if (normalize): |
532 533 | out = out / 0.5 |
533 534 | result.append(out) |
534 535 | return(result) |
535 536 | |
536 537 | def jinc(self, x, jwidth): |
537 - | """ |
538 + | """Compute JINC value. |
539 + | |
538 540 | The peak of this function is 0.5. |
539 541 | Migrated from AnalysisUtils (revision 1.2204, 2015/02/18) |
540 542 | """ |
541 543 | argument = np.pi * np.abs(x) / jwidth |
542 544 | np.seterr(invalid='ignore') # prevent warning for central point |
543 545 | result = scipy.special.j1(argument) / argument |
544 546 | np.seterr(invalid='warn') |
545 547 | for i in range(len(x)): |
546 548 | if (abs(x[i]) < 1e-8): |
547 549 | result[i] = 0.5 |
548 550 | return result |
549 551 | |
550 552 | def gjincGauss(self, x, gwidth): |
551 553 | return (np.exp(-np.log(2) * (x / float(gwidth))**2)) |
552 554 | |
553 555 | # |
554 556 | # Airly disk |
555 557 | # |
556 558 | def get_pb_kernel(self, axis, diam, ref_freq, epsilon=0.0): |
557 - | """ |
559 + | """Return Airy Disk array. |
560 + | |
558 561 | Return Airy Disk array defined by the axis, diameter, reference frequency |
559 562 | and ratio of central hole and antenna diameter |
560 563 | |
561 564 | axis: x-axis values |
562 565 | diameter: antenna diameter in unit of m |
563 566 | reference frequency: the reference frequency of the image |
564 567 | epsilon: ratio of central hole diameter to antenna diameter |
565 568 | """ |
566 569 | a = (self.rootAiryIntensity(axis, epsilon))**2 |
567 570 | airyfwhm = findFWHM(axis, a) |
568 571 | fwhm = primaryBeamArcsec(frequency=ref_freq, diameter=diam, |
569 572 | taper=self.taper, showEquation=False, |
570 573 | obscuration=diam * epsilon, |
571 574 | fwhmfactor=None) |
572 575 | ratio = fwhm / airyfwhm |
573 576 | tempaxis = axis / ratio |
574 577 | a = self.rootAiryIntensity(tempaxis, epsilon) |
575 578 | return a**2 |
576 579 | |
577 580 | def buildAiryDisk(self, fwhm, xaxisLimitInUnitsOfFwhm, convolutionPixelSize, |
578 581 | truncate=False, obscuration=0.75, diameter=12.0): |
579 - | """ |
582 + | """Build airy disk. |
583 + | |
580 584 | This function computes the Airy disk (with peak of 1.0) across a grid of points |
581 585 | specified in units of the FWHM of the disk. |
582 586 | fwhm: a value in arcsec |
583 587 | xaxisLimitInUnitsOfFwhm: an integer or floating point unitless value |
584 588 | truncate: if True, truncate the function at the first null (on both sides) |
585 589 | obscuration: central hole diameter (meters) |
586 590 | diameter: dish diameter (meters) |
587 591 | obscuration and diameter are used to compute the blockage ratio (epsilon) |
588 592 | and its effect on the pattern |
589 593 | Migrated from AnalysisUtils (revision 1.2204, 2015/02/18) |
601 605 | (xaxisLimitInUnitsOfFwhm * fwhm + 0.5 * convolutionPixelSize) / ratio, |
602 606 | convolutionPixelSize / ratio) |
603 607 | a = self.rootAiryIntensity(myxaxis, epsilon) |
604 608 | if (truncate): |
605 609 | a = self.trunc(a) |
606 610 | a = a**2 |
607 611 | myxaxis *= ratio |
608 612 | return(myxaxis, a) |
609 613 | |
610 614 | def rootAiryIntensity(self, myxaxis, epsilon=0.0, showplot=False): |
611 - | """ |
615 + | """Compute the value of 2*J1(x)/x. |
616 + | |
612 617 | This function computes 2*J1(x)/x, which can be squared to get an Airy disk. |
613 618 | myxaxis: the x-axis values to use |
614 619 | epsilon: radius of central hole in units of the dish diameter |
615 620 | Migrated from AnalysisUtils.py (revision 1.2204, 2015/02/18) |
616 621 | """ |
617 622 | if (epsilon > 0): |
618 623 | a = (2 * spspec.j1(myxaxis) / myxaxis - |
619 624 | epsilon**2 * 2 * spspec.j1(myxaxis * epsilon) / |
620 625 | (epsilon * myxaxis)) / (1 - epsilon**2) |
621 626 | else: |
622 627 | a = 2 * spspec.j1(myxaxis) / myxaxis # simpler formula for epsilon=0 |
623 628 | return(a) |
624 629 | |
625 630 | def trunc(self, result): |
626 - | """ |
631 + | """Truncate a list at the first null. |
632 + | |
627 633 | Truncates a list at the first null on both sides of the center, |
628 634 | starting at the center and moving outward in each direction. |
629 635 | Assumes the list is positive in the center, e.g. a Gaussian beam. |
630 636 | -Todd Hunter |
631 637 | Migrated from AnalysisUtils (revision 1.2204, 2015/02/18) |
632 638 | """ |
633 639 | # casa default truncate=-1, which means truncate at radius of first null |
634 640 | mask = np.zeros(len(result)) |
635 641 | truncateBefore = 0 |
636 642 | truncateBeyond = len(mask) |
643 649 | truncateBefore = r |
644 650 | break |
645 651 | mask[truncateBefore:truncateBeyond] = 1 |
646 652 | # casalog.post( |
647 653 | # "Truncating outside of pixels %d-%d (len=%d)" % |
648 654 | # (truncateBefore,truncateBeyond-1,len(mask))) |
649 655 | result *= mask |
650 656 | return result |
651 657 | |
652 658 | def gaussfit_errfunc(self, parameters, x, y): |
653 - | """Migrated from AnalysisUtils (revision 1.2204, 2015/02/18)""" |
659 + | """Migrated from AnalysisUtils (revision 1.2204, 2015/02/18).""" |
654 660 | return (y - self.gauss(x, parameters)) |
655 661 | |
656 662 | def gaussfit(self, x, y, showplot=False, minlevel=0, verbose=False, |
657 663 | title=None, truncate=False): |
658 - | """ |
664 + | """Perform 1D Gaussian fit. |
665 + | |
659 666 | Fits a 1D Gaussian assumed to be centered at x=0 with amp=1 to the |
660 667 | specified data, with an option to truncate it at some level. |
661 668 | Returns the FWHM and truncation point. |
662 669 | Migrated from AnalysisUtils (revision 1.2204, 2015/02/18) |
663 670 | """ |
664 671 | fwhm_guess = findFWHM(x, y) |
665 672 | if (truncate == False): |
666 673 | parameters = np.asarray([fwhm_guess], dtype=np.float64) |
667 674 | else: |
668 675 | parameters = np.asarray([fwhm_guess, truncate], dtype=np.float64) |
694 701 | casalog.post("fitted FWHM = %f" % (fwhm)) |
695 702 | if (truncate != False): |
696 703 | truncate = bestParameters[1] |
697 704 | casalog.post("optimized truncation = %f" % (truncate)) |
698 705 | else: |
699 706 | fwhm = bestParameters |
700 707 | return(fwhm, truncate) |
701 708 | |
702 709 | |
703 710 | def findFWHM(x, y, level=0.5, s=0): |
704 - | """ |
705 - | Measures the FWHM of the specified profile. This works |
706 - | well in a noise-free environment. The data are assumed to |
707 - | be sorted by the x variable. |
711 + | """Measure the FWHM of the specified profile. |
712 + | |
713 + | This works well in a noise-free environment. |
714 + | The data are assumed to be sorted by the x variable. |
708 715 | x: the position variable |
709 716 | y: the intensity variable |
710 717 | level: the signal level for which to find the full width |
711 718 | s: see help scipy.interpolate.UnivariateSpline |
712 719 | -Todd Hunter |
713 720 | Migrated from AnalysisUtils (revision 1.2204, 2015/02/18) |
714 721 | """ |
715 722 | halfmax = np.max(y) * level |
716 723 | spline = scipy.interpolate.UnivariateSpline(x, y - halfmax, s=s) |
717 724 | result = spline.roots() |
728 735 | # ) |
729 736 | # result = 2*findZeroCrossingBySlope(x, y-halfmax) |
730 737 | # return(result) |
731 738 | errmsg = "Unsupported FWHM search in CASA. " + \ |
732 739 | f"More than two corssings ({len(result)}) at level {halfmax} ({level} % of peak)." |
733 740 | raise Exception(errmsg) |
734 741 | |
735 742 | |
736 743 | def primaryBeamArcsec(frequency, diameter, obscuration, taper, |
737 744 | showEquation=True, use2007formula=True, fwhmfactor=None): |
738 - | """ |
739 - | Implements the Baars formula: b*lambda / D. |
745 + | """Implement the Baars formula: b*lambda / D. |
746 + | |
740 747 | if use2007formula==False, use the formula from ALMA Memo 456 |
741 748 | if use2007formula==True, use the formula from Baars 2007 book |
742 749 | (see au.baarsTaperFactor) |
743 750 | In either case, the taper value is expected to be entered as positive. |
744 751 | Note: if a negative value is entered, it is converted to positive. |
745 752 | The effect of the central obstruction on the pattern is also accounted for |
746 753 | by using a spline fit to Table 10.1 of Schroeder's Astronomical Optics. |
747 754 | fwhmfactor: if given, then ignore the taper |
748 755 | For further help and examples, see https://safe.nrao.edu/wiki/bin/view/ALMA/PrimaryBeamArcsec |
749 756 | -- Todd Hunter |
770 777 | else: |
771 778 | formula = "ALMA memo 456 Eq. 18" |
772 779 | casalog.post( |
773 780 | f"Coefficient from {formula} for a -{taper:.1f}dB edge taper " |
774 781 | f"and obscuration ratio={obscuration:g}/{diameter:g} = {b:.3f}*lambda/D") |
775 782 | return(b * lambdaMeters * 3600 * 180 / (diameter * np.pi)) |
776 783 | |
777 784 | |
778 785 | def effectiveTaper(fwhmFactor=1.16, diameter=12, obscuration=0.75, |
779 786 | use2007formula=True): |
780 - | """ |
787 + | """Compute effective taper from constant factor for illumination pattern. |
788 + | |
781 789 | The inverse of (Baars formula multiplied by the central |
782 790 | obstruction factor). Converts an observed value of the constant X in |
783 791 | the formula FWHM=X*lambda/D into a taper in dB (positive value). |
784 792 | if use2007formula == False, use Equation 18 from ALMA Memo 456 |
785 793 | if use2007formula == True, use Equation 4.13 from Baars 2007 book |
786 794 | -- Todd Hunter |
787 795 | Migrated from AnalysisUtils (revision 1.2204, 2015/02/18) |
788 796 | """ |
789 797 | cOF = centralObstructionFactor(diameter, obscuration) |
790 798 | if (fwhmFactor < 1.02 or fwhmFactor > 1.22): |
797 805 | break |
798 806 | else: |
799 807 | increment = -0.01 |
800 808 | for taper_dB in np.arange(10, 10 + increment * 1000, increment): |
801 809 | if (baarsTaperFactor(taper_dB, use2007formula) * cOF - fwhmFactor < 0): |
802 810 | break |
803 811 | return(taper_dB) |
804 812 | |
805 813 | |
806 814 | def baarsTaperFactor(taper_dB, use2007formula=True): |
807 - | """ |
815 + | """Convert a taper to the constant factor for illumination pattern. |
816 + | |
808 817 | Converts a taper in dB to the constant X |
809 818 | in the formula FWHM=X*lambda/D for the parabolic illumination pattern. |
810 819 | We assume that taper_dB comes in as a positive value. |
811 820 | use2007formula: False --> use Equation 18 from ALMA Memo 456. |
812 821 | True --> use Equation 4.13 from Baars 2007 book |
813 822 | - Todd Hunter |
814 823 | Migrated from AnalysisUtils (revision 1.2204, 2015/02/18) |
815 824 | """ |
816 825 | tau = 10**(-0.05 * taper_dB) |
817 826 | if (use2007formula): |
818 827 | return(1.269 - 0.566 * tau + 0.534 * (tau**2) - 0.208 * (tau**3)) |
819 828 | else: |
820 829 | return(1.243 - 0.343 * tau + 0.12 * (tau**2)) |
821 830 | |
822 831 | |
823 832 | def centralObstructionFactor(diameter=12.0, obscuration=0.75): |
824 - | """ |
833 + | """Compute the scale factor of an Airy pattern. |
834 + | |
825 835 | Computes the scale factor of an Airy pattern as a function of the |
826 836 | central obscuration, using Table 10.1 of Schroeder's 'Astronomical Optics'. |
827 837 | -- Todd Hunter |
828 838 | Migrated from AnalysisUtils (revision 1.2204, 2015/02/18) |
829 839 | """ |
830 840 | epsilon = obscuration / diameter |
831 841 | myspline = scipy.interpolate.UnivariateSpline( |
832 842 | [0, 0.1, 0.2, 0.33, 0.4], [1.22, 1.205, 1.167, 1.098, 1.058], s=0) |
833 843 | factor = myspline(epsilon) / 1.22 |
834 844 | if (type(factor) == np.float64): |