oxigrid 0.1.2

/// Wind farm wake models.
///
/// Provides:
/// - **Jensen (1983) top-hat wake model** — linear cone expansion, widely
///   used for layout optimization and energy yield assessment.
/// - **Frandsen (2006) model** — improved wake width model.
use serde::{Deserialize, Serialize};

/// Position of a turbine in a wind farm layout.
#[derive(Debug, Clone, Copy, Serialize, Deserialize)]
pub struct TurbinePosition {
    /// Downstream distance from reference origin `m`
    pub x: f64,
    /// Crosswind distance `m`
    pub y: f64,
}

/// Wake deficit parameters for a single upstream turbine.
#[derive(Debug, Clone, Copy, Serialize, Deserialize)]
pub struct WakeSource {
    /// Hub-height free-stream wind speed [m/s]
    pub u_inf: f64,
    /// Rotor diameter `m`
    pub d: f64,
    /// Thrust coefficient (−)
    pub ct: f64,
    /// Wake decay constant (Jensen: 0.04 offshore, 0.075 onshore)
    pub k: f64,
}

impl WakeSource {
    /// Compute wake wind speed at point (dx, dy) relative to source.
    ///
    /// Uses the Jensen top-hat model.
    /// Returns `u_inf` if outside the wake cone.
    pub fn jensen_wake_speed(&self, dx: f64, dy: f64) -> f64 {
        if dx <= 0.0 {
            return self.u_inf;
        }

        let r0 = self.d / 2.0;
        let r_wake = r0 + self.k * dx;

        // Wake deficit factor for a point fully inside the wake
        let deficit = (1.0 - (1.0 - self.ct).sqrt()) * (r0 / r_wake).powi(2);

        // Check overlap: is the evaluation point inside the wake cone?
        if dy.abs() <= r_wake {
            // Partial overlap area fraction (top-hat: binary inside/outside)
            let overlap = overlap_fraction(dy.abs(), r_wake, r0);
            self.u_inf * (1.0 - deficit * overlap)
        } else {
            self.u_inf
        }
    }

    /// Compute wake wind speed using the Frandsen (2006) improved model.
    pub fn frandsen_wake_speed(&self, dx: f64, dy: f64) -> f64 {
        if dx <= 0.0 {
            return self.u_inf;
        }

        let r0 = self.d / 2.0;
        // Frandsen wake expansion coefficient
        let beta = 0.5 * (1.0 + (1.0 - self.ct).sqrt()) / (1.0 - self.ct).sqrt();
        let r_wake = r0 * (beta + self.k * dx / r0).sqrt();

        if dy.abs() <= r_wake {
            let deficit = 0.5 * (1.0 - (1.0 - self.ct * (r0 / r_wake).powi(2)).sqrt());
            self.u_inf * (1.0 - 2.0 * deficit)
        } else {
            self.u_inf
        }
    }
}

/// Geometric overlap fraction (simplified linear approximation).
fn overlap_fraction(dy: f64, r_wake: f64, _r_rotor: f64) -> f64 {
    if dy > r_wake {
        0.0
    } else {
        1.0
    }
}

/// Compute effective wind speed at each turbine in a farm using the
/// square-sum wake superposition principle (Lissaman 1979).
///
/// `positions` — turbine positions
/// `wind_dir_deg` — wind direction (270° = westerly, blowing east)
/// `sources` — WakeSource parameters for each turbine
pub fn farm_wake_speeds(
    positions: &[TurbinePosition],
    wind_dir_deg: f64,
    sources: &[WakeSource],
) -> Vec<f64> {
    assert_eq!(positions.len(), sources.len());
    let n = positions.len();
    let u_inf = sources[0].u_inf;

    // Rotate coordinates to wind-aligned frame.
    // Meteorological convention: wind_dir_deg is the direction the wind comes FROM,
    // clockwise from north.  Downwind unit vector (x=east, y=north):
    //   u_dw = (-sin(θ), -cos(θ))
    // x_w = downstream distance (positive = further downwind)
    // y_w = crosswind distance
    let dir_rad = wind_dir_deg.to_radians();
    let (sin_d, cos_d) = (dir_rad.sin(), dir_rad.cos());

    let rotate = |p: &TurbinePosition| -> (f64, f64) {
        let x_w = -p.x * sin_d - p.y * cos_d;
        let y_w = -p.x * cos_d + p.y * sin_d;
        (x_w, y_w)
    };

    let rotated: Vec<(f64, f64)> = positions.iter().map(rotate).collect();

    let mut speeds = vec![u_inf; n];

    for i in 0..n {
        let mut deficit_sq_sum = 0.0_f64;
        for j in 0..n {
            if i == j {
                continue;
            }
            let dx = rotated[i].0 - rotated[j].0;
            let dy = rotated[i].1 - rotated[j].1;
            if dx <= 0.0 {
                continue;
            } // j is not upstream of i
            let u_wake = sources[j].jensen_wake_speed(dx, dy);
            let d_rel = (u_inf - u_wake) / u_inf;
            deficit_sq_sum += d_rel * d_rel;
        }
        speeds[i] = u_inf * (1.0 - deficit_sq_sum.sqrt());
    }

    speeds
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_jensen_no_deficit_upstream() {
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        let u = src.jensen_wake_speed(-100.0, 0.0);
        assert_eq!(u, 10.0);
    }

    #[test]
    fn test_jensen_deficit_on_centreline() {
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        let u = src.jensen_wake_speed(500.0, 0.0);
        assert!(u < 10.0 && u > 0.0, "u={:.2}", u);
    }

    #[test]
    fn test_jensen_deficit_outside_cone() {
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        // Very far crosswind — outside wake cone
        let u = src.jensen_wake_speed(500.0, 1000.0);
        assert_eq!(u, 10.0);
    }

    #[test]
    fn test_deficit_decreases_with_distance() {
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        let u1 = src.jensen_wake_speed(300.0, 0.0);
        let u2 = src.jensen_wake_speed(600.0, 0.0);
        // Deficit should recover with distance
        assert!(u2 > u1, "u1={:.2} u2={:.2}", u1, u2);
    }

    #[test]
    fn test_farm_two_turbines_in_line() {
        let positions = vec![
            TurbinePosition { x: 0.0, y: 0.0 },
            TurbinePosition { x: 500.0, y: 0.0 },
        ];
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        let sources = vec![src; 2];
        let speeds = farm_wake_speeds(&positions, 270.0, &sources); // westerly wind
                                                                    // Second turbine should be in wake of first
        assert!(speeds[1] < speeds[0], "speeds={:?}", speeds);
    }

    #[test]
    fn test_frandsen_upstream_returns_u_inf() {
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        let u = src.frandsen_wake_speed(-100.0, 0.0);
        assert_eq!(u, 10.0, "Frandsen upstream must return u_inf, got {u:.4}");
    }

    #[test]
    fn test_frandsen_deficit_on_centreline() {
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        let u = src.frandsen_wake_speed(500.0, 0.0);
        assert!(
            u < 10.0 && u > 0.0,
            "Frandsen centreline speed must be in (0, u_inf), got {u:.4}"
        );
    }

    #[test]
    fn test_frandsen_outside_wake_returns_u_inf() {
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        // 5 km crosswind — far outside any realistic wake cone
        let u = src.frandsen_wake_speed(500.0, 5000.0);
        assert_eq!(
            u, 10.0,
            "Frandsen far crosswind must return u_inf, got {u:.4}"
        );
    }

    #[test]
    fn test_jensen_zero_distance_gives_deficit() {
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        // Just a tiny step downstream: the rotor is effectively at the source,
        // so the evaluation point is inside the wake cone at dx → 0+.
        let u = src.jensen_wake_speed(0.001, 0.0);
        assert!(
            u <= 10.0,
            "Jensen at dx≈0 on centreline must be ≤ u_inf, got {u:.4}"
        );
    }

    #[test]
    fn test_farm_single_turbine() {
        let positions = vec![TurbinePosition { x: 0.0, y: 0.0 }];
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        let sources = vec![src];
        let speeds = farm_wake_speeds(&positions, 270.0, &sources);
        assert_eq!(speeds.len(), 1, "Expected exactly one speed value");
        assert_eq!(
            speeds[0], 10.0,
            "Single turbine must see free-stream speed, got {:.4}",
            speeds[0]
        );
    }

    #[test]
    fn test_farm_crosswind_no_wake() {
        // Two turbines side-by-side perpendicular to a westerly (270°) wind.
        // With wind blowing east (+x direction), turbines placed at large ±y
        // separation are purely crosswind — neither is in the other's wake.
        let positions = vec![
            TurbinePosition { x: 0.0, y: 5000.0 },
            TurbinePosition { x: 0.0, y: -5000.0 },
        ];
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        let sources = vec![src; 2];
        let speeds = farm_wake_speeds(&positions, 270.0, &sources);
        assert_eq!(
            speeds[0], 10.0,
            "Crosswind turbine 0 must see u_inf, got {:.4}",
            speeds[0]
        );
        assert_eq!(
            speeds[1], 10.0,
            "Crosswind turbine 1 must see u_inf, got {:.4}",
            speeds[1]
        );
    }

    #[test]
    fn test_farm_three_turbines_in_line_speed_decreasing() {
        // Three turbines in a row along the x-axis under a westerly (270°) wind.
        // Wind blows from west to east (+x), so turbine at x=0 is most upwind,
        // x=500 is mid, and x=1000 is most downwind.
        let positions = vec![
            TurbinePosition { x: 0.0, y: 0.0 },
            TurbinePosition { x: 500.0, y: 0.0 },
            TurbinePosition { x: 1000.0, y: 0.0 },
        ];
        let src = WakeSource {
            u_inf: 10.0,
            d: 90.0,
            ct: 0.8,
            k: 0.04,
        };
        let sources = vec![src; 3];
        let speeds = farm_wake_speeds(&positions, 270.0, &sources);
        assert_eq!(speeds.len(), 3, "Expected three speed values");
        assert!(
            speeds[0] > speeds[2],
            "Most upwind turbine must see higher speed than most downwind; speeds={speeds:?}"
        );
    }
}