From 893ebb9b2e75ac1dccfdb6778caa668c3544b6e9 Mon Sep 17 00:00:00 2001
From: Karl-Michael Schindler <karl-michael.schindler@web.de>
Date: Wed, 16 Oct 2024 13:45:25 +0200
Subject: [PATCH] Update seteps.f90
---
source/f90/seteps.f90 | 39 ++++++++++++++++++++++++++-------------
1 file changed, 26 insertions(+), 13 deletions(-)
diff --git a/source/f90/seteps.f90 b/source/f90/seteps.f90
index 8b6265a..80fd4f1 100644
--- a/source/f90/seteps.f90
+++ b/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
--
GitLab