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) < epsilon(1.0d0) / 2) 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