module reflections
  use const_and_precisions, only : wp_, comp_tiny, comp_eps, comp_huge, zero, one
  implicit none

  private
  public :: reflect, inters_linewall
  public :: linecone_coord, interssegm_coord, interssegm
  public :: wall_refl

contains

subroutine reflect(ki,nsurf,ko)
  implicit none
  real(wp_), intent(in), dimension(3) :: ki
  real(wp_), intent(in), dimension(3) :: nsurf
  real(wp_), intent(out), dimension(3) :: ko
  real(wp_) :: twokn,norm2
  norm2 = dot_product(nsurf,nsurf)
  if (norm2>zero) then
    twokn = 2.0_wp_*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(wp_), intent(in), dimension(3) :: xv,kv
  integer, intent(in) :: nw
  real(wp_), dimension(nw), intent(in) :: rw,zw
  real(wp_), intent(out) :: sint
  real(wp_), dimension(3), intent(out) :: normw
  integer :: i,j,ni,iint,first
  real(wp_), dimension(2) :: si,ti
  real(wp_) :: drw,dzw,xint,yint,rint,l,kxy
  sint=comp_huge
  iint=0
  normw=zero
  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)
    !discard solutions with s<=0
    first=ni+1
    do j=1,ni
      if (si(j)>zero) then
        first=j
        exit
      end if
    end do
    do j=first,ni
      if ((si(j)<sint .or. iint==0) .and. ti(j)>=zero .and. ti(j)<=one) 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 = hypot(xint, yint)
  l = hypot(drw, dzw)
  kxy = norm2(kv(1:2))
  normw(3) = -drw/l
  if (rint>zero) 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)>zero) normw=-normw
end subroutine inters_linewall



subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
  use utils, only : bubble
  implicit none
  real(wp_), intent(in), dimension(3) :: xv,kv
  real(wp_), intent(in), dimension(2) :: rs,zs
  real(wp_), dimension(2), intent(out) :: s,t
  integer, intent(out) :: n
  real(wp_) :: x0,y0,z0,kx,ky,kz
  real(wp_) :: 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)<comp_tiny) then
    !surface in horizontal plane
    if (abs(kz)<comp_tiny .or. abs(dr)<comp_tiny) 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)<comp_tiny) then
      !line parallel to cone generator
      if (abs(dr)<comp_tiny) 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<comp_tiny .and. abs(zrmin-zvertex)<comp_tiny) then
          !line passing by cone vertex
          !s(1) = srmin
          !t(1) = tvertex
          !n = 1
          n = 0
        else
          s(1) = -0.5_wp_*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(wp_), dimension(2), intent(in) :: xa,ya,xb,yb
  real(wp_), intent(out) :: s,t
  integer, intent(out) :: ierr
  real(wp_) :: 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)<comp_tiny) then
    s = zero
    t = zero
    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(wp_), dimension(2), intent(in) :: xa,ya,xb,yb
  logical :: interssegm
  real(wp_) :: s,t
  integer :: ierr
  interssegm = .false.
  call interssegm_coord(xa,ya,xb,yb,s,t,ierr)
  if (ierr==0 .and. s>=zero .and. s<=one .and. &
                    t>=zero .and. t<=one) interssegm = .true.
end function interssegm


subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,irfl)
  use limiter, only : rlim,zlim,nlim
  use utils,   only : inside
  implicit none
! arguments
  integer :: irfl
  real(wp_), dimension(3) :: xv,anv,xvrfl,anvrfl,walln
  complex(wp_) :: ext,eyt,extr,eytr
! local variables
  real(wp_) :: smax,rrm,zzm
  real(wp_), dimension(3) :: anv0,vv1,vv2,vv3
  complex(wp_) :: eztr
  complex(wp_), dimension(3) :: evin,evrfl 
!
  anv0=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2)
  rrm=1.0e-2_wp_*sqrt(xv(1)**2+xv(2)**2)
  zzm=1.0e-2_wp_*xv(3)
!
! computation of reflection coordinates and normal to the wall
  call inters_linewall(xv/1.0e2_wp_,anv0,rlim(1:nlim),zlim(1:nlim), &
                       nlim,smax,walln)
  smax=smax*1.0e2_wp_
  xvrfl=xv+smax*anv0
  irfl=1
  if (.not.inside(rlim,zlim,rrm,zzm)) then
    ! first wall interface is outside-inside
    if (dot_product(walln,walln)<tiny(walln)) then
      ! wall never hit
      xvrfl=xv
      anvrfl=anv0
      extr=ext
      eytr=eyt
      irfl=0
      return
    end if
    ! search second wall interface (inside-outside)
    call inters_linewall(xvrfl/1.0e2_wp_,anv0,rlim(1:nlim), &
                         zlim(1:nlim),nlim,smax,walln)
    smax=smax*1.0e2_wp_
    xvrfl=xvrfl+smax*anv0
    irfl=2
  end if
!
! rotation matrix from local to lab frame
  vv1(1)=anv0(2)
  vv1(2)=-anv0(1)
  vv1(3)=0.0_wp_
  vv2(1)=anv0(1)*anv0(3)
  vv2(2)=anv0(2)*anv0(3)
  vv2(3)=-anv0(1)*anv0(1)-anv0(2)*anv0(2)
  vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
  vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
  vv3=anv0
!
  evin=ext*vv1+eyt*vv2
! wave vector and electric field after reflection in lab frame
  anvrfl=anv0-2.0_wp_* &
         (anv0(1)*walln(1)+anv0(2)*walln(2)+anv0(3)*walln(3))*walln
  evrfl=-evin+2.0_wp_* &
         (evin(1)*walln(1)+evin(2)*walln(2)+evin(3)*walln(3))*walln
!
  vv1(1)=anvrfl(2)
  vv1(2)=-anvrfl(1)
  vv1(3)=0.0_wp_
  vv2(1)=anvrfl(1)*anvrfl(3)
  vv2(2)=anvrfl(2)*anvrfl(3)
  vv2(3)=-anvrfl(1)*anvrfl(1)-anvrfl(2)*anvrfl(2)
  vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
  vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
  vv3=anvrfl/sqrt(anvrfl(1)**2+anvrfl(2)**2+anvrfl(3)**2)
!
  extr=dot_product(vv1,evrfl)
  eytr=dot_product(vv2,evrfl)
  eztr=dot_product(vv3,evrfl)
end subroutine wall_refl

end module reflections