stormbird 0.9.0

// Copyright (C) 2024, NTNU
// Author: Jarle Vinje Kramer <jarlekramer@gmail.com; jarle.a.kramer@ntnu.no>
// License: GPL v3.0 (see separate file LICENSE or https://www.gnu.org/licenses/gpl-3.0.html)

//! Functionality to represent the constant variables in a wind environment, and methods necessary 
//! to compute velocities in spaces.s

use std::ops::Range;

use stormath::{
    type_aliases::Float,
    spatial_vector::SpatialVector
};
use serde::{Serialize, Deserialize};
use serde_json;

use crate::error::Error;
use crate::line_force_model::LineForceModel;
use super::inflow_corrections::InflowCorrections;
use super::wind_condition::WindCondition;

#[derive(Debug, Clone, Serialize, Deserialize)]
#[serde(deny_unknown_fields)]
/// Structure used to represent a wind environment, meaning all settings for the wind that is 
/// assumed to be constant for a given model. Wind parameters that vary is given in 
/// [`crate::wind::wind_condition::WindCondition`]. The wind environment defines the definition of 
/// different *directions*, and can be used to generate velocity vectors for different locations, 
/// based on the input wind condition and linear velocity
pub struct WindEnvironment {
    #[serde(default="WindEnvironment::default_up_direction")]
    /// The up direction specifies which direction is considered to be up relative to the ocean or 
    /// ground. This vector is predominantly used to compute the height for any given point, which 
    /// are later used to get the true wind speed.  
    pub up_direction: SpatialVector,
    #[serde(default="WindEnvironment::default_wind_rotation_axis")]
    /// The wind rotation axis which axis that is used when rotating the wind vector based on the 
    /// wind direction. This is kept separate from the *up-direction* to allow for arbitrary 
    /// definitions of wind direction
    pub wind_rotation_axis: SpatialVector,
    #[serde(default="WindEnvironment::default_zero_direction_vector")]
    /// The zero direction vector is the direction at which the true wind will point when the wind
    /// direction is zero. 
    pub zero_direction_vector: SpatialVector,
    #[serde(default)]
    /// The height of the water plane. Can be used to change the height that is used as input to 
    /// a model for the atmospheric boundary, for instance to allow the sails to be defined in a 
    /// different coordinate system than using the water plane as z=0.0
    pub water_plane_height: Float,
    #[serde(default)]
    /// Optional empirical inflow corrections to use when calculating the effective inflow to each 
    /// sail
    pub inflow_corrections: Option<InflowCorrections>,
}

impl Default for WindEnvironment {
    fn default() -> Self {
        Self {
            up_direction: Self::default_up_direction(),
            wind_rotation_axis: Self::default_wind_rotation_axis(),
            zero_direction_vector: Self::default_zero_direction_vector(),
            water_plane_height: 0.0,
            inflow_corrections: None
        }
    }
}

impl WindEnvironment {
    fn default_zero_direction_vector() -> SpatialVector {SpatialVector::from([1.0, 0.0, 0.0])}
    fn default_up_direction() -> SpatialVector {SpatialVector::from([0.0, 0.0, 1.0])}
    fn default_wind_rotation_axis() -> SpatialVector {SpatialVector::from([0.0, 0.0, 1.0])}

    /// Creates a wind environment from a JSON string
    pub fn from_json_string(json_string: &str) -> Result<Self, Error> {
        let serde_res = serde_json::from_str(json_string)?;

        Ok(serde_res)
    }

    /// Creates a wind environment from a JSON file
    pub fn from_json_file(file_path: &str) -> Result<Self, Error> {
        let json_string = std::fs::read_to_string(file_path)?;

        Self::from_json_string(&json_string)
    }
    
    #[inline(always)]
    /// Returns the height from a given vector, based on the definitions on the up-direction and 
    /// water plane height given in the structure.
    pub fn height_from_location(&self, location: SpatialVector) -> Float {
        (
            location.dot(self.up_direction) - self.water_plane_height
        ).max(0.0)
    }
    
    #[inline(always)]
    /// Returns the true wind direction as a vector
    pub fn true_wind_direction_vector(&self, direction_coming_from: Float) -> SpatialVector {
        self.zero_direction_vector.rotate_around_axis(
            direction_coming_from,
            self.wind_rotation_axis
        )
    }

    /// Returns the steady part of the true wind conditions at the specified location
    pub fn steady_true_wind_velocity_at_location(
        &self,
        condition: &WindCondition,
        location: SpatialVector,
    ) -> Float {
        let height = self.height_from_location(location);

        condition.steady_true_wind_velocity_at_height(height)
    }

    /// Computes unsteady true wind that is parallel to the steady wind. 
    pub fn unsteady_parallel_true_wind_velocity_at_location(
        &self,
        condition: &WindCondition,
        location: SpatialVector,
        time: Float
    ) -> Float {
        let height = self.height_from_location(location);

        condition.unsteady_parallel_true_wind_velocity_at_height(height, time)
    }

    /// Returns the unsteady true wind magnitude at the location given as input
    pub fn unsteady_true_wind_velocity_vector_at_location(
        &self,
        condition: &WindCondition,
        location: SpatialVector,
        time: Float
    ) -> SpatialVector {
        let parallel_vel = self.unsteady_parallel_true_wind_velocity_at_location(
            condition, 
            location, 
            time
        );
        let parallel_dir = self.true_wind_direction_vector(condition.direction_coming_from);
        
        let perpendicular_vel = condition.unsteady_perpendicular_true_wind_velocity(time);
        let perpendicular_dir = self.up_direction.cross(parallel_dir);
        
        let vertical_dir = self.up_direction;
        let vertical_vel = condition.unsteady_vertical_true_wind_velocity(time);

        parallel_vel * parallel_dir + 
        perpendicular_vel * perpendicular_dir + 
        vertical_vel * vertical_dir
    }
    
    /// Returns the true wind vector at the location given as input
    pub fn steady_true_wind_velocity_vector_at_location(
        &self,
        condition: &WindCondition,
        location: SpatialVector,
    ) -> SpatialVector {
        let parallel_vel = self.steady_true_wind_velocity_at_location(condition, location);
        let parallel_dir = self.true_wind_direction_vector(condition.direction_coming_from);

        parallel_vel * parallel_dir
    }

    /// Returns the steady apparent wind vector at the supplied location
    pub fn steady_apparent_wind_velocity_vector_at_location(
        &self,
        condition: &WindCondition,
        location: SpatialVector,
        linear_velocity: SpatialVector,
    ) -> SpatialVector {
        let true_wind = self.steady_true_wind_velocity_vector_at_location(
            condition, 
            location
        );
        
        true_wind + linear_velocity
    }

    /// Returns the unsteady apparent wind vector at the supplied location
    pub fn unsteady_apparent_wind_velocity_vector_at_location(
        &self,
        condition: &WindCondition,
        location: SpatialVector,
        linear_velocity: SpatialVector,
        time: Float
    ) -> SpatialVector {
        let true_wind = self.unsteady_true_wind_velocity_vector_at_location(
            condition, 
            location, 
            time
        );
        
        true_wind + linear_velocity
    }

    
    /// Computes the apparent wind velocity with corrections at the control points. 
    pub fn apparent_wind_velocity_vectors_at_ctrl_points_with_corrections_applied(
        &self,
        condition: &WindCondition,
        ctrl_points: &[SpatialVector],
        linear_velocity: SpatialVector,
        time: Float,
        wing_indices: &[Range<usize>]
    ) -> Vec<SpatialVector> {
        
        let nr_ctrl_points = ctrl_points.len();
        
        let mut wind_velocity = Vec::with_capacity(nr_ctrl_points);
        let mut average_height = 0.0;
        
        for i in 0..nr_ctrl_points {
            wind_velocity.push(
                self.unsteady_apparent_wind_velocity_vector_at_location(
                    condition, 
                    ctrl_points[i], 
                    linear_velocity, 
                    time
                )
            );
            
            average_height += ctrl_points[i].dot(self.up_direction);
        }
        
        average_height /= nr_ctrl_points as Float;
        
        let apparent_wind_direction = self.apparent_wind_direction_from_condition_and_linear_velocity(
            condition, linear_velocity, average_height
        );
        
        self.apply_inflow_corrections(
            apparent_wind_direction,
            &mut wind_velocity,
            ctrl_points,
            wing_indices,
        );

        wind_velocity
    }

    /// Applies inflow corrections to the first points in the input freestream velocity
    pub fn apply_inflow_corrections(
        &self,
        apparent_wind_direction: Float,
        freestream_velocity: &mut [SpatialVector],
        ctrl_points: &[SpatialVector],
        wing_indices: &[Range<usize>]
    ) {
        if let Some(corrections) = &self.inflow_corrections {
            let nr_ctrl_points = ctrl_points.len();
            let nr_wings = wing_indices.len();
            
            let mut height_values: Vec<Float> = Vec::with_capacity(
                nr_ctrl_points
            );
            
            for i in 0..nr_ctrl_points {
                height_values.push(
                    ctrl_points[i].dot(self.up_direction)
                );
            }
            
            for wing_index in 0..nr_wings {
                for i in wing_indices[wing_index].start..wing_indices[wing_index].end {
                    freestream_velocity[i] = corrections.correct_velocity_single_sail(
                        wing_index,
                        apparent_wind_direction,
                        height_values[i],
                        freestream_velocity[i],
                        self.up_direction
                    )
                }
            }
        }
    }

    /// Computes the apparent wind direction from the supplied condition and linear velocity
    pub fn apparent_wind_direction_from_condition_and_linear_velocity(
        &self,
        condition: &WindCondition,
        linear_velocity: SpatialVector,
        height: Float
    ) -> Float {
        let location = height * self.up_direction;
        
        let true_wind_vector = self.steady_true_wind_velocity_vector_at_location(
            condition,
            location
        );
        
        let apparent_velocity_vector = true_wind_vector + linear_velocity;

        self.zero_direction_vector.signed_angle_between(
            apparent_velocity_vector,
            self.wind_rotation_axis
        )
    }

    /// Measures the apparent wind direction based on the input velocity vectors, where the sign and
    /// magnitude is defined by the zero_direction_vector and the wind_rotation_axis.
    pub fn apparent_wind_direction_from_velocity_based_on_rotation_axis(
        &self,
        velocity: SpatialVector
    ) -> Float {
        self.zero_direction_vector.signed_angle_between(
            velocity,
            self.wind_rotation_axis
        )
    }

    /// Measures the apparent wind direction based on the input velocity vectors, where the sign is
    /// defined by the local, non-rotated, chord vector and rotation-axis of each wing in the line
    /// force model. This, then, gives the wind direction relative to the local coordinate system
    /// for each wing. A direction of zero means that the flow is aligned with the non-rotated chord
    pub fn apparent_wind_direction_from_velocity_and_line_force_model(
        &self,
        velocity: &[SpatialVector],
        line_force_model: &LineForceModel
    ) -> Vec<Float> {

        let nr_span_lines = line_force_model.nr_span_lines();

        let mut out = Vec::with_capacity(nr_span_lines);

        for i in 0..nr_span_lines {
            let wing_index = line_force_model.wing_index_from_global(i);

            let first_strip_index = line_force_model.wing_indices[wing_index].start;

            let rotation_axis = line_force_model.span_lines_global[first_strip_index].relative_vector().normalize();

            let chord_local_non_transformed = line_force_model.chord_vectors_local_not_rotated[first_strip_index];

            let chord_local_transformed = line_force_model.rigid_body_motion.transform_vector(chord_local_non_transformed);

            out.push(
                chord_local_transformed.signed_angle_between(
                    velocity[i],
                    rotation_axis
                )
            );
        }

        out
    }
}

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

    #[test]
    /// Performs basic tests of the rotation logic. 
    /// 
    /// The "north wind" is assumed to be in the direction of the positive x-axis by default 
    /// (e.g., positive x axis points backwards)
    fn test_true_wind_velocity_vectors() {
        
        let wind_environment = WindEnvironment::default();
        let location = SpatialVector::new(0.0, 0.0, 10.0);

        let wind_velocity = 8.2;

        let head_wind_condition = WindCondition::new_constant(
            0.0, 
            wind_velocity
        );
        let starboard_wind_condition = WindCondition::new_constant(
            Float::from(-90.0).to_radians(), 
            wind_velocity
        );

        let port_wind_condition = WindCondition::new_constant(
            Float::from(90.0).to_radians(), 
            wind_velocity
        );

        let tail_wind_condition = WindCondition::new_constant(
            Float::from(180.0).to_radians(), 
            wind_velocity
        );

        let head_vector = wind_environment.steady_true_wind_velocity_vector_at_location(
            &head_wind_condition, location
        );

        let starboard_vector = wind_environment.steady_true_wind_velocity_vector_at_location(
            &starboard_wind_condition, location
        );

        let port_vector = wind_environment.steady_true_wind_velocity_vector_at_location(
            &port_wind_condition, location
        );

        let tail_vector = wind_environment.steady_true_wind_velocity_vector_at_location(
            &tail_wind_condition, location
        );

        dbg!(head_vector);
        dbg!(starboard_vector);
        dbg!(port_vector);
        dbg!(tail_vector);

        assert!(head_vector[0] > 0.0); // x-axis points backwards, so head wind should be positive in x
        assert!(tail_vector[0] < 0.0); // Opposite of above
        
        assert!(starboard_vector[1] < 0.0); // y-axis points towards right, so wind from starboard should be negative in y
        assert!(port_vector[1] > 0.0); // Opposite of above
    }
}