anchor.diff:       C = C_in
        
      !  Find stiffness from inverse
  -     call inverse(6,6,C,S)
  +     S = CIJ_CtoS(C)
        
      !  Calculate Voigt moduli
        Kv = (1._rs/9._rs) * (C(1,1) + C(2,2) + C(3,3) + 2._rs*(C(1,2) + C(2,3) + C(3,1)))
  
positive.code:       enddo
        CIJ_is_stable = .true.
     end function CIJ_is_stable
  !-------------------------------------------------------------------------------
  
  !==============================================================================
  
negative.code:       
        voigt = C1*VF1 + C2*VF2
        reuss = C1_inv*VF1 + C2_inv*VF2
        reuss_inv = CIJ_StoC(reuss)
        rhave = rh1*VF1 + rh2*VF2
        
        Cave = (voigt + reuss_inv) /2.
  