prevent division by zero

source/f90/seteps.f90

@@ -71,7 +71,11 @@ subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
 ! because all cases of nos(l) > 1 are treated in case 1 above
       do k = 1, nos(l)
       j = j + 1
-      x = wn / osc(1, j)
+      if (osc(1, j) > 1d-3) then  ! prevent division by zero
+        x = wn / osc(1, j)
+      else
+        x = wn * 1d3
+      endif
       deno = x * dcmplx(x, osc(3, j))
       if (osc(2, j) >= 0.0d0) then
         deno = 1.0d0 - deno