Checking Energy Conservation

Its important that ray-tracers conserve energy. This demonstration gives an example case where we prove that the total power collected on the boundary equals to total power radiated from a simple emitting source. The principle can be extended to test cases with more complicated emission sources and boundaries.

As usual, we start by importing everything we need.

from math import pi
from raysect.core import translate, rotate
from raysect.primitive import Sphere
from raysect.optical import World, ConstantSF
from raysect.optical.observer import Pixel, PowerPipeline0D
from raysect.optical.material.emitter import UnityVolumeEmitter, UniformSurfaceEmitter

The emitting source will be a sphere at the origin, surrounded by an observing box. The number of samples requested will be relatively large so we can ensure good accuracy in the calculations. Note that this calculation isn’t spectrally resolved so we only need a narrow spectral range with a single wavelength bin. Initialise the constants and place the sphere in the scene.

samples = 100000

sphere_radius = 0.5
cube_size = 2

min_wl = 400
max_wl = 401

# set-up scenegraph
world = World()
emitter = Sphere(radius=sphere_radius, parent=world)

Note that due to symmetry we don’t need to observe each face of the cube. We only need to observe 1 side and multiply by 6. You could also model it with 6 individual pixel faces.

power = PowerPipeline0D(accumulate=False)
observing_plane = Pixel([power], x_width=cube_size, y_width=cube_size,
                        min_wavelength=min_wl, max_wavelength=max_wl,
                        spectral_bins=1, pixel_samples=samples,
                        parent=world, transform=rotate(0, 0, 0)*translate(0, 0, -cube_size / 2))

First, let’s consider the case where the sphere is a uniform volume emitter with emission of 1W/m^3/str/nm). The theoretical total emitted power would be given by 1W x the volume of the sphere x 4Pi steradians x the wavelength range sampled. Lets implement the calculation and compare it to the measured value.

print("Starting observations with volume emitter...")
calculated_volume_emission = 16 / 3 * pi**2 * sphere_radius**3 * (max_wl - min_wl)

emitter.material = UnityVolumeEmitter()
observing_plane.observe()
measured_volume_emission = 6 * power.value.mean
measured_volume_error = 6 * power.value.error()

print('Expected volume emission => {} W'.format(calculated_volume_emission))
print('Measured volume emission => {} +/- {} W'.format(measured_volume_emission, measured_volume_error))

For our test case, expected volume emission => 6.57 W, measured volume emission => 6.623 +/- 0.057 W.