ndspec 0.1.1

A package for ocean waves and energy spectra
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

//! The Ochi-Hubble energy density spectru:
//!
//!
//!
use crate::prelude::*;
use ndarray::{Array1, ArrayView1};
use special::Gamma;

use crate::waves::common::SpectralCommon;

const N_FREQ: usize = 256;

pub struct OchiHubble {
    pub hs: [f64; 2],
    pub tp: [f64; 2],
    pub lambda: [f64; 2],
    pub omega: Array1<f64>,
}

pub fn energy(omega: ArrayView1<f64>, hs: [f64; 2], tp: [f64; 2], lambda: [f64; 2]) -> Array1<f64> {
    let a = |w: f64, l: f64| w.powf(-4. * (l + 0.25)) * (-w.powf(-4.) * (l + 0.25)).exp();
    let b = |l: f64| 4. * (l + 0.25).powf(l) / Gamma::gamma(l);
    let c = |h: f64, t: f64| h.powf(2.) * t / (32. * PI);

    let energy: Array1<f64> = omega
        .iter()
        .map(|w_i| {
            let w_n = [w_i * tp[0] / TWO_PI, w_i * tp[1] / TWO_PI];
            let mut sum = 0.0;
            for i in 0..2 {
                sum += a(w_n[i], lambda[i]) * b(lambda[i]) * c(hs[i], tp[i]);
            }
            sum
        })
        .collect();

    energy
}

impl Default for OchiHubble {
    fn default() -> Self {
        Self {
            hs: [1.0, 0.5],
            tp: [20.0, 8.0],
            lambda: [0.5, 0.5],
            omega: Array1::<f64>::linspace(0.1, PI, N_FREQ),
        }
    }
}

impl SpectralCommon for OchiHubble {
    /// rms significant wave height approximation
    fn hs(&self) -> Result<f64, String> {
        match self.hs.len() {
            2 => Ok((self.hs[0].powf(2.) + self.hs[1].powf(2.)).sqrt()),
            _ => Err("hs is not length 2".to_string()),
        }
    }

    /// rms peak spectral wave period approximation
    fn tp(&self) -> Result<f64, String> {
        match self.tp.len() {
            2 => Ok((self.tp[0].powf(2.) + self.tp[1].powf(2.)).sqrt()),
            _ => Err("tp is not length 2".to_string()),
        }
    }

    fn omega(&self) -> &Array1<f64> {
        &self.omega
    }

    fn energy(&self) -> Array1<f64> {
        energy(self.omega.view(), self.hs, self.tp, self.lambda)
    }
}

impl OchiHubble {
    /// Create a new Ochi-Hubble spectrum
    pub fn new(hs: [f64; 2], tp: [f64; 2], lambda: [f64; 2]) -> Self {
        Self {
            hs,
            tp,
            lambda,
            ..Default::default()
        }
    }

    /// Set the significant wave height for each component
    pub fn set_hs(&mut self, hs: [f64; 2]) -> &mut Self {
        self.hs = hs;
        self
    }

    /// Set the peak spectral wave period for each component
    pub fn set_tp(&mut self, tp: [f64; 2]) -> &mut Self {
        self.tp = tp;
        self
    }

    /// Set the shape factor for each component
    pub fn set_lambda(&mut self, lambda: [f64; 2]) -> &mut Self {
        self.lambda = lambda;
        self
    }

    /// Return the peak spectral wave frequency for each component
    pub fn get_wp(&self) -> Result<[f64; 2], String> {
        match self.tp.len() {
            2 => Ok([self.tp[0] / TWO_PI, self.tp[1] / TWO_PI]),
            _ => Err("tp is not length 2".to_string()),
        }
    }
}

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

    #[test]
    fn test_ochi_hubble() {
        let ochi = OchiHubble::default();
        let energy = ochi.energy();
        assert_eq!(energy.len(), N_FREQ);
    }

    #[test]
    fn test_ochi_hubble_hs() {
        let ochi = OchiHubble::default();
        let hs = ochi.hs().unwrap();

        assert_eq!(ochi.hs[0], 1.0_f64);
        assert_eq!(ochi.hs[1], 0.5_f64);
        assert_eq!(hs, (1.25_f64).sqrt());
    }

    #[test]
    fn test_ochi_hubble_tp() {
        let ochi = OchiHubble::default();
        let tp = ochi.tp().unwrap();

        assert_eq!(ochi.tp[0], 20.0_f64);
        assert_eq!(ochi.tp[1], 8.0_f64);
        assert_eq!(tp, (20.0_f64.powf(2.) + 8.0_f64.powf(2.)).sqrt());
    }
}