diff --git a/source/f90/seteps.f90 b/source/f90/seteps.f90
index 2900bcd7fbc30eec4da2090782a6061973290b87..743ff75f721b926faf9371e65dc40147b984df7f 100644
--- a/source/f90/seteps.f90
+++ b/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