get setepsdriver from branch newSeteps

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

tests/seteps/seteps-old.f90

+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

tests/seteps/setepsdriver-f90-new-1.f90

+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

tests/seteps/setepsdriver-f90-new-2.f90

+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