Update myEels20-seteps.f90

make the code represent the usual formula more directly.

tests/SetEpsTestCasesFromScratch/WFW/myEels20-seteps.f90

@@ -28,7 +28,7 @@ subroutine seteps(nLayer, nos, osc, epsinf, wn, name, eps, layers, mode)
   !f2py depend(nLayer) nos, epsinf, eps
   
   double complex :: nomi, deno, addeps, addlorentz
-  double precision :: wn2, b
+  double precision :: wn2, b, gamma
   integer j, k, l, m
   logical debugFirstRun
 
@@ -52,7 +52,8 @@ subroutine seteps(nLayer, nos, osc, epsinf, wn, name, eps, layers, mode)
         deno =deno * (osc(1,j)**2 - wn * dcmplx( wn, osc(3,j) ) )
       
       else if (osc(1,j) < 0.) then! Negative TO mode means: _Additive_ Lorentz oscillator with Q
-        addlorentz = addlorentz + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 - wn2, -wn*abs(osc(1,j))*osc(3,j))  ! Sign of imaginary part changed (WFW)
+        gamma = abs(osc(1,j)) * osc(3,j)                                                          ! osc(3,j) is lambda. gamma = omega_TO * lambda
+        addlorentz = addlorentz + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 - wn2, -wn * gamma)  ! Sign of imaginary part changed (WFW)
       
       else                      ! osc(1,j) = 0   -> it is a Drude term
         addeps = addeps - dcmplx(osc(1,j+m)**2, wn*(osc(3,j)-osc(3,j+m))) /dcmplx(wn2, wn*osc(3,j))