From 640dfa63d95d53ad3c6a152af87ef294c82b54b4 Mon Sep 17 00:00:00 2001
From: Karl-Michael Schindler <karl-michael.schindler@web.de>
Date: Thu, 29 Aug 2024 16:43:22 +0200
Subject: [PATCH] add extended Drude

---
 source/f90/oscillatorType.f90 |  2 +-
 source/f90/seteps.f90         | 14 +++++++++-----
 2 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/source/f90/oscillatorType.f90 b/source/f90/oscillatorType.f90
index 48e802f..aa017ed 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 ea332ec..4b86b2c 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
-- 
GitLab