ray_draping
File Purpose: “field-of-view” grid of rays through a box draped across a sphere.
Useful for making synthetic observations based on a plane-parallel (uncurved) simulation.
The overarching idea is:
(1) Pick a FOV, like a grid of points in a rectangle on the sky.
(2) Anchor the simulation box to the surface of a sphere.
(3) “Drape” the box across the sphere, respecting curvature.
(4) Trace rays from observer, through FOV, through draped box.
There are two “coordinate” systems (though, chosen coords can vary…) to keep in mind:
(A) “SPH”: describes the situation “as understood by the sphere”.
Rays are all straight lines in this system.
Best coordinate systems for SPH:
- spherical coordinates centered at sphere center
- cartesian coordinates centered anywhere (but probably at sphere center)
(B) “BOX”: describes the situation “as understood by the (uncurved) simulation box”.
Rays are generally curved paths in this system.
Best coordinate systems for BOX:
- cartesian coordinates directly from simulation
- cartesian coordinates centered at anchor point
There are various miscellaneous intracies:
(1) Choosing FOV is simple enough. Using xarray makes it easy to have a grid of points.
(2) Choosing anchor is simple enough, just need to pick a point in both coord systems.
e.g. Abox=(x=X0+XEXT/2, y=Y0+YEXT/2, z=0); Asph=(r=R0, θ=(user choice), φ=(user choice)),
where R0 = radius of sphere,
XEXT, YEXT are x, y extents of the box,
and X0, Y0 are x, y coordinates of lower-left corner of the box, default (0,0).
(3) Draping the simulation is conceptually challenging but can be figured out,
and the resulting formulas “aren’t too bad”. Here’s what we want:
(a) (going up in BOX) == (going same distance radially outward in SPH).
–> zbox - Abox_z = rsph - Asph_r
–> zbox = rsph - R0 (assuming Abox_z=0 and Asph_r=R0 as done in this file.)
(b) (going horizontally in BOX at z=0) == (going same distance along sphere surface)
Combining (a) and (b) almost fully specifies the mapping constraints,
but there is one more choice to make: the direction for horizontal “motion”.
Basically, we need to pick two orthogonal unit vectors in the SPH system,
in the plane tangent to the sphere at Asph (“tangent” ensures no changes in rsph).
There are many ways to do this, but for this module let’s choose:
(c) (θhat_Asph <–> box x motion); (φhat_Asph <–> box y motion)
Where unit vectors are defined in the usual way:
rhat_Asph = sin(Asph_θ) cos(Asph_φ) xhat_sph
+ sin(Asph_θ) sin(Asph_φ) yhat_sph
+ cos(Asph_θ) zhat_sph
θhat_Asph = cos(Asph_θ) cos(Asph_φ) xhat_sph
+ cos(Asph_θ) sin(Asph_φ) yhat_sph
- sin(Asph_θ) zhat_sph
φhat_Asph = - sin(Asph_φ) xhat_sph
+ cos(Asph_φ) yhat_sph
(TODO: allow any e1hat and e2hat for x and y, mixing θhat_Asph and φhat_Asph.)
Now we have fully specified the mapping constraints.
(4) For tracing rays, we need to know when the ray overlaps with the “draped box” region.
Assuming we are allowed to periodically tile the box in x and y (and, around the sphere),
this is easy enough. For a ray expressed parametrically as w(s) (vector),
the ray overlaps the draped box region for all s such that:
R0 <= |w(s)| <= R0 + ZEXT,
where ZEXT is the extent of the simulation box above z=0.
(If desired, the lower bound can be decreased to allow boxes which extend “below the surface”.)
This expression might be solved analytically to choose s,
or it might just be imposed as a condition (e.g., mask values where it is False).
A bit more mapping logic:
For clarity, let’s first define Vsph:
Vsph(xbox,ybox) = (xbox - Abox_x) θhat_Asph + (ybox - Abox_y) φhat_Asph
Consider the tangent plane at Asph, before “draping” it onto the surface;
the point in that plane corresponding to (xbox,ybox) in the simulation is:
vsph(xbox,ybox,zbox=0) = Asph + Vsph(xbox,ybox)
This is a conceptually useful definition which will reappear in the formulas.
The mapping itself can now be determined algebraically by applying Rodrigues rotation formula,
with angle:
α = |Vsph(xbox,ybox)|/R0
around axis of rotation:
ghat_sph = rhat_Asph x Vsph_hat
(where Vsph_hat = Vsph/|Vsph|)
This determines the mapping because of (3b) above; “going a distance along the surface”
corresponds to “rotating around by the correct amount, around the axis which is perpendicular
to the plane containing the path as well as the sphere center.” That axis is determined by
rhat_Asph x Vhat_sph because those are two orthogonal vectors which lie in that plane.
Finding the inverse mapping just involves algebra, starting from the mapping itself.
The mapping itself:
ANCHOR POINT UNIT VECTORS (these are constants):
rhat_Asph = sin(Asph_θ) cos(Asph_φ) xhat_sph
+ sin(Asph_θ) sin(Asph_φ) yhat_sph
- cos(Asph_θ) zhat_sph
θhat_Asph = cos(Asph_θ) cos(Asph_φ) xhat_sph
+ cos(Asph_θ) sin(Asph_φ) yhat_sph
- sin(Asph_θ) zhat_sph
φhat_Asph = - sin(Asph_φ) xhat_sph
+ cos(Asph_φ) yhat_sph
BOX TO SPH: (xbox,ybox,zbox) –> wsph (vector):
rsph = R0 + zbox
Vsph = (xbox - Abox_x) θhat_Asph + (ybox - Abox_y) φhat_Asph
Vsph_hat = Vsph/|Vsph|
α = |Vsph|/R0
wsph = (cos(α) rhat_Asph + sin(α) Vsph_hat) rsph
SPH TO BOX: wsph (vector) –> (xbox,ybox,zbox):
rsph = |wsph|
|Vsph| = R0 cos^-1((wsph . rhat_Asph)/rsph)
wsph_perp = wsph - (wsph . rhat_Asph) rhat_Asph == rsph sin(α) Vsph_hat
Vsph_hat = wsph_perp/|wsph_perp|
Vsph = |Vsph| Vsph_hat
xbox = Abox_x + (Vsph . θhat_Asph)
ybox = Abox_y + (Vsph . φhat_Asph)
zbox = rsph - R0
TRICKY ANGLES:
When wsph_perp –> 0, it is impossible to get Vsph_hat;
this corresponds to the point xbox = Abox_x, ybox = Abox_y.
–> we need a test like “if |wsph_perp| < ε, impose Vsph = 0.
(maybe ε should be some small fraction of simulation grid cell size.)
If (vsph . rhat_Asph)/rsph > 1, the inverse cosine will fail.
this shouldn’t occur anywhere though, |wsph| = |rsph| by definition,
and rhat_Asph is a unit vector, so the whole thing should be 1 or less.
Numerically might be an issue? Can double-check, but low priority.
No other tricky angles, the rest should be fine!
RAY FORMULA:
Rays in sph system are straight lines:
w(s) = F + rayhat * s
where:
rayhat = (F - T)/|F - T|,
F is the point in FOV corresponding to this ray (there is one ray per FOV point),
T is the telescope/observer position. (F and T both expressed in sph system.)
The corresponding paths in box system are curved by the mapping formula above.
Assuming periodic box tiling is allowed, the rays overlap the draped box region wherever:
R0 <= |w(s)| <= R0 + ZEXT
and the periodic box tiling can be accomplished by taking the modulus:
xbox = xbox % XEXT
ybox = ybox % YEXT
zbox = zbox % ZEXT
RAY BOUNDARIES:
The values of s where rays enter and exit the draped box region can be computed from:
w(s) = F + rayhat * s
rayhat = (F - T)/|F - T|
R0 <= |w(s)| <= R0 + ZEXT
For boundary R = R0 or R0 + ZEXT:
R^2 = F^2 + 2 (F . rayhat) s + s^2
–> s^2 + 2 (F . rayhat) s + (F^2 - R^2) = 0
–> s = - (F . rayhat) ± sqrt((F . rayhat)^2 - F^2 + R^2)
The full ray intersects the boundary whenever the discriminant is non-negative.
In that case, is important to know which point the ray passes through first,
i.e. which is closest to T. That point is the “entry”.
If there are two intersections of a boundary, the other is the “exit”.
Thankfully, by definition, the smaller s value is always the entry point.
(Not necessarily smaller in magnitude; s might be negative. E.g. -7 is smaller than 3.)
The full decision tree looks like:
- find lower and upper boundary entry and exit points. (Up to 4 distinct points).
- if lower entry == lower exit, ray just grazes lower boundary; ignore lower boundary.
- if upper entry == upper exit, ray just grazes upper boundary; replace ray with NaNs.
- if no upper entry/exit, ray doesn’t overlap sphere; replace ray with NaNs.
- if no lower entry/exit, ray goes from upper entry to upper exit. (E.g., near disk edge.)
- else, upper entry & exit, and lower entry & exit are all distinct,
so ray goes from upper entry to lower entry (e.g., anywhere other than disk edge),
and the exits both correspond to the draped region on the other side of the sphere.
NOTE: when actually getting the standard ray_s values, the default is inclusive=(False, False),
which avoids points “on the edges”. This way all ray_s correspond to points INSIDE the box.
Without doing this, some interp results may give NaNs due to points on edges.
FOV POSITION AND SIZE:
The FOV can be any size and be placed anywhere,
as long as the sphere does not block it from T (not checked directly in this module…),
and as long as it is a plane (not actually strictly necessary but probably a good idea).
In practice, it is nice to choose FOV location based on box.
How big should FOV be? Ideally, big enough to cover box.
This defines “standard” FOV (assuming T along zsph axis):
get points on box surfaces in sph system, project onto zsph=constant plane,
then choose FOV x and y limits as min and max x and y across all projected points.
Projecting all 8 corners is a good start but insufficient;
also need to check top surface which can bulge out beyond the corners.
Internal algorithm checks all 8 corners and a low-res grid (21x21 points) for top surface.
(Currently, only implemented if T along zsph axis. Eventually, might handle arbitrary T.)
Classes
|
representation of Anchor point for ray draping problem. |
|
representation of Field-of-View grid for ray draping problem. |
|
representation of some "background" conditions for the ray draping problem. |
|
representation of the entire ray draping problem, including mapping calculations. |
|
core methods for the entire ray draping problem, including mapping calculations. |