.. _demo_energy_conservation:

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.