Coverage for frank/geometry.py: 92%

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

2# using non-parametric Gaussian Processes

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 contains methods for fitting a source's geometry and deprojecting

20 the visibilties by a fitted or a given geometry.

22 NOTE: The sign convention used here is east of north for position angle (PA)

23 and positive toward east for right ascension offset (dRA).

24"""

26import numpy as np

27from scipy.optimize import least_squares

28import logging

30from frank.constants import rad_to_arcsec, deg_to_rad

31from frank.radial_fitters import FourierBesselFitter

33def _fix_inc_and_PA_ranges(inc, PA):

34 """Make sure the inclination and PA are in the ranges [0,90] and [0-180]."""

35 inc = inc % 180

36 PA = PA % 180

37 if inc > 90:

38 inc = 180 - inc

39 return inc, PA

41def apply_phase_shift(u, v, V, dRA, dDec, inverse=False):

42 r"""

43 Apply a phase shift to the visibilities.

45 This is equivalent to moving the source in the image plane by the

46 vector (dRA, dDec).

48 Parameters

49 ----------

50 u : array of real, size = N, unit = :math:`\lambda`

51 u-points of the visibilities

52 v : array of real, size = N, unit = :math:`\lambda`

53 v-points of the visibilities

54 V : array of real, size = N, unit = Jy

55 Complex visibilites

56 dRA : float, unit = arcsec

57 Phase shift in right ascenion.

58 dDec : float, unit = arcsec

59 Phase shift in declination.

60 inverse : bool, default=False

61 If True, the phase shift is reversed (equivalent to

62 flipping the signs of dRA and dDec).

63 Returns

64 -------

65 shifted_vis : array of real, size = N, unit = Jy

66 Phase shifted visibilites

68 """

69 dRA *= 2. * np.pi / rad_to_arcsec

70 dDec *= 2. * np.pi / rad_to_arcsec

72 phi = u * dRA + v * dDec

74 if inverse:

75 shifted_vis = V / (np.cos(phi) + 1j * np.sin(phi))

76 else:

77 shifted_vis = V * (np.cos(phi) + 1j * np.sin(phi))

79 return shifted_vis

82def deproject(u, v, inc, PA, inverse=False):

83 r"""

84 Deproject the image in visibily space

86 Parameters

87 ----------

88 u : array of real, size = N, unit = :math:`\lambda`

89 u-points of the visibilities

90 v : array of real, size = N, unit = :math:`\lambda`

91 v-points of the visibilities

92 inc : float, unit = deg

93 Inclination

94 PA : float, unit = deg

95 Position angle, defined east of north.

96 inverse : bool, default=False

97 If True, the uv-points are reprojected rather than deprojected

99 Returns

100 -------

101 up : array, size = N, unit = :math:`\lambda`

102 Deprojected u-points

103 vp : array, size = N, unit = :math:`\lambda`

104 Deprojected v-points

105 wp : array of real, size = N, unit = :math:`\lambda`

106 Fourier w-points of the deprojected visibilities. Only returned if

107 deprojecting.

108

109 """

110

111 inc *= deg_to_rad

112 PA *= deg_to_rad

113

114 cos_t = np.cos(PA)

115 sin_t = np.sin(PA)

116

117 if inverse:

118 sin_t *= -1

119 u = u / np.cos(inc)

120

121 up = u * cos_t - v * sin_t

122 vp = u * sin_t + v * cos_t

123

124 if inverse:

125 return up, vp

126 else:

127 # Deproject

128 wp = up * np.sin(inc)

129 up = up * np.cos(inc)

130

131 return up, vp, wp

132

133def rescale_total_flux(V, weights, inc):

134 r"""

135 Scale the visibility amplitudes (and weights) according to the source

136 inclination.

137

138 Parameters

139 ----------

140 V : array of real, size = N, unit = Jy

141 Real component of the complex, deprojected visibilities

142 weights : array of real, size = N, unit = Jy

143 Weights on the visibilities

144 inc : float, unit = deg

145 Inclination of the disc

146

147 Returns

148 -------

149 V_scaled : array of real, size = N, unit = Jy

150 Rescaled real component of the complex visibilities

151 weights_scaled : array of real, size = N, unit = Jy

152 Rescaled weights on the visibilities

153

154 Notes

155 -----

156 This scaling accounts for the difference between the inclined (observed)

157 brightness and the assumed face-on brightness, assuming the

158 emission is optically thick. The source's integrated (2D) flux is assumed

159 to be

160 :math:`F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}`.

161 No rescaling would be appropriate in the optically thin limit.

162 """

163

164 # Ensure we're only altering the real component of the visibilities

165 V = V.real

166

167 V_scaled = V / np.cos(inc * deg_to_rad)

168 weights_scaled = weights * np.cos(inc * deg_to_rad) ** 2

169

170 return V_scaled, weights_scaled

171

172

173class SourceGeometry(object):

174 """

175 Base class for geometry corrections.

176

177 Centre and deproject the source to ensure axisymmetry

178

179 Parameters

180 ----------

181 inc : float, unit = deg

182 Inclination of the disc

183 PA : float, unit = deg

184 Position angle of the disc

185 dRA : float, unit = arcsec

186 Phase centre offset in right ascension.

187 dDec : float, units = arcsec

188 Phase centre offset in declination.

189

190 Notes

191 -----

192 The phase centre offsets, dRA and dDec, refer to the distance to the source

193 from the phase centre.

194 """

195

196 def __init__(self, inc=None, PA=None, dRA=None, dDec=None):

197 self._inc = inc

198 self._PA = PA

199 self._dRA = dRA

200 self._dDec = dDec

201

202 def apply_correction(self, u, v, V, use3D=False):

203 r"""

204 Correct the phase centre and deproject the visibilities

205

206 Parameters

207 ----------

208 u : array of real, size = N, unit = :math:`\lambda`

209 u-points of the visibilities

210 v : array of real, size = N, unit = :math:`\lambda`

211 v-points of the visibilities

212 V : array of real, size = N, units = Jy

213 Complex visibilites

214 use3D : bool, default=False

215 If True, also return the 3rd compoent of the

216 de-projected visibilities, wp.

217

218 Returns

219 -------

220 up : array of real, size = N, unit = :math:`\lambda`

221 Corrected u-points of the visibilities

222 vp : array of real, size = N, unit = :math:`\lambda`

223 Corrected v-points of the visibilities

224 wp : array of real, size = N, unit = :math:`\lambda`

225 [Optional] Corrected w-points of the visibilities

226 Vp : array of real, size = N, unit = Jy

227 Corrected complex visibilites

228

229 """

230 Vp = apply_phase_shift(u, v, V, self._dRA, self._dDec, inverse=True)

231 up, vp, wp = deproject(u, v, self._inc, self._PA)

232

233 if use3D:

234 return up, vp, wp, Vp

235 else:

236 return up, vp, Vp

237

238

239 def undo_correction(self, u, v, V):

240 r"""

241 Undo the phase centre correction and deprojection

242

243 Parameters

244 ----------

245 u : array of real, size = N, unit = :math:`\lambda`

246 u-points of the visibilities

247 v : array of real, size = N, unit = :math:`\lambda`

248 v-points of the visibilities

249 V : array of real, size = N, unit = Jy

250 Complex visibilites

251

252 Returns

253 -------

254 up : array of real, size = N, unit = :math:`\lambda`

255 Corrected u-points of the visibilities

256 vp : array of real, size = N, unit = :math:`\lambda`

257 Corrected v-points of the visibilities

258 Vp : array of real, size = N, unit = Jy

259 Corrected complex visibilites

260

261 """

262 up, vp = self.reproject(u, v)

263 Vp = apply_phase_shift(up, vp, V, self._dRA, self._dDec, inverse=False)

264

265 return up, vp, Vp

266

267 def deproject(self, u, v, use3D=False):

268 r"""Convert uv-points from sky-plane to deprojected space (u,v)

269

270 Parameters

271 ----------

272 u : array of real, size = N, unit = :math:`\lambda`

273 u-points of the visibilities

274 v : array of real, size = N, unit = :math:`\lambda`

275 v-points of the visibilities

276 use3D : bool, default=False

277 If True, also return the 3rd compoent of the

278 de-projected visibilities, wp.

279

280 Returns

281 -------

282 up : array of real, size = N, unit = :math:`\lambda`

283 Corrected u-points of the visibilities

284 vp : array of real, size = N, unit = :math:`\lambda`

285 Corrected v-points of the visibilities

286 wp : array of real, size = N, unit = :math:`\lambda`

287 [Optional] Corrected w-points of the visibilities

288

289 """

290 if use3D:

291 return deproject(u, v, self._inc, self._PA)

292 else:

293 return deproject(u, v, self._inc, self._PA)[:2]

294

295 def reproject(self, u, v):

296 r"""Convert uv-points from deprojected space to sky-plane

297

298 Parameters

299 ----------

300 u : array of real, size = N, unit = :math:`\lambda`

301 u-points of the visibilities

302 v : array of real, size = N, unit = :math:`\lambda`

303 v-points of the visibilities

304

305 Returns

306 -------

307 up : array of real, size = N, unit = :math:`\lambda`

308 Corrected u-points of the visibilities

309 vp : array of real, size = N, unit = :math:`\lambda`

310 Corrected v-points of the visibilities

311

312 """

313 return deproject(u, v, self._inc, self._PA, inverse=True)

314

315

316 def rescale_total_flux(self, V, weights):

317 return rescale_total_flux(V, weights, self._inc)

318

319 def fit(self, u, v, V, weights):

320 r"""

321 Determine geometry using the provided uv-data

322

323 Parameters

324 ----------

325 u : array of real, size = N, unit = :math:`\lambda`

326 u-points of the visibilities

327 v : array of real, size = N, unit = :math:`\lambda`

328 v-points of the visibilities

329 V : array of complex, size = N, unit = Jy

330 Complex visibilites

331 weights : array of real, size = N, unit = Jy

332 Weights on the visibilities

333 """

334

335 return

336

337 def clone(self):

338 """Save the geometry parameters in a seperate geometry object"""

339 return FixedGeometry(self.inc, self.PA, self.dRA, self.dDec)

340

341 @property

342 def dRA(self):

343 """Phase centre offset in right ascension, unit = arcsec"""

344 return self._dRA

345

346 @property

347 def dDec(self):

348 """Phase centre offset in declination, unit = arcsec"""

349 return self._dDec

350

351 @property

352 def PA(self):

353 """Position angle of the disc, unit = rad"""

354 return self._PA

355

356 @property

357 def inc(self):

358 """Inclination of the disc, unit = rad"""

359 return self._inc

360

361 @property

362 def rescale_factor(self):

363 """Factor used to rescale the visibility amplitudes, unit = 1 / rad"""

364 return 1.0 / np.cos(self._inc * deg_to_rad)

365

366 def __repr__(self):

367 return "SourceGeometry(inc={}, PA={}, dRA={}, dDec={})".format(

368 self.inc, self.PA, self.dRA, self.dDec

369 )

370

371

372class FixedGeometry(SourceGeometry):

373 """

374 Disc Geometry class using pre-determined parameters.

375

376 Centre and deproject the source to ensure axisymmetry

377

378 Parameters

379 ----------

380 inc : float, unit = deg

381 Disc inclination

382 PA : float, unit = deg

383 Disc positition angle.

384 dRA : float, default = 0, unit = arcsec

385 Phase centre offset in right ascension

386 dDec : float, default = 0, unit = arcsec

387 Phase centre offset in declination

388

389 Notes

390 -----

391 The phase centre offsets, dRA and dDec, refer to the distance to the source

392 from the phase centre.

393 """

394

395 def __init__(self, inc, PA, dRA=0.0, dDec=0.0):

396 super(FixedGeometry, self).__init__(inc, PA, dRA, dDec)

397

398 def __repr__(self):

399 return "FixedGeometry(inc={}, PA={}, dRA={}, dDEC={})".format(

400 self.inc, self.PA, self.dRA, self.dDec

401 )

402

403

404class FitGeometryGaussian(SourceGeometry):

405 """

406 Determine the disc geometry by fitting a Gaussian in Fourier space.

407

408 Centre and deproject the source to ensure axisymmetry

409

410 Parameters

411 ----------

412 inc_pa : tuple = (inclination, position angle) or None (default), unit = deg

413 Determine whether to fit for the source's inclination and position

414 angle. If inc_pa = None, the inclination and PA are fit for. Else

415 inc_pa should be provided as a tuple

416 phase_centre : tuple = (dRA, dDec) or None (default), unit = arcsec

417 Determine whether to fit for the source's phase centre. If

418 phase_centre = None, the phase centre is fit for. Else the phase

419 centre should be provided as a tuple

420 guess : list of len(4), default = None

421 Initial guess for the source's inclination [deg], position angle [deg],

422 right ascension offset [arcsec], declination offset [arcsec].

423

424 Notes

425 -----

426 The phase centre offsets, dRA and dDec, refer to the distance to the source

427 from the phase centre.

428 """

429

430 def __init__(self, inc_pa=None, phase_centre=None, guess=None):

431 super(FitGeometryGaussian, self).__init__()

432

433 self._inc_pa = inc_pa

434 self._phase_centre = phase_centre

435 self._guess = guess

436

437 if guess is None:

438 guess = [10.0, 10.0, 0.0, 0.0, 1.0, 1.0]

439 else:

440 guess.extend([1.0, 1.0])

441

442 if self._inc_pa is not None:

443 guess[0], guess[1] = self._inc_pa

444 if self._phase_centre is not None:

445 guess[2], guess[3] = self._phase_centre

446

447 self._guess = guess

448

449 def fit(self, u, v, V, weights):

450 r"""

451 Determine geometry using the provided uv-data

452

453 Parameters

454 ----------

455 u : array of real, size = N, unit = :math:`\lambda`

456 u-points of the visibilities

457 v : array of real, size = N, unit = :math:`\lambda`

458 v-points of the visibilities

459 V : array of complex, size = N, unit = Jy

460 Complex visibilites

461 weights : array of real, size=N, unit = Jy^-2

462 Weights on the visibilities

463 """

464

465 if self._inc_pa and self._phase_centre:

466 logging.info(' You requested a Gaussian fit to determine the geometry,'

467 ' but you provided values for inclination, PA, and the phase offset.'

468 ' --> Using your provided values (not fitting for the geometry)')

469 self._inc, self._PA = self._inc_pa

470 self._dRA, self._dDec = self._phase_centre

471

472 else:

473 if self._inc_pa:

474 logging.info(' Fitting Gaussian to determine geometry'

475 ' (not fitting for inc or PA)')

476

477 elif self._phase_centre:

478 logging.info(' Fitting Gaussian to determine geometry'

479 ' (not fitting for phase center)')

480

481 else:

482 logging.info(' Fitting Gaussian to determine geometry')

483

484 inc, PA, dRA, dDec = _fit_geometry_gaussian(

485 u, v, V, weights, guess=self._guess,

486 inc_pa=self._inc_pa,

487 phase_centre=self._phase_centre)

488

489 if not self._inc_pa:

490 inc, PA = _fix_inc_and_PA_ranges(inc, PA)

491

492 self._inc = inc

493 self._PA = PA

494 self._dRA = dRA

495 self._dDec = dDec

496

497

498def _fit_geometry_gaussian(u, v, V, weights, guess, inc_pa=None,

499 phase_centre=None):

500 r"""

501 Estimate the source geometry by fitting a Gaussian in uv-space

502

503 Parameters

504 ----------

505 u : array of real, size = N, unit = :math:`\lambda`

506 u-points of the visibilities

507 v : array of real, size = N, unit = :math:`\lambda`

508 v-points of the visibilities

509 V : array of complex, size = N, unit = Jy

510 Complex visibilites

511 weights : array of real, size = N, unit = Jy^-2

512 Weights on the visibilities

513 guess : list of len(6)

514 Initial guess for the source's inclination [deg], position angle [deg],

515 right ascension offset [arcsec], declination offset [arcsec],

516 the Gaussian's normalization, and its scaling. The latter 2 are forced

517 as 1.0

518 inc_pa: [inclination, position angle], optional, unit = deg

519 The inclination and position angle.

520 If not provided, these will be fit for

521 phase_centre: [dRA, dDec], optional, unit = arcsec

522 The phase centre offsets dRA and dDec.

523 If not provided, these will be fit for

524

525 Returns

526 -------

527 geometry : SourceGeometry object

528 Fitted geometry

529 """

530 fac = 2*np.pi / rad_to_arcsec

531 w = np.sqrt(weights)

532

533 if inc_pa is not None:

534 inc, PA = inc_pa

535 # Convert inc and PA from [deg] --> [rad]

536 inc *= deg_to_rad

537 PA *= deg_to_rad

538

539 if phase_centre is not None:

540 dRA, dDec = phase_centre

541 phi = dRA*fac * u + dDec*fac * v

542 V = V * (np.cos(phi) - 1j*np.sin(phi))

543

544 # Convert guess inc and PA from [deg] --> [rad]

545 guess[0] *= deg_to_rad

546 guess[1] *= deg_to_rad

547

548 def wrap(fun):

549 return np.concatenate([fun.real, fun.imag])

550

551 def _gauss_fun(params):

552 inc, PA, dRA, dDec, norm, scal = params

553

554 if phase_centre is None:

555 phi = dRA*fac * u + dDec*fac * v

556 Vp = V * (np.cos(phi) - 1j*np.sin(phi))

557 else:

558 Vp = V

559

560 c_t = np.cos(PA)

561 s_t = np.sin(PA)

562 c_i = np.cos(inc)

563 up = (u*c_t - v*s_t) * c_i / (scal*rad_to_arcsec)

564 vp = (u*s_t + v*c_t) / (scal*rad_to_arcsec)

565

566 fun = w*(norm * np.exp(-0.5*(up*up + vp*vp)) - Vp)

567

568 return wrap(fun)

569

570 def _gauss_jac(params):

571 inc, PA, dRA, dDec, norm, scal = params

572

573 jac = np.zeros([6, 2*len(w)])

574

575 if phase_centre is None:

576 phi = dRA*fac * u + dDec*fac * v

577 dVp = - w*V * (-np.sin(phi) - 1j*np.cos(phi)) * fac

578

579 jac[2] = wrap(dVp*u)

580 jac[3] = wrap(dVp*v)

581

582 c_t = np.cos(PA)

583 s_t = np.sin(PA)

584 c_i = np.cos(inc)

585 s_i = np.sin(inc)

586 up = (u*c_t - v*s_t)

587 vp = (u*s_t + v*c_t)

588

589 uv = (up*up*c_i*c_i + vp*vp)

590

591 G = w*np.exp(-0.5 * uv / (scal*rad_to_arcsec)**2)

592

593 norm = norm / (scal*rad_to_arcsec)**2

594

595 if inc_pa is None:

596 jac[0] = wrap(norm*G*up*up*c_i*s_i)

597 jac[1] = wrap(norm*G*up*vp*(c_i*c_i - 1)/2)

598

599 jac[4] = wrap(G)

600 jac[5] = wrap(norm*G*uv/scal)

601

602 return jac.T

603

604

605 res = least_squares(_gauss_fun, guess,

606 jac=_gauss_jac, method='lm')

607

608 inc, PA, dRA, dDec, _, _ = res.x

609

610 if inc_pa is not None:

611 inc, PA = inc_pa

612 else:

613 # convert back to [deg]

614 inc /= deg_to_rad

615 PA /= deg_to_rad

616

617 if phase_centre is not None:

618 dRA, dDec = phase_centre

619

620 geometry = inc, PA, dRA, dDec

621

622 return geometry

623

624

625class FitGeometryFourierBessel(SourceGeometry):

626 """

627 Determine the disc geometry by fitting a non-parametric brightness

628 profile in visibility space.

629

630 The best fit is obtained by finding the geometry that minimizes

631 the weighted chi^2 of the visibility fit.

632

633 The brightness profile is modelled using the FourierBesselFitter,

634 which is equivalent to a FrankFitter fit without the Gaussian

635 Process prior. For this reason, a small number of bins is

636 recommended for fit stability.

637

638 Parameters

639 ----------

640 Rmax : float, unit = arcsec

641 Radius of support for the functions to transform, i.e.,

642 f(r) = 0 for R >= Rmax

643 N : int

644 Number of collocation points

645 inc_pa : tuple = (inclination, position angle) or None (default), unit = deg

646 Determine whether to fit for the source's inclination and position

647 angle. If inc_pa = None, the inclination and PA are fit for. Else

648 inc_pa should be provided as a tuple

649 phase_centre : tuple = (dRA, dDec) or None (default), unit = arcsec

650 Determine whether to fit for the source's phase centre. If

651 phase_centre = None, the phase centre is fit for. Else the phase

652 centre should be provided as a tuple

653 guess : list of len(4), default = None

654 Initial guess for the source's inclination [deg], position angle [deg],

655 right ascension offset [arcsec], and declination offset [arcsec]

656 verbose : bool, default=False

657 Determines whether to print the iteration progress.

658 """

659

660 def __init__(self, Rmax, N, inc_pa=None, phase_centre=None, guess=None,

661 verbose=False):

662 self._N = N

663 self._R = Rmax

664 self._inc_pa = inc_pa

665 self._phase_centre = phase_centre

666

667 if guess is None:

668 guess = [10., 10., 0., 0.]

669 if self._inc_pa is not None:

670 guess[0], guess[1] = self._inc_pa

671 if self._phase_centre is not None:

672 guess[2], guess[3] = self._phase_centre

673

674 self._guess = guess

675

676 self._verbose = verbose

677

678 def _residual(self, params, uvdata=None):

679 inc, pa, dRA, dDec = params

680 if self._inc_pa is not None:

681 inc, pa = self._inc_pa

682 if self._phase_centre is not None:

683 dRA, dDec = self._phase_centre

684

685 geom = FixedGeometry(inc, pa, dRA, dDec)

686

687

688 FBF = FourierBesselFitter(self._R, self._N, geom, verbose=False)

689

690 u, v, vis, w_half = uvdata

691

692 sol = FBF.fit(u,v,vis, w_half*w_half)

693

694 error = w_half*(sol.predict(u,v) - vis)

695

696 if self._verbose:

697 Chi2 = 0.5 * np.sum(error.real**2 + error.imag**2) / len(w_half)

698 print('\n FitGeometryFourierBessel: Iteration {}, chi^2={:.8f}, inc={:.3f} PA={:.3f} dRA={:.5f} dDec={:.5f}'

699 ''.format(self._counter, Chi2, inc, pa, dRA, dDec),

700 end='', flush=True)

701 self._counter += 1

702

703 return np.concatenate([error.real, error.imag])

704

705 def fit(self, u, v, vis, w):

706 r"""

707 Determine geometry using the provided uv-data

708

709 Parameters

710 ----------

711 u : array of real, size = N, unit = :math:`\lambda`

712 u-points of the visibilities

713 v : array of real, size = N, unit = :math:`\lambda`

714 v-points of the visibilities

715 vis : array of complex, size = N, unit = Jy

716 Complex visibilites

717 w : array of real, size = N, unit = Jy

718 Weights on the visibilities

719 """

720 if self._inc_pa and self._phase_centre:

721 logging.info(' You requested a nonparametric fit to determine the geometry,'

722 ' but you provided values for inclination, PA, and the phase offset.'

723 ' --> Using your provided values (not fitting for the geometry)')

724 self._inc, self._PA = self._inc_pa

725 self._dRA, self._dDec = self._phase_centre

726

727 else:

728 if self._inc_pa:

729 logging.info(' Fitting nonparametric form to determine geometry'

730 ' (your supplied inclination and position angle will'

731 ' be applied at the end of the geometry fitting'

732 ' routine)')

733 elif self._phase_centre:

734 logging.info(' Fitting nonparametric form to determine geometry'

735 ' (your supplied phase center will be applied at the'

736 ' end of the geometry fitting routine)')

737

738 else:

739 logging.info(' Fitting nonparametric form to determine geometry')

740

741 uvdata= [u, v, vis, w**0.5]

742

743 self._counter = 0

744

745 result = least_squares(self._residual, self._guess, kwargs={'uvdata':uvdata},

746 method='lm')

747

748 if not result.success:

749 raise RuntimeError("FitGeometryFourierBessel failed to converge")

750

751 inc, pa, dRA, dDec = result.x

752 if self._inc_pa:

753 inc, pa = self._inc_pa

754 if self._phase_centre:

755 dRA, dDec = self._phase_centre

756

757 if not self._inc_pa:

758 inc, pa = _fix_inc_and_PA_ranges(inc, pa)

759

760 self._inc = inc

761 self._PA = pa

762 self._dRA = dRA

763 self._dDec = dDec