Working wave enrgy calculation

src/data/forecast_cbulletin_wave_data_record.rs

@@ -8,6 +8,7 @@ use serde::{Serialize, Deserialize};
 use crate::dimensional_data::DimensionalData;
 use crate::location::Location;
 use crate::swell::{Swell, SwellProvider, SwellSummary};
+use crate::tools::waves::wave_energy;
 use crate::units::{Direction, UnitConvertible, Unit, UnitSystem};
 
 use super::parseable_data_record::{DataRecordParsingError, ParseableDataRecord};
@@ -182,7 +183,7 @@ impl ParseableDataRecord for ForecastCBulletinWaveRecord {
                 period,
                 Direction::from_degrees(degrees),
                 None,
-                None,
+                Some(wave_energy(wave_height, period)),
                 None,
             ));
         }

src/data/gfs_wave_grib_point_data_record.rs

@@ -5,11 +5,7 @@ use gribberish::{message::Message, templates::product::tables::FixedSurfaceType}
 use serde::{Deserialize, Serialize};
 
 use crate::{
-    dimensional_data::DimensionalData,
-    location::Location,
-    model::NOAAModel,
-    swell::Swell,
-    units::{Direction, Unit, UnitConvertible, UnitSystem},
+    dimensional_data::DimensionalData, location::Location, model::NOAAModel, swell::Swell, tools::waves::wave_energy, units::{Direction, Unit, UnitConvertible, UnitSystem}
 };
 
 use super::parseable_data_record::DataRecordParsingError;
@@ -83,7 +79,7 @@ impl GFSWaveGribPointDataRecord {
             *period,
             Direction::from_degrees(*wave_direction as i32),
             None,
-            None,
+            Some(wave_energy(*wave_height, *period)),
             None,
         );
 
@@ -119,7 +115,7 @@ impl GFSWaveGribPointDataRecord {
                     data[&per_key],
                     Direction::from_degrees(data[&dir_key] as i32),
                     None,
-                    None,
+                    Some(wave_energy(data[&ht_key], data[&per_key])),
                     None,
                 );
 
@@ -140,7 +136,7 @@ impl GFSWaveGribPointDataRecord {
                 data["WVPER"],
                 Direction::from_degrees(data["WVDIR"] as i32),
                 None,
-                None,
+                Some(wave_energy(data["WVHGT"], data["WVPER"])),
                 None,
             );
 

src/swell.rs

@@ -15,6 +15,8 @@ pub struct Swell {
     pub partition: Option<usize>,
 }
 
+//      let energy = (1029.0 * ((9.81f64).powf(2.0))/(16.0 * PI)) * hs.powf(2.0) * peak_period.powf(2.0) / 1000.0;
+
 impl Swell {
     pub fn new(
         units: &UnitSystem,

src/tools/waves.rs

@@ -78,7 +78,7 @@ pub fn wavenu3(si: f64, h: f64) -> (f64, f64) {
     (k, cg)
 }
 
-// Chen and Thomson wavenumber approximation.
+/// Chen and Thomson wavenumber approximation.
 pub fn wavenuma(angle_freq: f64, water_depth: f64) -> f64 {
     let koh = 0.10194 * angle_freq * angle_freq * water_depth;
     let d = [0.0, 0.6522, 0.4622, 0.0, 0.0864, 0.0675];
@@ -90,8 +90,8 @@ pub fn wavenuma(angle_freq: f64, water_depth: f64) -> f64 {
     (koh * (1.0 + 1.0 / (koh * a)).sqrt()) / water_depth
 }
 
-// Wave celerity C
-// When depth is not supplied use deep water approximation
+/// Wave celerity C
+/// When depth is not supplied use deep water approximation
 pub fn celerity(freq: f64, depth: Option<f64>) -> f64 {
     if let Some(depth) = depth {
         let angle_freq = 2.0 * PI * freq;
@@ -101,8 +101,8 @@ pub fn celerity(freq: f64, depth: Option<f64>) -> f64 {
     }
 }
 
-// Wavelength L
-// When depth is not suppliked use deep water approximation
+/// Wavelength L
+/// When depth is not suppliked use deep water approximation
 pub fn wavelength(freq: f64, depth: Option<f64>) -> f64 {
     if let Some(depth) = depth {
         let angle_freq = 2.0 * PI * freq;
@@ -112,6 +112,11 @@ pub fn wavelength(freq: f64, depth: Option<f64>) -> f64 {
     }
 }
 
+/// Computes the wave energy for a given wave height and period. Units are metric, gravity is 9.81 m/s.
+pub fn wave_energy(hs: f64, tp: f64) -> f64 {
+    (1029.0 * ((9.81f64).powf(2.0)) / (16.0 * PI)) * hs.powf(2.0) * tp.powf(2.0) / 1000.0
+}
+
 /// Computes an estimate of the wave height for a given swell and beach conditions.
 pub fn estimate_breaking_wave_height(
     swell: &Swell,
@@ -480,7 +485,7 @@ pub fn pt_mean(
             }
         }
 
-        let energy = (1029.0 * (9.8f64).powf(2.0))/(4.0 * PI) * sume[ip] * peak_period.powf(2.0);
+        let energy = wave_energy(hs, peak_period);
 
         // let wind_sea_fraction = sumew[ip] / sume[ip];
 

src/units/mod.rs

@@ -292,6 +292,14 @@ impl UnitSystem {
             _ => 0.0,
         }
     }
+
+    pub fn density_of_seawater(&self) -> f64 {
+        match self {
+            UnitSystem::Metric => 1029.0,
+            UnitSystem::English => 64.0,
+            _ => 0.0,
+        }
+    }
 }
 
 impl Display for UnitSystem {

tests/buoy_data_tests.rs

@@ -137,7 +137,9 @@ fn read_wave_spectra_data() {
             component.to_units(&UnitSystem::English);
         });
 
-    let cart_e = record.spectra.project_cartesian(&record.spectra.energy, 50, Some(25.0), None);
+    let cart_e = record
+        .spectra
+        .project_cartesian(&record.spectra.energy, 50, Some(25.0), None);
     let (min_e, max_e) = record.spectra.energy_range();
     let _binned_cart_e = bin(&cart_e, &min_e, &max_e, &255);
 }
@@ -158,10 +160,21 @@ fn read_cbulletin_forecast_station_data() {
     }
 
     // Verify a random timestep
-    let swell_data = bulletin_records[7].swell_data().clone().unwrap().to_units(&UnitSystem::English).clone();
+    let swell_data = bulletin_records[7]
+        .swell_data()
+        .clone()
+        .unwrap()
+        .to_units(&UnitSystem::English)
+        .clone();
     assert_eq!(swell_data.summary.wave_height.get_value().ceil() as i32, 4);
-    assert_eq!(swell_data.components[0].wave_height.get_value().ceil() as i32, 3);
-    assert_eq!(swell_data.components[0].period.get_value().ceil() as i32, 13);
+    assert_eq!(
+        swell_data.components[0].wave_height.get_value().ceil() as i32,
+        3
+    );
+    assert_eq!(
+        swell_data.components[0].period.get_value().ceil() as i32,
+        13
+    );
 
     let raw_data = read_mock_data("tpc55.cbull");
     let mut data_collection = ForecastCBulletinWaveRecordCollection::from_data(raw_data.as_str());
@@ -177,10 +190,7 @@ fn read_spectral_forecast_station_data() {
     let spectral_records_iter = data_collection.records();
     assert!(spectral_records_iter.is_ok());
 
-    let spectral_records = spectral_records_iter
-        .unwrap()
-        .1
-        .collect::<Vec<_>>();
+    let spectral_records = spectral_records_iter.unwrap().1.collect::<Vec<_>>();
 
     assert_eq!(spectral_records[0].reference_date, spectral_records[0].date);
     assert_eq!(spectral_records[3].reference_date, spectral_records[0].date);