Ray casting in code (not in 3D viewer)

Hello community.

We actually implemented an algorithm to compute all hits performed by a vector (sun direction) to a satellite model converted in the OpenCascade format.

To do so, in fact, we produce a huge grid a rays (sized accordingly to the model bouding box). Then, for each CAD object, we retrieve all faces and use the GeomAPI_IntCS to compute intersection points between each ray and le face. We can then retrieve the normal of each hit point because some of our processes require the angle between the ray and hit points normal (for example, to compute disturbing torques of sun radiation).
With these information, we are able to calculate approximatively the percentage of shading on each surface.
See attachment to understand the process.

The problem is that we can have big models and we have to launch this treatments many time across the simulation (each 10 seconds during one complete day for example) and this can take very long time (4 to 8 hours depending on the model / simulation time).

Instead of analysing each face and each ray, I wonder if any class or module (we're open for buying such module) could help us to perform something like what is actually done but for a complete model and that could return a reference to every hit points and concerned faces with high performances (no need to do the treatment on each face of each shape).

Everything should be performed in code behind (not within a 3D viewer), but if there is any way to use the GPU, it would be really appreciated (we are capable of producing triangulation of the model).

Note: we also use this algorithm to compute hits between rays calculated from the Earth, creating a grid from the sphere surface and each ray direction going to the top of the bounding sphere (no screenshot actually...).

Thank you very much for your help.