really fix of the Lorentz term.

source/WfW seteps/myEels20-seteps.f90

@@ -28,7 +28,12 @@ subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps)
    
   j = 0
   do l = 1, nLayer      ! loop over different thin film layers
+    ! The numbering does not work with a mix of Lorentz and Kurosawa/extended Drude terms
+    ! Either Lorentz with no Drude or Kurosawa and extended Drude terms
     m = nos(l)/2      ! m number of TO modes = offset to reach the LO mode in the joint TO-LO list
+    if (osc(1,1) < 0.0) then
+      m = nos(l)
+    endif
     nomi = dcmplx(1.0d0, 0.0d0)
     deno = dcmplx(1.0d0, 0.0d0)
     addeps = dcmplx(0.0d0, 0.0d0)
@@ -37,21 +42,21 @@ subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps)
     do k = 1, m     ! loop over all TO modes
       j = j + 1
 
-      if (osc(1,j) > 0.) then     ! positive TO mode: 'Kurosawa' form: _Multiplicative_ phonon mode
+      if (osc(1,j) > 0.0) then     ! positive TO mode: 'Kurosawa' form: _Multiplicative_ phonon mode
         b = wn/osc(1, j+m)
         nomi = nomi * osc(1, j+m)**2 * (1.0 - b * dcmplx(b, osc(3,j+m)/osc(1, j+m)) )
         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
+        j = j + m     ! we have already looped over the LO modes, therefore increase the index
+      else if (osc(1,j) < 0.0) then! Negative TO mode means: _Additive_ Lorentz oscillator with Q
         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
+      else                      ! osc(1,j) = 0   -> it is a double damping extended Drude term
         addeps = addeps - dcmplx(osc(1,j+m)**2, wn*(osc(3,j)-osc(3,j+m))) /dcmplx(wn2, wn*osc(3,j))
+        j = j + m     ! we have already looped over the LO modes, therefore increase the index
       end if
 
     enddo
-    j = j+m     ! we have already looped over the LO modes, therefore increase the index
     eps(l) = epsinf(l) * (nomi / deno + addeps) + addlorentz ! brackets changed by HHe 230915
   enddo
   debugFirstRun = .false.

tests/seteps/setepsdriver-WFW.f90

@@ -58,6 +58,8 @@ program setepsdriver
         endif
         jos = jos + 1
         read(11, *) (osc(k, jos), k = 1, 3)
+! make the first parameter negative as a sign this is a Lorentz oscillator
+        osc(1,jos) = -osc(1,jos)
       enddo
     endif
   enddo