!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              COMPUTE_NORMALS    Apr. 2007                                    |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine compute_normals (ns,x,y,z,nt,icon,xn,yn,zn)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! given a set of points, and a connectivity array of their triangulation,
! this routine computes the normal to each point through cross-products:
! at first one computes the normal to each triangle and add the vector to the 
! three points that make the triangle
! then one loops over the points, and normalises the vectors.

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

!use mpi

implicit none

include 'mpif.h'

integer ns
double precision x(ns),y(ns),z(ns)
integer nt
integer icon(3,nt)
double precision xn(ns),yn(ns),zn(ns)

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer i,i1,i2,i3,ij,k
integer iproc,nproc,ierr
integer, dimension(:),allocatable :: nb
double precision x1,x2,x3
double precision y1,y2,y3
double precision z1,z2,z3
double precision xne,yne,zne
double precision xyzn

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

allocate(nb(ns))
nb=0

xn=0.d0
yn=0.d0
zn=0.d0
do i=1,nt
   i1=icon(1,i)
   i2=icon(2,i)
   i3=icon(3,i)
   x1=x(i1)
   x2=x(i2)
   x3=x(i3)
   y1=y(i1)
   y2=y(i2)
   y3=y(i3)
   z1=z(i1) 
   z2=z(i2) 
   z3=z(i3) 
   xne=(y2-y1)*(z3-z1)-(y3-y1)*(z2-z1)
   yne=(z2-z1)*(x3-x1)-(z3-z1)*(x2-x1)
   zne=(x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)
   xyzn=sqrt(xne**2+yne**2+zne**2)
   xne=xne/xyzn
   yne=yne/xyzn
   zne=zne/xyzn
   do k=1,3
      ij=icon(k,i)
      xn(ij)=xn(ij)+xne
      yn(ij)=yn(ij)+yne
      zn(ij)=zn(ij)+zne
      nb(ij)=nb(ij)+1
   enddo
enddo
do i=1,ns
   xyzn=sqrt(xn(i)**2+yn(i)**2+zn(i)**2)
   xn(i)=xn(i)/xyzn
   yn(i)=yn(i)/xyzn
   zn(i)=zn(i)/xyzn
enddo

if (minval(nb).eq.0) then
   do i=1,ns
      if (nb(i).eq.0.and.iproc.eq.0) then
         write (8,*) 'Particle ',i,' at ',x(i),y(i),z(i),' not connected'
      end if
   enddo
   call stop_run ('error in compute_normals_from_triangles$')
else
   if (iproc.eq.0) then
      write(8,*) 'in compute_normals_from_triangles:'
      write(8,*) 'minval(nb)=',minval(nb),' neighbours'
      write(8,*) 'maxval(nb)=',maxval(nb),' neighbours'
   end if            
end if

deallocate(nb)

return
end

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------