diff --git a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity.py b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity.py
new file mode 100644
index 000000000..9fe69d335
--- /dev/null
+++ b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity.py
@@ -0,0 +1,79 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+import floris.layout_visualization as layoutviz
+from floris import FlorisModel
+from floris.flow_visualization import visualize_cut_plane
+
+
+visualize = True
+
+# Get a test fi (single turbine at 0,0)
+fmodel = FlorisModel("../inputs/gch.yaml")
+fmodel.set(layout_x=[0, 0, 600, 600], layout_y=[0, -300, 0, -300])
+
+fmodel.set(
+    wind_directions=[270.0],
+    wind_speeds=[8.0],
+    turbulence_intensities=[0.06],
+    wind_shear=0.0
+)
+fmodel.run()
+P_wo_het = fmodel.get_turbine_powers()/1e6
+
+het_inflow_config = {
+    "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]),
+    "y": np.array([-500.0, 500.0, -500.0, 500.0]),
+    "speed_multipliers": np.array([[1.0, 1.0, 1.0, 1.0]]),
+    "u": np.array([[8.0, 8.0, 8.0, 8.0]]),
+    "v": np.array([[-2.0, 0.0, -2.0, 0.0]]),
+    #"v": np.array([[0.0, 0.0, 0.0, 0.0]]),
+}
+
+fmodel.set(heterogeneous_inflow_config=het_inflow_config)
+fmodel.run()
+P_w_het = fmodel.get_turbine_powers()/1e6
+
+print("Difference (MW):", P_w_het - P_wo_het)
+
+if visualize:
+    fig, ax = plt.subplots(2, 1, figsize=(4,8))
+    fmodel.set(
+        heterogeneous_inflow_config={
+            "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]),
+            "y": np.array([-1000.0, 1000.0, -1000.0, 1000.0]),
+            "speed_multipliers":np.array([[1.0, 1.0, 1.0, 1.0]])
+        }
+    )
+    fmodel.run()
+    horizontal_plane = fmodel.calculate_horizontal_plane(
+        x_resolution=200,
+        y_resolution=100,
+        height=90.0,
+        x_bounds=(-200, 1000),
+        y_bounds=(-500, 500),
+    )
+    visualize_cut_plane(
+        horizontal_plane,
+        ax=ax[0],
+        label_contours=False,
+        title="Without WD het"
+    )
+
+    fmodel.set(heterogeneous_inflow_config=het_inflow_config)
+    horizontal_plane = fmodel.calculate_horizontal_plane(
+        x_resolution=200,
+        y_resolution=100,
+        height=90.0,
+        x_bounds=(-200, 1000),
+        y_bounds=(-500, 500),
+    )
+    #import ipdb; ipdb.set_trace()
+    visualize_cut_plane(
+        horizontal_plane,
+        ax=ax[1],
+        label_contours=False,
+        title="With WD het"
+    )
+
+    plt.show()
diff --git a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py
new file mode 100644
index 000000000..f9a521bb0
--- /dev/null
+++ b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py
@@ -0,0 +1,71 @@
+"""Example: Visualize flow by sweeping turbines
+
+Demonstrate the use calculate_horizontal_plane_with_turbines
+
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import floris.flow_visualization as flowviz
+from floris import FlorisModel
+
+
+fmodel = FlorisModel("../inputs/gch.yaml")
+
+viz_by_sweep = False
+
+# # Some wake models may not yet have a visualization method included, for these cases can use
+# # a slower version which scans a turbine model to produce the horizontal flow
+x = np.array([-100, 0, 100, 200, 300, 400, 500, 600])
+v = -np.array([[0, 1, 2, 3, 4, 5, 6, 7]])
+
+het_inflow_config = {
+    "x": np.repeat(x, 2),
+    "y": np.tile(np.array([-100.0, 100.0]), 8),
+    #"z": np.array([90.0, 90.0, 90.0, 91.0]),
+    "speed_multipliers": 1.0*np.ones((1, 16)),
+    "u": 8.0*np.ones((1, v.shape[1]*2)),
+    "v": np.repeat(v, 2, axis=1)
+    #"v": np.array([[0.0, 0.0, 0.0, 0.0]]),
+}
+
+# Set a 2 turbine layout
+fmodel.set(
+    layout_x=[0, 300],
+    layout_y=[0, 0],
+    wind_directions=[270],
+    wind_speeds=[8],
+    turbulence_intensities=[0.06],
+    heterogeneous_inflow_config=het_inflow_config,
+)
+
+if viz_by_sweep:
+    horizontal_plane = flowviz.calculate_horizontal_plane_with_turbines(
+        fmodel,
+        x_resolution=20,
+        y_resolution=10,
+        x_bounds=(-100, 500),
+        y_bounds=(-100, 100),
+    )
+    title = "Coarse turbine scan method"
+else:
+    horizontal_plane = fmodel.calculate_horizontal_plane(
+        x_resolution=200,
+        y_resolution=100,
+        height=90.0,
+        x_bounds=(-100, 500),
+        y_bounds=(-100, 100),
+    )
+    title = "Standard flow visualization calculation"
+
+fig, ax = plt.subplots(figsize=(10, 4))
+flowviz.visualize_cut_plane(
+    horizontal_plane,
+    ax=ax,
+    label_contours=True,
+    title=title
+)
+
+
+plt.show()
diff --git a/floris/core/core.py b/floris/core/core.py
index d12005ba5..6e6dea703 100644
--- a/floris/core/core.py
+++ b/floris/core/core.py
@@ -96,6 +96,7 @@ def __attrs_post_init__(self) -> None:
                 turbine_diameters=self.farm.rotor_diameters,
                 wind_directions=self.flow_field.wind_directions,
                 grid_resolution=self.solver["turbine_grid_points"],
+                heterogeneous_inflow_config=self.flow_field.heterogeneous_inflow_config,
             )
         elif self.solver["type"] == "turbine_cubature_grid":
             self.grid = TurbineCubatureGrid(
@@ -103,6 +104,7 @@ def __attrs_post_init__(self) -> None:
                 turbine_diameters=self.farm.rotor_diameters,
                 wind_directions=self.flow_field.wind_directions,
                 grid_resolution=self.solver["turbine_grid_points"],
+                heterogeneous_inflow_config=self.flow_field.heterogeneous_inflow_config,
             )
         elif self.solver["type"] == "flow_field_grid":
             self.grid = FlowFieldGrid(
@@ -110,6 +112,7 @@ def __attrs_post_init__(self) -> None:
                 turbine_diameters=self.farm.rotor_diameters,
                 wind_directions=self.flow_field.wind_directions,
                 grid_resolution=self.solver["flow_field_grid_points"],
+                heterogeneous_inflow_config=self.flow_field.heterogeneous_inflow_config,
             )
         elif self.solver["type"] == "flow_field_planar_grid":
             self.grid = FlowFieldPlanarGrid(
@@ -119,6 +122,7 @@ def __attrs_post_init__(self) -> None:
                 normal_vector=self.solver["normal_vector"],
                 planar_coordinate=self.solver["planar_coordinate"],
                 grid_resolution=self.solver["flow_field_grid_points"],
+                heterogeneous_inflow_config=self.flow_field.heterogeneous_inflow_config,
                 x1_bounds=self.solver["flow_field_bounds"][0],
                 x2_bounds=self.solver["flow_field_bounds"][1],
             )
@@ -243,6 +247,7 @@ def solve_for_points(self, x, y, z):
             turbine_coordinates=self.farm.coordinates,
             turbine_diameters=self.farm.rotor_diameters,
             wind_directions=self.flow_field.wind_directions,
+            heterogeneous_inflow_config=self.flow_field.heterogeneous_inflow_config,
             grid_resolution=1,
             x_center_of_rotation=self.grid.x_center_of_rotation,
             y_center_of_rotation=self.grid.y_center_of_rotation
diff --git a/floris/core/flow_field.py b/floris/core/flow_field.py
index 19b488518..fc9fb4cc1 100644
--- a/floris/core/flow_field.py
+++ b/floris/core/flow_field.py
@@ -105,12 +105,18 @@ def heterogeneous_config_validator(self, instance: attrs.Attribute, value: dict
             return
 
         # Check that the correct keys are supplied for the heterogeneous_inflow_config dict
-        for k in ["speed_multipliers", "x", "y"]:
+        for k in ["x", "y"]:
             if k not in value.keys():
                 raise ValueError(
-                    "heterogeneous_inflow_config must contain entries for 'speed_multipliers',"
-                    f"'x', and 'y', with 'z' optional. Missing '{k}'."
+                    "heterogeneous_inflow_config must contain entries for "
+                    f"'x' and 'y', with 'z' optional. Missing '{k}'."
                 )
+        if ("speed_multipliers" not in value.keys() and
+            ("u" not in value.keys() or "v" not in value.keys())):
+            raise ValueError(
+                "heterogeneous_inflow_config must contain entries for either 'speed_multipliers'"
+                " or 'u' and 'v'."
+            )
         if "z" not in value:
             # If only a 2D case, add "None" for the z locations
             value["z"] = None
@@ -130,7 +136,7 @@ def het_map_validator(self, instance: attrs.Attribute, value: list | None) -> No
 
 
     def __attrs_post_init__(self) -> None:
-        if self.heterogeneous_inflow_config is not None:
+        if self._speed_heterogeneity:
             self.generate_heterogeneous_wind_map()
 
 
@@ -281,7 +287,13 @@ def generate_heterogeneous_wind_map(self):
                 - **y**: A list of y locations at which the speed up factors are defined.
                 - **z** (optional): A list of z locations at which the speed up factors are defined.
         """
-        speed_multipliers = self.heterogeneous_inflow_config['speed_multipliers']
+        if "speed_multipliers" in self.heterogeneous_inflow_config:
+            speed_multipliers = self.heterogeneous_inflow_config['speed_multipliers']
+        else:
+            speed_multipliers = np.sqrt(
+                np.array(self.heterogeneous_inflow_config['u'])**2
+                + np.array(self.heterogeneous_inflow_config['v'])**2
+            ) / self.wind_speeds
         x = self.heterogeneous_inflow_config['x']
         y = self.heterogeneous_inflow_config['y']
         z = self.heterogeneous_inflow_config['z']
@@ -305,6 +317,10 @@ def generate_heterogeneous_wind_map(self):
 
         self.het_map = in_region
 
+    @property
+    def _speed_heterogeneity(self):
+        return self.heterogeneous_inflow_config is not None
+
     @staticmethod
     def interpolate_multiplier_xy(x: NDArrayFloat,
                                   y: NDArrayFloat,
diff --git a/floris/core/grid.py b/floris/core/grid.py
index 5e7353cd7..95a84211a 100644
--- a/floris/core/grid.py
+++ b/floris/core/grid.py
@@ -16,6 +16,7 @@
     NDArrayInt,
 )
 from floris.utilities import (
+    apply_wind_direction_heterogeneity_warping,
     reverse_rotate_coordinates_rel_west,
     rotate_coordinates_rel_west,
 )
@@ -52,6 +53,7 @@ class Grid(ABC, BaseClass):
     turbine_diameters: NDArrayFloat = field(converter=floris_array_converter)
     wind_directions: NDArrayFloat = field(converter=floris_array_converter)
     grid_resolution: int | Iterable = field()
+    heterogeneous_inflow_config: dict | None = field(init=True, default=None)
 
     n_turbines: int = field(init=False)
     n_findex: int = field(init=False)
@@ -63,6 +65,11 @@ class Grid(ABC, BaseClass):
     z_sorted_inertial_frame: NDArrayFloat = field(init=False)
     cubature_weights: NDArrayFloat = field(init=False, default=None)
 
+    use_turbine_specific_layouts: NDArrayFloat = field(init=False, default=False)
+    x_sorted_per_turbine: NDArrayFloat = field(init=False, default=None)
+    y_sorted_per_turbine: NDArrayFloat = field(init=False, default=None)
+    z_sorted_per_turbine: NDArrayFloat = field(init=False, default=None)
+
     @turbine_coordinates.validator
     def check_coordinates(self, instance: attrs.Attribute, value: np.ndarray) -> None:
         """
@@ -104,7 +111,7 @@ def grid_resolution_validator(self, instance: attrs.Attribute, value: int | Iter
     def set_grid(self) -> None:
         raise NotImplementedError("Grid.set_grid")
 
-@define
+@define(kw_only=True)
 class TurbineGrid(Grid):
     """See `Grid` for more details.
 
@@ -245,6 +252,40 @@ def set_grid(self) -> None:
         self.y_sorted = np.take_along_axis(_y, self.sorted_indices, axis=1)
         self.z_sorted = np.take_along_axis(_z, self.sorted_indices, axis=1)
 
+        if (self.heterogeneous_inflow_config is not None
+            and "u" in self.heterogeneous_inflow_config):
+            self.use_turbine_specific_layouts = True
+
+            # Rotate the input x, y, z to match the wind direction
+            het_coords = np.concatenate(
+                (
+                    np.array(self.heterogeneous_inflow_config["x"])[:, None],
+                    np.array(self.heterogeneous_inflow_config["y"])[:, None],
+                    np.zeros((len(self.heterogeneous_inflow_config["x"]), 1))
+                ),
+                axis=1
+            )
+            x_het, y_het, _, _, _ = rotate_coordinates_rel_west(
+                self.wind_directions,
+                het_coords,
+                x_center_of_rotation=self.x_center_of_rotation,
+                y_center_of_rotation=self.y_center_of_rotation
+            )
+
+            self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \
+                apply_wind_direction_heterogeneity_warping(
+                    self.x_sorted.mean(axis=(2,3)),
+                    self.y_sorted.mean(axis=(2,3)),
+                    self.z_sorted.mean(axis=(2,3)),
+                    self.x_sorted,
+                    self.y_sorted,
+                    self.z_sorted,
+                    x_het,
+                    y_het,
+                    self.heterogeneous_inflow_config["u"],
+                    self.heterogeneous_inflow_config["v"],
+                )
+
         # Now calculate grid coordinates in original frame (from 270 deg perspective)
         self.x_sorted_inertial_frame, self.y_sorted_inertial_frame, self.z_sorted_inertial_frame = \
             reverse_rotate_coordinates_rel_west(
@@ -256,7 +297,7 @@ def set_grid(self) -> None:
                 y_center_of_rotation=self.y_center_of_rotation,
             )
 
-@define
+@define(kw_only=True)
 class TurbineCubatureGrid(Grid):
     """
     This grid type arranges points throughout the swept area of the rotor based on the cubature
@@ -437,7 +478,7 @@ def get_cubature_coefficients(cls, N: int):
             "B": np.pi/N,
         }
 
-@define
+@define(kw_only=True)
 class FlowFieldGrid(Grid):
     """
     Args:
@@ -506,7 +547,7 @@ def set_grid(self) -> None:
                 y_center_of_rotation=self.y_center_of_rotation,
             )
 
-@define
+@define(kw_only=True)
 class FlowFieldPlanarGrid(Grid):
     """
     Args:
@@ -569,10 +610,6 @@ def set_grid(self) -> None:
                 indexing="ij"
             )
 
-            self.x_sorted = x_points[None, :, :, :]
-            self.y_sorted = y_points[None, :, :, :]
-            self.z_sorted = z_points[None, :, :, :]
-
         elif self.normal_vector == "x":  # Rules of thumb for cross plane
             if self.x1_bounds is None:
                 self.x1_bounds = (np.min(y) - 2 * max_diameter, np.max(y) + 2 * max_diameter)
@@ -587,10 +624,6 @@ def set_grid(self) -> None:
                 indexing="ij"
             )
 
-            self.x_sorted = x_points[None, :, :, :]
-            self.y_sorted = y_points[None, :, :, :]
-            self.z_sorted = z_points[None, :, :, :]
-
         elif self.normal_vector == "y":  # Rules of thumb for y plane
             if self.x1_bounds is None:
                 self.x1_bounds = (np.min(x) - 2 * max_diameter, np.max(x) + 10 * max_diameter)
@@ -605,10 +638,44 @@ def set_grid(self) -> None:
                 indexing="ij"
             )
 
-            self.x_sorted = x_points[None, :, :, :]
-            self.y_sorted = y_points[None, :, :, :]
-            self.z_sorted = z_points[None, :, :, :]
+        self.x_sorted = x_points[None, :, :, :]
+        self.y_sorted = y_points[None, :, :, :]
+        self.z_sorted = z_points[None, :, :, :]
+
+        if (self.heterogeneous_inflow_config is not None
+            and "u" in self.heterogeneous_inflow_config):
+            self.use_turbine_specific_layouts = True
+
+            # Rotate the input x, y, z to match the wind direction
+            het_coords = np.concatenate(
+                (
+                    np.array(self.heterogeneous_inflow_config["x"])[:, None],
+                    np.array(self.heterogeneous_inflow_config["y"])[:, None],
+                    np.zeros((len(self.heterogeneous_inflow_config["x"]), 1))
+                ),
+                axis=1
+            )
+            x_het, y_het, _, _, _ = rotate_coordinates_rel_west(
+                self.wind_directions,
+                het_coords,
+                x_center_of_rotation=self.x_center_of_rotation,
+                y_center_of_rotation=self.y_center_of_rotation
+            )
 
+            # Apply wind direction heterogeneity to generate turbine-specific layouts
+            self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \
+                apply_wind_direction_heterogeneity_warping(
+                    np.take_along_axis(x, x.argsort(axis=1), axis=1),
+                    np.take_along_axis(y, x.argsort(axis=1), axis=1),
+                    np.take_along_axis(z, x.argsort(axis=1), axis=1),
+                    self.x_sorted,
+                    self.y_sorted,
+                    self.z_sorted,
+                    x_het,
+                    y_het,
+                    self.heterogeneous_inflow_config["u"],
+                    self.heterogeneous_inflow_config["v"],
+                )
         # Now calculate grid coordinates in original frame (from 270 deg perspective)
         self.x_sorted_inertial_frame, self.y_sorted_inertial_frame, self.z_sorted_inertial_frame = \
             reverse_rotate_coordinates_rel_west(
@@ -620,13 +687,13 @@ def set_grid(self) -> None:
                 y_center_of_rotation=self.y_center_of_rotation,
             )
 
-@define
+@define(kw_only=True)
 class PointsGrid(Grid):
     """
     Args:
         turbine_coordinates (:py:obj:`NDArrayFloat`): Not used for PointsGrid, but
             required for the `Grid` super-class.
-        turbine_diameters (:py:obj:`NDArrayFloat`):  Not used for PointsGrid, but
+        turbine_diameters (:py:obj:`NDArrayFloat`): Not used for PointsGrid, but
             required for the `Grid` super-class.
         wind_directions (:py:obj:`NDArrayFloat`): Wind directions supplied by the user.
         grid_resolution (:py:obj:`int` | :py:obj:`Iterable(int,)`): Not used for PointsGrid, but
diff --git a/floris/core/solver.py b/floris/core/solver.py
index a7c3d8796..2feded092 100644
--- a/floris/core/solver.py
+++ b/floris/core/solver.py
@@ -25,7 +25,7 @@
 )
 from floris.core.wake_velocity.empirical_gauss import awc_added_wake_mixing
 from floris.type_dec import NDArrayFloat
-from floris.utilities import cosd
+from floris.utilities import cosd, update_model_args
 
 
 def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ngrid):
@@ -89,6 +89,20 @@ def sequential_solver(
         u_i = flow_field.u_sorted[:, i:i+1]
         v_i = flow_field.v_sorted[:, i:i+1]
 
+        if grid.use_turbine_specific_layouts:
+            deficit_model_args = update_model_args(
+                deficit_model_args,
+                grid.x_sorted_per_turbine[:, :, :, :, i],
+                grid.y_sorted_per_turbine[:, :, :, :, i],
+                grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+            deflection_model_args = update_model_args(
+                deflection_model_args,
+                grid.x_sorted_per_turbine[:, :, :, :, i],
+                grid.y_sorted_per_turbine[:, :, :, :, i],
+                grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+
         ct_i = thrust_coefficient(
             velocities=flow_field.u_sorted,
             turbulence_intensities=flow_field.turbulence_intensity_field_sorted,
@@ -287,6 +301,7 @@ def full_flow_sequential_solver(
         turbine_coordinates=turbine_grid_farm.coordinates,
         turbine_diameters=turbine_grid_farm.rotor_diameters,
         wind_directions=turbine_grid_flow_field.wind_directions,
+        heterogeneous_inflow_config=turbine_grid_flow_field.heterogeneous_inflow_config,
         grid_resolution=3,
     )
     turbine_grid_farm.expand_farm_properties(
@@ -327,6 +342,20 @@ def full_flow_sequential_solver(
         u_i = turbine_grid_flow_field.u_sorted[:, i:i+1]
         v_i = turbine_grid_flow_field.v_sorted[:, i:i+1]
 
+        if flow_field_grid.use_turbine_specific_layouts:
+            deficit_model_args = update_model_args(
+                deficit_model_args,
+                flow_field_grid.x_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.y_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+            deflection_model_args = update_model_args(
+                deflection_model_args,
+                flow_field_grid.x_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.y_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+
         ct_i = thrust_coefficient(
             velocities=turbine_grid_flow_field.u_sorted,
             turbulence_intensities=turbine_grid_flow_field.turbulence_intensity_field_sorted,
@@ -494,6 +523,20 @@ def cc_solver(
         z_i = np.mean(grid.z_sorted[:, i:i+1], axis=(2, 3))
         z_i = z_i[:, :, None, None]
 
+        if grid.use_turbine_specific_layouts:
+            deficit_model_args = update_model_args(
+                deficit_model_args,
+                grid.x_sorted_per_turbine[:, :, :, :, i],
+                grid.y_sorted_per_turbine[:, :, :, :, i],
+                grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+            deflection_model_args = update_model_args(
+                deflection_model_args,
+                grid.x_sorted_per_turbine[:, :, :, :, i],
+                grid.y_sorted_per_turbine[:, :, :, :, i],
+                grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+
         rotor_diameter_i = farm.rotor_diameters_sorted[:, i:i+1, None, None]
 
         mask2 = (
@@ -723,6 +766,7 @@ def full_flow_cc_solver(
         turbine_coordinates=turbine_grid_farm.coordinates,
         turbine_diameters=turbine_grid_farm.rotor_diameters,
         wind_directions=turbine_grid_flow_field.wind_directions,
+        heterogeneous_inflow_config=turbine_grid_flow_field.heterogeneous_inflow_config,
         grid_resolution=3,
     )
     turbine_grid_farm.expand_farm_properties(
@@ -766,6 +810,20 @@ def full_flow_cc_solver(
         u_i = turbine_grid_flow_field.u_sorted[:, i:i+1]
         v_i = turbine_grid_flow_field.v_sorted[:, i:i+1]
 
+        if flow_field_grid.use_turbine_specific_layouts:
+            deficit_model_args = update_model_args(
+                deficit_model_args,
+                flow_field_grid.x_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.y_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+            deflection_model_args = update_model_args(
+                deflection_model_args,
+                flow_field_grid.x_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.y_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+
         turb_avg_vels = average_velocity(turbine_grid_flow_field.u_sorted)
         turb_Cts = thrust_coefficient(
             velocities=turb_avg_vels,
@@ -930,6 +988,20 @@ def turbopark_solver(
         z_i = np.mean(grid.z_sorted[:, i:i+1], axis=(2, 3))
         z_i = z_i[:, :, None, None]
 
+        if grid.use_turbine_specific_layouts:
+            deficit_model_args = update_model_args(
+                deficit_model_args,
+                grid.x_sorted_per_turbine[:, :, :, :, i],
+                grid.y_sorted_per_turbine[:, :, :, :, i],
+                grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+            deflection_model_args = update_model_args(
+                deflection_model_args,
+                grid.x_sorted_per_turbine[:, :, :, :, i],
+                grid.y_sorted_per_turbine[:, :, :, :, i],
+                grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+
         Cts = thrust_coefficient(
             velocities=flow_field.u_sorted,
             turbulence_intensities=flow_field.turbulence_intensity_field_sorted,
@@ -1198,6 +1270,20 @@ def empirical_gauss_solver(
         z_i = np.mean(grid.z_sorted[:, i:i+1], axis=(2, 3))
         z_i = z_i[:, :, None, None]
 
+        if grid.use_turbine_specific_layouts:
+            deficit_model_args = update_model_args(
+                deficit_model_args,
+                grid.x_sorted_per_turbine[:, :, :, :, i],
+                grid.y_sorted_per_turbine[:, :, :, :, i],
+                grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+            deflection_model_args = update_model_args(
+                deflection_model_args,
+                grid.x_sorted_per_turbine[:, :, :, :, i],
+                grid.y_sorted_per_turbine[:, :, :, :, i],
+                grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+
         ct_i = thrust_coefficient(
             velocities=flow_field.u_sorted,
             turbulence_intensities=flow_field.turbulence_intensity_field_sorted,
@@ -1383,6 +1469,7 @@ def full_flow_empirical_gauss_solver(
         turbine_coordinates=turbine_grid_farm.coordinates,
         turbine_diameters=turbine_grid_farm.rotor_diameters,
         wind_directions=turbine_grid_flow_field.wind_directions,
+        heterogeneous_inflow_config=turbine_grid_flow_field.heterogeneous_inflow_config,
         grid_resolution=3,
     )
     turbine_grid_farm.expand_farm_properties(
@@ -1421,6 +1508,20 @@ def full_flow_empirical_gauss_solver(
         z_i = np.mean(turbine_grid.z_sorted[:, i:i+1], axis=(2,3))
         z_i = z_i[:, :, None, None]
 
+        if flow_field_grid.use_turbine_specific_layouts:
+            deficit_model_args = update_model_args(
+                deficit_model_args,
+                flow_field_grid.x_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.y_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+            deflection_model_args = update_model_args(
+                deflection_model_args,
+                flow_field_grid.x_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.y_sorted_per_turbine[:, :, :, :, i],
+                flow_field_grid.z_sorted_per_turbine[:, :, :, :, i],
+            )
+
         ct_i = thrust_coefficient(
             velocities=turbine_grid_flow_field.u_sorted,
             turbulence_intensities=turbine_grid_flow_field.turbulence_intensity_field_sorted,
diff --git a/floris/floris_model.py b/floris/floris_model.py
index 10f410225..27b60649c 100644
--- a/floris/floris_model.py
+++ b/floris/floris_model.py
@@ -1036,17 +1036,25 @@ def set_for_viz(self, findex: int, solver_settings: dict) -> None:
         """
 
         # If not None, set the heterogeneous inflow configuration
+        # TODO: add wind direction info here
         if self.core.flow_field.heterogeneous_inflow_config is not None:
             heterogeneous_inflow_config = {
                 'x': self.core.flow_field.heterogeneous_inflow_config['x'],
                 'y': self.core.flow_field.heterogeneous_inflow_config['y'],
-                'speed_multipliers':
-                    self.core.flow_field.heterogeneous_inflow_config['speed_multipliers'][findex:findex+1],
             }
-            if 'z' in self.core.flow_field.heterogeneous_inflow_config:
-                heterogeneous_inflow_config['z'] = (
-                    self.core.flow_field.heterogeneous_inflow_config['z']
-                )
+            for k in ["z"]: # Optional, independent of findex
+                if k in self.core.flow_field.heterogeneous_inflow_config:
+                    heterogeneous_inflow_config[k] = (
+                        self.core.flow_field.heterogeneous_inflow_config[k]
+                    )
+            for k in ["u", "v", "speed_multipliers"]: # Optional, dependent on findex
+                if k in self.core.flow_field.heterogeneous_inflow_config:
+                    if self.core.flow_field.heterogeneous_inflow_config[k] is not None:
+                        heterogeneous_inflow_config[k] = (
+                            self.core.flow_field.heterogeneous_inflow_config[k][findex:findex+1,:]
+                        )
+                    else:
+                        heterogeneous_inflow_config[k] = None
         else:
             heterogeneous_inflow_config = None
 
@@ -1058,7 +1066,7 @@ def set_for_viz(self, findex: int, solver_settings: dict) -> None:
             power_setpoints=self.core.farm.power_setpoints[findex:findex+1,:],
             awc_modes=self.core.farm.awc_modes[findex:findex+1,:],
             awc_amplitudes=self.core.farm.awc_amplitudes[findex:findex+1,:],
-            heterogeneous_inflow_config = heterogeneous_inflow_config,
+            heterogeneous_inflow_config=heterogeneous_inflow_config,
             solver_settings=solver_settings,
         )
 
diff --git a/floris/utilities.py b/floris/utilities.py
index 4f539f371..ef4e966e3 100644
--- a/floris/utilities.py
+++ b/floris/utilities.py
@@ -11,9 +11,12 @@
     Tuple,
 )
 
+import numexpr as ne
 import numpy as np
 import yaml
 from attrs import define, field
+from scipy.integrate import solve_ivp
+from scipy.interpolate import LinearNDInterpolator
 
 from floris.type_dec import floris_array_converter, NDArrayFloat
 
@@ -240,7 +243,7 @@ def rotate_coordinates_rel_west(
 
     # Calculate the difference in given wind direction from 270 / West
     wind_deviation_from_west = wind_delta(wind_directions)
-    wind_deviation_from_west = np.reshape(wind_deviation_from_west, (len(wind_directions), 1))
+    wind_deviation_from_west = np.expand_dims(wind_deviation_from_west, axis=-1)
 
     # Construct the arrays storing the turbine locations
     x_coordinates, y_coordinates, z_coordinates = coordinates.T
@@ -325,6 +328,252 @@ def reverse_rotate_coordinates_rel_west(
 
     return grid_x_reversed, grid_y_reversed, grid_z_reversed
 
+def apply_wind_direction_heterogeneity_warping(
+        x_turbines,
+        y_turbines,
+        z_turbines,
+        x_points,
+        y_points,
+        z_points,
+        wind_direction_x_points,
+        wind_direction_y_points,
+        wind_direction_u_values,
+        wind_direction_v_values,
+    ):
+    """
+    Args:
+        x_turbines (np.ndarray): x-coordinates of the turbines (n_findex x n_turbines).
+        y_turbines (np.ndarray): y-coordinates of the turbines (n_findex x n_turbines).
+        z_turbines (np.ndarray): z-coordinates of the turbines (n_findex x n_turbines).
+        x_points (np.ndarray): x-coordinates to be warped (n_findex x n_turbines x n_grid x n_grid).
+        y_points (np.ndarray): y-coordinates to be warped (n_findex x n_turbines x n_grid x n_grid).
+        z_points (np.ndarray): z-coordinates to be warped (n_findex x n_turbines x n_grid x n_grid).
+        wind_direction_x_points (np.ndarray): x-coordinates of the specified flow points
+            (n_findex x n_points).
+        wind_direction_y_points (np.ndarray): y-coordinates of the specified flow points
+            (n_findex x n_points).
+        wind_direction_u_values (np.ndarray): u-velocity values at the specified flow points
+            (n_findex x n_points).
+        wind_direction_v_values (np.ndarray): v-velocity values at the specified flow points
+            (n_findex x n_points).
+        wind_direction_w_values (np.ndarray): w-velocity values at the specified flow points
+            (n_findex x n_points).
+    """
+
+    # TODO: create switch that uses 3D versions, possible, based on whether
+    # wind_direction_z_points is None
+    n_findex, n_turbines = x_turbines.shape
+
+    # Initialize storage
+    x_points_per_turbine = np.repeat(x_points[:,:,:,:,None], n_turbines, axis=4)
+    y_points_per_turbine = np.repeat(y_points[:,:,:,:,None], n_turbines, axis=4)
+    z_points_per_turbine = np.repeat(z_points[:,:,:,:,None], n_turbines, axis=4)
+
+    # Ensure all needed elements are arrays
+    wind_direction_u_values = np.array(wind_direction_u_values)
+    wind_direction_v_values = np.array(wind_direction_v_values)
+
+    # Compute new relative locations
+    # TODO: Can I avoid any of these looping operations? Do I need to?
+    for findex in range(n_findex):
+        for tindex in range(n_turbines):
+            print("computing streamline for tindex", tindex)
+            # 1. Compute the streamline
+            streamline_start = np.array([
+                x_turbines[findex,tindex],
+                y_turbines[findex,tindex]
+            ])
+            streamline = compute_2D_streamline(
+                wind_direction_x_points[findex,:],
+                wind_direction_y_points[findex,:],
+                wind_direction_u_values[findex,:],
+                wind_direction_v_values[findex,:],
+                streamline_start
+            )
+
+            # 2. Compute the point movements
+            x_shifted, y_shifted, z_shifted = shift_points_by_streamline(
+                streamline,
+                x_points[findex,:,:,:],
+                y_points[findex,:,:,:],
+                z_points[findex,:,:,:],
+            )
+
+            # 3. Apply to per_turbine values
+            x_points_per_turbine[findex,:,:,:,tindex] = x_shifted
+            y_points_per_turbine[findex,:,:,:,tindex] = y_shifted
+            z_points_per_turbine[findex,:,:,:,tindex] = z_shifted
+
+    # Return full set of locations
+    return x_points_per_turbine, y_points_per_turbine, z_points_per_turbine
+
+def compute_2D_streamline(x, y, u, v, xy_0, n_points=1000):
+    """
+    Compute 2D streamline (in x-y plane) starting at xy_0, according to (u,v)
+    velocity field.
+    """
+    s = np.linspace(0, 1, n_points) # parametric variable
+
+    scale = np.array([x.max() - x.min(), y.max() - y.min()])
+    offset = np.array([(x.max() + x.min())/2, (y.max() + y.min())/2])
+    # Need to compute the offset too
+
+    # Collect
+    xy = np.array([x, y]).T
+    uv = np.array([u,v]).T
+
+    # Normalize
+    xy_normalized = (xy - offset) / scale
+    xy_0_normalized = (xy_0 - offset) / scale
+    # TODO: Loop over findices! Maybe not necessary---can I compute all streamlines at once??
+    # Would it help? Currently assuming a single findex.
+
+    interp = LinearNDInterpolator(xy_normalized, uv, fill_value=0)
+
+    def velocity_field_interpolate_reg(s, xy_):
+        # if (xy_ < -0.5).any() or (xy_ > 0.5).any():
+        #     return np.array([0.1, 0.1])*np.sign(xy_)
+        # else:
+        return interp(xy_)[0]
+
+    ode_solution = solve_ivp(
+        velocity_field_interpolate_reg,
+        (0, 1),
+        xy_0_normalized,
+        t_eval=np.linspace(0,1,n_points)
+    )
+    # Handle the ODE solution
+    streamline = ode_solution.y.T
+
+    # Determine point of exit from specified domain
+    if (streamline < -0.5).any():
+        low_exit = np.argwhere(streamline < -0.5)[0,0]
+    else:
+        low_exit = n_points
+    if (streamline > 0.5).any():
+        high_exit = np.argwhere(streamline > 0.5)[0,0]
+    else:
+        high_exit = n_points
+    exit_index = np.minimum(low_exit, high_exit)
+    if exit_index == 0:
+        raise ValueError("Streamline could not be propagated inside the domain.")
+    if exit_index == n_points:
+        # TODO: recommend ensuring domain covers the buffer region
+        if (abs(streamline) > 0.8*0.5).any(): # 10% buffer given,
+            pass
+        else:
+            import ipdb; ipdb.set_trace()
+            raise ValueError("Streamline did not reach the edge of the domain.")
+
+    # Copy points after exit
+    #streamline[exit_index:] = streamline[exit_index-1]
+
+    # Denormalize
+    streamline_denormalized = streamline * scale + offset
+
+    return streamline_denormalized
+
+def shift_points_by_streamline(streamline, x, y, z):
+    ## I think should be dist by turbine only (not each rotor point); and then
+    # shift all rotor points by the same amount.
+    # TODO: How will this work over findices? Need to work that out.
+    xy = np.concatenate((x.mean(axis=(1,2))[None,:],y.mean(axis=(1,2))[None,:]), axis=0)
+
+    rotor_y = y - y.mean(axis=(1,2), keepdims=True) # TODO: still ok for viz solve?
+
+    # Storage for new points
+    x_new = x.copy()
+    y_new = y.copy()
+
+    disps_to_streamline = xy[None,:,:] - streamline[:, :,None]
+    dists_to_streamline = np.linalg.norm(disps_to_streamline, axis=1)
+    min_idx = np.argmin(dists_to_streamline, axis=0)
+    #points = streamline[min_idx, :]
+
+    # TODO: for viz/solve for points, there are many more points than turbines.
+    # It would be nice not to have a loop here.
+    for i in range(len(min_idx)):
+        streamline_diffs = np.diff(streamline[:min_idx[i], :], axis=0)
+        along_stream_dist = np.sum(np.linalg.norm(streamline_diffs, axis=1))
+        off_stream_dist = dists_to_streamline[min_idx[i], i]
+        off_stream_disp = disps_to_streamline[min_idx[i], :, i]
+        disp_x = off_stream_disp[0]
+        disp_y = off_stream_disp[1]
+
+        theta = np.arctan2(disp_y, disp_x)
+
+        x_new[i,:,:] = streamline[0,0] + along_stream_dist
+        # y_new[tindex,:,:] = streamline[0,1] + disp_x*np.sin(theta) + disp_y*np.cos(theta)\
+        # + rotor_y[tindex,:,:]
+        y_new[i,:,:] = streamline[0,1] + np.sign(theta)*off_stream_dist + rotor_y[i,:,:]
+        #y_new[tindex,:,:] = streamline[0,1] + off_stream_dist + rotor_y[tindex,:,:]
+        #if tindex == 1:
+    return x_new, y_new, z
+
+    #return x_new.reshape(x.shape), y_new.reshape(y.shape), z
+
+def shift_points_by_streamline2(streamline, x, y, z, rotor_points_only=True):
+    ## I think should be dist by turbine only (not each rotor point); and then
+    # shift all rotor points by the same amount.
+    # TODO: How will this work over findices? Need to work that out.
+
+    #xy = np.concatenate((x.mean(axis=(1,2))[None,:],y.mean(axis=(1,2))[None,:]), axis=0)
+    xy = np.concatenate((x.flatten()[None,:], y.flatten()[None,:]), axis=0)
+
+    #rotor_y = y - y.mean(axis=(1,2), keepdims=True) # TODO: still ok for viz solve?
+
+    # Storage for new points
+    # x_new = x.copy()
+    # y_new = y.copy()
+
+    disps_to_streamline = xy[None,:,:] - streamline[:, :,None]
+    # TODO: This line slow. What if norm is delayed?
+    dists_to_streamline = np.linalg.norm(disps_to_streamline, axis=1)
+    disps_x = disps_to_streamline[:, 0, :]
+    disps_y = disps_to_streamline[:, 1, :]
+    #dists_to_streamline = ne.evaluate("sqrt(disps_x ** 2 + disps_y ** 2)")
+    #min_idx = np.argmin(dists_to_streamline, axis=0)
+    dists_to_streamline_sq = disps_x ** 2 + disps_y ** 2
+    min_idx = np.argmin(dists_to_streamline_sq, axis=0)
+    #points = streamline[min_idx, :]
+
+    # TODO: for viz/solve for points, there are many more points than turbines.
+    # It would be nice not to have a loop here.
+    # x_new = np.zeros(len(min_idx))
+    # y_new = np.zeros(len(min_idx))
+
+    s_diffs = np.diff(streamline, axis=0)
+    as_dist_all = np.insert(np.cumsum(np.linalg.norm(s_diffs, axis=1)), 0, 0)
+    as_dists = as_dist_all[min_idx]
+
+    disps = disps_to_streamline[min_idx, :, np.arange(0, xy.shape[1], 1)]
+    theta = np.arctan2(disps[:,1], disps[:,0])
+    #os_dists = dists_to_streamline[min_idx, np.arange(0, xy.shape[1], 1)]
+    os_dists = np.sqrt(dists_to_streamline_sq[min_idx, np.arange(0, xy.shape[1], 1)])
+    # for i in range(len(min_idx)):
+    #     streamline_diffs = np.diff(streamline[:min_idx[i], :], axis=0)
+    #     along_stream_dist = np.sum(np.linalg.norm(streamline_diffs, axis=1))
+    #     off_stream_dist = dists_to_streamline[min_idx[i], i]
+    #     off_stream_disp = disps_to_streamline[min_idx[i], :, i]
+    #     disp_x = off_stream_disp[0]
+    #     disp_y = off_stream_disp[1]
+
+    #     theta = np.arctan2(disp_y, disp_x)
+
+    #     x_new[i] = streamline[0,0] + along_stream_dist
+    #     # y_new[tindex,:,:] = streamline[0,1] + disp_x*np.sin(theta) + disp_y*np.cos(theta)\
+    #     # + rotor_y[tindex,:,:]
+    #     y_new[i] = streamline[0,1] + np.sign(theta)*off_stream_dist
+    #     #y_new[tindex,:,:] = streamline[0,1] + off_stream_dist + rotor_y[tindex,:,:]
+    #     #if tindex == 1:
+    #     #    import ipdb; ipdb.set_trace()
+
+    x_new = (streamline[0,0] + as_dists).reshape(x.shape)
+    y_new = (streamline[0,1] + np.sign(theta)*os_dists).reshape(y.shape)
+
+    return x_new, y_new, z
+
+    #return x_new.reshape(x.shape), y_new.reshape(y.shape), z
 
 class Loader(yaml.SafeLoader):
 
@@ -439,3 +688,28 @@ def print_nested_dict(dictionary: Dict[str, Any], indent: int = 0) -> None:
             print_nested_dict(value, indent + 4)
         else:
             print(" " * (indent + 4) + str(value))
+
+def update_model_args(
+        model_args: Dict[str, Any],
+        x: np.ndarray,
+        y: np.ndarray,
+        z: np.ndarray
+    ):
+    """
+    Update the x_sorted, y_sorted, and z_sorted entries of the model_args dict.
+    """
+    # Check that new x, y, z are the same length as the old
+    if (("x" in model_args.keys() and x.shape != model_args["x"].shape) or \
+        ("y" in model_args.keys() and y.shape != model_args["y"].shape) or \
+        ("z" in model_args.keys() and z.shape != model_args["z"].shape)):
+        raise ValueError("Size of new x,y,z arrays must match existing ones.")
+
+    # Update the x, y, z values
+    if "x" in model_args.keys():
+        model_args["x"] = x
+    if "y" in model_args.keys():
+        model_args["y"] = y
+    if "z" in model_args.keys():
+        model_args["z"] = z
+
+    return model_args
diff --git a/tests/data/input_full.yaml b/tests/data/input_full.yaml
index 49c41273d..2ed76a01a 100644
--- a/tests/data/input_full.yaml
+++ b/tests/data/input_full.yaml
@@ -60,6 +60,12 @@ wake:
       ad: 0.0
       bd: 0.0
       kd: 0.05
+    empirical_gauss:
+      horizontal_deflection_gain_D: 3.0
+      vertical_deflection_gain_D: -1
+      deflection_rate: 22
+      mixing_gain_deflection: 0.0
+      yaw_added_mixing_gain: 0.0
 
   wake_velocity_parameters:
     cc:
@@ -81,6 +87,19 @@ wake:
     turboparkgauss:
       A: 0.04
       include_mirror_wake: True
+    turbopark:
+      A: 0.04
+      sigma_max_rel: 4.0
+    empirical_gauss:
+      wake_expansion_rates:
+      - 0.023
+      - 0.008
+      breakpoints_D:
+      - 10
+      sigma_0_D: 0.28
+      smoothing_length_D: 2.0
+      mixing_gain_velocity: 2.0
+
 
   wake_turbulence_parameters:
     crespo_hernandez:
@@ -88,3 +107,5 @@ wake:
       constant: 0.9
       ai: 0.83
       downstream: -0.25
+    wake_induced_mixing:
+      atmospheric_ti_gain: 0.0
diff --git a/tests/wind_direction_heterogeneity_unit_test.py b/tests/wind_direction_heterogeneity_unit_test.py
new file mode 100644
index 000000000..533e4e384
--- /dev/null
+++ b/tests/wind_direction_heterogeneity_unit_test.py
@@ -0,0 +1,296 @@
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from floris import FlorisModel
+from floris.heterogeneous_map import HeterogeneousMap
+
+
+TEST_DATA = Path(__file__).resolve().parent / "data"
+YAML_INPUT = TEST_DATA / "input_full.yaml"
+
+layout_x = np.array([0.0, 500.0])
+layout_y = np.array([0.0, 0.0])
+
+heterogeneous_wd_inflow_config_2D = {
+    "x": [-1000, -1000, 1000, 1000],
+    "y": [-1000, 1000, -1000, 1000],
+    "u": [[8.0, 8.0, 8.0, 8.0]],
+    "v": [[0.0, 0.0, 4.0, 4.0]],
+    "speed_multipliers": [[1.0, 1.0, 1.0, 1.0]] # Removes wind speed variations
+}
+
+def test_power():
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+
+    # First run without heterogeneity
+    fmodel.set(layout_x=layout_x, layout_y=layout_y)
+    fmodel.run()
+    P_without_het = fmodel.get_turbine_powers()[0,:]
+
+    # Add wind direction heterogeneity and rerun
+    fmodel.set(heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D)
+    fmodel.run()
+    P_with_het = fmodel.get_turbine_powers()[0,:]
+
+    # Confirm upstream power the same; downstream power increased
+    assert np.allclose(P_without_het[0], P_with_het[0])
+    assert P_with_het[1] > P_without_het[1]
+
+def test_symmetry():
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+
+    # Set layout and heterogeneity
+    fmodel.set(
+        layout_x=layout_x,
+        layout_y=layout_y,
+        heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D
+    )
+
+    # Run in this direction
+    fmodel.run()
+    P_positive_v = fmodel.get_turbine_powers()[0,:]
+
+    # Switch the sign of the inflow v component and rurun
+    fmodel.set(
+        heterogeneous_inflow_config={
+            "x": heterogeneous_wd_inflow_config_2D["x"],
+            "y": heterogeneous_wd_inflow_config_2D["y"],
+            "u": heterogeneous_wd_inflow_config_2D["u"],
+            "v": -np.array(heterogeneous_wd_inflow_config_2D["v"]),
+            "speed_multipliers": heterogeneous_wd_inflow_config_2D["speed_multipliers"]
+        }
+    )
+    fmodel.run()
+    P_negative_v = fmodel.get_turbine_powers()[0,:]
+
+    # Confirm symmetry
+    assert np.allclose(P_positive_v, P_negative_v)
+
+def test_input_types():
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+
+    # First run without heterogeneity
+    fmodel.set(
+        layout_x=layout_x,
+        layout_y=layout_y,
+        heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D
+    )
+    fmodel.run()
+    P_lists = fmodel.get_turbine_powers()
+
+    # Convert all elements of dictionary to arrays
+    for k in heterogeneous_wd_inflow_config_2D:
+        heterogeneous_wd_inflow_config_2D[k] = np.array(heterogeneous_wd_inflow_config_2D[k])
+    assert isinstance(heterogeneous_wd_inflow_config_2D["u"], np.ndarray)
+
+    fmodel.set(heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D)
+    fmodel.run()
+    P_arrays = fmodel.get_turbine_powers()
+
+    # Confirm upstream power the same
+    assert np.allclose(P_lists, P_arrays)
+
+def test_multiple_findices():
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+
+    # Set layout and heterogeneity
+    fmodel.set(
+        layout_x=layout_x,
+        layout_y=layout_y,
+        heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D
+    )
+    fmodel.run()
+    P_1findex = fmodel.get_turbine_powers()[0,:]
+
+    # Duplicate the inflow conditons
+    fmodel.set(
+        heterogeneous_inflow_config={
+            "x": heterogeneous_wd_inflow_config_2D["x"],
+            "y": heterogeneous_wd_inflow_config_2D["y"],
+            "u": np.repeat(heterogeneous_wd_inflow_config_2D["u"], 2, axis=0),
+            "v": np.repeat(heterogeneous_wd_inflow_config_2D["v"], 2, axis=0),
+            "speed_multipliers": np.repeat(
+                heterogeneous_wd_inflow_config_2D["speed_multipliers"], 2, axis=0
+            )
+        },
+        wind_directions = [270.0, 270.0],
+        wind_speeds=[8.0, 8.0],
+        turbulence_intensities=[0.06, 0.06]
+    )
+    fmodel.run()
+    P_2findex = fmodel.get_turbine_powers()
+
+    # Confirm the same results
+    assert np.allclose(P_1findex, P_2findex)
+
+def test_flipped_direction():
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+
+    # Set up a case with 2 findices, central turbines, aligned in
+    fmodel.set(
+        wind_directions=[270, 90],
+        wind_speeds=[8.0, 8.0],
+        turbulence_intensities=[0.06, 0.06],
+        layout_x=np.array([-250.0, 250.0]),
+        layout_y=np.array([0.0, 0.0]),
+        heterogeneous_inflow_config={
+            "x": heterogeneous_wd_inflow_config_2D["x"],
+            "y": heterogeneous_wd_inflow_config_2D["y"],
+            "u": np.concatenate(
+                (
+                    heterogeneous_wd_inflow_config_2D["u"],
+                    np.flip(heterogeneous_wd_inflow_config_2D["u"], axis=1)
+                ),
+                axis=0
+            ),
+            "v": np.concatenate(
+                (
+                    heterogeneous_wd_inflow_config_2D["v"],
+                    np.flip(heterogeneous_wd_inflow_config_2D["v"], axis=1)
+                ),
+                axis=0
+            ),
+            "speed_multipliers": np.repeat(
+                heterogeneous_wd_inflow_config_2D["speed_multipliers"], 2, axis=0
+            )
+        }
+    )
+    fmodel.run()
+    powers = fmodel.get_turbine_powers()
+    assert np.allclose(powers[0,:], np.flip(powers[:,1]))
+
+def test_rotational_symmetry():
+    layout_x_symmetry = np.array([-250.0, 250.0])
+    layout_y_symmetry = np.array([0.0, 0.0])
+
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+    fmodel.set(
+        layout_x=layout_x_symmetry,
+        layout_y=layout_y_symmetry,
+        heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D
+    )
+    fmodel.run()
+    P_base = fmodel.get_turbine_powers()
+
+    # Rotate to a series of new positions, and rotate layout also
+    wind_directions_test = np.linspace(0, 360, 10)
+
+    for wd in wind_directions_test:
+        het_wd_inflow_test = heterogeneous_wd_inflow_config_2D.copy()
+        het_wd_inflow_test["x"] = \
+            np.cos(np.deg2rad(270-wd))*np.array(heterogeneous_wd_inflow_config_2D["x"]) \
+            - np.sin(np.deg2rad(270-wd))*np.array(heterogeneous_wd_inflow_config_2D["y"])
+        het_wd_inflow_test["y"] = \
+            np.sin(np.deg2rad(270-wd))*np.array(heterogeneous_wd_inflow_config_2D["x"]) \
+            + np.cos(np.deg2rad(270-wd))*np.array(heterogeneous_wd_inflow_config_2D["y"])
+
+        layout_x_test = (
+            np.cos(np.deg2rad(270-wd))*layout_x_symmetry
+            - np.sin(np.deg2rad(270-wd))*layout_y_symmetry
+        )
+        layout_y_test = (
+            np.sin(np.deg2rad(270-wd))*layout_x_symmetry
+            + np.cos(np.deg2rad(270-wd))*layout_y_symmetry
+        )
+
+        fmodel.set(
+            wind_directions=[wd],
+            layout_x=layout_x_test,
+            layout_y=layout_y_test,
+            heterogeneous_inflow_config=het_wd_inflow_test
+        )
+        fmodel.run()
+        P_new = fmodel.get_turbine_powers()
+
+        assert np.allclose(P_base, P_new, rtol=1e-4)
+
+def test_wake_models():
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+    fmodel_dict = fmodel.core.as_dict()
+    velocity_models = fmodel_dict["wake"]["wake_velocity_parameters"].keys()
+
+    # Switch off secondary effects
+    fmodel_dict["wake"]["enable_secondary_steering"] = False
+    fmodel_dict["wake"]["enable_yaw_added_recovery"] = False
+    fmodel_dict["wake"]["enable_transverse_velocities"] = False
+
+    # Loop through velocity models
+    for m in velocity_models:
+        fmodel_dict["wake"]["model_strings"]["velocity_model"] = m
+        if m == "empirical_gauss":
+            fmodel_dict["wake"]["model_strings"]["turbulence_model"] = "wake_induced_mixing"
+            fmodel_dict["wake"]["model_strings"]["deflection_model"] = "empirical_gauss"
+        elif m == "jensen":
+            fmodel_dict["wake"]["model_strings"]["turbulence_model"] = "crespo_hernandez"
+            fmodel_dict["wake"]["model_strings"]["deflection_model"] = "jimenez"
+        else:
+            fmodel_dict["wake"]["model_strings"]["turbulence_model"] = "crespo_hernandez"
+            fmodel_dict["wake"]["model_strings"]["deflection_model"] = "gauss"
+
+        fmodel = FlorisModel(fmodel_dict)
+
+        # First run without heterogeneity
+        fmodel.set(layout_x=layout_x, layout_y=layout_y)
+        fmodel.run()
+        P_without_het = fmodel.get_turbine_powers()[0,:]
+
+        # Add wind direction heterogeneity and rerun
+        fmodel.set(heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D)
+        fmodel.run()
+        P_with_het = fmodel.get_turbine_powers()[0,:]
+
+        # Confirm upstream power the same; downstream power increased
+        assert np.allclose(P_without_het[0], P_with_het[0])
+        assert P_with_het[1] > P_without_het[1]
+
+def test_speed_multipliers_option():
+    # Check that setting speed heterogeneity by u, v, works as expected
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+    fmodel.set(
+        heterogeneous_inflow_config={
+            "x": [0, 0, 1, 1],
+            "y": [0, 1, 0, 1],
+            "u": np.ones((1, 4)),
+            "v": np.ones((1, 4)),
+        },
+        wind_shear=0.0
+    )
+    assert (fmodel.core.flow_field.u_initial_sorted == np.sqrt(2)).all()
+
+    # Look at more realistic case
+    heterogeneous_wd_inflow_config_test = {
+        "x": [-1000, -1000, 1000, 1000],
+        "y": [-1000, 1000, -1000, 1000],
+        "u": [[8.0, 8.0, 8.0, 8.0]],
+        "v": [[0.0, 0.0, 4.0, 4.0]],
+        "speed_multipliers": [[1.0, 1.0, 1.0, 1.0]]
+    }
+
+    # First run with speed_multipliers specified
+    fmodel.set(
+        layout_x=layout_x,
+        layout_y=layout_y,
+        heterogeneous_inflow_config=heterogeneous_wd_inflow_config_test
+    )
+    fmodel.run()
+    P_with_sm = fmodel.get_turbine_powers()[0,:]
+
+    # Same speed at second turbine
+    assert (fmodel.core.flow_field.u_initial_sorted[0,1,:,:]
+            == fmodel.core.flow_field.u_initial_sorted[0,0,:,:]).all()
+
+    # Set without specified speed multipliers and rerun
+    # u, v will result in higher speeds at the back turbine
+    del heterogeneous_wd_inflow_config_test["speed_multipliers"]
+    fmodel.set(heterogeneous_inflow_config=heterogeneous_wd_inflow_config_test)
+    fmodel.run()
+    P_without_sm = fmodel.get_turbine_powers()[0,:]
+
+    # Higher speeds at second turbine
+    assert (fmodel.core.flow_field.u_initial_sorted[0,1,:,:]
+            > fmodel.core.flow_field.u_initial_sorted[0,0,:,:]).all()
+
+    # Also higher power than case without
+    assert P_without_sm[1] > P_with_sm[1]