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:       real(rs),intent(in)  :: VF_in(n),C_in(n,6,6),rh_in(n)
        real(rs),intent(out) :: Cave(6,6),rhave
        real(rs)             :: VF(n),C(n,6,6),C_inv(n,6,6),rh(n),&
                                voigt(6,6),reuss(6,6),reuss_inv(6,6)
        
  !  Normalise the volume fractions to sum to unity
        VF = VF_in / sum(VF_in)
  
negative.code:          endif
           if (average_type == 'r' .or. average_type == 'R' .or. &
              average_type == 'vrh' .or. average_type == 'VRH') then
              Srot = CIJ_CtoS(Crot)
              Sout = Sout + Srot/real(nrotations,rs)
           endif
           a = a + da
        enddo
  
        if (average_type == 'vrh' .or. average_type == 'VRH') then
           Crot = CIJ_StoC(Sout)
           Cout = (Cout + Crot)/2._rs
        endif
  
  