diff --git a/source/f90/CHANGELOG b/source/f90/CHANGELOG
index d3d00f0c15273c752c4d71154935adb3336577ab..58f984750a4335d7fd81ad3a70a6dafc1bc40289 100644
--- a/source/f90/CHANGELOG
+++ b/source/f90/CHANGELOG
@@ -1,5 +1,7 @@
History
+2025 07 09 use new seteps function in doeels, eels, eels-boson.
+
2025 04 17 Makefile also creates eels-boson
2014 07 22 separate file IO and create subroutine calls.
diff --git a/source/f90/Makefile b/source/f90/Makefile
index baab7f77d8a8603d6207867f2a1d7c924d4cfcd1..70f44d156dabdddcf54eb906a1e4a22a4dbdafb1 100644
--- a/source/f90/Makefile
+++ b/source/f90/Makefile
@@ -16,6 +16,9 @@ endif
# python3 interpreter
PYTHON3 = python3
+# Platform, i.e. operating system
+OS := $(shell echo `uname -s`)
+
# implicit rules for compiling fortran files
%.o: %.f90 ; $(FC) $(FFLAGS) -c $<
%_mod.mod: %.f90 ; $(FC) $(FFLAGS) -c $<
@@ -31,7 +34,7 @@ bosonsubs = doboson.o respon.o sicot.o sintr.o rcffi.o o.o
boson: boson.f90 change_working_dir.o $(bosonsubs) $(bosonmods)
$(FC) $(FFLAGS) -o bosonf90 boson.f90 change_working_dir.o $(bosonsubs)
-eelsmods = quanc8_mod.mod queels_mod.mod seteps_mod.mod param_mod.mod
+eelsmods = quanc8_mod.mod queels_mod.mod seteps_mod.mod param_mod.mod oscillatorType_mod.mod
doeels.o: doeels.f90 $(eelsmods)
queels_mod.mod: quanc8_mod.mod fint1_mod.mod fint2_mod.mod fint3_mod.mod
@@ -40,6 +43,7 @@ fint1_mod.mod: surlos_mod.mod param_mod.mod
fint2_mod.mod: surlos_mod.mod param_mod.mod
fint3_mod.mod: surlos_mod.mod param_mod.mod
fun_mod.mod: param_mod.mod
+seteps_mod.mod: oscillatorType_mod.mod
getoptsubs = getopt.o text.o system.o constants.o date_and_time.o kinds.o dummy_variables.o
@@ -66,7 +70,7 @@ sufr_kinds.mod:
sufr_dummy.mod:
$(FC) $(FFLAGS) -c -o dummy_variables.o ../getopt-libs/libsufr-0.7.7/src/dummy_variables.f90
-eelssubs = doeels.o usurlo.o quanc8.o fun.o queels.o fint1.o fint2.o fint3.o surlos.o seteps.o phint.o qrat.o param.o
+eelssubs = doeels.o usurlo.o quanc8.o fun.o queels.o fint1.o fint2.o fint3.o surlos.o seteps.o phint.o qrat.o param.o oscillatorType.o
eels: eels.f90 change_working_dir.o $(eelssubs) $(eelsmods)
$(FC) $(FFLAGS) -o eelsf90 eels.f90 change_working_dir.o $(eelssubs)
@@ -78,11 +82,27 @@ eels-boson: eels-boson.f90 get_commandline_options.o $(getoptsubs) $(getoptmods)
pylibs: doeels-py doboson-py
-doeels-py: doeels.f90 quanc8.f90 surlos.f90 fint1.f90 fint2.f90 fint3.f90 queels.f90 phint.f90 qrat.f90 usurlo.f90
- $(PYTHON3) -m numpy.f2py -c seteps.f90 doeels.f90 quanc8.f90 surlos.f90 fint1.f90 fint2.f90 fint3.f90 queels.f90 phint.f90 qrat.f90 usurlo.f90 -m doeels
+doeels-py: doeels.f90 quanc8.f90 surlos.f90 fint1.f90 fint2.f90 fint3.f90 queels.f90 phint.f90 qrat.f90 usurlo.f90 param.f90 oscillatorType.f90 fun.f90
+ifeq ($(OS),Darwin)
+ $(PYTHON3) -m numpy.f2py -c seteps.f90 doeels.f90 quanc8.f90 surlos.f90 fint1.f90 fint2.f90 fint3.f90 queels.f90 phint.f90 qrat.f90 usurlo.f90 param.f90 oscillatorType.f90 fun.f90 -m doeels
+ install_name_tool -change @rpath/libgfortran.5.dylib /opt/local/lib/libgcc/libgfortran.5.dylib doeels.cpython-312-darwin.so
+else
+ $(PYTHON3) -m numpy.f2py -c seteps.f90 doeels.f90 quanc8.f90 surlos.f90 fint1.f90 fint2.f90 fint3.f90 queels.f90 phint.f90 qrat.f90 usurlo.f90 param.f90 oscillatorType.f90 fun.f90 -m doeels
+endif
doboson-py: doboson.f90 sicot.f90 sintr.f90 rcffi.f90 respon.f90 o.f90
- $(PYTHON3) -m numpy.f2py -c o.90 sintr.f90 rcffi.f90 respon.f90 sicot.f90 doboson.f90 -m doboson
+ifeq ($(OS),Darwin)
+ $(PYTHON3) -m numpy.f2py -c o.f90 sintr.f90 rcffi.f90 respon.f90 sicot.f90 doboson.f90 -m doboson
+ install_name_tool -change @rpath/libgfortran.5.dylib /opt/local/lib/libgcc/libgfortran.5.dylib doboson.cpython-312-darwin.so
+else
+ $(PYTHON3) -m numpy.f2py -c o.f90 sintr.f90 rcffi.f90 respon.f90 sicot.f90 doboson.f90 -m doboson
+endif
+
+mytest:
+ echo $(OS)
+ifeq ( $(OS), Darwin )
+ echo $(OS)
+endif
clean:
rm -f *.o
diff --git a/source/f90/TODO b/source/f90/TODO
index 23c6832b5ffa5890ab0c760e1f000d62a43b4fa2..8385f973cd75836d4d95e81bc83a77584b0b21d7 100644
--- a/source/f90/TODO
+++ b/source/f90/TODO
@@ -1,3 +1 @@
-- usage output when called without parameter or with options -h --help
-- options for specifying input and output files.
- Makefile: make tests
\ No newline at end of file
diff --git a/source/f90/doeels.f90 b/source/f90/doeels.f90
index d0f03b8e85c444ab5c5960bbdf3a75767a265505..56f28dafc49654e1d85607d4b051ad17d8362910 100644
--- a/source/f90/doeels.f90
+++ b/source/f90/doeels.f90
@@ -1,8 +1,8 @@
module doeels_mod
contains
subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size, &
- layers, neps, nper, name, name_size, thick, epsinf, nos, osc, osc_size,&
- contrl, mode, wn_array, debug, f_array, wn_array_size)
+ layers, neps, nper, name, name_size, thick, epsinf, nos, oscType, osc, osc_size, &
+ contrl, wn_array, debug, f_array, wn_array_size)
! ******************************************************************
! * *
@@ -24,10 +24,10 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
! * thick: layer thickness *
! * epsinf: epsilon at infinite frequency *
! * nos: number of oscillators *
+! * oscType: oscillator type: Lorentz, Kurosawa, drude *
! * osc: oscillator parameters, wTO, Q, lambda *
! * osc_size: number of oscillators *
! * contrl: 'image' for image-charge screening *
-! * mode: 'kurosawa' for kurosawa model *
! * wn_array: frequencies *
! * debug: logical *
! * f_array: spectrum *
@@ -47,10 +47,10 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
double precision, intent(in) :: e0, theta, phia, phib, wmin, wmax, dw
character (len = 72) :: comment(comment_size)
character (len = 10) :: name(name_size)
- double precision, intent(in) :: thick(name_size), epsinf(name_size), osc(3, osc_size)
- character (len = 10) :: contrl, mode
- integer, intent(in) :: comment_size, name_size, osc_size, wn_array_size
- integer, intent(in out) :: layers, nper, nos(name_size)
+ double precision, intent(in) :: thick(name_size), epsinf(name_size), osc(4, osc_size)
+ character (len = 10) :: contrl
+ integer, intent(in) :: comment_size, name_size, osc_size, wn_array_size, oscType(:)
+ integer, intent(in) :: layers, nper, nos(name_size)
double precision, intent(in out) :: wn_array(wn_array_size), f_array(wn_array_size)
logical, intent(in), optional :: debug
@@ -113,10 +113,15 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
nout = 1 + nw / 20
ener = 8065 * e0
psia = phia / 180 * pi
+! *** kms
psii = theta / 180 * pi
cospsi = dcos(psii)
sinpsi = dsin(psii)
tanpsi = dtan(psii)
+! *** f2023 replacement suggestion
+! cospsi = dcosd(phia)
+! sinpsi = dsind(phia)
+! tanpsi = dtand(phia)
prefac = dsqrt(255500 / e0)/(137 * cospsi)
facru = psia / cospsi * dsqrt(0.2624664d0 * e0)
elleps = (1.0d0 - phia / phib) * (1.0d0 + phia / phib)
@@ -192,7 +197,7 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
! if (debug) write (*,*) 'wn: ', wn
if (wn >= 0.0d0) then
if (wn /= 0.0d0) then
- if (.not. user) call seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
+ if (.not. user) call seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
x = wn / (2 * ener * psia)
if (rational) then
diff --git a/source/f90/eels-boson.f90 b/source/f90/eels-boson.f90
index 1bff0549ab08eea4df19c418ab69106a18c9d0bf..778aa7db9cff9692582d2f4cf2b8098e6cb6243e 100644
--- a/source/f90/eels-boson.f90
+++ b/source/f90/eels-boson.f90
@@ -10,6 +10,7 @@ program eels_boson
use get_commandline_options_mod
use doeels_mod
+ use oscillator_mod
implicit none
@@ -19,16 +20,16 @@ program eels_boson
integer :: i, j, jos, k, l, layers, neps, nper, nw
logical :: user
character (len = 72) :: comment(2)
- character (len = 10) :: contrl, mode
+ character (len = 10) :: contrl
double precision, allocatable :: xout(:), yout(:)
integer :: nout
double precision, allocatable :: thick(:), epsinf(:), osc(:, :)
- integer, allocatable :: nos(:)
+ integer, allocatable :: nos(:), oscType(:), tmp_oscType(:)
character (len = 10), allocatable :: name(:)
double precision, allocatable :: wn_array(:), f(:)
- integer :: old_size_1, old_size_2
+ integer :: old_size, old_size_1, old_size_2
double precision, allocatable :: tmp_osc(:, :)
double precision :: asym, emax, emin, gauss, t, width
@@ -79,14 +80,14 @@ program eels_boson
! *** read target specifications
- read(11, *) layers, nper, mode
+ read(11, *) layers, nper
user = layers == 0
if (.not. user) then
neps = layers
if (nper == -1) then
neps = layers + 1
nper = 1
- write(*,*) 'the substrate is a anisotropic uniaxial material'
+ write(*,*) 'The substrate is an anisotropic uniaxial material'
endif
if (layers < 0 .or. nper < 1 .or. nper > layers) then
stop '*** invalid target specifications ***'
@@ -106,7 +107,7 @@ program eels_boson
else
do j = 1, nos(l)
if (.not. allocated(osc)) then
- allocate(osc(3, nos(l)))
+ allocate(osc(4, nos(l)))
else if (j == 1) then
! call extend3(osc, nos(j))
old_size_1 = size(osc, 1)
@@ -116,20 +117,35 @@ program eels_boson
deallocate(osc)
call move_alloc(tmp_osc, osc)
endif
+ if (.not. allocated(oscType)) then
+ allocate(oscType(nos(l)))
+ else if (j == 1) then
+ old_size = size(oscType)
+ allocate(tmp_oscType(old_size + nos(l)))
+ tmp_oscType(1:old_size) = oscType
+ deallocate(oscType)
+ call move_alloc(tmp_oscType, oscType)
+ endif
jos = jos + 1
read(11, *) (osc(k, jos), k = 1, 3)
+ oscType(jos) = Lorentz
+ if (osc(2, jos) == -1.0) then ! Drude term
+ osc(2, jos) = osc(3, jos)
+ osc(3, jos) = 0.0
+ oscType(jos) = Drude
+ endif
if ((j == nos(l) / 2 + 1) .and. (nos(l) > 1)) then
write(6, 106)
write(6, 107)
endif
if (j == 1) then
if (l <= layers) then
- write(6, 104) l, name(l), thick(l), epsinf(l), (osc(i, jos), i = 1, 3)
+ write(6, 104) l, name(l), thick(l), epsinf(l), (osc(i, jos), i = 1, 4)
else
- write(6, 105) epsinf(l), (osc(i, jos), i = 1, 3)
+ write(6, 105) epsinf(l), (osc(i, jos), i = 1, 4)
endif
else
- write(6, 108) (osc(i, jos), i = 1, 3)
+ write(6, 108) (osc(i, jos), i = 1, 4)
endif
enddo
endif
@@ -148,8 +164,8 @@ program eels_boson
debug = .false.
call doeels(e0, theta, phia, phib, wmin, wmax, dw, comment, size(comment), &
- layers, neps, nper, name, size(name), thick, epsinf, nos, osc, size(osc, 2),&
- contrl, mode, wn_array, debug, f, size(wn_array))
+ layers, neps, nper, name, size(name), thick, epsinf, nos, oscType, osc, size(osc, 2),&
+ contrl, wn_array, debug, f, size(wn_array))
open(unit = 13, file = bosin_name)
! target temperature (Kelvin)
@@ -202,9 +218,9 @@ program eels_boson
'thickness', 5x, 'epsinf', 4x, 'wto , wp', 5x, 'q', 7x, 'gam/wto')
102 format(a10, d15.5)
103 format(1x, 72('-'))
-104 format(1x, i3, 2x, a10, g15.3, f10.4, f12.4, f10.4, f9.4)
-105 format(31x, f10.4, f12.4, f10.4, f9.4)
+104 format(1x, i3, 2x, a10, g15.3, f10.4, f12.4, f10.4, f9.4, f9.4)
+105 format(31x, f10.4, f12.4, f10.4, f9.4, f9.4)
106 format(45x, 'wlo , wp', 5x, 'q', 7x, 'gam/wlo')
107 format(45x, 28('-'))
-108 format(41x, f12.4, f10.4, f9.4)
+108 format(41x, f12.4, f10.4, f9.4, f9.4)
end program eels_boson
diff --git a/source/f90/eels.f90 b/source/f90/eels.f90
index 3eb1277e01072b48e1ee02fa199b009c24b473ad..99345eff336ef3809a613997ecc5513694456d0a 100644
--- a/source/f90/eels.f90
+++ b/source/f90/eels.f90
@@ -9,6 +9,7 @@ program eels
! ******************************************************************
use doeels_mod
+ use oscillator_mod
implicit none
@@ -16,14 +17,14 @@ program eels
integer :: i, j, jos, k, l, layers, neps, nper, nw
logical :: user
character (len = 72) :: comment(2)
- character (len = 10) :: contrl, mode
+ character (len = 10) :: contrl
double precision, allocatable :: thick(:), epsinf(:), osc(:, :)
- integer, allocatable :: nos(:)
+ integer, allocatable :: nos(:), oscType(:), tmp_oscType(:)
character (len = 10), allocatable :: name(:)
double precision, allocatable :: wn_array(:), f(:)
- integer :: old_size_1, old_size_2
+ integer :: old_size, old_size_1, old_size_2
double precision, allocatable :: tmp_osc(:, :)
logical :: debug
integer :: ioStatus
@@ -65,14 +66,14 @@ program eels
! *** read target specifications
- read(11, *) layers, nper, mode
+ read(11, *) layers, nper
user = layers == 0
if (.not. user) then
neps = layers
if (nper == -1) then
neps = layers + 1
nper = 1
- write(*,*) 'the substrate is a anisotropic uniaxial material'
+ write(*,*) 'The substrate is an anisotropic uniaxial material'
endif
if (layers < 0 .or. nper < 1 .or. nper > layers) then
stop '*** invalid target specifications ***'
@@ -92,7 +93,7 @@ program eels
else
do j = 1, nos(l)
if (.not. allocated(osc)) then
- allocate(osc(3, nos(l)))
+ allocate(osc(4, nos(l)))
else if (j == 1) then
! call extend3(osc, nos(j))
old_size_1 = size(osc, 1)
@@ -102,20 +103,35 @@ program eels
deallocate(osc)
call move_alloc(tmp_osc, osc)
endif
+ if (.not. allocated(oscType)) then
+ allocate(oscType(nos(l)))
+ else if (j == 1) then
+ old_size = size(oscType)
+ allocate(tmp_oscType(old_size + nos(l)))
+ tmp_oscType(1:old_size) = oscType
+ deallocate(oscType)
+ call move_alloc(tmp_oscType, oscType)
+ endif
jos = jos + 1
read(11, *) (osc(k, jos), k = 1, 3)
+ oscType(jos) = Lorentz
+ if (osc(2, jos) == -1.0) then ! Drude term
+ osc(2, jos) = osc(3, jos)
+ osc(3, jos) = 0.0
+ oscType(jos) = Drude
+ endif
if ((j == nos(l) / 2 + 1) .and. (nos(l) > 1)) then
write(6, 106)
write(6, 107)
endif
if (j == 1) then
if (l <= layers) then
- write(6, 104) l, name(l), thick(l), epsinf(l), (osc(i, jos), i = 1, 3)
+ write(6, 104) l, name(l), thick(l), epsinf(l), (osc(i, jos), i = 1, 4)
else
- write(6, 105) epsinf(l), (osc(i, jos), i = 1, 3)
+ write(6, 105) epsinf(l), (osc(i, jos), i = 1, 4)
endif
else
- write(6, 108) (osc(i, jos), i = 1, 3)
+ write(6, 108) (osc(i, jos), i = 1, 4)
endif
enddo
endif
@@ -134,8 +150,8 @@ program eels
debug = .false.
call doeels(e0, theta, phia, phib, wmin, wmax, dw, comment, size(comment), &
- layers, neps, nper, name, size(name), thick, epsinf, nos, osc, size(osc, 2),&
- contrl, mode, wn_array, debug, f, size(wn_array))
+ layers, neps, nper, name, size(name), thick, epsinf, nos, oscType, osc, size(osc, 2),&
+ contrl, wn_array, debug, f, size(wn_array))
open (unit = 12, file = 'eelsou')
write (12, 207) e0, theta, phia, phib
@@ -151,10 +167,10 @@ program eels
'thickness', 5x, 'epsinf', 4x, 'wto , wp', 5x, 'q', 7x, 'gam/wto')
102 format(a10, d15.5)
103 format(1x, 72('-'))
-104 format(1x, i3, 2x, a10, g15.3, f10.4, f12.4, f10.4, f9.4)
-105 format(31x, f10.4, f12.4, f10.4, f9.4)
+104 format(1x, i3, 2x, a10, g15.3, f10.4, f12.4, f10.4, f9.4, f9.4)
+105 format(31x, f10.4, f12.4, f10.4, f9.4, f9.4)
106 format(45x, 'wlo , wp', 5x, 'q', 7x, 'gam/wlo')
107 format(45x, 28('-'))
-108 format(41x, f12.4, f10.4, f9.4)
+108 format(41x, f12.4, f10.4, f9.4, f9.4)
207 format('e0 =', f6.2, ' theta =', f5.1, ' phia =', f5.2, ' phib =', f5.2)
end program eels
diff --git a/source/f90/oscillatorType.f90 b/source/f90/oscillatorType.f90
new file mode 100644
index 0000000000000000000000000000000000000000..48e802f4304a06da2c09724187fa900470f762ec
--- /dev/null
+++ b/source/f90/oscillatorType.f90
@@ -0,0 +1,4 @@
+module oscillator_mod
+ implicit none
+ integer, parameter :: Lorentz = 1, Kurosawa = 2, Drude = 3
+end module oscillator_mod
diff --git a/source/f90/seteps.f90 b/source/f90/seteps.f90
index 6ae2fc11c9762187db5b512e97d3b6f70b772de4..ea332ec625e185512fd984cfc50347251d5ea64d 100644
--- a/source/f90/seteps.f90
+++ b/source/f90/seteps.f90
@@ -1,6 +1,6 @@
module seteps_mod
contains
-subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
+subroutine seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
! ******************************************************************
! * *
@@ -9,83 +9,52 @@ subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
! * *
! * neps: number of epsilon values *
! * nos: number of oscillators *
+! * oscType: oscillator type: Lorentz, Kurosawa, drude *
! * osc: oscillator parameters, wTO, Q, lambda *
! * epsinf: epsilon at infinite frequency *
-! * wn: frequency *
-! * name: layer names *
+! * wn: frequency *
! * eps: epsilon *
! * layers: number of layers *
-! * mode: 'kurosawa' for kurosawa model *
! * *
! ******************************************************************
- implicit none
- integer, intent(in) :: neps
- integer, intent(in) :: nos(:)
+ use oscillator_mod
+
+ integer, intent(in) :: neps, nos(:), oscType(:)
double precision, intent(in) :: osc(:, :), epsinf(:), wn
- character (len=10), intent(in) :: name(:)
double complex, intent(in out) :: eps(:)
integer, intent(in) :: layers
- character (len=10), intent(in) :: mode
-
- double precision :: x
- double complex :: deno, nomi
- integer :: j, k, l, m
-! write (*,*) 'seteps:'
-! write (*,*) 'nos: ', size(nos)
-! write (*,*) 'osc: ', size(osc)
-! write (*,*) 'epsinf: ', size(epsinf)
-! write (*,*) 'name: ', size(name)
-! write (*,*) 'eps: ', size(eps)
-! write (*,*) 'thick: ', size(thick)
+ integer :: i, j, k
+ double complex :: one, nomi, deno, addLorentz, addDrude
+ double precision :: wn2, b
- j = 0
- do l = 1, neps
- m = nos(l)/2
- nomi = dcmplx(1.0d0, 0.0d0)
- deno = dcmplx(1.0d0, 0.0d0)
- if (mode == 'kurosawa') then
- do k = 1, m
- j = j + 1
-! since osc(1,*) and wn are real, the following should be equivalent
-! Check required.
-! Furthermore the first term can be rewritten as
-! (osc(1, j + m) - wn) * (osc(1, j + m) + wn)
-! nomi = nomi * cmplx(osc(1, j + m)**2 - wn**2, - wn * osc(3, j + m))
- nomi = nomi * (osc(1, j + m)**2 - wn**2 - dcmplx(0.0d0, wn * osc(3, j + m)))
- deno = deno * (osc(1, j )**2 - wn**2 - dcmplx(0.0d0, wn * osc(3, j )))
- enddo
- eps(l) = epsinf(l) * nomi / deno
- elseif (name(l) == 'metal') then
- j = j + 1
-! suggestion for replacement
-! eps(l) = -osc(1, j)**2 / cmplx(wn**2, wn * osc(3, j))
- eps(l) = -osc(1, j)**2 / ( wn**2 + dcmplx(0.0d0, wn * osc(3, j)) )
- else
- eps(l) = epsinf(l)
-! The original version had this additional loop. It seems, it has been removed
-! because all cases of nos(l) > 1 are treated in case 1 above
- do k = 1, nos(l)
- j = j + 1
- if (osc(1, j) > 1d-3) then ! prevent division by zero
- x = wn / osc(1, j)
- else
- x = wn * 1d3
- endif
- deno = x * dcmplx(x, osc(3, j))
- if (osc(2, j) >= 0.0d0) then
- deno = 1.0d0 - deno
- endif
- if (cdabs(deno) == 0.0d0) then ! replace 0 by machine epsilon
- ! if deno is always > 0 then this would do it:
- ! deno = cdmax(deno, epsilon(1.0d0) / 2)
- deno = epsilon(1.0d0) / 2
- endif
- eps(l) = eps(l) + osc(2, j) / deno
- enddo
- endif
+ one = dcmplx(1.0, 0.0)
+ wn2 = wn**2
+ k = 1
+ do i = 1, neps
+ nomi = one
+ deno = one
+ addLorentz = dcmplx(0.0, 0.0)
+ 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 (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 default
+ write (*,'(A)') '*** Error in oscillator type! ***'
+ end select
+ k = k + 1
+ eps(i) = epsinf(i) * (nomi / deno + addDrude) + addLorentz
+ enddo
enddo
+
if (neps == layers + 1) then
! the substrate is a anisotropic uniaxial material
eps(layers) = cdsqrt(eps(layers) * eps(layers + 1))