limit power

0.0

}

}

// Calculates power generated based on head and release values.

/// Performs interpolation on a table of (head, release, power) triplets.

pub fn interpolate_power_from_release(actual_release: f64, current_head: f64, 

    head_values: &[f64], release_values: &[f64], 

    power_values: &[f64]) -> f64 {

// Handle edge cases

if head_values.is_empty() || release_values.is_empty() || power_values.is_empty() {

return 0.0; // Default value if tables are empty

}

// If no release or negative release, no power generated

if actual_release <= 0.0 {

return 0.0;

}

// Find rows with similar head values

let head_tolerance = 1.0; // 1 meter tolerance

// Collect indices of rows with similar head values

let mut similar_head_indices = Vec::new();

for i in 0..head_values.len() {

if (head_values[i] - current_head).abs() <= head_tolerance {

similar_head_indices.push(i);

}

}

// If no similar head values found, find the closest head values

if similar_head_indices.is_empty() {

// Find the two closest head values (one below, one above)

let mut lower_head_idx = 0;

let mut upper_head_idx = 0;

let mut min_distance_lower = f64::INFINITY;

let mut min_distance_upper = f64::INFINITY;

for i in 0..head_values.len() {

let distance = head_values[i] - current_head;

if distance < 0.0 && distance.abs() < min_distance_lower {

min_distance_lower = distance.abs();

lower_head_idx = i;

} else if distance >= 0.0 && distance < min_distance_upper {

min_distance_upper = distance;

upper_head_idx = i;

}

}

similar_head_indices.push(lower_head_idx);

if upper_head_idx != lower_head_idx {

similar_head_indices.push(upper_head_idx);

}

}

// Find the best power value for each similar head

let mut head_power_pairs = Vec::new();

for &idx in &similar_head_indices {

let head = head_values[idx];

// Find indices where release is closest to actual_release

let release = release_values[idx];

// If release at this head/power combo is close to actual, use this power

if (release - actual_release).abs() < 0.1 {

head_power_pairs.push((head, power_values[idx]));

continue;

}

// Otherwise, find adjacent releases to interpolate between

let mut lower_release_idx = idx;

let mut upper_release_idx = idx;

// Search for adjacent release values in table

for i in 0..release_values.len() {

if head_values[i] == head {

if release_values[i] < actual_release && release_values[i] > release_values[lower_release_idx] {

lower_release_idx = i;

}

if release_values[i] > actual_release && release_values[i] < release_values[upper_release_idx] {

upper_release_idx = i;

}

}

}

// Interpolate power between release values

if lower_release_idx != upper_release_idx {

let lower_release = release_values[lower_release_idx];

let upper_release = release_values[upper_release_idx];

let lower_power = power_values[lower_release_idx];

let upper_power = power_values[upper_release_idx];

// Linear interpolation between power values

let release_factor = (actual_release - lower_release) / (upper_release - lower_release);

let power = lower_power + (upper_power - lower_power) * release_factor;

head_power_pairs.push((head, power));

} else {

// If we can't find a good interpolation, just use this point

head_power_pairs.push((head, power_values[idx]));

}

}

// Now interpolate between the power values at different heads

if head_power_pairs.len() > 1 {

// Sort by head

head_power_pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());

// Find bounding head values

let mut lower_idx = 0;

let mut upper_idx = 0;

for i in 0..head_power_pairs.len() {

if head_power_pairs[i].0 <= current_head {

lower_idx = i;

}

if head_power_pairs[i].0 >= current_head && upper_idx == 0 {

upper_idx = i;

}

}

let lower_head = head_power_pairs[lower_idx].0;

let upper_head = head_power_pairs[upper_idx].0;

let lower_power = head_power_pairs[lower_idx].1;

let upper_power = head_power_pairs[upper_idx].1;

// Avoid division by zero

if (upper_head - lower_head).abs() < 1e-10 {

return lower_power;

}

// Linear interpolation for power based on head

let head_factor = (current_head - lower_head) / (upper_head - lower_head);

let power = lower_power + (upper_power - lower_power) * head_factor;

power

} else if head_power_pairs.len() == 1 {

// Only one head value, use its power

head_power_pairs[0].1

} else {

// No valid head values found, return zero power

0.0

}

}

 No newline at end of file

mod helpers;

use helpers::interpolate_elevation;

use helpers::interpolate_release_for_power;

use helpers::interpolate_power_from_release;

#[derive(Debug)]

enum ObjectType {

    head: f64,

    release: f64,

    spill: f64,

    actual_power: f64

}

// Define results structures for each object type

    total_inflow: Vec<f64>,

    target_release: Vec<f64>,

    target_power: Vec<f64>,

    actual_power: Vec<f64>,

}

#[derive(FromPyObject, IntoPyObject)]

        .min(reservoir.max_release)

        .min(available_water);

    let potential_storage = available_water - release;

    // Calculate actual power based on the actual release and head

    let actual_power = interpolate_power_from_release(

        release,

        head,

        &reservoir.hrpt_head,

        &reservoir.hrpt_release,

        &reservoir.hrpt_power

    );

    let (storage, spill) = if potential_storage > reservoir.capacity {

        (reservoir.capacity, potential_storage - reservoir.capacity)

    } else {

        head,

        release,

        spill,

        actual_power,

    }

}

                    total_inflow: vec![0.0; n],

                    target_release: reservoir.target_release.clone(),

                    target_power: reservoir.target_power.clone(),

                    actual_power: vec![0.0; n],

                }));

            },

            ObjectType::River => {

                        current_results.release[t] = state.release;

                        current_results.spill[t] = state.spill;

                        current_results.total_inflow[t] = current_inflow;

                        current_results.actual_power[t] = state.actual_power;

                    }

                },

                ObjectType::River => {