fix lorentz in myyEels20-seteps.f90

source/WfW seteps/myEels20-seteps.f90

@@ -10,17 +10,17 @@ subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps)
   implicit none
 
   integer, intent(in) :: nLayer
-  integer, dimension(nLayer),intent(in) :: nos
+  integer, intent(in) :: nos(:)
   integer, intent(in) :: osc_size
-  double precision, dimension(3,osc_size),intent(in) :: osc
+  double precision, intent(in) :: osc(:,:)
   !f2py depend(osc_size) osc
   double precision, dimension(nLayer),intent(in) :: epsinf
   double precision,  intent(in) :: wn
   double complex, dimension(nLayer), intent(out) :: eps
   !f2py depend(nLayer) nos, epsinf, eps
   
-  double complex :: nomi, deno, addeps
-  double precision :: wn2, b
+  double complex :: nomi, deno, addeps, addlorentz
+  double precision :: wn2, b, gamma
   integer j, k, l, m
   logical debugFirstRun
 
@@ -32,6 +32,7 @@ subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps)
     nomi = dcmplx(1.0d0, 0.0d0)
     deno = dcmplx(1.0d0, 0.0d0)
     addeps = dcmplx(0.0d0, 0.0d0)
+    addlorentz = dcmplx(0.0d0, 0.0d0)
     wn2 = wn**2
     do k = 1, m     ! loop over all TO modes
       j = j + 1
@@ -42,7 +43,8 @@ subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps)
         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
-        addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 - wn2, -1*wn*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))
@@ -50,7 +52,7 @@ subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps)
 
     enddo
     j = j+m     ! we have already looped over the LO modes, therefore increase the index
-    eps(l) = epsinf(l) * (nomi / deno + addeps) ! brackets changed by HHe 230915
+    eps(l) = epsinf(l) * (nomi / deno + addeps) + addlorentz ! brackets changed by HHe 230915
   enddo
   debugFirstRun = .false.
   return

tests/seteps/Makefile

@@ -25,9 +25,9 @@ build-f90/seteps_mod.mod: ../../source/f90/seteps.f90
 
 MYEELS20_SETEPS = ../../source/WFW seteps/myEels20-seteps.f90
 	
-build-WFW/seteps.o:       $(wildcard ../../source/WFW\\ seteps/myEels20-seteps.f90)
+build-WFW/seteps.o:       ../../source/WFW\ seteps/myEels20-seteps.f90
 	$(FC) $(FFLAGS) -c -Jbuild-WFW -o $@ '$(MYEELS20_SETEPS)'
-build-WFW/seteps_mod.mod: $(wildcard ../../source/WFW\\ seteps/myEels20-seteps.f90)
+build-WFW/seteps_mod.mod: ../../source/WFW\ seteps/myEels20-seteps.f90
 	$(FC) $(FFLAGS) -c -Jbuild-WFW -o build-WFW/seteps.o '$(MYEELS20_SETEPS)'
 
 build-f77/setepsdriver: setepsdriver-f77.f90 build-f77/seteps.o