An option to limit the propagation in resonant cavities was added (#163)

pyOpTools failed to simulate resonand cavities. This was solved by adding 2 limits to the possible number of iterations in the propagation of a ray. One is to limit the maximum number o parents, a ray can have, another is the intensity threshold for a ray to be propagated. This 2 options are provided to the System.

solves #158

pyoptools/raytrace/ray/ray.pxd

@@ -14,9 +14,15 @@ cdef class Ray:
 
     cdef public list orig_surf  # path of the originating surface
     cdef public double intensity, wavelength, pop
-    cdef list __childs,
+    cdef list __childs
     cdef np.ndarray _dir
 
+    # amount of parents of this ray in the optical path. It counts the number of
+    # times  a ray have been propagated through surfaces. It is used to stop
+    # the propagation in resonant cavities.
+
+    cdef public int _parent_cnt
+
     cpdef Ray ch_coord_sys(self, np.ndarray no, np. ndarray ae)
     cpdef Ray ch_coord_sys_inv_f(self, np.ndarray no , np.ndarray ae, bool childs)
 
@@ -31,4 +37,4 @@ cdef class Ray:
     # ~ def optical_path_parent(self):
     # ~ def optical_path(self):
 cdef Ray Rayf(np.ndarray pos, np.ndarray dir, double intensity, double wavelength,
-              n, label, draw_color, parent, double pop, orig_surf, int order)
+              n, label, draw_color, parent, double pop, orig_surf, int order, int _parent_cnt)

pyoptools/raytrace/ray/ray.pyx

@@ -41,7 +41,7 @@ np.import_array()
 
 
 cdef Ray Rayf(np.ndarray pos, np.ndarray dir, double intensity, double wavelength,
-              n, label, draw_color, parent, double pop, orig_surf, int order):
+              n, label, draw_color, parent, double pop, orig_surf, int order, int parent_cnt):
     """Function to create and initialize Ray instances Fast from cython
     when changing the code, always check that the __init__, and the
     Rayf initialization do the same.
@@ -65,6 +65,7 @@ cdef Ray Rayf(np.ndarray pos, np.ndarray dir, double intensity, double wavelengt
     instance.orig_surf=orig_surf
     instance.__childs=[]
     instance.order=order
+    instance._parent_cnt=parent_cnt
     return instance
 
 
@@ -147,6 +148,7 @@ cdef class Ray:
             self.orig_surf=orig_surf
             self.__childs=[]
             self.order=order
+            self._parent_cnt = 0
 
         # HasPrivateTraits.__init__(self,**traits)
     def __reduce__(self):
@@ -235,7 +237,8 @@ cdef class Ray:
         #        wavelength=self.wavelength,n=self.n, label=self.label,
         #        orig_surf=self.orig_surf)
 
-        cdef Ray parent=Rayf(npos, ndir, self.intensity, self.wavelength, self.n, self.label, self.draw_color, None , 0, self.orig_surf, self.order)
+        cdef Ray parent=Rayf(npos, ndir, self.intensity, self.wavelength, self.n, self.label, self.draw_color, None , 0, self.orig_surf, self.order, self._parent_cnt)
+
 
         # Calculate the transform of the childs and link them
 
@@ -286,7 +289,7 @@ cdef class Ray:
 
         ndir=dot(tm, self.dir)
 
-        cdef Ray parent=Rayf(npos, ndir, self.intensity, self.wavelength, self.n, self.label, self.draw_color, None , 0, self.orig_surf, self.order)
+        cdef Ray parent=Rayf(npos, ndir, self.intensity, self.wavelength, self.n, self.label, self.draw_color, None , 0, self.orig_surf, self.order, self._parent_cnt)
 
         # Calculate the transform of the childs and link them
 
@@ -349,7 +352,7 @@ cdef class Ray:
         # npos=mvdot1(tm,t)
         # ndir=mvdot1(tm,self.dir)
 
-        return Rayf(npos, ndir, self.intensity, self.wavelength, self.n, self.label, self.draw_color, None , 0, self.orig_surf, self.order)
+        return Rayf(npos, ndir, self.intensity, self.wavelength, self.n, self.label, self.draw_color, None , 0, self.orig_surf, self.order, self._parent_cnt)
 
     def get_final_rays(self, inc_zeros=True):
         '''Find the final rays of the raytrace
@@ -410,6 +413,7 @@ cdef class Ray:
         assert(isinstance(cr, Ray)), "A ray child must be a ray"
         cr.parent=self
         cr.order=len(self.__childs)
+        # cr._parent_cnt = self._parent_cnt + 1
         self.__childs.append(cr)
 
     def optical_path_parent(self):

pyoptools/raytrace/surface/surface.pyx

@@ -440,7 +440,7 @@ cdef class Surface(Picklable):
             # return [Ray(pos=PI,dir=S2,intensity=ri.intensity,
             #           wavelength=ri.wavelength,n=nr,
             #           label=ri.label, orig_surf=self.id)]
-            return [Rayf(PI, S2, ri.intensity, ri.wavelength, nr, ri.label, ri.draw_color, None, 0, self.id, 0)]
+            return [Rayf(PI, S2, ri.intensity, ri.wavelength, nr, ri.label, ri.draw_color, None, 0, self.id, 0, ri._parent_cnt+1)]
         elif sometrue(npisnan(S2)):
             # Total internal refraction case
             gamma1 = -2.*ni*cos(I)
@@ -465,7 +465,7 @@ cdef class Surface(Picklable):
             #            intensity=ri.intensity,
             #            wavelength=ri.wavelength,n=ni,
             #            label=ri.label, orig_surf=self.id)]
-            return [Rayf(PI, S3, ri.intensity, ri.wavelength, ni, ri.label, ri.draw_color, None, 0, self.id, 0)]
+            return [Rayf(PI, S3, ri.intensity, ri.wavelength, ni, ri.label, ri.draw_color, None, 0, self.id, 0, ri._parent_cnt+1)]
 
         else:
             # BeamSplitter case
@@ -490,18 +490,18 @@ cdef class Surface(Picklable):
                     #    intensity=ri.intensity*(1.-self.reflectivity),
                     #    wavelength=ri.wavelength,n=nr, label=ri.label, orig_surf=self.id),
                     Rayf(PI, S2, ri.intensity*(1.-reflect), ri.wavelength,
-                         nr, ri.label, ri.draw_color, None, 0, self.id, 0),
+                         nr, ri.label, ri.draw_color, None, 0, self.id, 0,ri._parent_cnt+1),
                     # Ray(pos=PI,dir=S3,
                     #    intensity=ri.intensity*self.reflectivity,
                     #    wavelength=ri.wavelength,n=ni,label=ri.label, orig_surf=self.id)
                     Rayf(PI, S3, ri.intensity*reflect, ri.wavelength, ni,
-                         ri.label, ri.draw_color, None, 0, self.id, 0)
+                         ri.label, ri.draw_color, None, 0, self.id, 0, ri._parent_cnt+1)
                 ]
             else:
                 # return [Ray(pos=PI,dir=S3,
                 #            intensity=ri.intensity*self.reflectivity,
                 #            wavelength=ri.wavelength,n=ni,label=ri.label, orig_surf=self.id)]
-                return [Rayf(PI, S3, ri.intensity, ri.wavelength, ni, ri.label, ri.draw_color, None, 0, self.id, 0)]
+                return [Rayf(PI, S3, ri.intensity, ri.wavelength, ni, ri.label, ri.draw_color, None, 0, self.id, 0, ri._parent_cnt+1)]
 
     cpdef pw_propagate1(self, Ray ri, ni, nr, rsamples, isamples, knots):
         '''Method to calculate wavefront emerging from the surface

pyoptools/raytrace/system/system.pxd

@@ -7,6 +7,10 @@ cdef class System(Picklable):
     cdef public double n  # This could be changed to material
     cdef public list _np_rays
     cdef public list _p_rays
+    cdef public int _max_ray_parent_cnt
+    cdef public float intensity_threshold
+
+    cdef public int _exit_status_flag # 0: no error, 1: error
 
     cpdef distance(self, Ray ri)
     cpdef propagate_ray(self, Ray ri)