diff --git a/tests/seteps/Makefile b/tests/seteps/Makefile
index f6d12e24689fc8142505690d550076dee8fc721f..d461c7ed5e813428500dffe1e660c46843845d91 100644
--- a/tests/seteps/Makefile
+++ b/tests/seteps/Makefile
@@ -12,22 +12,32 @@ ifeq "$(GFORTAN_VERSION_GE_5)" "1"
FFLAGS += -fdiagnostics-color=auto
endif
-all: build-f77/setepsdriver build-f90/setepsdriver build-WFW/setepsdriver
+all: build-f77/setepsdriver build-f90/setepsdriver build-f90-new/setepsdriver-1 build-f90-new/setepsdriver-2 build-WFW/setepsdriver
# implicit rules for compiling fortran files
build-f77/seteps.o: ../../source/split-f77/seteps.f
$(FC) $(FFLAGS) -Wno-conversion -Wno-compare-reals -c -o $@ $<
-build-f90/seteps.o: ../../source/f90/seteps.f90
+build-f90/seteps.o: seteps-old.f90
$(FC) $(FFLAGS) -c -Jbuild-f90 -o $@ $<
-build-f90/seteps_mod.mod: ../../source/f90/seteps.f90
+build-f90/seteps_mod.mod: seteps-old.f90
$(FC) $(FFLAGS) -c -Jbuild-f90 -o build-f90/seteps.o $<
+build-f90-new/seteps.o: ../../source/f90/seteps.f90
+ $(FC) $(FFLAGS) -c -Jbuild-f90-new -o $@ $<
+build-f90-new/seteps_mod.mod: ../../source/f90/seteps.f90
+ $(FC) $(FFLAGS) -c -Jbuild-f90-new -o build-f90/seteps.o $<
+
+build-f90-new/oscillatorType.o: ../../source/f90/oscillatorType.f90
+ $(FC) $(FFLAGS) -c -Jbuild-f90-new -o $@ $<
+build-f90-new/oscillatorType_mod.mod: ../../source/f90/oscillatorType.f90
+ $(FC) $(FFLAGS) -c -Jbuild-f90-new -o build-f90/oscillatorType.o $<
+
MYEELS20_SETEPS = ../../source/WFW seteps/myEels20-seteps.f90
-build-WFW/seteps.o: ../../source/WFW\ seteps/myEels20-seteps.f90
+build-WFW/seteps.o: $(wildcard ../../source/WFW\\ seteps/myEels20-seteps.f90)
$(FC) $(FFLAGS) -c -Jbuild-WFW -o $@ '$(MYEELS20_SETEPS)'
-build-WFW/seteps_mod.mod: ../../source/WFW\ seteps/myEels20-seteps.f90
+build-WFW/seteps_mod.mod: $(wildcard ../../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
@@ -36,11 +46,17 @@ build-f77/setepsdriver: setepsdriver-f77.f90 build-f77/seteps.o
build-f90/setepsdriver: setepsdriver-f90.f90 build-f90/seteps.o build-f90/seteps_mod.mod
$(FC) $(FFLAGS) -Jbuild-f90 -o $@ $< build-f90/seteps.o
+build-f90-new/setepsdriver-1: setepsdriver-f90-new-1.f90 build-f90-new/seteps.o build-f90-new/seteps_mod.mod build-f90-new/oscillatorType.o build-f90-new/oscillatorType_mod.mod
+ $(FC) $(FFLAGS) -Jbuild-f90-new -o $@ $< build-f90-new/seteps.o build-f90-new/oscillatorType.o
+
+build-f90-new/setepsdriver-2: setepsdriver-f90-new-2.f90 build-f90-new/seteps.o build-f90-new/seteps_mod.mod build-f90-new/oscillatorType.o build-f90-new/oscillatorType_mod.mod
+ $(FC) $(FFLAGS) -Jbuild-f90-new -o $@ $< build-f90-new/seteps.o build-f90-new/oscillatorType.o
+
build-WFW/setepsdriver: setepsdriver-WFW.f90 build-WFW/seteps.o build-WFW/seteps_mod.mod
$(FC) $(FFLAGS) -Jbuild-WFW -o $@ $< build-WFW/seteps.o
-test: build-f90/setepsdriver build-f77/setepsdriver build-WFW/setepsdriver
- for NUMBER in 000 001 002 003 004 005 006 007 008 009 010 011 012 013 014 015 016 017 018 019 020 021 022 023 024 025 ; do \
+test: build-f90/setepsdriver build-f90-new/setepsdriver-1 build-f90-new/setepsdriver-2 build-f77/setepsdriver build-WFW/setepsdriver
+ for NUMBER in 001 002 003 004 005 006 007 008 009 010 011 012 013 014 015 016 017 018 019 020 021 022 023 024 025 ; do \
echo " eels input file number: "$$NUMBER ; \
cp ../../tests/inputFiles/eelsin$$NUMBER build-f77/eelsin ; \
build-f77/setepsdriver ; \
@@ -50,22 +66,35 @@ test: build-f90/setepsdriver build-f77/setepsdriver build-WFW/setepsdriver
build-f90/setepsdriver ; \
mv build-f90/epsout results-f90/epsout$$NUMBER ; \
rm build-f90/eelsin ; \
+ cp ../../tests/inputFiles/eelsin$$NUMBER build-f90-new/eelsin ; \
+ build-f90-new/setepsdriver-1 ; \
+ mv build-f90-new/epsout results-f90/epsout$$NUMBER ; \
+ rm build-f90-new/eelsin ; \
cp ../../tests/inputFiles/eelsin$$NUMBER build-WFW/eelsin ; \
build-WFW/setepsdriver ; \
mv build-WFW/epsout results-WFW/epsout$$NUMBER ; \
rm build-WFW/eelsin ; \
done
+ for NUMBER in 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 ; do \
+ echo " eels input file number: "$$NUMBER ; \
+ cp ../../tests/seteps/setepsIn/setepsIn$$NUMBER.txt build-f90-new/setepsIn.txt ; \
+ build-f90-new/setepsdriver-2 ; \
+ mv build-f90-new/epsout results-f90-new/epsout$$NUMBER ; \
+ rm build-f90-new/setepsIn.txt ; \
+ done
clean: cleanbuild cleanresults
cleanbuild:
rm -rf build-f77/*
rm -rf build-f90/*
+ rm -rf build-f90-new/*
rm -rf build-WFW/*
cleanresults:
rm -rf results-f77/*
rm -rf results-f90/*
+ rm -rf results-f90-new/*
rm -rf results-WFW/*
.PHONY: all test clean
diff --git a/tests/seteps/seteps-old.f90 b/tests/seteps/seteps-old.f90
new file mode 100644
index 0000000000000000000000000000000000000000..6ae2fc11c9762187db5b512e97d3b6f70b772de4
--- /dev/null
+++ b/tests/seteps/seteps-old.f90
@@ -0,0 +1,98 @@
+module seteps_mod
+contains
+subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
+
+! ******************************************************************
+! * *
+! * set up long-wavelength dielectric functions of the layers for *
+! * the present frequency wn (in cm**-1) *
+! * *
+! * neps: number of epsilon values *
+! * nos: number of oscillators *
+! * osc: oscillator parameters, wTO, Q, lambda *
+! * epsinf: epsilon at infinite frequency *
+! * wn: frequency *
+! * name: layer names *
+! * eps: epsilon *
+! * layers: number of layers *
+! * mode: 'kurosawa' for kurosawa model *
+! * *
+! ******************************************************************
+
+ implicit none
+ integer, intent(in) :: neps
+ integer, intent(in) :: nos(:)
+ 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)
+
+ 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
+ enddo
+ if (neps == layers + 1) then
+! the substrate is a anisotropic uniaxial material
+ eps(layers) = cdsqrt(eps(layers) * eps(layers + 1))
+ if (dimag(eps(layers)) < 0.0d0) then
+ eps(layers) = -eps(layers)
+ endif
+ endif
+ return
+end subroutine seteps
+end module seteps_mod
diff --git a/tests/seteps/setepsdriver-f90-new-1.f90 b/tests/seteps/setepsdriver-f90-new-1.f90
new file mode 100644
index 0000000000000000000000000000000000000000..add95cc82ad84f50ba542fcb1ca34ca1580bb991
--- /dev/null
+++ b/tests/seteps/setepsdriver-f90-new-1.f90
@@ -0,0 +1,93 @@
+program setepsdriver
+
+ use seteps_mod
+ use oscillator_mod
+
+ implicit none
+
+ double precision :: wnmin, wnmax, dwn, wn
+
+ double precision, allocatable :: osc(:,:), epsinf(:)
+ double complex, allocatable :: eps(:)
+ integer, allocatable :: nos(:), oscType(:)
+ integer :: jos, layers, neps, nper
+ character (len = 10), allocatable :: name(:)
+ character (len = 10) mode
+ integer :: old_size
+ double precision, allocatable :: tmp_osc(:,:)
+ integer, allocatable :: tmp_oscType(:)
+
+ integer :: i, j, k, l, noPoints
+
+ write(*,*) 'setepsdriver-f90-new-1 start'
+
+ open(unit = 11, file = 'build-f90-new/eelsin')
+ read(11, *)
+ read(11, *)
+ read(11, *)
+ read(11, *)
+ read(11, *) wnmin
+ read(11, *) wnmax
+ read(11, *) dwn
+ read(11, *)
+ read(11, *)
+ read(11, *) layers, nper, mode
+ neps = layers
+ if (nper == -1) then
+ neps = layers + 1
+ nper = 1
+ endif
+ jos = 0
+ allocate (name(neps), epsinf(neps), nos(neps))
+ do l = 1, neps
+ if (l <= layers) then
+ read(11, *) name(l)
+ endif
+ read(11, *) epsinf(l), nos(l)
+ if (nos(l) > 0) then
+ do j = 1, nos(l)
+ 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
+ if (.not. allocated(osc)) then
+ allocate(osc(4, nos(l)))
+ else if (j == 1) then
+ old_size = size(osc, 2)
+ allocate(tmp_osc(4, old_size + nos(l)))
+ tmp_osc(1:3, 1:old_size) = osc
+ deallocate(osc)
+ call move_alloc(tmp_osc, osc)
+ endif
+ jos = jos + 1
+ read(11, *) (osc(k, jos), k = 1, 3)
+ oscType(jos) = Lorentz
+ if (osc(2, jos) + 1 < 1e-7) then
+ oscType(jos) = Drude
+ osc(2, jos) = osc(3, jos)
+ osc(3, jos) = 0
+ endif
+ enddo
+ endif
+ enddo
+ close (unit = 11)
+
+ noPoints = 1 + int((wnmax - wnmin) / dwn)
+ wn = wnmin
+ allocate (eps(layers))
+
+ open(unit = 12, file = 'build-f90-new/epsout')
+ do i = 0, noPoints
+ call seteps(layers, nos, oscType, osc, epsinf, wn, eps, layers)
+ write (12, '(f8.2, 30g16.8)') wn, (eps(j), j = 1, layers)
+ wn = wn + dwn
+ end do
+ close (unit = 12)
+
+ write(*,*) 'setepsdriver-f90-new-1 stop'
+end program setepsdriver
diff --git a/tests/seteps/setepsdriver-f90-new-2.f90 b/tests/seteps/setepsdriver-f90-new-2.f90
new file mode 100644
index 0000000000000000000000000000000000000000..b228a92ad18920f63ebec2628efcb5cccf0bbe1c
--- /dev/null
+++ b/tests/seteps/setepsdriver-f90-new-2.f90
@@ -0,0 +1,96 @@
+program setepsdriver
+
+ use seteps_mod
+ use oscillator_mod
+
+ implicit none
+
+ double precision :: wnmin, wnmax, dwn, wn
+
+ double precision, allocatable :: osc(:,:), epsinf(:)
+ double complex, allocatable :: eps(:), epsReference(:,:)
+ integer, allocatable :: nos(:), oscType(:)
+ integer :: jos, layers
+ integer :: old_size
+ double precision, allocatable :: tmp_osc(:,:)
+ integer, allocatable :: tmp_oscType(:)
+
+ integer :: i, j, k, l, noPoints
+ integer :: IOstatus
+ logical :: testresult
+
+ write(*,*) 'setepsdriver-f90-new 2 start'
+
+ open(unit = 11, file = 'build-f90-new/setepsIn.txt')
+ read(11, *) wnmin, wnmax, dwn
+ read(11, *) layers
+
+ allocate (epsinf(layers), nos(layers))
+ allocate (eps(layers))
+ noPoints = 1 + int((wnmax - wnmin) / dwn)
+ allocate (epsReference(layers, noPoints))
+
+ jos = 0
+ do l = 1, layers
+ read(11, *) epsinf(l)
+ read(11, *) nos(l)
+ do j = 1, nos(l)
+ 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
+ if (.not. allocated(osc)) then
+ allocate(osc(4, nos(l)))
+ else if (j == 1) then
+ old_size = size(osc, 2)
+ allocate(tmp_osc(4, old_size + nos(l)))
+ tmp_osc(1:4, 1:old_size) = osc
+ deallocate(osc)
+ call move_alloc(tmp_osc, osc)
+ endif
+ jos = jos + 1
+ read(11, *) oscType(jos)
+ select case (oscType(jos))
+ case (Lorentz)
+ read(11, *) (osc(k, jos), k = 1, 3)
+ case (Kurosawa)
+ read(11, *) (osc(k, jos), k = 1, 4)
+ case (Drude)
+ read(11, *) (osc(k, jos), k = 1, 2)
+ case default
+ write (*,'(A)') '*** Error in oscillator type! ***'
+ end select
+ enddo
+ enddo
+ read(11, *, IOSTAT = IOstatus)
+ if (IOstatus == 0) then
+ do i = 1, noPoints
+ read(11, *) (epsReference(l, i), l = 1, layers)
+ enddo
+ endif
+ close (unit = 11)
+
+ wn = wnmin
+
+ testresult = .true.
+ open(unit = 12, file = 'build-f90-new/epsout')
+ do i = 1, noPoints
+ call seteps(layers, nos, oscType, osc, epsinf, wn, eps, layers)
+ write (12, '(30("(", g16.8, ",", g16.8, ") "))') (eps(j), j = 1, layers)
+ do j = 1, layers
+ if (abs((eps(j) - epsReference(j, i))/(eps(j) + epsReference(j, i))) > 1e-5) then
+! write (*,*) 'wn: ', wn, 'eps: ', eps(j), 'epsReference: ', epsReference(j, i)
+ testresult = .false.
+ endif
+ enddo
+ wn = wn + dwn
+ end do
+ close (unit = 12)
+
+ write(*,*) 'setepsdriver-f90-new 2 stop with test result: ', testresult
+end program setepsdriver