diff --git a/src/reflections.f90 b/src/reflections.f90
new file mode 100644
index 0000000..ffe59bd
--- /dev/null
+++ b/src/reflections.f90
@@ -0,0 +1,287 @@
+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
+