Update seteps.f90

remove sumnos and add the anisotropic uniaxial substrate layer

testprograms/NewSeteps/seteps.f90

 module seteps_mod
 contains
 subroutine seteps(neps, nos, oscType, osc, epsinf, wn, name, eps, layers)
+
   use oscillator_mod
+
   implicit none
+
   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
   
-  integer :: i, j, k, sumnos
+  integer :: i, j, k
   double complex :: one, nomi, deno, addeps
   double precision :: wn2, b
 
-  sumnos = sum(nos)
   write (*,*) 'neps: ', neps
   write (*,*) 'nos:  ', nos
   write (*,*) 'Oscillators:'
-  do i = 1, sumnos
+  do i = 1, sum(nos)
     select case (oscType(i))
       case (Lorentz)
         write (*,'(A)') 'Lorentz'
@@ -55,13 +57,20 @@ subroutine seteps(neps, nos, oscType, osc, epsinf, wn, name, eps, layers)
          case (Drude)
           addeps = addeps - dcmplx(osc(1,k)**2, wn*(osc(2,k)-osc(2,k))) / dcmplx(wn2, wn*osc(2,j))
         case default
-          write (*,'(A)') '*** Error in oscillator type!'
+          write (*,'(A)') '*** Error in oscillator type! ***'
       end select
       k = k + 1
       eps(i) = epsinf(i) * (nomi / deno + addeps)
     enddo
   enddo
   
+  if (neps == layers + 1) then
+! The substrate is an anisotropic uniaxial material
+    eps(layers) = sqrt(eps(layers) * eps(layers + 1))
+    if (aimag(eps(layers)) < 0.0) then
+      eps(layers) = -eps(layers)
+    endif
+  endif
   return
 end subroutine seteps
 end module seteps_mod