anchor.diff:      
       N_U = unit_vector(U)
       N_V = unit_vector(V)
  +!    write(*,*) n_u
  +!    write(*,*) n_v
       ANGLE = scalar_product_3(N_U,N_V)
       ANGLE = max(-1.0_dp,ANGLE)
       ANGLE = min(1.0_dp,ANGLE)
  -    ANGLE = cos(ANGLE)
  +    ANGLE = acos(ANGLE)
       
       angle_btwn_vectors=ANGLE
       
  
positive.code: 
  !!!###############################################################
  
    function point_internal_to_surface(triangles,point_xyz,vertex_xyz)
  !!! Cast a line in positive x-direction from each data point and 
  !!! then work out how many triangular elements it crosses. If even it is in the 
  !!! shape and if odd it is outside the shape
  
      integer,intent(in) :: triangles(:,:)
      real(dp),intent(in) :: point_xyz(3),vertex_xyz(:,:)
      logical :: point_internal_to_surface
  
  !!! Local Variables
      integer :: ntri,num_triangles
      real(dp) :: area,area_triangle,cofm_surfaces(3),denominator,&
           norm_v(3),point(3),P1(3),P2(3),P3(3),u
      real(dp),parameter :: dist_tol = 1.0e-4_dp, user_tol = 1.0e-14_dp
  
negative.code: 
  !!! list of functions
  
    ! angle_btwn_points
    ! .... returns the angle between three points
   
    ! angle_btwn_vectors
    ! .... returns the angle between two vectors
   
  