In [ ]:
import batoid
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
The most fundamental object in batoid is the Ray, which you can roughly think of as a photon. Rays are defined primarily by their position r and velocity vector v. Positions in batoid are always specified in meters, and velocities are always in units of the speed of light in vacuum. This makes the magnitude of the velocity vector equal to $1/n$, where $n$ is the refractive index of the medium in which the ray is currently propagating. If you're just planning on geometrically tracing rays through reflective optics, then positions and velocities are the only required Ray parameters (and are the only required elements in the Ray constructors).
The next most important Ray attribute after position and velocity is wavelength, which is (almost always) required for tracing through refractive optics. (The exception is if you define a refractive medium that doesn't depend on wavelength. Real media do have wavelength-dependent refractive indices of course.) Ray wavelengths are always specified in meters, and always in vacuum -- even when the Ray being created isn't currently in vacuum. If you don't specify a value for wavelength explicitly when constructing a Ray, a (somewhat nonsensical) value of 0.0 meters will be used.
Next in the list of Ray attributes is time t. Batoid doesn't track time directly, but rather tracks the product of time and the speed of light in vacuum. That product has dimensions of length, so time in batoid is measured in meters. Note that this convention is complementary to the convention of tracking velocity in units of the speed of light in vacuum -- we still have the familiar relation $\Delta \vec{r} = \vec{v} \Delta t$. The default value of time t for a Ray is 0.0 meters.
The final fundamental attributes of a Ray are a flux in arbitrary units and status indicators for whether a Ray has become vignetted during tracing, or if batoid failed to successfully trace a Ray for some reason. These and a number of additional derived attributes will be discussed more below.
In [ ]:
# Create a Ray from position and velocity 3-tuples (or 3-elements lists, arrays, ...)
x, y, z = 0.1, 0.2, 0.3 # meters
vx, vy, vz = 0.1, 0.2, np.sqrt(1 - 0.1**2 - 0.2**2) # in units of c
r = (x, y, z)
v = (vx, vy, vz)
ray = batoid.Ray(r, v)
In [ ]:
# Printing the ray reveals its attributes.
print(ray)
In [ ]:
# You can access the position and velocity attributes either as scalars or as vectors.
print(ray.x)
print(ray.vz)
print(ray.r)
print(ray.v)
assert ray.r[0] == ray.x
assert ray.r[1] == ray.y
assert ray.r[2] == ray.z
assert ray.v[0] == ray.vx
assert ray.v[1] == ray.vy
assert ray.v[2] == ray.vz
In [ ]:
# Rays have a number of optional fields too
print(ray.t) # The time in meters (see explanation above)
print(ray.wavelength) # Vacuum wavelength in meters.
print(ray.flux) # Flux in arbitrary units.
print(ray.vignetted) # Whether or not the ray has been vignetted
print(ray.failed) # Whether or not the ray is failed.
# Usually this means an intersection between the Ray and a Surface could not be found.
In [ ]:
# You can create rays with specific values for these optional fields too
# (except for failed, which always starts out as False)
ray = batoid.Ray(r, v, t=0.3, wavelength=400e-9, flux=2.3, vignetted=True)
print(ray)
Two additional attributes derivable from the fundamental attributes detailed above are the wavevector $\vec{k}$ and the temporal frequency $\omega$. The wavevector is defined by
$\vec{k} = \frac{2 \pi \vec{v}}{\lambda |\vec{v}|^2}$
(with $\lambda$ being the vacuum wavelength) and has units of radians per meter.
The normal physics definition of angular temporal frequency is $\omega = \frac{2 \pi c}{\lambda}$, but in batoid, we divide out the speed of light in vacuum to get
$\omega = \frac{2 \pi}{\lambda}$
The units are again radians per meter (but this time in a scalar variable). This choice for the frequency means that $\omega t$ is a phase angle in batoid (just as $\vec{k} \cdot \vec{r}$ is).
In [ ]:
print(ray.omega) # angular frequency
print(ray.k) # wave vector
print(ray.kx) # individual components of the wave vector are also accessible directly
assert ray.k[0] == ray.kx
assert ray.k[1] == ray.ky
assert ray.k[2] == ray.kz
The first set of methods available to Ray is for propagation of the photon through time. The positionAtTime method accepts a time (still in meters) and returns the position of the Ray at that past, present, or future time. The propagate method similarly propagates the ray to a given time. Note that the return value of propagate is a reference to the original Ray (not a copy).
For the second set of Ray methods, we change our conception of what a ray is: from a propagating infinitesimal point to a propagating plane wave. I.e., we imagine that the velocity v attribute is that of a plane wave, and that the position r and time t together indicate a point in space-time where the wave amplitude is a maximum. Note that these conditions will remain true as we propagate the wave in time.
The phase method accepts a position and time and returns the phase $\phi$ of this plane wave in radians. The amplitude method returns the complex amplitude $\exp(i \phi)$ .
In [ ]:
ray = batoid.Ray([0,0,0], [0,0,1])
print(ray)
print(ray.positionAtTime(1.0)) # propagated from (0,0,0) to (0,0,1)
ray2 = ray.propagate(1.0) # propagate ray forward to t=1
print(ray) # propagate works in-place
print(ray2) # the return value is a reference to the original
print(ray == ray2) # so they're equal
print(ray is ray2) # in fact, they're the same object
# If you want to preserver the original ray when performing a propagation, first make a copy.
ray3 = ray.copy().propagate(2.0)
print(ray)
print(ray3)
In [ ]:
# Plot one period of the plane wave with wavelength 800e-9 m
wavelength = 800e-9
ray = batoid.Ray([0,0,0], [0,0,1], t=0, wavelength=wavelength)
ts = np.linspace(0, wavelength)
plt.plot(ts, [np.cos(ray.phase([0,0,0], t)) for t in ts])
plt.axhline(0, c='k')
plt.xlabel("c t (meters)")
plt.ylabel("amplitude")
In [ ]:
# We can also look at a fixed time over a small region of space
# Since the wave is travelling in the +z direction, its amplitude
# varies with z, but not with x or y
wavelength = 800e-9
ray = batoid.Ray([0,0,0], [0,0,1], t=0, wavelength=wavelength)
xs = ys = zs = np.linspace(0, wavelength)
plt.plot(xs, [np.cos(ray.phase([x,0,0], 0)) for x in xs], label='x')
plt.plot(ys, [np.cos(ray.phase([0,y,0], 0))+0.1 for y in ys], label='y+0.1')
plt.plot(zs, [np.cos(ray.phase([0,0,z], 0)) for z in zs], label='z')
plt.xlabel("position (meters)")
plt.ylabel("amplitude")
plt.legend()
There are many times when you may want to work with more than one ray simultaneously. While this could be done using, say, a regular python list of Rays, it's more efficient to use a batoid RayVector. The advantage is that this allows more of the looping over Rays to occur in the c++ layer as opposed to the python layer.
RayVectors behave similarly to python lists, and can even be constructed from python lists of Rays.
In [ ]:
ray1 = batoid.Ray((0,0,0),(0,0,1),0,500e-9,1.1)
ray2 = batoid.Ray((1,0,0),(0,0,1),0,600e-9,1.2)
ray3 = batoid.Ray((0,1,0),(0,0,1),0,700e-9,1.3)
ray4 = batoid.Ray((0,0,1),(0,0,1),0,800e-9,1.4)
rayVector = batoid.RayVector([ray1, ray2, ray3, ray4])
assert rayVector[0] == ray1
assert rayVector[1] == ray2
assert rayVector[2] == ray3
assert rayVector[3] == ray4
RayVectors have the same attributes as Rays, but return numpy arrays with the first dimension corresponding to their different constituent Rays.
In [ ]:
print(rayVector.vx) # x component of each Ray's velocity
print(rayVector.r) # 2D array with shape (nray, 3) containing each Ray's position.
print(rayVector.r.shape) # First index ranges over the different Rays.
print(rayVector.omega) # Derived attributes work too.
RayVectors have the same methods as Rays. Note that they still only accept single position and time arguments though, which are then used for all of the constituent Rays.
In [ ]:
rayVector = batoid.RayVector([ray1, ray2, ray3, ray4])
print(rayVector.positionAtTime(1))
rayVector2 = rayVector.copy().propagate(1)
assert rayVector != rayVector2
rayVector.propagate(1)
assert rayVector == rayVector2
In [ ]:
print(rayVector.phase([0,0,0], 0))
print(rayVector.amplitude([0,0,0], 1))
RayVectors also have a few additional methods not present for Ray.
The sumAmplitude method is similar to the amplitude method, but adds the complex amplitudes of each Ray's plane wave together. This could be done in python, of course, but is significantly faster in some cases to do on the c++ side.
The trimVignetted and trimVignettedInPlace methods work to remove Rays that have .vignetted==True. We'll examine these methods in a future tutorial.
In [ ]:
assert sum(rayVector.amplitude([0,0,0],0)) == rayVector.sumAmplitude([0,0,0], 0)
It can be slow to construct large grids of Rays in python just to then assemble them into RayVectors. To help speed up this construction, batoid includes a few RayVector factory functions that can be used to rapidly create commonly desired grids of rays. These include RayVector.asGrid, RayVector.asPolar, and RayVector.asSpokes. In each case, the goal is to sample either a plane wave or spherical wave incident upon some optical system by
creating a grid of either parallel or co-radial co-phased Rays. For RayVector.asGrid, the Rays are assembled onto a parallelogram (commonly just a square), for RayVector.asPolar the pattern is more circularly symmetric. The RayVector.asSpokes function is a bit special, in that it does not uniformly sample the region of interest. It does have a use, however, for integrating functions on a circle or annulus using Gaussian quadrature. There are a number of different options for each of these factory functions to support different use cases.
In [ ]:
# An example of a square grid of rays:
rv1 = batoid.RayVector.asGrid(
backDist=10, # roughly how far back from the origin the rays are created
nx=15, # number of rays on a side
lx=10.0, # length of one side of the grid in meters
theta_x=np.deg2rad(0.0), # Field angle of rays in x-direction
theta_y=np.deg2rad(0.0), # Field angle of rays in y-direction
wavelength=500e-9
)
# Polar RayVector with similar extent
rv2 = batoid.RayVector.asPolar(
backDist=10,
nrad=20, # number of radii
naz=50, # number of azimuths on outermost ring.
outer=5.0, # outer radius. Default inner radius is 0.0
theta_x=np.deg2rad(0.0),
theta_y=np.deg2rad(0.0),
wavelength=500e-9
)
# Spokes RayVector with similar extent
rv3 = batoid.RayVector.asSpokes(
backDist=10,
rings=10, # number of rings
spokes=10, # number of spokes
outer=5.0, # outer radius
inner=2.5, # inner radius
theta_x=np.deg2rad(0.0),
theta_y=np.deg2rad(0.0),
wavelength=500e-9
)
plt.scatter(rv1.x, rv1.y, c='b', s=25)
plt.scatter(rv2.x, rv2.y, c='r', s=25)
plt.scatter(rv3.x, rv3.y, c='c', s=25)
plt.xlabel("x")
plt.ylabel("y")
plt.xlim(-6, 6)
plt.ylim(-6, 6)
plt.gca().set_aspect("equal")
plt.show()
plt.scatter(rv1.x, rv1.z, c='b', s=25)
plt.scatter(rv2.x, rv2.z, c='r', s=10)
plt.scatter(rv3.x, rv3.z, c='c', s=5)
plt.xlabel("x")
plt.ylabel("z")
plt.xlim(-6, 6)
plt.ylim(-1, 11)
plt.gca().set_aspect("equal")
plt.show()
In [ ]:
# Repeat with a non-zero field angle.
# An example of a square grid of rays:
rv1 = batoid.RayVector.asGrid(
backDist=10, # roughly how far back from the origin the rays are created
nx=15, # number of rays on a side
lx=10.0, # length of one side of the grid in meters
theta_x=np.deg2rad(5), # Field angle of rays in x-direction
theta_y=np.deg2rad(10), # Field angle of rays in y-direction
wavelength=500e-9
)
# Polar RayVector with similar extent
rv2 = batoid.RayVector.asPolar(
backDist=10,
nrad=20, # number of radii
naz=50, # number of azimuths on outermost ring.
outer=5.0, # outer radius. Default inner radius is 0.0
theta_x=np.deg2rad(5),
theta_y=np.deg2rad(10),
wavelength=500e-9
)
# Spokes RayVector with similar extent
rv3 = batoid.RayVector.asSpokes(
backDist=10,
rings=10, # number of rings
spokes=10, # number of spokes
outer=5.0, # outer radius
inner=2.5, # inner radius
theta_x=np.deg2rad(5),
theta_y=np.deg2rad(10),
wavelength=500e-9
)
plt.scatter(rv1.x, rv1.y, c='b', s=25)
plt.scatter(rv2.x, rv2.y, c='r', s=25)
plt.scatter(rv3.x, rv3.y, c='c', s=25)
plt.xlabel("x")
plt.ylabel("y")
plt.xlim(-6, 6)
plt.ylim(-6, 6)
plt.gca().set_aspect("equal")
plt.show()
plt.scatter(rv1.x, rv1.z, c='b', s=25)
plt.scatter(rv2.x, rv2.z, c='r', s=10)
plt.scatter(rv3.x, rv3.z, c='c', s=5)
plt.xlabel("x")
plt.ylabel("z")
plt.xlim(-6, 6)
plt.ylim(-1, 11)
plt.gca().set_aspect("equal")
plt.show()
In [ ]:
# If you have ipyvolume installed, you can visualize the ray grids in 3D.
has_ipv = True
try:
import ipyvolume as ipv
except ImportError:
has_ipv = False
In [ ]:
if has_ipv:
ipv.figure()
ipv.quiver(rv1.x, rv1.y, rv1.z, rv1.vx, rv1.vy, rv1.vz, size=5, color='red')
ipv.quiver(rv2.x, rv2.y, rv2.z, rv2.vx, rv2.vy, rv2.vz, size=5, color='blue')
ipv.quiver(rv3.x, rv3.y, rv3.z, rv3.vx, rv3.vy, rv3.vz, size=5, color='cyan')
ipv.xyzlim(10)
ipv.view(0,90)
ipv.show()
In [ ]: