getCD Function

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

    alphaIndx = getInterval(this%AD,alphaQuery)
    machIndx = getInterval(this%MaD,machQuery)

    if (alphaIndx(1) .eq. alphaIndx(2)) then
      if (machIndx(1) .eq. machIndx(2)) then
        getCD = this%CD(alphaIndx(1),machIndx(1))
      else
        getCD = this%CD(alphaIndx(1),machIndx(1))+ &
          (machQuery-this%MaD(machIndx(1))) * &
          (this%CD(alphaIndx(1),machIndx(2))- &
          this%CD(alphaIndx(1),machIndx(1)))/ &
          (this%MaD(machIndx(2))-this%MaD(machIndx(1)))
      endif
    elseif (machIndx(1) .eq. machIndx(2)) then
      getCD = this%CD(alphaIndx(1),machIndx(1))+ &
        (alphaQuery-this%AD(alphaIndx(1))) * &
        (this%CD(alphaIndx(2),machIndx(1))- &
        this%CD(alphaIndx(1),machIndx(1)))/ &
        (this%AD(alphaIndx(2))-this%AD(alphaIndx(1)))
    else
      getCD = getBilinearInterp(alphaQuery,machQuery, &
        (/this%AD(alphaIndx(1)),this%AD(alphaIndx(2))/), &
        (/this%MaD(machIndx(1)),this%MaD(machIndx(2))/), &
        this%CD(alphaIndx(1),machIndx(1)), &
        this%CD(alphaIndx(1),machIndx(2)), &
        this%CD(alphaIndx(2),machIndx(1)), &
        this%CD(alphaIndx(2),machIndx(2)))
    endif
  end function getCD