diff --git a/source/f90/oscillatorType.f90 b/source/f90/oscillatorType.f90
index 48e802f4304a06da2c09724187fa900470f762ec..aa017edef72a322009e7fb787f498404bf3e0335 100644
--- a/source/f90/oscillatorType.f90
+++ b/source/f90/oscillatorType.f90
@@ -1,4 +1,4 @@
 module oscillator_mod
   implicit none
-  integer, parameter :: Lorentz = 1, Kurosawa = 2, Drude = 3
+  integer, parameter :: Lorentz = 1, Kurosawa = 2, Drude = 3, extendedDrude = 4
 end module oscillator_mod
diff --git a/source/f90/seteps.f90 b/source/f90/seteps.f90
index ea332ec625e185512fd984cfc50347251d5ea64d..4b86b2cf9153d677dfc6ca39c0ae46a29583bfec 100644
--- a/source/f90/seteps.f90
+++ b/source/f90/seteps.f90
@@ -27,7 +27,7 @@ subroutine seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
 
   integer :: i, j, k
   double complex :: one, nomi, deno, addLorentz, addDrude
-  double precision :: wn2, b
+  double precision :: wn2, b, gamma
 
   one = dcmplx(1.0, 0.0)
   wn2 = wn**2
@@ -39,14 +39,18 @@ subroutine seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
     addDrude = dcmplx(0.0, 0.0)
     do j = 1, nos(i)
       select case (oscType(i))
-        case (Lorentz)
-          addLorentz = addLorentz + osc(1,k)**2 * osc(2,k) / dcmplx(osc(1,k)**2 - wn2, -wn*osc(1,k)*osc(3,k))
+        case (Lorentz) ! osc3 is lambda and not gammma
+          gamma = osc(1,k) * osc(3,k)
+          addLorentz = addLorentz + osc(1,k)**2 * osc(2,k) / dcmplx(osc(1,k)**2 - wn2, -wn * gamma)
         case (Kurosawa)
           b = wn / osc(1,k)
           nomi = nomi *  osc(3,k)**2 * (1.0 - b * dcmplx( b, osc(4,k)/osc(3,k)) )
           deno = deno * (osc(1,k)**2 -       wn * dcmplx(wn, osc(2,k)         ) )
-        case (Drude)
-          addDrude = addDrude - osc(1,k)**2 / dcmplx(wn2, wn*osc(1,j)*osc(2,j))
+        case (Drude)            ! As with lorentz, osc2 is lambda and not gammma
+          gamma = osc(1,k) * osc(2,k)
+          addDrude = addDrude - osc(1,k)**2 / dcmplx(wn2, wn * gamma)
+        case (extendedDrude)    ! osc2 and osc3 are gamma_Plasma and Gamma_0
+          addDrude = addDrude - dcmplx(osc(1,k)**2, wn*(osc(2,j)-osc(3,j))) / dcmplx(wn2, wn * osc(2,j))
         case default
           write (*,'(A)') '*** Error in oscillator type! ***'
       end select