anchor.diff:  c === ./dirac.f      
             real*8 dirac_value
             integer nearest_node_index
  +          logical calculate_dirac
   c === ./lininterp.f
             real*8 linterpvalue
   c ==========      
  
positive.code:           parameter (thermo_coefficient = 1d0)
  c === probably should fix it later      
  
            iSYS(1) = 3
            iSYS(2) = 3
  
            Do i = 1, 3
               Do j = 1, 3
                  Coef(i, j) = 0D0
               End do
            End do
  
            Coef(1, 1) = thermo_coefficient
            Coef(2, 2) = thermo_coefficient
            Coef(3, 3) = thermo_coefficient
  
            Ddiff = TENSOR_SYMMETRIC
  
            Return
        End
  
negative.code:               if (tet_nodes(i) .eq. nearest_node_index) then
                    calculate_dirac = .true.
                    
                    swap_temporary = tet_nodes(i)
                    tet_nodes(i) = tet_nodes(4)
                    tet_nodes(4) = swap_temporary
                    goto 10
                end if
            end do
  