Coverage for starry/_core/core.py : 94%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
2from .. import config
3from .._constants import *
4from .. import _c_ops
5from .ops import (
6 sTOp,
7 rTReflectedOp,
8 dotROp,
9 tensordotRzOp,
10 FOp,
11 tensordotDOp,
12 spotYlmOp,
13 pTOp,
14 minimizeOp,
15 LDPhysicalOp,
16 LimbDarkOp,
17 GetClOp,
18 RaiseValueErrorOp,
19 RaiseValueErrorIfOp,
20)
21from .utils import logger, autocompile
22from .math import math
23import theano
24import theano.tensor as tt
25import theano.sparse as ts
26from theano.ifelse import ifelse
27import numpy as np
28from astropy import units
30try: # pragma: no cover
31 # starry requires exoplanet >= v0.2.0
32 from packaging import version
33 import exoplanet
35 if version.parse(exoplanet.__version__) < version.parse("0.2.0"):
36 exoplanet = None
37except ModuleNotFoundError: # pragma: no cover
38 exoplanet = None
41__all__ = ["OpsYlm", "OpsLD", "OpsReflected", "OpsRV", "OpsSystem"]
44class OpsYlm(object):
45 """Class housing Theano operations for spherical harmonics maps."""
47 def __init__(
48 self, ydeg, udeg, fdeg, drorder, nw, reflected=False, **kwargs
49 ):
50 # Ingest kwargs
51 self.ydeg = ydeg
52 self.udeg = udeg
53 self.fdeg = fdeg
54 self.deg = ydeg + udeg + fdeg
55 self.filter = (fdeg > 0) or (udeg > 0)
56 self.drorder = drorder
57 self.diffrot = drorder > 0
58 self.nw = nw
59 self._reflected = reflected
61 # Instantiate the C++ Ops
62 config.rootHandler.terminator = ""
63 logger.info("Pre-computing some matrices... ")
64 self._c_ops = _c_ops.Ops(ydeg, udeg, fdeg, drorder)
65 config.rootHandler.terminator = "\n"
66 logger.info("Done.")
68 # Solution vectors
69 self._sT = sTOp(self._c_ops.sT, self._c_ops.N)
70 self._rT = tt.shape_padleft(tt.as_tensor_variable(self._c_ops.rT))
71 self._rTA1 = tt.shape_padleft(tt.as_tensor_variable(self._c_ops.rTA1))
73 # Change of basis matrices
74 self._A = ts.as_sparse_variable(self._c_ops.A)
75 self._A1 = ts.as_sparse_variable(self._c_ops.A1)
76 self._A1Inv = ts.as_sparse_variable(self._c_ops.A1Inv)
78 # Rotation operations
79 self._tensordotRz = tensordotRzOp(self._c_ops.tensordotRz)
80 self._dotR = dotROp(self._c_ops.dotR)
82 # Filter
83 # TODO: Make the filter operator sparse
84 self._F = FOp(self._c_ops.F, self._c_ops.N, self._c_ops.Ny)
86 # Differential rotation
87 self._tensordotD = tensordotDOp(self._c_ops.tensordotD)
89 # Misc
90 self._spotYlm = spotYlmOp(self._c_ops.spotYlm, self.ydeg, self.nw)
91 self._pT = pTOp(self._c_ops.pT, self.deg)
92 if self.nw is None:
93 if self._reflected:
94 self._minimize = minimizeOp(
95 self.unweighted_intensity,
96 self.P,
97 self.ydeg,
98 self.udeg,
99 self.fdeg,
100 )
101 else:
102 self._minimize = minimizeOp(
103 self.intensity, self.P, self.ydeg, self.udeg, self.fdeg
104 )
105 else:
106 # TODO: Implement minimization for spectral maps?
107 self._minimize = None
108 self._LimbDarkIsPhysical = LDPhysicalOp(_c_ops.nroots)
110 @property
111 def rT(self):
112 return self._rT
114 @property
115 def rTA1(self):
116 return self._rTA1
118 @property
119 def A(self):
120 return self._A
122 @property
123 def A1(self):
124 return self._A1
126 @property
127 def A1Inv(self):
128 return self._A1Inv
130 @autocompile
131 def sT(self, b, r):
132 return self._sT(b, r)
134 @autocompile
135 def tensordotRz(self, matrix, theta):
136 return self._tensordotRz(matrix, theta)
138 @autocompile
139 def dotR(self, matrix, ux, uy, uz, theta):
140 return self._dotR(matrix, ux, uy, uz, theta)
142 @autocompile
143 def F(self, u, f):
144 return self._F(u, f)
146 @autocompile
147 def tensordotD(self, matrix, wta):
148 return self._tensordotD(matrix, wta)
150 @autocompile
151 def spotYlm(self, amp, sigma, lat, lon):
152 return self._spotYlm(amp, sigma, lat, lon)
154 @autocompile
155 def pT(self, x, y, z):
156 return self._pT(x, y, z)
158 @autocompile
159 def limbdark_is_physical(self, u):
160 """Return True if the limb darkening profile is physical."""
161 return self._LimbDarkIsPhysical(u)
163 @autocompile
164 def get_minimum(self, y):
165 """Compute the location and value of the intensity minimum."""
166 return self._minimize(y)
168 @autocompile
169 def X(self, theta, xo, yo, zo, ro, inc, obl, u, f, alpha):
170 """Compute the light curve design matrix."""
171 # Determine shapes
172 rows = theta.shape[0]
173 cols = self.rTA1.shape[1]
174 X = tt.zeros((rows, cols))
176 # Compute the occultation mask
177 b = tt.sqrt(xo ** 2 + yo ** 2)
178 b_rot = tt.ge(b, 1.0 + ro) | tt.le(zo, 0.0) | tt.eq(ro, 0.0)
179 b_occ = tt.invert(b_rot)
180 i_rot = tt.arange(b.size)[b_rot]
181 i_occ = tt.arange(b.size)[b_occ]
183 # Compute filter operator
184 if self.filter:
185 F = self.F(u, f)
187 # Rotation operator
188 if self.filter:
189 rTA1 = ts.dot(tt.dot(self.rT, F), self.A1)
190 else:
191 rTA1 = self.rTA1
192 X = tt.set_subtensor(
193 X[i_rot], self.right_project(rTA1, inc, obl, theta[i_rot], alpha)
194 )
196 # Occultation + rotation operator
197 sT = self.sT(b[i_occ], ro)
198 sTA = ts.dot(sT, self.A)
199 theta_z = tt.arctan2(xo[i_occ], yo[i_occ])
200 sTAR = self.tensordotRz(sTA, theta_z)
201 if self.filter:
202 A1InvFA1 = ts.dot(ts.dot(self.A1Inv, F), self.A1)
203 sTAR = tt.dot(sTAR, A1InvFA1)
204 X = tt.set_subtensor(
205 X[i_occ], self.right_project(sTAR, inc, obl, theta[i_occ], alpha)
206 )
208 return X
210 @autocompile
211 def flux(self, theta, xo, yo, zo, ro, inc, obl, y, u, f, alpha):
212 """Compute the light curve."""
213 return tt.dot(self.X(theta, xo, yo, zo, ro, inc, obl, u, f, alpha), y)
215 @autocompile
216 def P(self, lat, lon):
217 """Compute the pixelization matrix, no filters or illumination."""
218 # Get the Cartesian points
219 xpt, ypt, zpt = self.latlon_to_xyz(lat, lon)
221 # Compute the polynomial basis at the point
222 pT = self.pT(xpt, ypt, zpt)[:, : (self.ydeg + 1) ** 2]
224 # Transform to the Ylm basis
225 pTA1 = ts.dot(pT, self.A1)
227 # We're done
228 return pTA1
230 @autocompile
231 def intensity(self, lat, lon, y, u, f, wta, ld):
232 """Compute the intensity at a point or a set of points."""
233 # Get the Cartesian points
234 xpt, ypt, zpt = self.latlon_to_xyz(lat, lon)
236 # Compute the polynomial basis at the point
237 pT = self.pT(xpt, ypt, zpt)
239 # Apply the differential rotation operator
240 if self.diffrot:
241 if self.nw is None:
242 y = tt.reshape(
243 self.tensordotD(tt.reshape(y, (1, -1)), [wta]), (-1,)
244 )
245 else:
246 y = tt.transpose(
247 self.tensordotD(tt.transpose(y), tt.ones(self.nw) * wta)
248 )
250 # Transform the map to the polynomial basis
251 A1y = ts.dot(self.A1, y)
253 # Apply the filter
254 if self.filter:
255 u0 = tt.zeros_like(u)
256 u0 = tt.set_subtensor(u0[0], -1.0)
257 A1y = ifelse(
258 ld, tt.dot(self.F(u, f), A1y), tt.dot(self.F(u0, f), A1y)
259 )
261 # Dot the polynomial into the basis
262 return tt.dot(pT, A1y)
264 @autocompile
265 def render(self, res, projection, theta, inc, obl, y, u, f, alpha):
266 """Render the map on a Cartesian grid."""
267 # Compute the Cartesian grid
268 xyz = ifelse(
269 tt.eq(projection, STARRY_RECTANGULAR_PROJECTION),
270 self.compute_rect_grid(res),
271 self.compute_ortho_grid(res),
272 )
274 # Compute the polynomial basis
275 pT = self.pT(xyz[0], xyz[1], xyz[2])
277 # If orthographic, rotate the map to the correct frame
278 if self.nw is None:
279 Ry = ifelse(
280 tt.eq(projection, STARRY_ORTHOGRAPHIC_PROJECTION),
281 self.left_project(
282 tt.transpose(tt.tile(y, [theta.shape[0], 1])),
283 inc,
284 obl,
285 theta,
286 alpha,
287 ),
288 tt.transpose(
289 self.tensordotD(
290 tt.tile(y, [theta.shape[0], 1]), theta * alpha
291 )
292 )
293 if self.diffrot
294 else tt.transpose(tt.tile(y, [theta.shape[0], 1])),
295 )
296 else:
297 Ry = ifelse(
298 tt.eq(projection, STARRY_ORTHOGRAPHIC_PROJECTION),
299 self.left_project(
300 y, inc, obl, tt.tile(theta[0], self.nw), alpha
301 ),
302 y,
303 )
305 # Change basis to polynomials
306 A1Ry = ts.dot(self.A1, Ry)
308 # Apply the filter *only if orthographic*
309 if self.filter:
310 f0 = tt.zeros_like(f)
311 f0 = tt.set_subtensor(f0[0], np.pi)
312 u0 = tt.zeros_like(u)
313 u0 = tt.set_subtensor(u0[0], -1.0)
314 A1Ry = ifelse(
315 tt.eq(projection, STARRY_ORTHOGRAPHIC_PROJECTION),
316 tt.dot(self.F(u, f), A1Ry),
317 tt.dot(self.F(u0, f0), A1Ry),
318 )
320 # Dot the polynomial into the basis
321 res = tt.reshape(tt.dot(pT, A1Ry), [res, res, -1])
323 # We need the shape to be (nframes, npix, npix)
324 return res.dimshuffle(2, 0, 1)
326 @autocompile
327 def expand_spot(self, amp, sigma, lat, lon):
328 """Return the spherical harmonic expansion of a Gaussian spot."""
329 return self.spotYlm(amp, sigma, lat, lon)
331 @autocompile
332 def compute_ortho_grid(self, res):
333 """Compute the polynomial basis on the plane of the sky."""
334 dx = 2.0 / res
335 y, x = tt.mgrid[-1:1:dx, -1:1:dx]
336 x = tt.reshape(x, [1, -1])
337 y = tt.reshape(y, [1, -1])
338 z = tt.sqrt(1 - x ** 2 - y ** 2)
339 return tt.concatenate((x, y, z))
341 @autocompile
342 def compute_rect_grid(self, res):
343 """Compute the polynomial basis on a rectangular lat/lon grid."""
344 dx = np.pi / res
345 lat, lon = tt.mgrid[
346 -np.pi / 2 : np.pi / 2 : dx, -3 * np.pi / 2 : np.pi / 2 : 2 * dx
347 ]
348 x = tt.reshape(tt.cos(lat) * tt.cos(lon), [1, -1])
349 y = tt.reshape(tt.cos(lat) * tt.sin(lon), [1, -1])
350 z = tt.reshape(tt.sin(lat), [1, -1])
351 R = self.RAxisAngle(tt.as_tensor_variable([1.0, 0.0, 0.0]), -np.pi / 2)
352 return tt.dot(R, tt.concatenate((x, y, z)))
354 @autocompile
355 def right_project(self, M, inc, obl, theta, alpha):
356 r"""Apply the projection operator on the right.
358 Specifically, this method returns the dot product :math:`M \cdot R`,
359 where ``M`` is an input matrix and ``R`` is the Wigner rotation matrix
360 that transforms a spherical harmonic coefficient vector in the
361 input frame to a vector in the observer's frame.
362 """
363 # Trivial case
364 if self.ydeg == 0:
365 return M
367 # Rotate to the sky frame
368 # TODO: Do this in a single compound rotation
369 M = self.dotR(
370 self.dotR(
371 self.dotR(
372 M,
373 -tt.cos(obl),
374 -tt.sin(obl),
375 math.to_tensor(0.0),
376 -(0.5 * np.pi - inc),
377 ),
378 math.to_tensor(0.0),
379 math.to_tensor(0.0),
380 math.to_tensor(1.0),
381 obl,
382 ),
383 math.to_tensor(1.0),
384 math.to_tensor(0.0),
385 math.to_tensor(0.0),
386 -0.5 * np.pi,
387 )
389 # Rotate to the correct phase
390 if theta.ndim > 0:
391 M = self.tensordotRz(M, theta)
392 else:
393 M = self.dotR(
394 M,
395 math.to_tensor(0.0),
396 math.to_tensor(0.0),
397 math.to_tensor(1.0),
398 theta,
399 )
401 # Rotate to the polar frame
402 M = self.dotR(
403 M,
404 math.to_tensor(1.0),
405 math.to_tensor(0.0),
406 math.to_tensor(0.0),
407 0.5 * np.pi,
408 )
410 # Apply the differential rotation
411 if self.diffrot:
412 if theta.ndim > 0:
413 M = self.tensordotD(M, -theta * alpha)
414 else:
415 M = RaiseValueErrorOp(
416 "Code this branch up if needed.", M.shape
417 )
419 return M
421 @autocompile
422 def left_project(self, M, inc, obl, theta, alpha):
423 r"""Apply the projection operator on the left.
425 Specifically, this method returns the dot product :math:`R \cdot M`,
426 where ``M`` is an input matrix and ``R`` is the Wigner rotation matrix
427 that transforms a spherical harmonic coefficient vector in the
428 input frame to a vector in the observer's frame.
430 """
431 # Trivial case
432 if self.ydeg == 0:
433 return M
435 # Note that here we are using the fact that R . M = (M^T . R^T)^T
436 MT = tt.transpose(M)
438 # Apply the differential rotation
439 if self.diffrot:
440 if theta.ndim > 0:
441 MT = self.tensordotD(MT, theta * alpha)
442 else:
443 MT = RaiseValueErrorOp(
444 "Code this branch up if needed.", MT.shape
445 )
447 # Rotate to the polar frame
448 MT = self.dotR(
449 MT,
450 math.to_tensor(1.0),
451 math.to_tensor(0.0),
452 math.to_tensor(0.0),
453 -0.5 * np.pi,
454 )
456 # Rotate to the correct phase
457 if theta.ndim > 0:
458 MT = self.tensordotRz(MT, -theta)
459 else:
460 MT = self.dotR(
461 MT,
462 math.to_tensor(0.0),
463 math.to_tensor(0.0),
464 math.to_tensor(1.0),
465 -theta,
466 )
468 # Rotate to the sky frame
469 # TODO: Do this in a single compound rotation
470 MT = self.dotR(
471 self.dotR(
472 self.dotR(
473 MT,
474 math.to_tensor(1.0),
475 math.to_tensor(0.0),
476 math.to_tensor(0.0),
477 0.5 * np.pi,
478 ),
479 math.to_tensor(0.0),
480 math.to_tensor(0.0),
481 math.to_tensor(1.0),
482 -obl,
483 ),
484 -tt.cos(obl),
485 -tt.sin(obl),
486 math.to_tensor(0.0),
487 (0.5 * np.pi - inc),
488 )
490 return tt.transpose(MT)
492 @autocompile
493 def set_map_vector(self, vector, inds, vals):
494 """Set the elements of the theano map coefficient tensor."""
495 res = tt.set_subtensor(vector[inds], vals * tt.ones_like(vector[inds]))
496 return res
498 @autocompile
499 def latlon_to_xyz(self, lat, lon):
500 """Convert lat-lon points to Cartesian points."""
501 if lat.ndim == 0:
502 lat = tt.shape_padleft(lat, 1)
503 if lon.ndim == 0:
504 lon = tt.shape_padleft(lon, 1)
505 R1 = self.RAxisAngle([1.0, 0.0, 0.0], -lat)
506 R2 = self.RAxisAngle([0.0, 1.0, 0.0], lon)
507 R = tt.batched_dot(R2, R1)
508 xyz = tt.transpose(tt.dot(R, [0.0, 0.0, 1.0]))
509 return xyz[0], xyz[1], xyz[2]
511 @autocompile
512 def RAxisAngle(self, axis=[0, 1, 0], theta=0):
513 """Wigner axis-angle rotation matrix."""
515 def compute(axis=[0, 1, 0], theta=0):
516 axis = tt.as_tensor_variable(axis)
517 axis /= axis.norm(2)
518 cost = tt.cos(theta)
519 sint = tt.sin(theta)
521 return tt.reshape(
522 tt.as_tensor_variable(
523 [
524 cost + axis[0] * axis[0] * (1 - cost),
525 axis[0] * axis[1] * (1 - cost) - axis[2] * sint,
526 axis[0] * axis[2] * (1 - cost) + axis[1] * sint,
527 axis[1] * axis[0] * (1 - cost) + axis[2] * sint,
528 cost + axis[1] * axis[1] * (1 - cost),
529 axis[1] * axis[2] * (1 - cost) - axis[0] * sint,
530 axis[2] * axis[0] * (1 - cost) - axis[1] * sint,
531 axis[2] * axis[1] * (1 - cost) + axis[0] * sint,
532 cost + axis[2] * axis[2] * (1 - cost),
533 ]
534 ),
535 [3, 3],
536 )
538 # If theta is a vector, this is a tensor!
539 if hasattr(theta, "ndim") and theta.ndim > 0:
540 fn = lambda theta, axis: compute(axis=axis, theta=theta)
541 R, _ = theano.scan(fn=fn, sequences=[theta], non_sequences=[axis])
542 return R
543 else:
544 return compute(axis=axis, theta=theta)
547class OpsLD(object):
548 """Class housing Theano operations for limb-darkened maps."""
550 def __init__(
551 self, ydeg, udeg, fdeg, drorder, nw, reflected=False, **kwargs
552 ):
553 # Sanity checks
554 assert ydeg == fdeg == drorder == 0
555 assert reflected is False
557 # Ingest kwargs
558 self.udeg = udeg
559 self.nw = nw
561 # Set up the ops
562 self._get_cl = GetClOp()
563 self._limbdark = LimbDarkOp()
564 self._LimbDarkIsPhysical = LDPhysicalOp(_c_ops.nroots)
566 @autocompile
567 def limbdark_is_physical(self, u):
568 """Return True if the limb darkening profile is physical."""
569 return self._LimbDarkIsPhysical(u)
571 @autocompile
572 def intensity(self, mu, u):
573 """Compute the intensity at a set of points."""
574 basis = tt.reshape(1.0 - mu, (-1, 1)) ** np.arange(self.udeg + 1)
575 return -tt.dot(basis, u)
577 @autocompile
578 def flux(self, xo, yo, zo, ro, u):
579 """Compute the light curve."""
580 # Initialize flat light curve
581 flux = tt.ones_like(xo)
583 # Compute the occultation mask
584 b = tt.sqrt(xo ** 2 + yo ** 2)
585 b_occ = tt.invert(tt.ge(b, 1.0 + ro) | tt.le(zo, 0.0) | tt.eq(ro, 0.0))
586 i_occ = tt.arange(b.size)[b_occ]
588 # Get the Agol `c` coefficients
589 c = self._get_cl(u)
590 if self.udeg == 0:
591 c_norm = c / (np.pi * c[0])
592 else:
593 c_norm = c / (np.pi * (c[0] + 2 * c[1] / 3))
595 # Compute the occultation flux
596 los = zo[i_occ]
597 r = ro * tt.ones_like(los)
598 flux = tt.set_subtensor(
599 flux[i_occ], self._limbdark(c_norm, b[i_occ], r, los)[0]
600 )
601 return flux
603 @autocompile
604 def X(self, theta, xo, yo, zo, ro, inc, obl, u, f, alpha):
605 """
606 Convenience function for integration of limb-darkened maps
607 with the ``System`` class. The design matrix for limb-darkened
608 maps is just a column vector equal to the total flux, since the
609 spherical harmonic coefficient vector is ``[1.0]``.
611 """
612 flux = self.flux(xo, yo, zo, ro, u)
613 X = tt.reshape(flux, (-1, 1))
614 return X
616 @autocompile
617 def render(self, res, projection, theta, inc, obl, y, u, f, alpha):
618 """Render the map on a Cartesian grid."""
619 nframes = tt.shape(theta)[0]
620 image = self.render_ld(res, u)
621 return tt.tile(image, (nframes, 1, 1))
623 @autocompile
624 def render_ld(self, res, u):
625 """Simplified version of `render` w/o the extra params.
627 The method `render` requires a bunch of dummy params for
628 compatibility with the `System` class. This method is a
629 convenience method for use in the `Map` class.
630 """
631 # TODO: There may be a bug in Theano related to
632 # tt.mgrid; I get different results depending on whether the
633 # function is compiled using `theano.function()` or if it
634 # is evaluated using `.eval()`. The small perturbation to `res`
635 # is a temporary fix that ensures that `y` and `x` are of the
636 # correct length in all cases I've tested.
638 # Compute the Cartesian grid
639 dx = 2.0 / (res - 0.01)
640 y, x = tt.mgrid[-1:1:dx, -1:1:dx]
642 x = tt.reshape(x, [1, -1])
643 y = tt.reshape(y, [1, -1])
644 mu = tt.sqrt(1 - x ** 2 - y ** 2)
646 # Compute the intensity
647 intensity = self.intensity(mu, u)
649 # We need the shape to be (nframes, npix, npix)
650 return tt.reshape(intensity, (1, res, res))
652 @autocompile
653 def set_map_vector(self, vector, inds, vals):
654 """Set the elements of the theano map coefficient tensor."""
655 res = tt.set_subtensor(vector[inds], vals * tt.ones_like(vector[inds]))
656 return res
659class OpsRV(OpsYlm):
660 """Class housing Theano operations for radial velocity maps."""
662 @autocompile
663 def compute_rv_filter(self, inc, obl, veq, alpha):
664 """Compute the radial velocity field Ylm multiplicative filter."""
665 # Define some angular quantities
666 cosi = tt.cos(inc)
667 sini = tt.sin(inc)
668 cosl = tt.cos(obl)
669 sinl = tt.sin(obl)
670 A = sini * cosl
671 B = sini * sinl
672 C = cosi
674 # Compute the Ylm expansion of the RV field
675 return (
676 tt.reshape(
677 [
678 0,
679 veq
680 * np.sqrt(3)
681 * B
682 * (-(A ** 2) * alpha - B ** 2 * alpha - C ** 2 * alpha + 5)
683 / 15,
684 0,
685 veq
686 * np.sqrt(3)
687 * A
688 * (-(A ** 2) * alpha - B ** 2 * alpha - C ** 2 * alpha + 5)
689 / 15,
690 0,
691 0,
692 0,
693 0,
694 0,
695 veq * alpha * np.sqrt(70) * B * (3 * A ** 2 - B ** 2) / 70,
696 veq
697 * alpha
698 * 2
699 * np.sqrt(105)
700 * C
701 * (-(A ** 2) + B ** 2)
702 / 105,
703 veq
704 * alpha
705 * np.sqrt(42)
706 * B
707 * (A ** 2 + B ** 2 - 4 * C ** 2)
708 / 210,
709 0,
710 veq
711 * alpha
712 * np.sqrt(42)
713 * A
714 * (A ** 2 + B ** 2 - 4 * C ** 2)
715 / 210,
716 veq * alpha * 4 * np.sqrt(105) * A * B * C / 105,
717 veq * alpha * np.sqrt(70) * A * (A ** 2 - 3 * B ** 2) / 70,
718 ],
719 [-1],
720 )
721 * np.pi
722 )
724 @autocompile
725 def rv(self, theta, xo, yo, zo, ro, inc, obl, y, u, veq, alpha):
726 """Compute the observed radial velocity anomaly."""
727 # Compute the velocity-weighted intensity
728 f = self.compute_rv_filter(inc, obl, veq, alpha)
729 Iv = self.flux(theta, xo, yo, zo, ro, inc, obl, y, u, f, alpha)
731 # Compute the inverse of the intensity
732 f0 = tt.zeros_like(f)
733 f0 = tt.set_subtensor(f0[0], np.pi)
734 I = self.flux(theta, xo, yo, zo, ro, inc, obl, y, u, f0, alpha)
735 invI = tt.ones((1,)) / I
736 invI = tt.where(tt.isinf(invI), 0.0, invI)
738 # The RV signal is just the product
739 return Iv * invI
742class OpsReflected(OpsYlm):
743 """Class housing Theano operations for reflected light maps."""
745 def __init__(self, *args, **kwargs):
746 super(OpsReflected, self).__init__(*args, reflected=True, **kwargs)
747 self._rT = rTReflectedOp(self._c_ops.rTReflected, self._c_ops.N)
748 self._A1Big = ts.as_sparse_variable(self._c_ops.A1Big)
750 @property
751 def A1Big(self):
752 return self._A1Big
754 @autocompile
755 def rT(self, b):
756 return self._rT(b)
758 @autocompile
759 def intensity(self, lat, lon, y, u, f, xs, ys, zs, wta, ld):
760 """Compute the intensity at a series of lat-lon points on the surface."""
761 # Get the Cartesian points
762 xpt, ypt, zpt = self.latlon_to_xyz(lat, lon)
764 # Compute the polynomial basis at the point
765 pT = self.pT(xpt, ypt, zpt)
767 # Apply the differential rotation operator
768 if self.diffrot:
769 if self.nw is None:
770 y = tt.reshape(
771 self.tensordotD(tt.reshape(y, (1, -1)), [wta]), (-1,)
772 )
773 else:
774 y = tt.transpose(
775 self.tensordotD(tt.transpose(y), tt.ones(self.nw) * wta)
776 )
778 # Transform the map to the polynomial basis
779 A1y = ts.dot(self.A1, y)
781 # Apply the filter
782 if self.filter:
783 u0 = tt.zeros_like(u)
784 u0 = tt.set_subtensor(u0[0], -1.0)
785 A1y = ifelse(
786 ld, tt.dot(self.F(u, f), A1y), tt.dot(self.F(u0, f), A1y)
787 )
789 # Dot the polynomial into the basis
790 intensity = tt.shape_padright(tt.dot(pT, A1y))
792 # Weight the intensity by the illumination
793 xyz = tt.concatenate(
794 (
795 tt.reshape(xpt, [1, -1]),
796 tt.reshape(ypt, [1, -1]),
797 tt.reshape(zpt, [1, -1]),
798 )
799 )
800 I = self.compute_illumination(xyz, xs, ys, zs)
802 # Add an extra dimension for the wavelength
803 if self.nw is not None:
804 I = tt.shape_padaxis(I, 1)
806 intensity = tt.switch(tt.isnan(intensity), intensity, intensity * I)
807 return intensity
809 @autocompile
810 def unweighted_intensity(self, lat, lon, y, u, f, ld):
811 """Compute the intensity in the absence of an illumination source."""
812 # Get the Cartesian points
813 xpt, ypt, zpt = self.latlon_to_xyz(lat, lon)
815 # Compute the polynomial basis at the point
816 pT = self.pT(xpt, ypt, zpt)
818 # Transform the map to the polynomial basis
819 A1y = ts.dot(self.A1, y)
821 # Apply the filter
822 if self.filter:
823 u0 = tt.zeros_like(u)
824 u0 = tt.set_subtensor(u0[0], -1.0)
825 A1y = ifelse(
826 ld, tt.dot(self.F(u, f), A1y), tt.dot(self.F(u0, f), A1y)
827 )
829 # Dot the polynomial into the basis
830 return tt.dot(pT, A1y)
832 @autocompile
833 def X(self, theta, xs, ys, zs, xo, yo, zo, ro, inc, obl, u, f, alpha):
834 """Compute the light curve design matrix."""
835 # Determine shapes
836 rows = theta.shape[0]
837 cols = (self.ydeg + 1) ** 2
838 X = tt.zeros((rows, cols))
840 # Compute the occultation mask
841 b = tt.sqrt(xo ** 2 + yo ** 2)
842 b_rot = tt.ge(b, 1.0 + ro) | tt.le(zo, 0.0) | tt.eq(ro, 0.0)
843 b_occ = tt.invert(b_rot)
844 i_rot = tt.arange(b.size)[b_rot]
845 i_occ = tt.arange(b.size)[b_occ]
847 # Rotation operator
848 # -----------------
850 # Compute the semi-minor axis of the terminator
851 # and the reflectance integrals
852 r2 = xs[i_rot] ** 2 + ys[i_rot] ** 2 + zs[i_rot] ** 2
853 bterm = -zs[i_rot] / tt.sqrt(r2)
854 rT = self.rT(bterm)
856 # Transform to Ylms and rotate on the sky plane
857 rTA1 = ts.dot(rT, self.A1Big)
858 theta_z = tt.arctan2(xs[i_rot], ys[i_rot])
859 rTA1Rz = self.tensordotRz(rTA1, theta_z)
861 # Apply limb darkening?
862 F = self.F(u, f)
863 A1InvFA1 = ts.dot(ts.dot(self.A1Inv, F), self.A1)
864 rTA1Rz = tt.dot(rTA1Rz, A1InvFA1)
866 # Rotate to the correct phase and weight by the distance to the source
867 # The factor of 2/3 ensures that the flux from a uniform map
868 # with unit amplitude seen at noon is unity.
869 X = tt.set_subtensor(
870 X[i_rot],
871 self.right_project(rTA1Rz, inc, obl, theta[i_rot], alpha)
872 / (2.0 / 3.0 * tt.shape_padright(r2)),
873 )
875 # Occultation operator
876 # --------------------
878 X = tt.set_subtensor(
879 X[i_occ],
880 RaiseValueErrorIfOp(
881 "Occultations in reflected light not yet implemented."
882 )(b_occ.any()),
883 )
885 # We're done
886 return X
888 @autocompile
889 def flux(
890 self, theta, xs, ys, zs, xo, yo, zo, ro, inc, obl, y, u, f, alpha
891 ):
892 """Compute the reflected light curve."""
893 return tt.dot(
894 self.X(theta, xs, ys, zs, xo, yo, zo, ro, inc, obl, u, f, alpha), y
895 )
897 @autocompile
898 def render(
899 self,
900 res,
901 projection,
902 illuminate,
903 theta,
904 inc,
905 obl,
906 y,
907 u,
908 f,
909 alpha,
910 xs,
911 ys,
912 zs,
913 ):
914 """Render the map on a Cartesian grid."""
915 # Compute the Cartesian grid
916 xyz = ifelse(
917 tt.eq(projection, STARRY_RECTANGULAR_PROJECTION),
918 self.compute_rect_grid(res),
919 self.compute_ortho_grid(res),
920 )
922 # Compute the polynomial basis
923 pT = self.pT(xyz[0], xyz[1], xyz[2])
925 # If orthographic, rotate the map to the correct frame
926 if self.nw is None:
927 Ry = ifelse(
928 tt.eq(projection, STARRY_ORTHOGRAPHIC_PROJECTION),
929 self.left_project(
930 tt.transpose(tt.tile(y, [theta.shape[0], 1])),
931 inc,
932 obl,
933 theta,
934 alpha,
935 ),
936 tt.transpose(
937 self.tensordotD(
938 tt.tile(y, [theta.shape[0], 1]), theta * alpha
939 )
940 )
941 if self.diffrot
942 else tt.transpose(tt.tile(y, [theta.shape[0], 1])),
943 )
944 else:
945 Ry = ifelse(
946 tt.eq(projection, STARRY_ORTHOGRAPHIC_PROJECTION),
947 self.left_project(
948 y, inc, obl, tt.tile(theta[0], self.nw), alpha
949 ),
950 y,
951 )
953 # Transform to polynomials
954 A1Ry = ts.dot(self.A1, Ry)
956 # Apply the filter *only if orthographic*
957 f0 = tt.zeros_like(f)
958 f0 = tt.set_subtensor(f0[0], np.pi)
959 u0 = tt.zeros_like(u)
960 u0 = tt.set_subtensor(u0[0], -1.0)
961 A1Ry = ifelse(
962 tt.eq(projection, STARRY_ORTHOGRAPHIC_PROJECTION),
963 tt.dot(self.F(u, f), A1Ry),
964 tt.dot(self.F(u0, f0), A1Ry),
965 )
967 # Dot the polynomial into the basis
968 image = tt.dot(pT, A1Ry)
970 # Compute the illumination profile
971 I = self.compute_illumination(xyz, xs, ys, zs)
973 # Add an extra dimension for the wavelength
974 if self.nw is not None:
975 I = tt.repeat(I, self.nw, 1)
977 # Weight the image by the illumination
978 image = ifelse(
979 illuminate, tt.switch(tt.isnan(image), image, image * I), image
980 )
982 # We need the shape to be (nframes, npix, npix)
983 return tt.reshape(image, [res, res, -1]).dimshuffle(2, 0, 1)
985 @autocompile
986 def compute_illumination(self, xyz, xs, ys, zs):
987 """Compute the illumination profile when rendering maps."""
988 r2 = xs ** 2 + ys ** 2 + zs ** 2
989 b = -zs / tt.sqrt(r2) # semi-minor axis of terminator
990 invsr = 1.0 / tt.sqrt(xs ** 2 + ys ** 2)
991 cosw = ys * invsr
992 sinw = -xs * invsr
993 xrot = (
994 tt.shape_padright(xyz[0]) * cosw + tt.shape_padright(xyz[1]) * sinw
995 )
996 yrot = (
997 -tt.shape_padright(xyz[0]) * sinw
998 + tt.shape_padright(xyz[1]) * cosw
999 )
1000 I = tt.sqrt(1.0 - b ** 2) * yrot - b * tt.shape_padright(xyz[2])
1001 I = tt.switch(
1002 tt.eq(tt.abs_(b), 1.0),
1003 tt.switch(
1004 tt.eq(b, 1.0),
1005 tt.zeros_like(I), # midnight
1006 tt.shape_padright(xyz[2]), # noon
1007 ),
1008 I,
1009 )
1010 I = tt.switch(tt.gt(I, 0.0), I, tt.zeros_like(I)) # set night to zero
1012 # Weight by the distance to the source
1013 # The factor of 2/3 ensures that the flux from a uniform map
1014 # with unit amplitude seen at noon is unity.
1015 I /= (2.0 / 3.0) * tt.shape_padleft(r2)
1016 return I
1019class OpsSystem(object):
1020 """Class housing ops for modeling Keplerian systems."""
1022 def __init__(
1023 self,
1024 primary,
1025 secondaries,
1026 reflected=False,
1027 rv=False,
1028 light_delay=False,
1029 texp=None,
1030 oversample=7,
1031 order=0,
1032 ):
1033 # System members
1034 self.primary = primary
1035 self.secondaries = secondaries
1036 self._reflected = reflected
1037 self._rv = rv
1038 self.nw = self.primary._map.nw
1039 self.light_delay = light_delay
1040 self.texp = texp
1041 self.oversample = oversample
1042 self.order = order
1044 # Require exoplanet
1045 assert exoplanet is not None, "This class requires exoplanet >= 0.2.0."
1047 @autocompile
1048 def position(
1049 self,
1050 t,
1051 pri_m,
1052 pri_t0,
1053 sec_m,
1054 sec_t0,
1055 sec_porb,
1056 sec_ecc,
1057 sec_w,
1058 sec_Omega,
1059 sec_iorb,
1060 ):
1061 """Compute the Cartesian positions of all bodies."""
1062 orbit = exoplanet.orbits.KeplerianOrbit(
1063 period=sec_porb,
1064 t0=sec_t0,
1065 incl=sec_iorb,
1066 ecc=sec_ecc,
1067 omega=sec_w,
1068 Omega=sec_Omega,
1069 m_planet=sec_m,
1070 m_star=pri_m,
1071 r_star=1.0, # doesn't matter
1072 )
1073 # Position of the primary
1074 if len(self.secondaries) > 1:
1075 x_pri, y_pri, z_pri = tt.sum(orbit.get_star_position(t), axis=-1)
1076 else:
1077 x_pri, y_pri, z_pri = orbit.get_star_position(t)
1079 # Positions of the secondaries
1080 try:
1081 x_sec, y_sec, z_sec = orbit.get_planet_position(
1082 t, light_delay=self.light_delay
1083 )
1084 except TypeError:
1085 if self.light_delay:
1086 logger.warn(
1087 "This version of `exoplanet` does not model light delays."
1088 )
1089 x_sec, y_sec, z_sec = orbit.get_planet_position(t)
1091 # Concatenate them
1092 x = tt.transpose(tt.concatenate((x_pri, x_sec), axis=-1))
1093 y = tt.transpose(tt.concatenate((y_pri, y_sec), axis=-1))
1094 z = tt.transpose(tt.concatenate((z_pri, z_sec), axis=-1))
1096 return x, y, z
1098 @autocompile
1099 def X(
1100 self,
1101 t,
1102 pri_r,
1103 pri_m,
1104 pri_prot,
1105 pri_t0,
1106 pri_theta0,
1107 pri_L,
1108 pri_inc,
1109 pri_obl,
1110 pri_u,
1111 pri_f,
1112 pri_alpha,
1113 sec_r,
1114 sec_m,
1115 sec_prot,
1116 sec_t0,
1117 sec_theta0,
1118 sec_porb,
1119 sec_ecc,
1120 sec_w,
1121 sec_Omega,
1122 sec_iorb,
1123 sec_L,
1124 sec_inc,
1125 sec_obl,
1126 sec_u,
1127 sec_f,
1128 sec_alpha,
1129 ):
1130 """Compute the system light curve design matrix."""
1131 # Exposure time integration?
1132 if self.texp != 0.0:
1134 texp = tt.as_tensor_variable(self.texp)
1135 oversample = int(self.oversample)
1136 oversample += 1 - oversample % 2
1137 stencil = np.ones(oversample)
1139 # Construct the exposure time integration stencil
1140 if self.order == 0:
1141 dt = np.linspace(-0.5, 0.5, 2 * oversample + 1)[1:-1:2]
1142 elif self.order == 1:
1143 dt = np.linspace(-0.5, 0.5, oversample)
1144 stencil[1:-1] = 2
1145 elif self.order == 2:
1146 dt = np.linspace(-0.5, 0.5, oversample)
1147 stencil[1:-1:2] = 4
1148 stencil[2:-1:2] = 2
1149 else:
1150 raise ValueError("Parameter `order` must be <= 2")
1151 stencil /= np.sum(stencil)
1153 if texp.ndim == 0:
1154 dt = texp * dt
1155 else:
1156 dt = tt.shape_padright(texp) * dt
1157 t = tt.shape_padright(t) + dt
1158 t = tt.reshape(t, (-1,))
1160 # Compute the relative positions of all bodies
1161 orbit = exoplanet.orbits.KeplerianOrbit(
1162 period=sec_porb,
1163 t0=sec_t0,
1164 incl=sec_iorb,
1165 ecc=sec_ecc,
1166 omega=sec_w,
1167 Omega=sec_Omega,
1168 m_planet=sec_m,
1169 m_star=pri_m,
1170 r_star=pri_r,
1171 )
1172 try:
1173 x, y, z = orbit.get_relative_position(
1174 t, light_delay=self.light_delay
1175 )
1176 except TypeError:
1177 if self.light_delay:
1178 logger.warn(
1179 "This version of `exoplanet` does not model light delays."
1180 )
1181 x, y, z = orbit.get_relative_position(t)
1183 # Get all rotational phases
1184 pri_prot = ifelse(
1185 tt.eq(pri_prot, 0.0), math.to_tensor(np.inf), pri_prot
1186 )
1187 theta_pri = (2 * np.pi) / pri_prot * (t - pri_t0) + pri_theta0
1188 sec_prot = tt.switch(
1189 tt.eq(sec_prot, 0.0), math.to_tensor(np.inf), sec_prot
1190 )
1191 theta_sec = (2 * np.pi) / tt.shape_padright(sec_prot) * (
1192 tt.shape_padleft(t) - tt.shape_padright(sec_t0)
1193 ) + tt.shape_padright(sec_theta0)
1195 # Compute all the phase curves
1196 phase_pri = pri_L * self.primary.map.ops.X(
1197 theta_pri,
1198 tt.zeros_like(t),
1199 tt.zeros_like(t),
1200 tt.zeros_like(t),
1201 math.to_tensor(0.0),
1202 pri_inc,
1203 pri_obl,
1204 pri_u,
1205 pri_f,
1206 pri_alpha,
1207 )
1208 if self._reflected:
1209 phase_sec = [
1210 pri_L
1211 * sec_L[i]
1212 * sec.map.ops.X(
1213 theta_sec[i],
1214 -x[:, i],
1215 -y[:, i],
1216 -z[:, i],
1217 -x[:, i], # not used
1218 -y[:, i], # not used
1219 -z[:, i], # not used, since...
1220 math.to_tensor(0.0), # occultor of zero radius
1221 sec_inc[i],
1222 sec_obl[i],
1223 sec_u[i],
1224 sec_f[i],
1225 sec_alpha[i],
1226 )
1227 for i, sec in enumerate(self.secondaries)
1228 ]
1229 else:
1230 phase_sec = [
1231 sec_L[i]
1232 * sec.map.ops.X(
1233 theta_sec[i],
1234 -x[:, i],
1235 -y[:, i],
1236 -z[:, i],
1237 math.to_tensor(0.0), # occultor of zero radius
1238 sec_inc[i],
1239 sec_obl[i],
1240 sec_u[i],
1241 sec_f[i],
1242 sec_alpha[i],
1243 )
1244 for i, sec in enumerate(self.secondaries)
1245 ]
1247 # Compute any occultations
1248 occ_pri = tt.zeros_like(phase_pri)
1249 occ_sec = [tt.zeros_like(ps) for ps in phase_sec]
1251 # Compute the period if we were given a semi-major axis
1252 sec_porb = tt.switch(
1253 tt.eq(sec_porb, 0.0),
1254 (G_grav * (pri_m + sec_m) * sec_porb ** 2 / (4 * np.pi ** 2))
1255 ** (1.0 / 3),
1256 sec_porb,
1257 )
1259 # Compute transits across the primary
1260 for i, _ in enumerate(self.secondaries):
1261 xo = x[:, i] / pri_r
1262 yo = y[:, i] / pri_r
1263 zo = z[:, i] / pri_r
1264 ro = sec_r[i] / pri_r
1265 b = tt.sqrt(xo ** 2 + yo ** 2)
1266 b_occ = tt.invert(
1267 tt.ge(b, 1.0 + ro) | tt.le(zo, 0.0) | tt.eq(ro, 0.0)
1268 )
1269 idx = tt.arange(b.shape[0])[b_occ]
1270 occ_pri = tt.set_subtensor(
1271 occ_pri[idx],
1272 occ_pri[idx]
1273 + pri_L
1274 * self.primary.map.ops.X(
1275 theta_pri[idx],
1276 xo[idx],
1277 yo[idx],
1278 zo[idx],
1279 ro,
1280 pri_inc,
1281 pri_obl,
1282 pri_u,
1283 pri_f,
1284 pri_alpha,
1285 )
1286 - phase_pri[idx],
1287 )
1289 # Compute occultations by the primary
1290 for i, sec in enumerate(self.secondaries):
1291 xo = -x[:, i] / sec_r[i]
1292 yo = -y[:, i] / sec_r[i]
1293 zo = -z[:, i] / sec_r[i]
1294 ro = pri_r / sec_r[i]
1295 b = tt.sqrt(xo ** 2 + yo ** 2)
1296 b_occ = tt.invert(
1297 tt.ge(b, 1.0 + ro) | tt.le(zo, 0.0) | tt.eq(ro, 0.0)
1298 )
1299 idx = tt.arange(b.shape[0])[b_occ]
1300 if self._reflected:
1301 occ_sec[i] = tt.set_subtensor(
1302 occ_sec[i][idx],
1303 occ_sec[i][idx]
1304 + pri_L
1305 * sec_L[i]
1306 * sec.map.ops.X(
1307 theta_sec[i, idx],
1308 xo[idx], # the primary is both the source...
1309 yo[idx],
1310 zo[idx],
1311 xo[idx], # ... and the occultor
1312 yo[idx],
1313 zo[idx],
1314 ro,
1315 sec_inc[i],
1316 sec_obl[i],
1317 sec_u[i],
1318 sec_f[i],
1319 sec_alpha[i],
1320 )
1321 - phase_sec[i][idx],
1322 )
1323 else:
1324 occ_sec[i] = tt.set_subtensor(
1325 occ_sec[i][idx],
1326 occ_sec[i][idx]
1327 + sec_L[i]
1328 * sec.map.ops.X(
1329 theta_sec[i, idx],
1330 xo[idx],
1331 yo[idx],
1332 zo[idx],
1333 ro,
1334 sec_inc[i],
1335 sec_obl[i],
1336 sec_u[i],
1337 sec_f[i],
1338 sec_alpha[i],
1339 )
1340 - phase_sec[i][idx],
1341 )
1343 # Compute secondary-secondary occultations
1344 for i, sec in enumerate(self.secondaries):
1345 for j, _ in enumerate(self.secondaries):
1346 if i == j:
1347 continue
1348 xo = (-x[:, i] + x[:, j]) / sec_r[i]
1349 yo = (-y[:, i] + y[:, j]) / sec_r[i]
1350 zo = (-z[:, i] + z[:, j]) / sec_r[i]
1351 ro = sec_r[j] / sec_r[i]
1352 b = tt.sqrt(xo ** 2 + yo ** 2)
1353 b_occ = tt.invert(
1354 tt.ge(b, 1.0 + ro) | tt.le(zo, 0.0) | tt.eq(ro, 0.0)
1355 )
1356 idx = tt.arange(b.shape[0])[b_occ]
1357 if self._reflected:
1358 xs = -x[:, i] / sec_r[i]
1359 ys = -y[:, i] / sec_r[i]
1360 zs = -z[:, i] / sec_r[i]
1361 occ_sec[i] = tt.set_subtensor(
1362 occ_sec[i][idx],
1363 occ_sec[i][idx]
1364 + sec_L[i]
1365 * pri_L
1366 * sec.map.ops.X(
1367 theta_sec[i, idx],
1368 xs[idx], # the primary is the source
1369 ys[idx],
1370 zs[idx],
1371 xo[idx], # another secondary is the occultor
1372 yo[idx],
1373 zo[idx],
1374 ro,
1375 sec_inc[i],
1376 sec_obl[i],
1377 sec_u[i],
1378 sec_f[i],
1379 sec_alpha[i],
1380 )
1381 - phase_sec[i][idx],
1382 )
1383 else:
1384 occ_sec[i] = tt.set_subtensor(
1385 occ_sec[i][idx],
1386 occ_sec[i][idx]
1387 + sec_L[i]
1388 * sec.map.ops.X(
1389 theta_sec[i, idx],
1390 xo[idx],
1391 yo[idx],
1392 zo[idx],
1393 ro,
1394 sec_inc[i],
1395 sec_obl[i],
1396 sec_u[i],
1397 sec_f[i],
1398 sec_alpha[i],
1399 )
1400 - phase_sec[i][idx],
1401 )
1403 # Concatenate the design matrices
1404 X_pri = phase_pri + occ_pri
1405 X_sec = [ps + os for ps, os in zip(phase_sec, occ_sec)]
1406 X = tt.horizontal_stack(X_pri, *X_sec)
1408 # Sum and return
1409 if self.texp == 0.0:
1411 return X
1413 else:
1415 stencil = tt.shape_padright(tt.shape_padleft(stencil, 1), 1)
1416 return tt.sum(
1417 stencil * tt.reshape(X, (-1, self.oversample, X.shape[1])),
1418 axis=1,
1419 )
1421 @autocompile
1422 def rv(
1423 self,
1424 t,
1425 pri_r,
1426 pri_m,
1427 pri_prot,
1428 pri_t0,
1429 pri_theta0,
1430 pri_L,
1431 pri_inc,
1432 pri_obl,
1433 pri_y,
1434 pri_u,
1435 pri_alpha,
1436 pri_veq,
1437 sec_r,
1438 sec_m,
1439 sec_prot,
1440 sec_t0,
1441 sec_theta0,
1442 sec_porb,
1443 sec_ecc,
1444 sec_w,
1445 sec_Omega,
1446 sec_iorb,
1447 sec_L,
1448 sec_inc,
1449 sec_obl,
1450 sec_y,
1451 sec_u,
1452 sec_alpha,
1453 sec_veq,
1454 keplerian,
1455 ):
1456 """Compute the observed system radial velocity (RV maps only)."""
1457 # TODO: This method is currently very inefficient, as it
1458 # calls `X` twice per call and instantiates an `orbit`
1459 # instance up to three separate times per call. We should
1460 # re-code the logic from `X()` in here to optimize it.
1462 # Compute the RV filter
1463 pri_f = self.primary.map.ops.compute_rv_filter(
1464 pri_inc, pri_obl, pri_veq, pri_alpha
1465 )
1466 sec_f = tt.as_tensor_variable(
1467 [
1468 sec.map.ops.compute_rv_filter(
1469 sec_inc[k], sec_obl[k], sec_veq[k], sec_alpha[k]
1470 )
1471 for k, sec in enumerate(self.secondaries)
1472 ]
1473 )
1475 # Compute the identity filter
1476 pri_f0 = tt.zeros_like(pri_f)
1477 pri_f0 = tt.set_subtensor(pri_f0[0], np.pi)
1478 sec_f0 = tt.as_tensor_variable([pri_f0 for sec in self.secondaries])
1480 # Compute the two design matrices
1481 X = self.X(
1482 t,
1483 pri_r,
1484 pri_m,
1485 pri_prot,
1486 pri_t0,
1487 pri_theta0,
1488 pri_L,
1489 pri_inc,
1490 pri_obl,
1491 pri_u,
1492 pri_f,
1493 pri_alpha,
1494 sec_r,
1495 sec_m,
1496 sec_prot,
1497 sec_t0,
1498 sec_theta0,
1499 sec_porb,
1500 sec_ecc,
1501 sec_w,
1502 sec_Omega,
1503 sec_iorb,
1504 sec_L,
1505 sec_inc,
1506 sec_obl,
1507 sec_u,
1508 sec_f,
1509 sec_alpha,
1510 )
1512 X0 = self.X(
1513 t,
1514 pri_r,
1515 pri_m,
1516 pri_prot,
1517 pri_t0,
1518 pri_theta0,
1519 pri_L,
1520 pri_inc,
1521 pri_obl,
1522 pri_u,
1523 pri_f0,
1524 pri_alpha,
1525 sec_r,
1526 sec_m,
1527 sec_prot,
1528 sec_t0,
1529 sec_theta0,
1530 sec_porb,
1531 sec_ecc,
1532 sec_w,
1533 sec_Omega,
1534 sec_iorb,
1535 sec_L,
1536 sec_inc,
1537 sec_obl,
1538 sec_u,
1539 sec_f0,
1540 sec_alpha,
1541 )
1543 # Get the indices of X corresponding to each body
1544 pri_inds = np.arange(0, self.primary.map.Ny, dtype=int)
1545 sec_inds = [None for sec in self.secondaries]
1546 n = self.primary.map.Ny
1547 for i, sec in enumerate(self.secondaries):
1548 sec_inds[i] = np.arange(n, n + sec.map.Ny, dtype=int)
1549 n += sec.map.Ny
1551 # Compute the integral of the velocity-weighted intensity
1552 Iv = tt.as_tensor_variable(
1553 [tt.dot(X[:, pri_inds], pri_y)]
1554 + [
1555 tt.dot(X[:, sec_inds[n]], sec_y[n])
1556 for n in range(len(self.secondaries))
1557 ]
1558 )
1560 # Compute the inverse of the integral of the intensity
1561 invI = tt.as_tensor_variable(
1562 [tt.ones((1,)) / tt.dot(X0[:, pri_inds], pri_y)]
1563 + [
1564 tt.ones((1,)) / tt.dot(X0[:, sec_inds[n]], sec_y[n])
1565 for n in range(len(self.secondaries))
1566 ]
1567 )
1568 invI = tt.where(tt.isinf(invI), 0.0, invI)
1570 # The RV anomaly is just the product
1571 rv = Iv * invI
1573 # Compute the Keplerian RV
1574 orbit = exoplanet.orbits.KeplerianOrbit(
1575 period=sec_porb,
1576 t0=sec_t0,
1577 incl=sec_iorb,
1578 ecc=sec_ecc,
1579 omega=sec_w,
1580 Omega=sec_Omega,
1581 m_planet=sec_m,
1582 m_star=pri_m,
1583 r_star=pri_r,
1584 )
1585 return ifelse(
1586 keplerian,
1587 tt.inc_subtensor(
1588 rv[0],
1589 tt.reshape(
1590 orbit.get_radial_velocity(
1591 t, output_units=units.m / units.s
1592 ),
1593 (-1,),
1594 ),
1595 ),
1596 rv,
1597 )
1599 @autocompile
1600 def render(
1601 self,
1602 t,
1603 res,
1604 pri_r,
1605 pri_m,
1606 pri_prot,
1607 pri_t0,
1608 pri_theta0,
1609 pri_L,
1610 pri_inc,
1611 pri_obl,
1612 pri_y,
1613 pri_u,
1614 pri_f,
1615 pri_alpha,
1616 sec_r,
1617 sec_m,
1618 sec_prot,
1619 sec_t0,
1620 sec_theta0,
1621 sec_porb,
1622 sec_ecc,
1623 sec_w,
1624 sec_Omega,
1625 sec_iorb,
1626 sec_L,
1627 sec_inc,
1628 sec_obl,
1629 sec_y,
1630 sec_u,
1631 sec_f,
1632 sec_alpha,
1633 ):
1634 """Render all of the bodies in the system."""
1635 # Compute the relative positions of all bodies
1636 orbit = exoplanet.orbits.KeplerianOrbit(
1637 period=sec_porb,
1638 t0=sec_t0,
1639 incl=sec_iorb,
1640 ecc=sec_ecc,
1641 omega=sec_w,
1642 Omega=sec_Omega,
1643 m_planet=sec_m,
1644 m_star=pri_m,
1645 r_star=pri_r,
1646 )
1647 try:
1648 x, y, z = orbit.get_relative_position(
1649 t, light_delay=self.light_delay
1650 )
1651 except TypeError:
1652 if self.light_delay:
1653 logger.warn(
1654 "This version of `exoplanet` does not model light delays."
1655 )
1656 x, y, z = orbit.get_relative_position(t)
1658 # Get all rotational phases
1659 pri_prot = ifelse(
1660 tt.eq(pri_prot, 0.0), math.to_tensor(np.inf), pri_prot
1661 )
1662 theta_pri = (2 * np.pi) / pri_prot * (t - pri_t0) + pri_theta0
1663 sec_prot = tt.switch(
1664 tt.eq(sec_prot, 0.0), math.to_tensor(np.inf), sec_prot
1665 )
1666 theta_sec = (2 * np.pi) / tt.shape_padright(sec_prot) * (
1667 tt.shape_padleft(t) - tt.shape_padright(sec_t0)
1668 ) + tt.shape_padright(sec_theta0)
1670 # Compute all the maps
1671 img_pri = self.primary.map.ops.render(
1672 res,
1673 STARRY_ORTHOGRAPHIC_PROJECTION,
1674 theta_pri,
1675 pri_inc,
1676 pri_obl,
1677 pri_y,
1678 pri_u,
1679 pri_f,
1680 pri_alpha,
1681 )
1682 if self._reflected:
1683 img_sec = tt.as_tensor_variable(
1684 [
1685 sec.map.ops.render(
1686 res,
1687 STARRY_ORTHOGRAPHIC_PROJECTION,
1688 theta_sec[i],
1689 sec_inc[i],
1690 sec_obl[i],
1691 sec_y[i],
1692 sec_u[i],
1693 sec_f[i],
1694 sec_alpha[i],
1695 -x[:, i],
1696 -y[:, i],
1697 -z[:, i],
1698 )
1699 for i, sec in enumerate(self.secondaries)
1700 ]
1701 )
1702 else:
1703 img_sec = tt.as_tensor_variable(
1704 [
1705 sec.map.ops.render(
1706 res,
1707 STARRY_ORTHOGRAPHIC_PROJECTION,
1708 theta_sec[i],
1709 sec_inc[i],
1710 sec_obl[i],
1711 sec_y[i],
1712 sec_u[i],
1713 sec_f[i],
1714 sec_alpha[i],
1715 )
1716 for i, sec in enumerate(self.secondaries)
1717 ]
1718 )
1720 # Return the images and secondary orbital positions
1721 return img_pri, img_sec, x, y, z