anchor.diff:      do nd=1,num_data
          nsp=map_seed_to_elem(nd) !space # that random point belongs to
          if(nsp.eq.ne1)then !random point belongs to this element space
  -          NPOINTS=NPOINTS+1
             dist = -scalar_product_3(norml,data_xyz(1,nd)) - norml(4) ! distance between two planes
             if(dist.ge.0.0_dp)then
                if(dat1.eq.0) nd1_1st = nd
  
positive.code: 
  !!! check whether enough seed points remaining in BOTH seed groups for branching to be done
  !!! (note: this could be improved to continue branching in one set of seeds)
               if(enough_points)then
                  do N = 1,2 !for each of the two new branches
                     ! Set up arrays for new element and node
                     ! after create_new_node the current element == ne and current node == np
  
negative.code:     COUNT=0
      do while(.NOT.INTERNAL.and.COUNT.lt.2)
         candidate_xyz(1:3) = node_xyz(1:3,np-1)
         internal=.true.
         if(.not.internal)then ! halve the length, and halve the angle from parent
            length = distance_between_points(node_xyz(1,np1),node_xyz(1,np-1))
            node_xyz(1:3,np-1) = node_xyz(1:3,np1) + 0.5_dp*length*elem_direction(1:3,ne-1)
            candidate_xyz(1:3) = node_xyz(1:3,np-1)
            call reduce_branch_angle(np1,np2,np-1,candidate_xyz,0.5_dp) ! reduces the branch angle by a half
  