Coverage for frank/utilities.py: 93%

1# Frankenstein: 1D disc brightness profile reconstruction from Fourier data 

2# using non-parametric Gaussian Processes 

3# 

4# Copyright (C) 2019-2020 R. Booth, J. Jennings, M. Tazzari 

5# 

6# This program is free software: you can redistribute it and/or modify 

7# it under the terms of the GNU General Public License as published by 

8# the Free Software Foundation, either version 3 of the License, or 

9# (at your option) any later version. 

10# 

11# This program is distributed in the hope that it will be useful, 

12# but WITHOUT ANY WARRANTY; without even the implied warranty of 

13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

14# GNU General Public License for more details. 

15# 

16# You should have received a copy of the GNU General Public License 

17# along with this program. If not, see <https://www.gnu.org/licenses/> 

18# 

19"""This module has various functions useful for plotting and analyzing fit 

20results. 

21""" 

22import logging 

23import numpy as np 

24from scipy.interpolate import interp1d 

25 

26from frank.constants import deg_to_rad, sterad_to_arcsec, rad_to_arcsec 

27from frank.geometry import FixedGeometry 

28from frank.hankel import DiscreteHankelTransform 

29from frank.statistical_models import VisibilityMapping 

30 

31def arcsec_baseline(x): 

32 """ 

33 Provide x as a radial scale [arcsec] to return the corresponding baseline 

34 [lambda], or vice-versa 

35 

36 Parameters 

37 ---------- 

38 x : float 

39 Radial scale [arcsec] or baseline [lambda] 

40 

41 Returns 

42 ------- 

43 converted : float 

44 Baseline [lambda] or radial scale [arcsec] 

45 

46 """ 

47 

48 converted = 1 / (x / 60 / 60 * np.pi / 180) 

49 

50 return converted 

51 

52 

53def radius_convert(x, dist, conversion='arcsec_au'): 

54 """ 

55 Provide x as a radius/radii in [arcsec] to convert to [au] (or vice-versa), 

56 assuming a distance in [pc] 

57 

58 Parameters 

59 ---------- 

60 x : float or array 

61 Radius or radii ([arcsec] or [au]) 

62 dist : float 

63 Distance to source [pc] 

64 conversion : {'arcsec_au', 'au_arcsec'}, default = 'arcsec_au' 

65 The unit conversion to perform, e.g. 'arcsec_au' converts from [arcsec] 

66 to [au] 

67 

68 Returns 

69 ------- 

70 converted : float or array 

71 Radius or radii ([au] or [arcsec]) 

72 

73 """ 

74 

75 if conversion == 'arcsec_au': 

76 converted = x * dist 

77 elif conversion == 'au_arcsec': 

78 converted = x / dist 

79 else: 

80 raise AttributeError("conversion must be one of {}" 

81 "".format(['arcsec_au', 'au_arcsec'])) 

82 

83 return converted 

84 

85 

86def jy_convert(x, conversion, bmaj=None, bmin=None): 

87 """ 

88 Provide x as a brightness in one of the units [Jy / beam], [Jy / arcsec^2], 

89 [Jy / sterad] to convert x to another of these units 

90 

91 Parameters 

92 ---------- 

93 x : array 

94 Brightness in one of: [Jy / beam], [Jy / arcsec^2], [Jy / sterad] 

95 conversion : { 'beam_sterad', 'beam_arcsec2', 'arcsec2_beam', 

96 'arcsec2_sterad', 'sterad_beam', 'sterad_arcsec2'} 

97 The unit conversion to perform, e.g., 'beam_sterad' converts 

98 [Jy / beam] to [Jy / sterad] 

99 bmaj : float, optional 

100 Beam FWHM along the major axis [arcsec] 

101 bmin : float, optional 

102 Beam FWHM along the minor axis [arcsec] 

103 

104 Returns 

105 ------- 

106 converted : float 

107 Brightness in converted units 

108 

109 """ 

110 if (bmaj is None or bmin is None) and conversion in ['beam_sterad', 

111 'beam_arcsec2', 

112 'arcsec2_beam', 

113 'sterad_beam']: 

114 raise ValueError('bmaj and bmin must be specified to perform the' 

115 ' conversion {}'.format(conversion)) 

116 

117 if bmaj is not None and bmin is not None: 

118 beam_solid_angle = np.pi * bmaj * bmin / (4 * np.log(2)) 

119 

120 if conversion == 'beam_arcsec2': 

121 converted = x / beam_solid_angle 

122 elif conversion == 'arcsec2_beam': 

123 converted = x * beam_solid_angle 

124 elif conversion == 'arcsec2_sterad': 

125 converted = x * sterad_to_arcsec 

126 elif conversion == 'sterad_arcsec2': 

127 converted = x / sterad_to_arcsec 

128 elif conversion == 'beam_sterad': 

129 converted = x / beam_solid_angle * sterad_to_arcsec 

130 elif conversion == 'sterad_beam': 

131 converted = x * beam_solid_angle / sterad_to_arcsec 

132 else: 

133 raise AttributeError("conversion must be one of {}" 

134 "".format(['beam_sterad', 'beam_arcsec2', 

135 'arcsec2_beam', 'arcsec2_sterad', 

136 'sterad_beam', 'sterad_arcsec2'])) 

137 

138 return converted 

139 

140 

141def get_fit_stat_uncer(fit): 

142 r""" 

143 Retrieve a model brightness profile's statistical uncertainty (that due to  

144 the uncertainty on the observed baselines). This is only a lower bound on  

145 the total uncertainty (which also has a systematic contribution due to  

146 sparse sampling of the (u, v) plane). 

147 

148 Parameters 

149 ---------- 

150 fit : FrankFitter result object 

151 Fitted profile to obtain statistical uncertainty of  

152 

153 Returns 

154 ------- 

155 sigma: array 

156 Statistical uncertainty on the fitted brightness, at the fitted radial 

157 points. 

158 """ 

159 

160 if 'method' not in fit._info.keys(): 

161 raise AttributeError("'fit' object lacks '_info.method' key.") 

162 

163 if fit._info["method"] == "Normal": 

164 sigma = np.sqrt(np.diag(fit.covariance)) 

165 

166 elif fit._info["method"] == "LogNormal": 

167 # convert stored variance from var(log(brightness)) to var(brightness). 

168 # see Appendix A.2 in https://ui.adsabs.harvard.edu/abs/2016A%26A...586A..76J/abstract; 

169 # https://en.wikipedia.org/wiki/Log-normal_distribution 

170 log_variance = np.diag(fit.covariance) 

171 lin_variance = (np.exp(log_variance) - 1) * np.exp(2 * np.log(fit.I)) 

172 sigma = np.sqrt(lin_variance) 

173 

174 else: 

175 raise ValueError("fit._info['method'] {} must be one of " 

176 "['Normal', 'LogNormal']".format(fit._info['method'])) 

177 

178 return sigma 

179 

180 

181class UVDataBinner(object): 

182 r""" 

183 Average uv-data into bins of equal size. 

184 

185 Compute the weighted mean of the visibilities in each bin 

186 

187 Parameters 

188 ---------- 

189 uv : array, unit = :math:`\lambda` 

190 Baselines of the data to bin 

191 V : array, unit = Jy 

192 Observed visibility. If complex, both the real and imaginary 

193 components will be binned. Else only the real part will be binned. 

194 weights : array, unit = Jy^-2 

195 Weights on the visibility points 

196 bin_width : float, unit = :math:`\lambda` 

197 Width of the uv-bins 

198 

199 Notes 

200 ----- 

201 Uses numpy masked arrays to mask bins with no uv points. 

202 

203 """ 

204 

205 def __init__(self, uv, V, weights, bin_width): 

206 nbins = np.ceil(uv.max() / bin_width).astype('int') 

207 #Protect against rounding 

208 if nbins * bin_width < uv.max(): 

209 nbins += 1 

210 bins = np.arange(nbins+1, dtype='float64') * bin_width 

211 

212 self._bins = bins 

213 self._nbins = nbins 

214 self._norm = 1 / bin_width 

215 

216 bin_uv, bin_wgt, bin_vis, bin_n = \ 

217 self.bin_quantities(uv, weights, 

218 uv, np.ones_like(uv), V, 

219 bin_counts=True) 

220 

221 # Normalize 

222 idx = bin_n > 0 

223 w = bin_wgt[idx] 

224 bin_uv[idx] /= w 

225 bin_vis[idx] /= w 

226 

227 # Store the binned data, masking empty bins: 

228 mask = (bin_n == 0) 

229 self._uv = np.ma.masked_where(mask, bin_uv) 

230 self._V = np.ma.masked_where(mask, bin_vis) 

231 self._w = np.ma.masked_where(mask, bin_wgt) 

232 self._count = np.ma.masked_where(mask, bin_n) 

233 

234 self._uv_left = np.ma.masked_where(mask, bins[:-1]) 

235 self._uv_right = np.ma.masked_where(mask, bins[1:]) 

236 

237 # Compute the uncertainty on the means: 

238 bin_vis_err = np.full(nbins, np.nan, dtype=V.dtype) 

239 

240 # 1) Get the binned mean for each vis: 

241 mu = self._V[self.determine_uv_bin(uv)] 

242 

243 # 2) Compute weighted error for bins with n > 1 

244 quantities = (V-mu).real**2 

245 if np.iscomplexobj(V): 

246 quantities = quantities + 1j * (V-mu).imag**2 

247 err = self.bin_quantities(uv, weights**2, quantities) 

248 

249 idx2 = bin_n > 1 

250 err[idx2] /= bin_wgt[idx2]**2 * (1 - 1 / bin_n[idx2]) 

251 

252 temp_error = np.sqrt(err.real[idx2]) 

253 if np.iscomplexobj(V): 

254 temp_error = temp_error + 1.j * np.sqrt(err.imag[idx2]) 

255 bin_vis_err[idx2] = temp_error 

256 

257 # 3) Use a sensible error for bins with one baseline 

258 idx1 = bin_n == 1 

259 bin_vis_err[idx1].real = 1 / np.sqrt(bin_wgt[idx1]) 

260 if np.iscomplexobj(V): 

261 bin_vis_err[idx1].imag = bin_vis_err[idx1].real 

262 bin_vis_err[mask] = np.nan 

263 

264 # 4) Store the error 

265 self._Verr = np.ma.masked_where(mask, bin_vis_err) 

266 

267 

268 def determine_uv_bin(self, uv): 

269 r"""Determine the bin that the given uv points belong too. 

270 

271 Parameters 

272 ---------- 

273 uv : array, unit = :math:`\lambda` 

274 Baselines to determine the bins of 

275 

276 Returns 

277 ------- 

278 idx : array, 

279 Bins that the uv point belongs to. Will be -1 if the bin does not 

280 exist. 

281 

282 """ 

283 

284 bins = self._bins 

285 nbins = self._nbins 

286 

287 idx = np.floor(uv * self._norm).astype('int32') 

288 idx[uv < bins[idx]] -= 1 # Fix rounding error 

289 idx[uv == bins[nbins]] -= 1 # Handle point exaclty on outer boundary 

290 

291 too_high = idx >= nbins 

292 idx[too_high] = -1 

293 

294 # Increase points on the RH boundary, unless it is the outer one 

295 idx_tmp = idx[~too_high] 

296 increment = (uv[~too_high] >= bins[idx_tmp+1]) & (idx_tmp+1 < nbins) 

297 idx_tmp[increment] += 1 

298 

299 return idx 

300 

301 def bin_quantities(self, uv, w, *quantities, bin_counts=False): 

302 r"""Bin the given quantities according the to uv points and weights. 

303 

304 Parameters 

305 ---------- 

306 uv : array, unit = :math:`\lambda` 

307 Baselines of the data to bin 

308 weights : array, unit = Jy^-2 

309 Weights on the visibility points 

310 quantities : arrays, 

311 Quantities evaluated at the uv points to bin. 

312 bin_counts : bool, default = False 

313 Determines whether to count the number of uv points per bin. 

314 

315 Returns 

316 ------- 

317 results : arrays, same type as quantities 

318 Binned data 

319 bin_counts : array, int64, optional 

320 If bin_counts=True, then this array contains the number of uv 

321 points in each bin. Otherwise, it is not returned. 

322 

323 """ 

324 

325 bins = self._bins 

326 nbins = self._nbins 

327 norm = self._norm 

328 

329 results = [np.zeros(nbins, dtype=x.dtype) for x in quantities] 

330 

331 if bin_counts: 

332 counts = np.zeros(nbins, dtype='int64') 

333 

334 BLOCK = 65536 

335 for i in range(0, len(uv), BLOCK): 

336 tmp_uv = uv[i:i+BLOCK] 

337 tmp_wgt = w[i:i+BLOCK] 

338 

339 idx = np.floor(tmp_uv * norm).astype('int32') 

340 

341 # Fix rounding: 

342 idx[tmp_uv < bins[idx]] -= 1 

343 # Move points exactly on the outer boundary inwards: 

344 idx[idx == nbins] -= 1 

345 

346 # Only the last bin includes its edge 

347 increment = (tmp_uv >= bins[idx+1]) & (idx+1 != nbins) 

348 idx[increment] += 1 

349 

350 if bin_counts: 

351 counts += np.bincount(idx, minlength=nbins) 

352 

353 for qty, res in zip(quantities, results): 

354 tmp_qty = tmp_wgt*qty[i:i+BLOCK] 

355 

356 if np.iscomplexobj(qty): 

357 res.real += np.bincount(idx, weights=tmp_qty.real, 

358 minlength=nbins) 

359 res.imag += np.bincount(idx, weights=tmp_qty.imag, 

360 minlength=nbins) 

361 else: 

362 res += np.bincount(idx, weights=tmp_qty, minlength=nbins) 

363 

364 if bin_counts: 

365 return results + [counts] 

366 if len(results) == 1: 

367 return results[0] 

368 return results 

369 

370 def __len__(self): 

371 return len(self._uv) 

372 

373 @property 

374 def uv(self): 

375 r"""Binned uv points, unit = :math:`\lambda`""" 

376 return self._uv 

377 

378 @property 

379 def V(self): 

380 """Binned visibility, unit = Jy""" 

381 return self._V 

382 

383 @property 

384 def weights(self): 

385 """Binned weights, unit = Jy^-2""" 

386 return self._w 

387 

388 @property 

389 def error(self): 

390 """Uncertainty on the binned visibilities, unit = Jy""" 

391 return self._Verr 

392 

393 @property 

394 def bin_counts(self): 

395 """Number of points in each bin""" 

396 return self._count 

397 

398 @property 

399 def bin_edges(self): 

400 """Edges of the histogram bins""" 

401 return [self._uv_left, self._uv_right] 

402 

403 

404def normalize_uv(u, v, wle): 

405 r""" 

406 Normalize data u and v coordinates by the observing wavelength 

407 

408 Parameters 

409 ---------- 

410 u, v : array, unit = [m] 

411 u and v coordinates of observations 

412 wle : float or array, unit = [m] 

413 Observing wavelength of observations. If an array, it should be the 

414 pointwise wavelength for each (u,v) point 

415 

416 Returns 

417 ------- 

418 u_normed, v_normed : array, unit = :math:`\lambda` 

419 u and v coordinates normalized by observing wavelength 

420 

421 """ 

422 

423 logging.info(' Normalizing u and v coordinates by provided' 

424 ' observing wavelength of {} m'.format(wle)) 

425 

426 wle = np.atleast_1d(wle).astype('f8') 

427 if len(wle) != 1 and len(wle) != len(u): 

428 raise ValueError("len(wle) = {}. It should be equal to " 

429 "len(u) = {} (or 1 if all wavelengths are the same)".format(len(wle), len(u))) 

430 u_normed = u / wle 

431 v_normed = v / wle 

432 

433 return u_normed, v_normed 

434 

435 

436def cut_data_by_baseline(u, v, vis, weights, cut_range, geometry=None): 

437 r""" 

438 Truncate the data to be within a chosen baseline range. 

439 

440 The cut will be done in deprojected baseline space if the geometry is 

441 provided. 

442 

443 Parameters 

444 ---------- 

445 u, v : array, unit = :math:`\lambda` 

446 u and v coordinates of observations 

447 vis : array, unit = Jy 

448 Observed visibilities (complex: real + imag * 1j) 

449 weights : array, unit = Jy^-2 

450 Weights assigned to observed visibilities, of the form 

451 :math:`1 / \sigma^2` 

452 cut_range : list of float, length = 2, unit = [\lambda] 

453 Lower and upper baseline bounds outside of which visibilities are 

454 truncated 

455 geometry : SourceGeometry object, optional 

456 Fitted geometry (see frank.geometry.SourceGeometry). 

457 

458 

459 Returns 

460 ------- 

461 u_cut, v_cut : array, unit = :math:`\lambda` 

462 u and v coordinates in the chosen baseline range 

463 vis_cut : array, unit = Jy 

464 Visibilities in the chosen baseline range 

465 weights_cut : array, unit = Jy^-2 

466 Weights in the chosen baseline range 

467 

468 """ 

469 

470 logging.info(' Cutting data outside of the minimum and maximum baselines' 

471 ' of {} and {}' 

472 ' klambda'.format(cut_range[0] / 1e3, 

473 cut_range[1] / 1e3)) 

474 if geometry is not None: 

475 up, vp = geometry.deproject(u, v) 

476 else: 

477 up, vp = u, v 

478 

479 baselines = np.hypot(up, vp) 

480 above_lo = baselines >= cut_range[0] 

481 below_hi = baselines <= cut_range[1] 

482 in_range = above_lo & below_hi 

483 u_cut, v_cut, vis_cut, weights_cut = [x[in_range] for x in [u, v, vis, weights]] 

484 

485 return u_cut, v_cut, vis_cut, weights_cut 

486 

487 

488def estimate_weights(u, v=None, V=None, nbins=300, log=True, use_median=False, 

489 verbose=True): 

490 r""" 

491 Estimate the weights using the variance of the binned visibilities. 

492 

493 The estimation is done assuming that the variation in each bin is dominated 

494 by the noise. This will be true if: 

495 1) The source is axi-symmetric, 

496 2) The uv-points have been deprojected, 

497 3) The bins are not too wide, 

498 Otherwise the variance may be dominated by intrinsic variations in the 

499 visibilities. 

500 

501 Parameters 

502 ---------- 

503 u, v : array, unit = :math:`\lambda` 

504 u and v coordinates of observations (deprojected). Data will be binned 

505 by baseline. If v is not None, np.hpot(u,v) will be used instead. Note 

506 that if V is None the argument v will be intepreted as V instead 

507 V : array, unit = Jy, default = None 

508 Observed visibility. If complex, the weights will be computed from the 

509 average of the variance of the real and imaginary components, as in 

510 CASA's statwt. Otherwise the variance of the real part is used. 

511 nbins : int, default = 300 

512 Number of bins used. 

513 log : bool, default = True 

514 If True, the uv bins will be constructed in log space, otherwise linear 

515 spaced bins will be used. 

516 use_median : bool, default = False 

517 If True all of the weights will be set to the median of the variance 

518 estimated across the bins. Otherwise, the baseline dependent variance 

519 will be used. 

520 verbose : bool, default = True 

521 If true, the logger will record calls to this function, along with 

522 whether the median estimate was used. 

523 

524 Returns 

525 ------- 

526 weights : array 

527 Estimate of the weight for each uv point. 

528 

529 Notes 

530 ----- 

531 - This function does not use the original weights in the estimation. 

532 - Bins with only one uv point do not have a variance estimate. Thus 

533 the mean of the variance in the two adjacent bins is used instead. 

534 

535 Examples 

536 -------- 

537 All of the following calls will work as expected: 

538 `estimate_weights(u, v, V) ` 

539 `estimate_weights(u, V)` 

540 `estimate_weights(u, V=V)` 

541 In each case the variance of V in the uv-bins is used to estimate the 

542 weights. The first call will use q = np.hypot(u, v) in the uv-bins. The 

543 second and third calls are equivalent to the first with u=0. 

544 

545 """ 

546 

547 if verbose: 

548 logging.info(' Estimating visibility weights') 

549 

550 if V is None: 

551 if v is not None: 

552 V = v 

553 q = np.abs(u) 

554 else: 

555 raise ValueError("The visibilities, V, must be supplied") 

556 elif v is not None: 

557 q = np.hypot(u,v) 

558 else: 

559 q = np.abs(u) 

560 

561 if log: 

562 q = np.log(q) 

563 q -= q.min() 

564 

565 bin_width = (q.max()-q.min()) / nbins 

566 

567 uvBin = UVDataBinner(q, V, np.ones_like(q), bin_width) 

568 

569 if uvBin.bin_counts.max() == 1: 

570 raise ValueError("No bin contains more than one uv point, can't" 

571 " estimate the variance. Use fewer bins.") 

572 

573 if np.iscomplex(V.dtype): 

574 var = 0.5*(uvBin.error.real**2 + uvBin.error.imag**2) * uvBin.bin_counts 

575 else: 

576 var = uvBin.error.real**2 * uvBin.bin_counts 

577 

578 if use_median: 

579 if verbose: 

580 logging.info(' Setting all weights as median binned visibility ' 

581 'variance') 

582 return np.full(len(u), 1/np.ma.median(var[uvBin.bin_counts > 1])) 

583 else: 

584 if verbose: 

585 logging.info(' Setting weights according to baseline-dependent ' 

586 'binned visibility variance') 

587 # For bins with 1 uv point, use the average of the adjacent bins 

588 no_var = np.argwhere(uvBin.bin_counts == 1).reshape(-1) 

589 if len(no_var) > 0: 

590 # Find the location `loc` of the bad points in the array of good points 

591 good_var = np.argwhere(uvBin.bin_counts > 1).reshape(-1) 

592 loc = np.searchsorted(good_var, no_var, side='right') 

593 

594 # Set the variance to the average of the two adjacent bins 

595 im = good_var[np.maximum(loc-1, 0)] 

596 ip = good_var[np.minimum(loc, len(good_var)-1)] 

597 var[no_var] = 0.5*(var[im] + var[ip]) 

598 

599 bin_id = uvBin.determine_uv_bin(q) 

600 assert np.all(bin_id != -1), "Error in binning" # Should never occur 

601 

602 weights = 1/var[bin_id] 

603 

604 return weights 

605 

606 

607def draw_bootstrap_sample(u, v, vis, weights): 

608 r""" 

609 Obtain the sample for a bootstrap, drawing, with replacement, N samples from 

610 a length N dataset 

611 

612 Parameters 

613 ---------- 

614 u, v : array, unit = :math:`\lambda` 

615 u and v coordinates of observations 

616 vis : array, unit = Jy 

617 Observed visibilities (complex: real + imag * 1j) 

618 weights : array, unit = Jy^-2 

619 Weights on the visibilities, of the form 

620 :math:`1 / \sigma^2` 

621 

622 Returns 

623 ------- 

624 u_boot, v_boot : array, unit = :math:`\lambda` 

625 Bootstrap sampled u and v coordinates 

626 vis_boot : array, unit = Jy 

627 Bootstrap sampled visibilities 

628 weights_boot : array, unit = Jy^-2 

629 Boostrap sampled weights on the visibilities 

630 

631 """ 

632 idxs = np.random.randint(low=0, high=len(u), size=len(u)) 

633 

634 u_boot = u[idxs] 

635 v_boot = v[idxs] 

636 vis_boot = vis[idxs] 

637 weights_boot = weights[idxs] 

638 

639 return u_boot, v_boot, vis_boot, weights_boot 

640 

641 

642def sweep_profile(r, I, project=False, phase_shift=False, geom=None, axis=0, 

643 xmax=None, ymax=None, dr=None): 

644 r""" 

645 Sweep a 1D radial brightness profile over :math:`2 \pi` to yield a 2D 

646 brightness distribution. Optionally project this sweep by a supplied 

647 geometry. 

648 

649 Parameters 

650 ---------- 

651 r : array 

652 Radial coordinates at which the 1D brightness profile is defined 

653 I : array 

654 Brightness values at r 

655 project : bool, default = False 

656 Whether to project the swept profile by the supplied geom 

657 phase_shift : bool, default = False 

658 Whether to phase shift the projected profile by the supplied geom. 

659 If False, the source will be centered in the image 

660 geom : SourceGeometry object, default=None 

661 Fitted geometry (see frank.geometry.SourceGeometry). Here we use 

662 geom.inc [deg], geom.PA [deg], geom.dRA [arcsec], geom.dDec [arcsec] if 

663 project=True 

664 axis : int, default = 0 

665 Axis over which to interpolate the 1D profile 

666 xmax, ymax : float, optional, default = None 

667 Value setting the x- and y-bounds of the image (same units as r). The 

668 positive and negative bounds are both set to this value (modulo sign). 

669 If not provided, these will be set to r.max() 

670 dr : float, optional, default = None 

671 Pixel size (same units as r). If not provided, it will be set at the 

672 same spatial scale as r 

673 

674 Returns 

675 ------- 

676 I2D : array, shape = (len(r), len(r)) 

677 2D brightness distribution (projected if project=True) 

678 xmax : float 

679 Maximum x-value of the 2D grid 

680 ymax : float 

681 Maximum y-value of the 2D grid 

682 

683 Notes 

684 ----- 

685 Sign convention: a negative geom.dRA shifts the source to the right 

686 in the image 

687 

688 """ 

689 if project: 

690 inc, pa = geom.inc, geom.PA 

691 inc *= deg_to_rad 

692 pa *= deg_to_rad 

693 

694 cos_i = np.cos(inc) 

695 cos_pa, sin_pa = np.cos(pa), np.sin(pa) 

696 

697 if xmax is None: 

698 xmax = r.max() 

699 if ymax is None: 

700 ymax = r.max() 

701 

702 if dr is None: 

703 dr = np.mean(np.diff(r)) 

704 

705 x = np.linspace(-xmax, xmax, int(xmax/dr) + 1) 

706 y = np.linspace(-ymax, ymax, int(ymax/dr) + 1) 

707 

708 if phase_shift: 

709 xi, yi = np.meshgrid(x + geom.dRA, y - geom.dDec) 

710 else: 

711 xi, yi = np.meshgrid(x, y) 

712 

713 if project: 

714 xp = xi * cos_pa + yi * sin_pa 

715 yp = -xi * sin_pa + yi * cos_pa 

716 xp /= cos_i 

717 r1D = np.hypot(xp, yp) 

718 else: 

719 r1D = np.hypot(xi, yi) 

720 

721 im_shape = r1D.shape + I.shape[1:] 

722 

723 interp = interp1d(r, I, bounds_error=False, fill_value=0., axis=axis) 

724 I2D = interp(r1D.ravel()).reshape(*im_shape) 

725 

726 return I2D, xmax, ymax 

727 

728 

729def make_image(fit, Npix, xmax=None, ymax=None, project=True): 

730 """Make an image of a model fit. 

731  

732 Parameters 

733 ---------- 

734 fit : FrankFitter result object 

735 Fitted profile to make an image of 

736 Npix : int or list 

737 Number of pixels in the x-direction, or [x-,y-] direction 

738 xmax: float or None, unit=arcsec 

739 Size of the image is [-xmax, xmax]. By default this is twice 

740 fit.Rmax to avoid aliasing. 

741 ymax: float or None, unit=arcsec 

742 Size of the image is [-ymax,ymax]. Defaults to xmax if ymax=None 

743 project: bool, default=True 

744 Whether to produce a projected image. 

745 

746 Returns 

747 ------- 

748 x : array, 1D; unit=arcsec 

749 Locations of the x-points in the image. 

750 y : array, 1D; unit=arcsec 

751 Locations of the y-points in the image. 

752 I : array, 2D; unit=Jy/Sr 

753 Image of the surface brightness. 

754 """ 

755 if xmax is None: 

756 xmax = 2*fit.Rmax 

757 if ymax is None: 

758 ymax = xmax 

759 

760 try: 

761 Nx, Ny = Npix 

762 except TypeError: 

763 Nx = Npix 

764 Ny = int(Nx*(ymax/xmax)) 

765 

766 dx = 2*xmax/(Nx*rad_to_arcsec) 

767 dy = 2*ymax/(Ny*rad_to_arcsec) 

768 

769 

770 # The easiest way to produce an image is to predict the visibilities 

771 # on a regular grid in the Fourier plane and then transform it back. 

772 # All frank models must be able to compute the visibilities so this 

773 # method is completely general. 

774 u = np.fft.fftfreq(Nx)/dx 

775 v = np.fft.fftfreq(Ny)/dy 

776 u, v = np.meshgrid(u,v, indexing='ij') 

777 

778 # Get the visibilities: 

779 if project: 

780 Vis = fit.predict(u.reshape(-1), v.reshape(-1)).reshape(*u.shape) 

781 else: 

782 q = np.hypot(u,v) 

783 Vis = fit.predict_deprojected(q.reshape(-1)).reshape(*u.shape) 

784 

785 # Convert to the image plane 

786 I = np.fft.ifft(np.fft.ifft(Vis, axis=0), axis=1).real 

787 I /= dx*dy 

788 

789 # numpy's fft has zero in the corner. We want it in the middle so we need 

790 # to wrap: 

791 tmp = I.copy() 

792 tmp[:Nx//2,], tmp[Nx//2:] = I[Nx//2:], I[:Nx//2] 

793 I[:,:Ny//2], I[:,Ny//2:] = tmp[:,Ny//2:], tmp[:,:Ny//2] 

794 

795 xe = np.linspace(-xmax, xmax, Nx+1) 

796 x = 0.5*(xe[1:] + xe[:-1]) 

797 ye = np.linspace(-ymax, ymax, Ny+1) 

798 y = 0.5*(ye[1:] + ye[:-1]) 

799 

800 return x, y, I 

801 

802 

803def convolve_profile(r, I, disc_i, disc_pa, clean_beam, 

804 n_per_sigma=5, axis=0): 

805 r""" 

806 Convolve a 1D radial brightness profile with a 2D Gaussian beam, degrading 

807 the profile's resolution 

808 

809 Parameters 

810 ---------- 

811 r : array 

812 Radial coordinates at which the 1D brightness profile is defined 

813 I : array 

814 Brightness values at r 

815 disc_i : float, unit = deg 

816 Disc inclination 

817 disc_pa : float, unit = deg 

818 Disc position angle 

819 clean_beam : dict 

820 Dictionary with beam `bmaj` (FWHM of beam along its major axis) [arcsec], 

821 `bmin` (FWHM of beam along its minor axis) [arcsec], 

822 `pa` (beam position angle) [deg] 

823 n_per_sigma : int, default = 50 

824 Number of points per standard deviation of the Gaussian kernel (used 

825 for gridding) 

826 axis : int, default = 0 

827 Axis over which to interpolate the 1D profile 

828 

829 Returns 

830 ------- 

831 I_smooth : array, shape = (len(r), len(r)) 

832 Convolved brightness profile I at coordinates r 

833 

834 """ 

835 

836 from scipy.constants import c 

837 from scipy.interpolate import interp1d 

838 from scipy.ndimage import gaussian_filter 

839 

840 # Set up the geometry for the smoothing mesh. 

841 # We align the beam with the grid (major axis aligned) and rotate the 

842 # image accordingly 

843 

844 # Convert beam FWHM to sigma 

845 bmaj = clean_beam['bmaj'] / np.sqrt(8 * np.log(2)) 

846 bmin = clean_beam['bmin'] / np.sqrt(8 * np.log(2)) 

847 

848 PA = (disc_pa - clean_beam['beam_pa']) * np.pi / 180. 

849 

850 cos_i = np.cos(disc_i * np.pi/180.) 

851 cos_PA = np.cos(PA) 

852 sin_PA = np.sin(PA) 

853 

854 # Pixel size in terms of bmin 

855 rmax = r.max() 

856 dx = bmin / n_per_sigma 

857 

858 xmax = rmax * (np.abs(cos_i * cos_PA) + np.abs(sin_PA)) 

859 ymax = rmax * (np.abs(cos_i * sin_PA) + np.abs(cos_PA)) 

860 

861 xmax = int(xmax / dx + 1) * dx 

862 ymax = int(ymax / dx + 1) * dx 

863 

864 x = np.arange(-xmax, xmax + dx / 2, dx) 

865 y = np.arange(-ymax, ymax + dx / 2, dx) 

866 

867 xi, yi = np.meshgrid(x, y) 

868 

869 xp = xi * cos_PA + yi * sin_PA 

870 yp = -xi * sin_PA + yi * cos_PA 

871 xp /= cos_i 

872 

873 r1D = np.hypot(xi, yi) 

874 

875 im_shape = r1D.shape + I.shape[1:] 

876 

877 # Interpolate to grid and apply smoothing 

878 interp = interp1d(r, I, bounds_error=False, fill_value=0., axis=axis) 

879 

880 I2D = interp(r1D.ravel()).reshape(*im_shape) 

881 sigma = [float(n_per_sigma), (bmaj / bmin) * n_per_sigma] 

882 I2D = gaussian_filter(I2D, sigma, mode='nearest', cval=0.) 

883 

884 # Now convert back to a 1D profile 

885 edges = np.concatenate(([r[0] * r[0] / r[1]], r, [r[-1] * r[-1] / r[-2]])) 

886 edges = 0.5 * (edges[:-1] + edges[1:]) 

887 

888 I_smooth = np.empty_like(I) 

889 I_smooth = np.histogram(r1D.ravel(), weights=I2D.ravel(), bins=edges)[0]