diff --git a/examples/examples_turbine/004_compare_yaw_loss.py b/examples/examples_turbine/004_compare_yaw_loss.py
new file mode 100644
index 000000000..1510689cc
--- /dev/null
+++ b/examples/examples_turbine/004_compare_yaw_loss.py
@@ -0,0 +1,87 @@
+"""
+Example: Change operation model and compare power loss in yaw.
+
+This example illustrates how to define different operational models and compares
+the power loss resulting from yaw misalignment across these various models.
+"""
+
+import itertools
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from floris import FlorisModel, TimeSeries
+
+
+# Parameters
+N = 101  # How many steps to cover yaw range in
+yaw_max = 30  # Maximum yaw angle to test
+
+# Set up the yaw angle sweep
+yaw_angles = np.zeros((N, 1))
+yaw_angles[:, 0] = np.linspace(-yaw_max, yaw_max, N)
+print(yaw_angles.shape)
+
+
+def evaluate_yawed_power(wsp: float, op_model: str) -> float:
+    print(f"Evaluating model: {op_model}   wind speed: {wsp} m/s")
+
+    # Grab model of FLORIS
+    fmodel = FlorisModel("../inputs/gch.yaml")
+
+    # Run N cases by setting up a TimeSeries (which is just several independent simulations)
+    wind_directions = np.ones(N) * 270.0
+    fmodel.set(
+        wind_data=TimeSeries(
+            wind_speeds=wsp,
+            wind_directions=wind_directions,
+            turbulence_intensities=0.06,
+        )
+    )
+
+    yaw_angles = np.array(
+        [(yaw, 0.0, 0.0) for yaw in np.linspace(-yaw_max, yaw_max, N)]
+    )
+    fmodel.set_operation_model(op_model)
+    fmodel.set(yaw_angles=yaw_angles)
+    fmodel.run()
+
+    # Save the power output results in kW
+    return fmodel.get_turbine_powers()[:, 0] / 1000
+
+
+# Loop over the operational models and wind speeds to compare
+op_models = ["simple", "cosine-loss", "mit-loss"]
+wind_speeds = [11.0, 11.5, 15.0]
+results = {}
+for op_model, wsp in itertools.product(op_models, wind_speeds):
+
+    # Save the power output results in kW
+    results[(op_model, wsp)] = evaluate_yawed_power(wsp, op_model)
+# Plot the results
+fig, axes = plt.subplots(1, len(wind_speeds), sharey=True)
+
+colors = ["C0", "k", "r"]
+linestyles = ["solid", "dashed", "dotted"]
+for wsp, ax in zip(wind_speeds, axes):
+    ax.set_title(f"wsp: {wsp} m/s")
+    ax.set_xlabel("Yaw angle [deg]")
+    ax.grid(True)
+    for op_model, c, ls in zip(op_models, colors, linestyles):
+
+        upstream_yaw_angle = yaw_angles[:, 0]
+        central_power = results[(op_model, wsp)][upstream_yaw_angle == 0]
+        ax.plot(
+            upstream_yaw_angle,
+            results[(op_model, wsp)] / central_power,
+            label=op_model,
+            color=c,
+            linestyle=ls,
+        )
+
+ax.grid(True)
+ax.legend()
+axes[0].set_xlabel("Yaw angle [deg]")
+axes[0].set_ylabel("Normalized turbine power [-]")
+
+plt.show()
diff --git a/floris/core/rotor_velocity.py b/floris/core/rotor_velocity.py
index 43d4e3077..c167f7f13 100644
--- a/floris/core/rotor_velocity.py
+++ b/floris/core/rotor_velocity.py
@@ -15,6 +15,7 @@
     NDArrayObject,
 )
 from floris.utilities import cosd
+from floris.core.turbine.mit_models.momentum import Heck
 
 
 def rotor_velocity_yaw_cosine_correction(
@@ -239,3 +240,49 @@ def rotor_velocity_air_density_correction(
     # Produce equivalent velocities at the reference air density
 
     return (air_density/ref_air_density)**(1/3) * velocities
+
+
+def mit_rotor_axial_induction(Cts: NDArrayFloat, yaw_angles: NDArrayFloat)-> NDArrayFloat:
+    """
+    Computes the axial induction of a yawed rotor given the yaw-aligned thrust
+    coefficient and yaw angles using the yawed actuator disk model developed at
+    MIT as described in Heck et al. 2023. Assumes the modified thrust
+    coefficient, C_T', is invariant to yaw misalignment angle.
+
+    Args
+        Cts (NDArrayFloat): Yaw-aligned thrust coefficient(s).
+        yaw_angles (NDArrayFloat): Rotor yaw angle(s) in degrees.
+
+    Returns: NDArrayFloat: Axial induction factor(s) of the yawed rotor.
+    """
+    Ctprime = - 4 * (Cts + 2 * np.sqrt(1 - Cts) - 2) / Cts
+    sol = Heck()(Ctprime, np.deg2rad(yaw_angles))
+
+    return sol.an
+
+
+def mit_rotor_velocity_yaw_correction(
+        Cts: NDArrayFloat,
+        yaw_angles: NDArrayFloat,
+        axial_inductions: NDArrayFloat,
+        rotor_effective_velocities: NDArrayFloat,
+) -> NDArrayFloat:
+    """
+    Computes adjusted rotor wind speeds given the yaw-aligned thrust
+    coefficient, yaw angles, and axial induction values using the yawed actuator
+    disk model developed at MIT as described in Heck et al. 2023. Assumes the
+    modified thrust coefficient, C_T', is invariant to yaw misalignment angle.
+
+    Args
+        Cts (NDArrayFloat): Yaw-aligned thrust coefficient(s).
+        yaw_angles (NDArrayFloat): Rotor yaw angle(s) in degrees.
+        axial_induction (NDArrayFloat): Rotor axial induction(s).
+        rotor_effective_velocities (NDArrayFloat) rotor effective wind speed(s) at the rotor. .
+
+    Returns: NDArrayFloat: corrected rotor effective wind speed(s) of the yawed rotor.
+    """
+    Ctprime = - 4 * (Cts + 2 * np.sqrt(1 - Cts) - 2) / Cts
+
+    ratio = (1 + 0.25 * Ctprime) * (1 - axial_inductions) * np.cos(np.deg2rad(yaw_angles))
+
+    return ratio * rotor_effective_velocities
diff --git a/floris/core/turbine/__init__.py b/floris/core/turbine/__init__.py
index a7cde822a..7c27ba1df 100644
--- a/floris/core/turbine/__init__.py
+++ b/floris/core/turbine/__init__.py
@@ -6,4 +6,5 @@
     PeakShavingTurbine,
     SimpleDeratingTurbine,
     SimpleTurbine,
+    MITTurbine,
 )
diff --git a/floris/core/turbine/mit_models/FixedPointIteration.py b/floris/core/turbine/mit_models/FixedPointIteration.py
new file mode 100644
index 000000000..e895adeac
--- /dev/null
+++ b/floris/core/turbine/mit_models/FixedPointIteration.py
@@ -0,0 +1,165 @@
+from dataclasses import dataclass
+from typing import Any, Callable, List, Protocol, Tuple
+
+import numpy as np
+from numpy.typing import ArrayLike
+
+
+class FixedPointIterationCompatible(Protocol):
+    def residual(self, *args, **kwargs) -> Tuple[ArrayLike]: ...
+
+    def initial_guess(self, *args, **kwargs) -> Tuple[ArrayLike]: ...
+
+
+@dataclass
+class FixedPointIterationResult:
+    converged: bool
+    niter: int
+    relax: float
+    max_resid: float
+    x: ArrayLike
+
+
+def _fixedpointiteration(
+    f: Callable[[ArrayLike, Any], np.ndarray],
+    x0: np.ndarray,
+    args=(),
+    kwargs={},
+    eps=0.00001,
+    maxiter=100,
+    relax=0,
+    callback=None,
+) -> FixedPointIterationResult:
+    """
+    Performs fixed-point iteration on function f until residuals converge or max
+    iterations is reached.
+
+    Args:
+        f (Callable): residual function of form f(x, *args, **kwargs) -> np.ndarray
+        x0 (np.ndarray): Initial guess
+        args (tuple): arguments to pass to residual function. Defaults to ().
+        kwargs (dict): keyword arguments to pass to residual function. Defaults to {}.
+        eps (float): Convergence tolerance. Defaults to 0.000001.
+        maxiter (int): Maximum number of iterations. Defaults to 100.
+        callback (Callable): optional callback function at each iteration of the form f(x0) -> None
+
+    Returns:
+        FixedPointIterationResult: Solution to residual function.
+    """
+
+    for c in range(maxiter):
+        residuals = f(x0, *args, **kwargs)
+
+        x0 = [_x0 + (1 - relax) * _r for _x0, _r in zip(x0, residuals)]
+        max_resid = [np.nanmax(np.abs(_r)) for _r in residuals]
+
+        if callback:
+            callback(x0)
+
+        if all(_r < eps for _r in max_resid):
+            converged = True
+            break
+    else:
+        converged = False
+
+    if maxiter == 0:
+        return FixedPointIterationResult(False, 0, np.nan, np.nan, x0)
+    return FixedPointIterationResult(converged, c, relax, max_resid, x0)
+
+
+def fixedpointiteration(
+    max_iter: int = 100, tolerance: float = 1e-6, relaxation: float = 0.0
+) -> FixedPointIterationCompatible:
+    """
+    Class decorator which adds a __call__ method to the class which performs
+    fixed-point iteration.
+
+    Args:
+        max_iter (int): Maximum number of iterations (default: 100)
+        tolerance (float): Convergence criteria (default: 1e-6)
+        relaxation (float): Relaxation factor between 0 and 1 (default: 0.0)
+
+    The class must contain 2 mandatory methods and 3
+    optional method:
+
+    mandatory:
+    initial_guess(self, *args, **kwargs)
+    residual(self, x, *args, **kwargs)
+
+    optional:
+    pre_process(self, *args, **kwargs) # Optional
+    post_process(self, result:FixedPointIterationResult) # Optional
+    callback(self, x) # Optional
+
+    """
+
+    def decorator(cls: FixedPointIterationCompatible) -> Callable:
+        def call(self, *args, **kwargs):
+            if hasattr(self, "pre_process"):
+                self.pre_process(*args, **kwargs)
+
+            callback = self.callback if hasattr(self, "callback") else None
+
+            x0 = self.initial_guess(*args, **kwargs)
+            result = _fixedpointiteration(
+                self.residual,
+                x0,
+                args=args,
+                kwargs=kwargs,
+                eps=tolerance,
+                maxiter=max_iter,
+                relax=relaxation,
+                callback=callback,
+            )
+
+            if hasattr(self, "post_process"):
+                return self.post_process(result, *args, **kwargs)
+            else:
+                return result
+
+        setattr(cls, "__call__", call)
+        return cls
+
+    return decorator
+
+
+def adaptivefixedpointiteration(
+    max_iter: int = 100, tolerance: float = 1e-6, relaxations: List[float] = [0.0]
+) -> Callable:
+    """
+    Class decorator which adds a __call__ method to the class which performs
+    fixed-point iteration. Same as `fixedpointiteration`, but takes a list of
+    relaxation factors, and iterates over all of them in order until convergence
+    is reached.
+    """
+
+    def decorator(cls: FixedPointIterationCompatible) -> Callable:
+        def call(self, *args, **kwargs):
+            if hasattr(self, "pre_process"):
+                self.pre_process(*args, **kwargs)
+            callback = self.callback if hasattr(self, "callback") else None
+
+            for relaxation in relaxations:
+                x0 = self.initial_guess(*args, **kwargs)
+                result = _fixedpointiteration(
+                    self.residual,
+                    x0,
+                    args=args,
+                    kwargs=kwargs,
+                    eps=tolerance,
+                    maxiter=max_iter,
+                    relax=relaxation,
+                    callback=callback,
+                )
+                if result.converged:
+                    break
+
+            if hasattr(self, "post_process"):
+                return self.post_process(result, *args, **kwargs)
+            else:
+                return result
+
+        setattr(cls, "__call__", call)
+        return cls
+
+    return decorator
diff --git a/floris/core/turbine/mit_models/__init__.py b/floris/core/turbine/mit_models/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/floris/core/turbine/mit_models/momentum.py b/floris/core/turbine/mit_models/momentum.py
new file mode 100644
index 000000000..245f07d86
--- /dev/null
+++ b/floris/core/turbine/mit_models/momentum.py
@@ -0,0 +1,142 @@
+from abc import ABCMeta
+from dataclasses import dataclass
+from typing import Optional, Union
+
+import numpy as np
+import numpy.typing as npt
+
+from .FixedPointIteration import fixedpointiteration
+
+
+@dataclass
+class MomentumSolution:
+    """Stores the results of the Unified Momentum model solution."""
+
+    Ctprime: float
+    yaw: float
+    an: Union[float, npt.ArrayLike]
+    u4: Union[float, npt.ArrayLike]
+    v4: Union[float, npt.ArrayLike]
+    x0: Union[float, npt.ArrayLike]
+    dp: Union[float, npt.ArrayLike]
+    dp_NL: Optional[Union[float, npt.ArrayLike]] = 0.0
+    niter: Optional[int] = 1
+    converged: Optional[bool] = True
+    beta: Optional[float] = 0.0
+
+    @property
+    def Ct(self):
+        """Returns the thrust coefficient Ct."""
+        return self.Ctprime * (1 - self.an) ** 2 * np.cos(self.yaw) ** 2
+
+    @property
+    def Cp(self):
+        """Returns the power coefficient Cp."""
+        return self.Ctprime * ((1 - self.an) * np.cos(self.yaw)) ** 3
+
+
+class MomentumBase(metaclass=ABCMeta):
+    pass
+
+
+class LimitedHeck(MomentumBase):
+    """
+    Solves the limiting case when v_4 << u_4. (Eq. 2.19, 2.20). Also takes Numpy
+    array arguments.
+    """
+
+    def __call__(self, Ctprime: float, yaw: float, **kwargs) -> MomentumSolution:
+        """
+        Args:
+            Ctprime (float): Rotor thrust coefficient.
+            yaw (float): Rotor yaw angle (radians).
+
+        Returns:
+            Tuple[float, float, float]: induction and outlet velocities.
+        """
+
+        a = Ctprime * np.cos(yaw) ** 2 / (4 + Ctprime * np.cos(yaw) ** 2)
+        u4 = (4 - Ctprime * np.cos(yaw) ** 2) / (4 + Ctprime * np.cos(yaw) ** 2)
+        v4 = (
+            -(4 * Ctprime * np.sin(yaw) * np.cos(yaw) ** 2)
+            / (4 + Ctprime * np.cos(yaw) ** 2) ** 2
+        )
+        dp = np.zeros_like(a)
+        x0 = np.inf * np.ones_like(a)
+        return MomentumSolution(Ctprime, yaw, a, u4, v4, x0, dp)
+
+
+@fixedpointiteration(max_iter=500, tolerance=0.00001, relaxation=0.1)
+class Heck(MomentumBase):
+    """
+    Solves the iterative momentum equation for an actuator disk model.
+    """
+
+    def __init__(self, v4_correction: float = 1.0):
+        """
+        Initialize the HeckModel instance.
+
+        Args:
+            v4_correction (float, optional): The premultiplier of v4 in the Heck
+            model. A correction factor applied to v4, with a default value of
+            1.0, indicating no correction. Lu (2023) suggests an empirical correction
+            of 1.5.
+
+        Example:
+            >>> model = HeckModel(v4_correction=1.5)
+        """
+        self.v4_correction = v4_correction
+
+    def initial_guess(self, Ctprime, yaw):
+        sol = LimitedHeck()(Ctprime, yaw)
+        return sol.an, sol.u4, sol.v4
+
+    def residual(self, x: np.ndarray, Ctprime: float, yaw: float) -> np.ndarray:
+        """
+        Residual function of yawed-actuator disk model in Eq. 2.15.
+
+        Args:
+            x (np.ndarray): (a, u4, v4)
+            Ctprime (float): Rotor thrust coefficient.
+            yaw (float): Rotor yaw angle (radians).
+
+        Returns:
+            np.ndarray: residuals of induction and outlet velocities.
+        """
+
+        a, u4, v4 = x
+        e_a = 1 - np.sqrt(1 - u4**2 - v4**2) / (np.sqrt(Ctprime) * np.cos(yaw)) - a
+
+        e_u4 = (1 - 0.5 * Ctprime * (1 - a) * np.cos(yaw) ** 2) - u4
+
+        e_v4 = (
+            -self.v4_correction
+            * 0.25
+            * Ctprime
+            * (1 - a) ** 2
+            * np.sin(yaw)
+            * np.cos(yaw) ** 2
+            - v4
+        )
+
+        return np.array([e_a, e_u4, e_v4])
+
+    def post_process(self, result, Ctprime: float, yaw: float):
+        if result.converged:
+            a, u4, v4 = result.x
+        else:
+            a, u4, v4 = np.nan * np.zeros_like([Ctprime, Ctprime, Ctprime])
+        dp = np.zeros_like(a)
+        x0 = np.inf * np.ones_like(a)
+        return MomentumSolution(
+            Ctprime,
+            yaw,
+            a,
+            u4,
+            v4,
+            x0,
+            dp,
+            niter=result.niter,
+            converged=result.converged,
+        )
+
diff --git a/floris/core/turbine/operation_models.py b/floris/core/turbine/operation_models.py
index a6c1ff160..9f461bc0d 100644
--- a/floris/core/turbine/operation_models.py
+++ b/floris/core/turbine/operation_models.py
@@ -20,6 +20,8 @@
     rotor_velocity_air_density_correction,
     rotor_velocity_tilt_cosine_correction,
     rotor_velocity_yaw_cosine_correction,
+    mit_rotor_axial_induction,
+    mit_rotor_velocity_yaw_correction,
 )
 from floris.type_dec import (
     NDArrayFloat,
@@ -698,3 +700,172 @@ def axial_induction(
         )
 
         return (1 - np.sqrt(1 - thrust_coefficient))/2
+
+
+class MITTurbine(BaseOperationModel):
+    """
+    Base class for turbine operation models. All turbine operation models must implement static
+    power(), thrust_coefficient(), and axial_induction() methods, which are called by power() and
+    thrust_coefficient() through the interface in the turbine.py module.
+
+    Args:
+        BaseClass (_type_): _description_
+
+    Raises:
+        NotImplementedError: _description_
+        NotImplementedError: _description_
+    """
+
+    @staticmethod
+    @abstractmethod
+    def power(
+        power_thrust_table: dict,
+        velocities: NDArrayFloat,
+        air_density: float,
+        yaw_angles: NDArrayFloat,
+        average_method: str = "cubic-mean",
+        cubature_weights: NDArrayFloat | None = None,
+        **kwargs,
+    ) -> None:
+
+
+        # Construct thrust coefficient interpolant
+        thrust_coefficient_interpolator = interp1d(
+            power_thrust_table["wind_speed"],
+            power_thrust_table["thrust_coefficient"],
+            fill_value=0.0001,
+            bounds_error=False,
+        )
+
+        # Compute the power-effective wind speed across the rotor
+        rotor_average_velocities = average_velocity(
+            velocities=velocities,
+            method=average_method,
+            cubature_weights=cubature_weights,
+        )
+
+        rotor_effective_velocities = rotor_velocity_air_density_correction(
+            velocities=rotor_average_velocities,
+            air_density=air_density,
+            ref_air_density=power_thrust_table["ref_air_density"]
+        )
+
+        thrust_coefficients = thrust_coefficient_interpolator(rotor_effective_velocities)
+
+
+        # mit_rotor_velocity_yaw_correction
+        axial_inductions = mit_rotor_axial_induction(thrust_coefficients, yaw_angles)
+
+
+        corrected_rotor_effective_velocities = mit_rotor_velocity_yaw_correction(
+           thrust_coefficients, yaw_angles, axial_inductions, rotor_effective_velocities
+        )
+
+        # Tilt correction?
+
+
+        # Construct power interpolant
+        power_interpolator = interp1d(
+            power_thrust_table["wind_speed"],
+            power_thrust_table["power"],
+            fill_value=0.0,
+            bounds_error=False,
+        )
+
+        # Compute power
+        power = power_interpolator(corrected_rotor_effective_velocities) * 1e3 # Convert to W
+
+        return power
+
+
+    @staticmethod
+    @abstractmethod
+    def thrust_coefficient(
+        power_thrust_table: dict,
+        velocities: NDArrayFloat,
+        air_density: float,
+        yaw_angles: NDArrayFloat,
+        average_method: str = "cubic-mean",
+        cubature_weights: NDArrayFloat | None = None,
+        **kwargs,
+    ) -> None:
+
+        # Construct thrust coefficient interpolant
+        thrust_coefficient_interpolator = interp1d(
+            power_thrust_table["wind_speed"],
+            power_thrust_table["thrust_coefficient"],
+            fill_value=0.0001,
+            bounds_error=False,
+        )
+
+        # Compute the power-effective wind speed across the rotor
+        rotor_average_velocities = average_velocity(
+            velocities=velocities,
+            method=average_method,
+            cubature_weights=cubature_weights,
+        )
+
+        rotor_effective_velocities = rotor_velocity_air_density_correction(
+            velocities=rotor_average_velocities,
+            air_density=air_density,
+            ref_air_density=power_thrust_table["ref_air_density"]
+        )
+
+        thrust_coefficients = thrust_coefficient_interpolator(rotor_effective_velocities)
+
+
+        # mit_rotor_velocity_yaw_correction
+        axial_inductions = mit_rotor_axial_induction(thrust_coefficients, yaw_angles)
+
+
+        corrected_rotor_effective_velocities = mit_rotor_velocity_yaw_correction(
+           thrust_coefficients, yaw_angles, axial_inductions, rotor_effective_velocities
+        )
+
+        # Tilt correction?
+
+        # Compute thrust coefficient
+        yawed_thrust_coefficients = thrust_coefficient_interpolator(corrected_rotor_effective_velocities) 
+
+        return yawed_thrust_coefficients
+
+    @staticmethod
+    @abstractmethod
+    def axial_induction(
+        power_thrust_table: dict,
+        velocities: NDArrayFloat,
+        air_density: float,
+        yaw_angles: NDArrayFloat,
+        average_method: str = "cubic-mean",
+        cubature_weights: NDArrayFloat | None = None,
+        **kwargs,
+    ):
+
+        # Construct thrust coefficient interpolant
+        thrust_coefficient_interpolator = interp1d(
+            power_thrust_table["wind_speed"],
+            power_thrust_table["thrust_coefficient"],
+            fill_value=0.0001,
+            bounds_error=False,
+        )
+
+        # Compute the power-effective wind speed across the rotor
+        rotor_average_velocities = average_velocity(
+            velocities=velocities,
+            method=average_method,
+            cubature_weights=cubature_weights,
+        )
+
+        rotor_effective_velocities = rotor_velocity_air_density_correction(
+            velocities=rotor_average_velocities,
+            air_density=air_density,
+            ref_air_density=power_thrust_table["ref_air_density"]
+        )
+
+        thrust_coefficients = thrust_coefficient_interpolator(rotor_effective_velocities)
+
+
+        # mit_rotor_velocity_yaw_correction
+        axial_inductions = mit_rotor_axial_induction(thrust_coefficients, yaw_angles)
+
+        return axial_inductions
diff --git a/floris/core/turbine/turbine.py b/floris/core/turbine/turbine.py
index 2f98c45ea..e586ffef7 100644
--- a/floris/core/turbine/turbine.py
+++ b/floris/core/turbine/turbine.py
@@ -19,6 +19,7 @@
     PeakShavingTurbine,
     SimpleDeratingTurbine,
     SimpleTurbine,
+    MITTurbine,
 )
 from floris.type_dec import (
     convert_to_path,
@@ -41,6 +42,7 @@
         "mixed": MixedOperationTurbine,
         "awc": AWCTurbine,
         "peak-shaving": PeakShavingTurbine,
+        "mit-loss": MITTurbine,
     },
 }
 
diff --git a/tests/turbine_operation_models_unit_test.py b/tests/turbine_operation_models_unit_test.py
index b50aab54b..c00fdc4b0 100644
--- a/tests/turbine_operation_models_unit_test.py
+++ b/tests/turbine_operation_models_unit_test.py
@@ -9,6 +9,7 @@
     POWER_SETPOINT_DEFAULT,
     SimpleDeratingTurbine,
     SimpleTurbine,
+    MITTurbine
 )
 from floris.utilities import cosd
 from tests.conftest import SampleInputs, WIND_SPEEDS
@@ -40,6 +41,10 @@ def test_submodel_attributes():
     assert hasattr(PeakShavingTurbine, "thrust_coefficient")
     assert hasattr(PeakShavingTurbine, "axial_induction")
 
+    assert hasattr(MITTurbine, "power")
+    assert hasattr(MITTurbine, "thrust_coefficient")
+    assert hasattr(MITTurbine, "axial_induction")
+
 def test_SimpleTurbine():
 
     n_turbines = 1
@@ -659,3 +664,55 @@ def test_PeakShavingTurbine():
     assert (test_power <= base_power).all()
     assert test_power[0,0] == base_power[0,0]
     assert test_power[-1,0] == base_power[-1,0]
+
+
+def test_MITTurbine():
+
+    # NOTE: These tests should be updated to reflect actual expected behavior
+    # of the MITTurbine model. Currently, match the CosineLossTurbine model.
+
+    n_turbines = 1
+    wind_speed = 10.0
+    turbine_data = SampleInputs().turbine
+
+    yaw_angles_nom = 0 * np.ones((1, n_turbines))
+    tilt_angles_nom = turbine_data["power_thrust_table"]["ref_tilt"] * np.ones((1, n_turbines))
+    yaw_angles_test = 20 * np.ones((1, n_turbines))
+    tilt_angles_test = 0 * np.ones((1, n_turbines))
+
+
+    # Check that power works as expected
+    test_power = MITTurbine.power(
+        power_thrust_table=turbine_data["power_thrust_table"],
+        velocities=wind_speed * np.ones((1, n_turbines, 3, 3)), # 1 findex, 1 turbine, 3x3 grid
+        air_density=turbine_data["power_thrust_table"]["ref_air_density"], # Matches ref_air_density
+        yaw_angles=yaw_angles_nom,
+        tilt_angles=tilt_angles_nom,
+        tilt_interp=None
+    )
+    truth_index = turbine_data["power_thrust_table"]["wind_speed"].index(wind_speed)
+    baseline_power = turbine_data["power_thrust_table"]["power"][truth_index] * 1000
+    assert np.allclose(baseline_power, test_power)
+
+    # Check that yaw and tilt angle have an effect
+    test_power = MITTurbine.power(
+        power_thrust_table=turbine_data["power_thrust_table"],
+        velocities=wind_speed * np.ones((1, n_turbines, 3, 3)), # 1 findex, 1 turbine, 3x3 grid
+        air_density=turbine_data["power_thrust_table"]["ref_air_density"], # Matches ref_air_density
+        yaw_angles=yaw_angles_test,
+        tilt_angles=tilt_angles_test,
+        tilt_interp=None
+    )
+    assert test_power < baseline_power
+
+    # Check that a lower air density decreases power appropriately
+    test_power = MITTurbine.power(
+        power_thrust_table=turbine_data["power_thrust_table"],
+        velocities=wind_speed * np.ones((1, n_turbines, 3, 3)), # 1 findex, 1 turbine, 3x3 grid
+        air_density=1.1,
+        yaw_angles=yaw_angles_nom,
+        tilt_angles=tilt_angles_nom,
+        tilt_interp=None
+    )
+    assert test_power < baseline_power
+