Update seteps.f90

source/f90/seteps.f90

@@ -19,6 +19,8 @@ subroutine seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
 ! ******************************************************************
 
   use oscillatorType_mod
+  
+  implicit none
 
   integer, intent(in) :: neps, nos(:), oscType(:)
   double precision, intent(in) :: osc(:, :), epsinf(:), wn
@@ -27,7 +29,8 @@ 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, gamma
+  double precision :: wn2, b
+  double precision :: wTO, wLO, wP, Q, yTO, yLO, y, y0, yp
 
   one = dcmplx(1.0, 0.0)
   wn2 = wn**2
@@ -39,22 +42,32 @@ subroutine seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
     addDrude = dcmplx(0.0, 0.0)
     do j = 1, nos(i)
       select case (oscType(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 (Lorentz)
+          wTO = osc(1,k)
+          Q   = osc(2,k)
+          yTO = osc(3,k)
+          addLorentz = addLorentz + wTO**2 * Q / dcmplx(wTO**2 - wn2, -wn * yTO)
         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)            ! 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,k)-osc(3,k))) / dcmplx(wn2, wn * osc(2,k))
+          wTO = osc(1,k)
+          yTO = osc(2,k)
+          wLO = osc(3,k)
+          yLO = osc(4,k)
+          b = wn / wTO
+          nomi = nomi *  wLO**2 * (1.0 - b * dcmplx( b, yLO/wLO) )
+          deno = deno * (wTO**2 -       wn * dcmplx(wn, yTO      ) )
+        case (Drude)
+          wp = osc(1,k)
+          y  = osc(2,k)
+          addDrude = addDrude - wp**2 / dcmplx(wn2, wn * y)
+        case (extendedDrude)
+          wp = osc(1,k)
+          yp = osc(2,k)
+          y0 = osc(3,k)
+          addDrude = addDrude - dcmplx(wp**2, wn*(y0 - yp)) / dcmplx(wn2, wn * y0)
         case default
           write (*,'(A)') '*** Error in oscillator type! ***'
       end select
-      eps(i) = epsinf(i) * (nomi / deno) + addDrude + addLorentz
+      eps(i) = epsinf(i) * (nomi / deno + addDrude) + addLorentz
       k = k + 1
     enddo
   enddo