module reflections
  implicit none
  private
  integer, parameter :: r8=selected_real_kind(15,300)
  real(r8), parameter :: tinyr8=tiny(1._r8)
  public :: reflect,inters_linewall,inside
  public :: linecone_coord,interssegm_coord,interssegm
contains

subroutine reflect(ki,nsurf,ko)
  implicit none
  real(r8), intent(in), dimension(3) :: ki
  real(r8), intent(in), dimension(3) :: nsurf
  real(r8), intent(out), dimension(3) :: ko
  real(r8) :: twokn,norm2
  norm2 = dot_product(nsurf,nsurf)
  if (norm2>0.0_r8) then
    twokn = 2.0_r8*dot_product(ki,nsurf)/norm2
    ko=ki-twokn*nsurf
  else
    ko=ki
  end if
end subroutine reflect

subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
  implicit none
  real(r8), intent(in), dimension(3) :: xv,kv
  integer, intent(in) :: nw
  real(r8), dimension(nw), intent(in) :: rw,zw
  real(r8), intent(out) :: sint
  real(r8), dimension(3), intent(out) :: normw
  integer :: i,j,ni,iint
  real(r8), dimension(2) :: si,ti
  real(r8) :: drw,dzw,xint,yint,rint,l,kxy
  real(r8) :: tol=sqrt(epsilon(1.0_r8))
  sint=huge(sint)
  iint=0
  normw=0.0_r8
  do i=1,nw-1
    !search intersections with i-th wall segment
    call linecone_coord(xv,kv,rw(i:i+1),zw(i:i+1),si,ti,ni)
    do while (ni>0 .and. si(1)<=tol)
      !remove solutions with s<=0
      ni = ni-1
      si(1) = si(2)
      ti(1) = ti(2)
    end do
    do j=1,ni
      if ((si(j)<sint .or. iint==0) .and. ti(j)>=0._r8 .and. ti(j)<=1._r8) then
        !check intersection is in r,z range and keep the closest
        sint = si(j)
        iint = i
      end if
    end do
  end do
  if (iint==0) return
  !calculate wall normal at intersection point
  drw = rw(iint+1)-rw(iint)
  dzw = zw(iint+1)-zw(iint)
  xint = xv(1) + sint*kv(1)
  yint = xv(2) + sint*kv(2)
  rint = sqrt(xint**2+yint**2)
  l = sqrt(drw**2+dzw**2)
  kxy = sqrt(kv(1)**2+kv(2)**2)
  normw(3) = -drw/l
  if (rint>0.0_r8) then
    normw(1) = xint/rint*dzw/l
    normw(2) = yint/rint*dzw/l
  else
    normw(1) = kv(1)/kxy*dzw/l
    normw(2) = kv(2)/kxy*dzw/l
  end if
  !reverse normal if k.n>0
  if (dot_product(normw,kv)>0.0_r8) normw=-normw
end subroutine inters_linewall

subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
  implicit none
  real(r8), intent(in), dimension(3) :: xv,kv
  real(r8), intent(in), dimension(2) :: rs,zs
  real(r8), dimension(2), intent(out) :: s,t
  integer, intent(out) :: n
  real(r8) :: x0,y0,z0,kx,ky,kz
  real(r8) :: dr,dz,r,a,bhalf,c,delta,tvertex,zvertex,srmin,rmin,zrmin
  x0=xv(1)
  y0=xv(2)
  z0=xv(3)
  kx=kv(1)
  ky=kv(2)
  kz=kv(3)
  dr = rs(2)-rs(1)
  dz = zs(2)-zs(1)
  s = 0
  t = 0
  if (abs(dz)<tinyr8) then
    !surface in horizontal plane
    if (abs(kz)<tinyr8 .or. abs(dr)<tinyr8) then
      n = 0
    else
      s(1) = (zs(1)-z0)/kz
      r = sqrt((x0+s(1)*kx)**2+(y0+s(1)*ky)**2)
      t(1) = (r-rs(1))/dr
      n = 1
    end if
  else
    a = (kx**2+ky**2) - (dr/dz*kz)**2
    bhalf = -dr/dz*kz*rs(1) + (kx*x0 + ky*y0) - (dr/dz)**2*kz*(z0-zs(1))
    c = (x0**2+y0**2) - (rs(1) + dr/dz*(z0-zs(1)))**2
    if (abs(a)<tinyr8) then
      !line parallel to cone generator
      if (abs(dr)<tinyr8) then
        !cylinder and vertical line
        n = 0
      else
        tvertex = -rs(1)/dr
        zvertex = zs(1) + tvertex*dz
        srmin = -(kx*x0 + ky*y0)/(kx**2+ky**2)
        rmin = sqrt((x0+srmin*kx)**2+(y0+srmin*ky)**2)
        zrmin = z0 + srmin*kz
        if (rmin<tinyr8 .and. abs(zrmin-zvertex)<tinyr8) then
          !line passing by cone vertex
          !s(1) = srmin
          !t(1) = tvertex
          !n = 1
          n = 0
        else
          s(1) = -0.5_r8*c/bhalf
          t(1) = (kz*s(1)+(z0-zs(1)))/dz
          n = 1
        end if
      end if
    else
      delta = bhalf**2 - a*c
      if (delta<0) then
        n = 0
      else
        s(1) = (-bhalf+sqrt(delta))/a
        s(2) = (-bhalf-sqrt(delta))/a
        call bubble(s,2)
        t(:) = (kz*s(:)+(z0-zs(1)))/dz
        n = 2
      end if
    end if
  end if
end subroutine linecone_coord

subroutine interssegm_coord(xa,ya,xb,yb,s,t,ierr)
  implicit none
  real(r8), dimension(2), intent(in) :: xa,ya,xb,yb
  real(r8), intent(out) :: s,t
  integer, intent(out) :: ierr
  real(r8) :: crossprod,dxa,dya,dxb,dyb
  dxa = xa(2)-xa(1)
  dya = ya(2)-ya(1)
  dxb = xb(2)-xb(1)
  dyb = yb(2)-yb(1)
  crossprod = dxb*dya - dxa*dyb
  if (abs(crossprod)<tiny(crossprod)) then
    s = 0.0_r8
    t = 0.0_r8
    ierr = 1
  else
    s = (dyb*(xa(1)-xb(1)) - dxb*(ya(1)-yb(1)))/crossprod
    t = (dya*(xa(1)-xb(1)) - dxa*(ya(1)-yb(1)))/crossprod
    ierr = 0
  end if
end subroutine interssegm_coord

function interssegm(xa,ya,xb,yb)
  implicit none
  real(r8), dimension(2), intent(in) :: xa,ya,xb,yb
  logical :: interssegm
  real(r8) :: s,t
  integer :: ierr
  interssegm = .false.
  call interssegm_coord(xa,ya,xb,yb,s,t,ierr)
  if (ierr==0 .and. s>=0._r8 .and. s<=1._r8 .and. &
                    t>=0._r8 .and. t<=1._r8) interssegm = .true.
end function interssegm

function inside(xc,yc,n,x,y)
  implicit none
  integer, intent(in) :: n
  real(r8), dimension(n), intent(in) :: xc,yc
  real(r8), intent(in) :: x,y
  logical :: inside
  integer, dimension(n) :: jint
  real(r8), dimension(n) :: xint
  real(r8), dimension(n+1) :: xclosed,yclosed
  integer :: i,nj
  xclosed(1:n)=xc(1:n)
  yclosed(1:n)=yc(1:n)
  xclosed(n+1)=xc(1)
  yclosed(n+1)=yc(1)
  call locate_unord(yclosed,n+1,y,jint,n,nj)
  inside=.false.
  if (nj==0) return
  do i=1,nj
    xint(i)=intlin(yclosed(jint(i)),xclosed(jint(i)), &
                   yclosed(jint(i)+1),xclosed(jint(i)+1),y)
  end do
  call bubble(xint,nj)
  inside=(mod(locate(xint,nj,x),2)==1)
end function inside

function intlin(x1,y1,x2,y2,x) result(y)
  !linear interpolation
  !must be x1 != x2
  implicit none
  real(r8),intent(in) :: x1,y1,x2,y2,x
  real(r8) :: y
  real(r8) :: a
  a=(x2-x)/(x2-x1)
  y=a*y1+(1._r8-a)*y2
end function intlin

subroutine locate_unord(a,n,x,j,m,nj)
  implicit none
  integer, intent(in) :: n,m
  integer, intent(out) :: nj
  real(r8), dimension(n), intent(in) :: a
  real(r8), intent(in) :: x
  integer, dimension(m), intent(inout) :: j
  integer :: i
  nj=0
  do i=1,n-1
    if (x>a(i).neqv.x>a(i+1)) then
      nj=nj+1
      if (nj<=m) j(nj)=i
    end if
  end do
end subroutine locate_unord

function locate(a,n,x) result(j)
  !Given an array a(n), and a value x, with a(n) monotonic, either
  !increasing or decreasing, returns a value j such that
  !a(j) < x <= a(j+1) for a increasing, and such that
  !a(j+1) < x <= a(j) for a decreasing.
  !j=0 or j=n  indicate that x is out of range  (Numerical Recipes)
  implicit none
  integer, intent(in) :: n
  real(r8), dimension(n), intent(in) :: a
  real(r8), intent(in) :: x
  integer :: j
  integer :: jl,ju,jm
  logical :: incr
  jl=0
  ju=n+1
  incr=a(n)>a(1)
  do while ((ju-jl)>1)
    jm=(ju+jl)/2
    if(incr.eqv.(x>a(jm))) then
      jl=jm
    else
      ju=jm
    endif
  end do
  j=jl
end function locate

subroutine order(p,q)
  !returns p,q in ascending order
  implicit none
  real(r8), intent(inout) :: p,q
  real(r8) :: temp
  if (p>q) then
    temp=p
    p=q
    q=temp
  end if
end subroutine order

subroutine bubble(a,n)
  !bubble sorting of array a
  implicit none
  integer, intent(in) :: n
  real(r8), dimension(n), intent(inout) :: a
  integer :: i, j
  do i=1,n
    do j=n,i+1,-1
      call order(a(j-1), a(j))
    end do
  end do
end subroutine bubble

end module reflections