fix make of setepsdriver-f90.f90

tests/seteps/Makefile

@@ -18,9 +18,9 @@ all: build-f77/setepsdriver build-f90/setepsdriver build-WFW/setepsdriver
 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

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