chbergmann / OpticsWorkbench

GNU Lesser General Public License v3.0
69 stars 23 forks source link

Computing Cone Spherical Beam fails for small number of `ConeAngle` and target Number of Rays #44

Closed DAsh0244 closed 7 months ago

DAsh0244 commented 7 months ago

Hello, I have been experimenting with your workbench and noticed that, as you set the ConeAngle property for a spherical beam to a small enough number relative to the number of rays used to render the cone, you can run into divide by zero errors, where no rays are rendered. It is not very rigourous, but any cone angle less than a degree with any amount of rays less than having a few hundred seems to reliably throw this error.

a typical traceback looks as below, when asking to render the attached file cone_test.zip:

22:54:19  Traceback (most recent call last):
  File "C:\Users\...\AppData\Roaming\FreeCAD\Mod\OpticsWorkbench\.\Ray.py", line 59, in execute
    self.redrawRay(fp)
  File "C:\Users\...\AppData\Roaming\FreeCAD\Mod\OpticsWorkbench\.\Ray.py", line 101, in redrawRay
    d_angle1 = math.radians(coneAngle/2)/M_angle1 # calculate the distance between the circles of the latitude
<class 'ZeroDivisionError'>: float division by zero

22:54:19  Recompute failed!

After poking through the source files, I think I have discovered the issue here , where the patch area a is estimated as 2*math.radians(coneAngle)/N, whereas instead it should be 2*math.pi*(1-math.cos(math.radians(coneAngle/2)))/N (see wiki for derivation). I tried this test out and it seems to work fine with the limited testing around in the above attached file. The angle is easily adjusted by changing Sketch.Constraints.cone_angle which is propagated to the other sketch and the beam coneAngle parameter, while the BeamNumber params are controlled within the Beam instance itself. This wouldn't break rendering spherical sources either from this code, as the expression reduces to the expected 4*pi*r^2 when coneAngle==360.

This approximation has the issue of causing a zero value to show up for small enough cone angles and low ray counts as shown in the below link see here for a plot of the differences between the figures, which causes the error. There is a very small (few milliradians) region that ends up having this zero region, and this region gets smaller and smaller the more rays you render with, as the patch areas a get smaller, causing the characteristic target dimension d for the patch size to decrease even more.

I have not done any extensive testing aside from this one instance, but I might have time to throw together a PR for this tomorrow or next weekend if desired. I wanted to try extending the workbench to support elllpitic beams and beams with non-uniform divergence angles, but might not find time for that until a few weeks from now. That said, this is probably a small enough fix I can squeeze the time for it in somewhere sooner.

P.S. I think using N = int(fp.BeamNrColumns * fp.BeamNrRows) instead N = int(fp.BeamNrColumns) of as shown here may also be a better choice for the target number of rays to render?

Regards, -D

chbergmann commented 7 months ago

I fixed it. Thank you for reporting