anchor.diff:        
         voigt = C1*VF1 + C2*VF2
         reuss = C1_inv*VF1 + C2_inv*VF2
  -      call inverse(6,6,reuss,reuss_inv)
  +      reuss_inv = CIJ_StoC(reuss)
         rhave = rh1*VF1 + rh2*VF2
         
         Cave = (voigt + reuss_inv) /2.
  
positive.code:       C = C_in
        rh = rh_in
        do i=1,n
           call inverse(6,6,C(i,:,:),C_inv(i,:,:))
        enddo
  
  !  Initialise matrices to 0s
  
negative.code:       do i=1,n
           !  Find compliance from input stiffness
           C_temp = C_in(i,:,:)
           S_in = CIJ_CtoS(C_temp)
           S = S + VF_in(i)*S_in/sum(VF_in)
           rhave = rhave + VF_in(i)*rh_in(i)/sum(VF_in)
        enddo
        
        
  !  Find stiffness from compliance
        Cave = CIJ_StoC(S)
              
     end subroutine CIJ_Reuss_av
  !-------------------------------------------------------------------------------
  