getCM Function

  function getCM(this,alphaQuery,machQuery)
    !! Returns value of 2-d linearly interpolated CM
    !! for given alphaQuery and machQuery values
  class(C81_class) :: this
    real, intent(in) :: alphaQuery, machQuery
    real :: getCM
    integer, dimension(2) :: alphaIndx, machIndx

    alphaIndx = getInterval(this%AM,alphaQuery)
    machIndx = getInterval(this%MaM,machQuery)

    if (alphaIndx(1) .eq. alphaIndx(2)) then
      if (machIndx(1) .eq. machIndx(2)) then
        getCM = this%CM(alphaIndx(1),machIndx(1))
      else
        getCM = this%CM(alphaIndx(1),machIndx(1))+ &
          (machQuery-this%MaM(machIndx(1))) * &
          (this%CM(alphaIndx(1),machIndx(2))- &
          this%CM(alphaIndx(1),machIndx(1)))/ &
          (this%MaM(machIndx(2))-this%MaM(machIndx(1)))
      endif
    elseif (machIndx(1) .eq. machIndx(2)) then
      getCM = this%CM(alphaIndx(1),machIndx(1))+ &
        (alphaQuery-this%AM(alphaIndx(1))) * &
        (this%CM(alphaIndx(2),machIndx(1))- &
        this%CM(alphaIndx(1),machIndx(1)))/ &
        (this%AM(alphaIndx(2))-this%AM(alphaIndx(1)))
    else
      getCM = getBilinearInterp(alphaQuery,machQuery, &
        (/this%AM(alphaIndx(1)),this%AM(alphaIndx(2))/), &
        (/this%MaM(machIndx(1)),this%MaM(machIndx(2))/), &
        this%CM(alphaIndx(1),machIndx(1)), &
        this%CM(alphaIndx(1),machIndx(2)), &
        this%CM(alphaIndx(2),machIndx(1)), &
        this%CM(alphaIndx(2),machIndx(2)))
    endif
  end function getCM