Solid angle for multiple apertures - arbitrary points

[Didou09]

Framework:

• tofu needs to be able to compute efficiently the solid angle subtended:

  • by a set of nd closed 3d polygons representing nd detectors (planar)
  • observed through a series of na closed 3d polygons representing na successive apertures
  • as seen from any npts points in space

• The nd closed 3d polygons are:

  • representing nd detectors sharing the same apertures
  • 3d, defined in (X, Y, Z) coordinates
  • 3d, but simple (no self-crossing)
  • planar and oriented (provided with a normal vector to indicate which side is valid) (i.e.: photons can only be detected from one side)
  • typically of the order of nd = 10 - 1000

• The na closed 3d polygons are:

  • representing na apertures, shared by all detectors
  • 3d, defined in (X, Y, Z) coordinates
  • 3d, but simple (no self-crossing)
  • not necessarily planar
  • typically of the order of na = 1-10

• the npts points are:

  • arbitrary and user-provided
  • provided in (X, Y, Z) coordinates
  • provided as 3 (X, Y, Z) arrays of the same shape
  • the shape is also arbitrary and the routine should return the result as a single array of the same shape as the input
  • typically many (npts = 10^4-6)

• It will be necessary, for each observation point, to compute the solid angle subtended by the intersection between all apertures and the detector, for each detector.

  • It will require first to compute the intersection of all polygons as seen by the point
  • Then to compute the solid angle subtended by that unique intersection polygon

• Additionally, there should be an option to return not only the solid angle, but also a unit vector pointed from the observation point, towards the center of mass of the intersection polygon (to indicate, for each points, in which main direction are emitted the photons, useful for anisotropic radiation like hard X-rays)

• Additionally, a test should be run to check whether an observation point can actually see the intersection polygon by running the ray-tracing routine that can tell whether an obstacle lies on the path (i.e. a structural element, like a PFC).

  • the code for this already exists
  • this check should be activated by default, by it should be possible to de-activate it

• The branch already exists, the user-interface is almost done and can be found in: tofu/geom/_comp_solidangles.py It is the calc_solidangle_apertures() routine

  Don't hesitate to start working on it when you feel ready

Suggestions & code structure:

• I will create a new dedicated branch from devel called Issue_SolidAnglesArbitrary on which I will prepare the basics and you can take over then

• I will create the pure-python routine that will serve as entry-point for the user

  • this entry point will deal with the inputs pre-formatting and output post-formatting
  • it the middle it will call a cython routine that you will create

• I will prepare basic unit tests that you can modify / complement as you see fit

• The typical loop would look like this (to be discussed based on your expertise):

# input options: 
    visibility = True / False
    return_vector = True / False

# initialize solid angle array with zeros
# (optional, if return_vector) initialize unit vector array with nans

# loop 1: on npts (observation points)

    # loop 2: on n1 (detectors)

        # test 1: sides
           check if point lies on the good side of the detector
           => stop if not

        # loop 3: on na (apertures)

            # `computation 1`:
               compute intersection between detector and all apertures as seen from point
               stop as soon as null is found (no intersection)

        # `computation 2`:
            compute the solid angle subtended by the intersection (fill the array that was initialized if non-zero)

       # computation 3 (optional, necessary if visibility or return_vector):
           compute the direction of photon emission: 
           the unit vector from the observation point towards the center of mass of the intersection polygon

        #` test 2`: visibility (optional, requires computation 3)
             check ray-tracing to make sure there is no obstacle between the observation point and the projection, on the detector surface, of the direction of photon emission (unit vector)
             => if no visibility, set solid angle to 0 and unit vector to nan

Code re-use:

• test 2: a visibility test has already been coded by Laura for single-point, called from a cython routine only (inlined)

• computation 2: a routine already exists to compute the solid angle subtended by a closed 3d polygon as seen from any point

Possible difficulties / discussion points:

• computation 1: the intersection polygon

  • There are (at least) 2 ways of doing this:
    • The most general way : the polygons are projected on a sphere around the observation point, and we compute their intersection polygon in 2d but in spherical coordinates
    • The easiest but less generic way: the polygons are projected on a common plane located on the center of mass of the one which is closest to the observation point, and we compute their intersection polygon in 2d in local planar coordinates
  • The easiest way is just a simple 2d polygon intersection algorithms (there are C libraries doing this, and algorithms online). But it will fail if one of the aperture polygons is very large non-planar such that we can't find a plane to project it on
  • The most generic way (2d spherical coordinates) will work in all configurations but I'm guessing it is significantly more difficult to code
  • The easiest way is a proxy for the generic way: in reality we are looking for a direction of observation (does it exist ?), which is more faithfully represented by spherical coordinates
  • We should discuss this choice when you are ready and decide which way to go

• I am guessing the best is to parallelize the major loop (loop 1), because the computation is independent for each point, but you might have a better idea.

Priorities:

• It is probably best to make 3 versions of the code, with the following order of priority:

---------	----------------------------------------------------	----------	---------------------		Priority	name	visibility	return_unit_vector
1	compute_solid_angle_apertures_unitvectors()	False	True
2	compute_solid_angle_apertures_light()	False	False
3	compute_solid_angle_apertures_visibility()	True	False

The last combination (visibility = True and return_unit_vector=True) can be handled by re-using the first one and combining it with already-coded visibility checks.

Illustrations:

In the example below, we define three arbitrary 3d polygons, the square represents the detector, the other 2 are its associated apertures. The red lines indicate normal vectors.

And after computing the solid angle for all points on a fine sampling of the tokamak volume, we plot the area where the solid angle is non-zero as a grey area. We used visibility = True to detect the fact that the internal wall hides a fraction of the viewing volume of the detector.

This was done in 2014 using a very ill-coded version of the code that I think you should not use (the computation took ~2 days....).

[Didou09]

Hi @obrejank , now the Issue653_SoliAngleArbitrary branch is up-to-date, you can start on it when you feel like it

The code you can work on is in tofu/geom/_GG.pyx (bottom of the file, routine compute_solid_angle_apertures_unitvectors())

I have prepared basci unit tests that are running and you can use then as you make progress to debug. They can be launched in 2 ways: from the terminal or from a python console.

From the terminal (with verbose and debugger option)

$ pytest tofu/tests/tests01_geom/test_02_compute.py -v -x --pdb

From inside a python console:

In [1]: import tofu as tf

In [2]: test = tf.tests.tests01_geom.test_02_compute.Test01_SolidAngles()

In [3]: test.setup()

In [4]: test.test01_solid_angle_apertures()

Thanks for your help!