 SUBROUTINE raycasting(sol,x,y,z,chi)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    NUM/raycasting *
!
!    NOM
! raycasting(sol,x,y,z,chi)
!
!    DESCRIPTION  
!*
!  Algorithme de raycasting permettant de determiner si la cellule courante est dans le "solide" ou non.
!  Raycasting algorithm to determine whether the current cell is in the ‘solid’ or not.
!*
! 
!    ENTREES / INPUT
!    sol   : solide / solid
!    x,y,z : coordonnees du point a tester / coordinates of the point to be tested
!
!    SORTIES / OUTPUT
!    chi   : indicatrice de presence (0: non, 1: oui) / presence indicator (0: no, 1: yes)
!
!    SOURCE
! Ce logiciel est regi par la licence [CeCILL-B]
! This software is governed by the [CeCILL-B] license
!=========================== DEBUT DES DECLARATIONS ==================== 
!.1-----  Implicit, Use
  use num_typ
  IMPLICIT NONE
!.2-----  Declaration
  real(kind=kind(0.d0)), intent(in)  :: x,y,z 
!! coordonnees du centre de la cellule a tester / coordinates of the centre of the cell to be tested
  type(solide), intent(inout)        :: sol !! Objet solide / Solid object
  real(kind=kind(0.d0)), intent(out) :: chi !! Indicatrice de presence (0: non, 1: oui) / Presence indicator (0: no, 1: yes)
  
  integer                            :: k,ind
  real(kind=kind(0.d0))              :: xn
  
  type(point)                        :: ori,v0,v1,v2
  integer                            :: cpt,intersect
  type(vect_3)                       :: dir
!=========================== DEBUT DU CODE EXECUTABLE ==================

  ind=0

  !~!  
  !~!  ! 3D : voir "Fast, minimum storage Ray/Triangle intersection" (1997) - Moller et Trumbore
  !~!  
    ori%x=x ! Origine du rayon
    ori%y=y ! Origine du rayon
    ori%z=z ! Origine du rayon
    dir%v(1)=sol%xg%v(1) - x                ! Vecteur directeur du rayon  
    dir%v(2)=sol%xg%v(2) - y
    dir%v(3)=sol%xg%v(3) - z
    dir%v(1)=1.d0                ! Vecteur directeur du rayon  
    dir%v(2)=0.d0
    dir%v(3)=0.d0
    xn=dsqrt(dir%v(1)**2+dir%v(2)**2+dir%v(3)**2)
    cpt = 0
    if (xn.ne.0.d0)dir%v=dir%v/xn
        !call display(ori)
        !write(*,*)dir,sol%xg

    DO k=1,sol%nb_face ! Boucle sur les faces/triangles
        
        ! Construction des points du triangle
        v0=sol%list_vertex(sol%list_facette(k)%vertex(1))
        v1=sol%list_vertex(sol%list_facette(k)%vertex(2))
        v2=sol%list_vertex(sol%list_facette(k)%vertex(3))
       
        ! Compteur d intersection
        cpt = cpt + intersect(V0,V1,V2,ori,dir)
        !write(*,*)intersect(V0,V1,V2,ori,dir)
        !call display(v0)
        !call display(v1)
        !call display(v2)
    END DO

    ind=0
    IF (mod(cpt,2).EQ.1) ind=1 ! Si le nombre d intersection est impair, alors la maille est dans le polygone.
  
  IF (ind.NE.0) chi=sol%id*1.d0
!===========================   FIN DE LA ROUTINE    ====================
END SUBROUTINE raycasting
!DOCFIN
!================================================================================================================
  
  

FUNCTION signe(det) ! Determine le signe d un reel
  integer :: signe
  real(kind=kind(0.d0)) :: det
  
  if (det>0) signe = 1
  if (det==0) signe = 0
  if (det<0) signe = -1
    
END FUNCTION

FUNCTION determinant2D(x1,y1,x2,y2) ! Calcul un déterminant 2x2
  real(kind=kind(0.d0)) :: determinant2D,x1,y1,x2,y2
  
  determinant2D = x1*y2 - x2*y1
    
END FUNCTION

FUNCTION SUB(P1,P2) ! Construit un vecteur à partir de 2 points en 3D
  use geo_typ
  type(point) :: P1,P2
  type(vect_3) :: SUB
  
  SUB%v(1)=P1%x-P2%x
  SUB%v(2)=P1%y-P2%y
  SUB%v(3)=P1%z-P2%z
END FUNCTION SUB

FUNCTION CROSS(V1,V2) ! Effectue le produit vectoriel entre 2 vecteurs 3D
  use geo_typ
  type(vect_3) :: CROSS,V1,V2
  
  CROSS%v(1)= V1%v(2)*V2%v(3)-V1%v(3)*V2%v(2)
  CROSS%v(2)= V1%v(3)*V2%v(1)-V1%v(1)*V2%v(3)
  CROSS%v(3)= V1%v(1)*V2%v(2)-V1%v(2)*V2%v(1)
     
END FUNCTION CROSS

FUNCTION DOT(V1,V2) ! Effectue le produit scalaire entre 2 vecteurs 3D
  use geo_typ
  type(vect_3) :: V1,V2 
  real(kind=kind(0.d0)) :: DOT
  
  DOT = V1%v(1)*V2%v(1)+V1%v(2)*V2%v(2)+V1%v(3)*V2%v(3)  
   
END FUNCTION DOT


FUNCTION intersect(V0,V1,V2,ori,dir) ! Determine s il existe une intersection entre un rayon et un triangle (voir "Fast, minimum storage Ray/Triangle intersection" (1997) - Moller et Trumbore)

  use geo_typ, only : point, vect_3
  type(point) :: V0,V1,V2,ori
  type(vect_3) :: E1,E2,dir,SUB,CROSS,P,T,Q
  real(kind=kind(0.d0)), parameter:: eps=1.d-10
  real(kind=kind(0.d0)) :: DOT,det,inv_det,u,v,tray
  integer :: intersect
  
  intersect = 0
  E1 = SUB(V1,V0)
  E2 = SUB(V2,V0)
  P = CROSS(dir,E2)
  det = DOT(E1,P)
  
  IF (abs(det) < eps) return ! Si le determinant est nul, le rayon et le triangle sont dans le meme plan et ne s intersectent donc pas => on s arrete
  
  inv_det=1.d0/det 
  T = SUB(ori,V0)  
  u = DOT(T,P)*inv_det
  
  IF (u<0.d0 .OR. u>1.d0) return ! La coordonnée barycentrique doit 0 < u < 1 sinon on est en dehors du triangle et on s arrete 
  
  Q=CROSS(T,E1)
  v=DOT(dir,Q)*inv_det
  IF (v<0.d0 .OR. (u+v)>1.d0) return ! La coordonnee barycentrique doit 0 < v < 1 et 0 < u + v < 1 sinon on est en dehors du triangle et on s arrete 
  
  tray = DOT(E2,Q)*inv_det
  IF (tray<0) return ! Il y a intersection que si t > 0 sinon on s arrete
  
  intersect = 1 ! Si aucune des conditions précedentes n a ete remplie alors le rayon a bien intersecte le triangle

END FUNCTION intersect
!================================================================================================================