Update seteps.f90

1 to 1 correspondence of formulas to literature.

source/f90/seteps.f90

@@ -29,8 +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
-  double precision :: wTO, wLO, wP, Q, yTO, yLO, y, y0, yp
+  double precision :: wn2
+  double precision :: wTO, wLO, wP, Qk, yk, yTO, yLO, y, y0, yp
 
   one = dcmplx(1.0, 0.0)
   wn2 = wn**2
@@ -42,28 +42,27 @@ 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)
+        case (Lorentz) ! Lambin 1990, equ. 7
           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)
+          Qk  = osc(2,k)
+          yk  = osc(3,k)
+          addLorentz = addLorentz + Qk * wTO**2 / dcmplx(wTO**2 - wn2, -yk * wn)
+        case (Kurosawa) ! Gervais 2002, equ. 2.16
           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)
+          nomi = nomi * dcmplx(wLO**2 - wn2, yLO * wn)
+          deno = deno * dcmplx(wTO**2 - wn2, yTO * wn)
+        case (Drude)! Lambin 1990, equ. 8
           wp = osc(1,k)
           y  = osc(2,k)
-          addDrude = addDrude - wp**2 / dcmplx(wn2, wn * y)
-        case (extendedDrude)
+          addDrude = addDrude - wp**2 / dcmplx(wn2, y * wn)
+        case (extendedDrude) ! Gervais 2002, equ. 2.14
           wp = osc(1,k)
           yp = osc(2,k)
           y0 = osc(3,k)
-          addDrude = addDrude - dcmplx(wp**2, wn*(y0 - yp)) / dcmplx(wn2, wn * y0)
+          addDrude = addDrude - dcmplx(wp**2, (y0 - yp) * wn) / dcmplx(wn2, wn * y0)
         case default
           write (*,'(A)') '*** Error in oscillator type! ***'
       end select