Update JS and PM to match IEC (#128)

* Update JS and PM to match IEC

* Fix bug and adjust tests.

* remove BS and update JS (#2)

* remove BS and update JS

* updated function reference to IEC TS62600-2

* replace BS calls with JS in tests

Co-authored-by: Kelley Ruehl <[email protected]>

Author: ssolson

mhkit/tests/test_wave.py

@@ -41,15 +41,18 @@ def tearDownClass(self):
         pass
 
     def test_pierson_moskowitz_spectrum(self):
-        S = wave.resource.pierson_moskowitz_spectrum(self.f,self.Tp)
+        S = wave.resource.pierson_moskowitz_spectrum(self.f,self.Tp,self.Hs)
+        Hm0 = wave.resource.significant_wave_height(S).iloc[0,0]
         Tp0 = wave.resource.peak_period(S).iloc[0,0]
 
-        error = np.abs(self.Tp - Tp0)/self.Tp
-        
-        self.assertLess(error, 0.01)
+        errorHm0 = np.abs(self.Tp - Tp0)/self.Tp
+        errorTp0 = np.abs(self.Hs - Hm0)/self.Hs
 
-    def test_bretschneider_spectrum(self):
-        S = wave.resource.bretschneider_spectrum(self.f,self.Tp,self.Hs)
+        self.assertLess(errorHm0, 0.01)
+        self.assertLess(errorTp0, 0.01)
+
+    def test_jonswap_spectrum(self):
+        S = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)
         Hm0 = wave.resource.significant_wave_height(S).iloc[0,0]
         Tp0 = wave.resource.peak_period(S).iloc[0,0]
 
@@ -60,7 +63,7 @@ def test_bretschneider_spectrum(self):
         self.assertLess(errorTp0, 0.01)
 
     def test_surface_elevation_seed(self):
-        S = wave.resource.bretschneider_spectrum(self.f,self.Tp,self.Hs)
+        S = wave.resource.jonswap_spectrum(self.f,self.Tp,self.Hs)
 
         sig = inspect.signature(wave.resource.surface_elevation)
         seednum = sig.parameters['seed'].default
@@ -71,7 +74,7 @@ def test_surface_elevation_seed(self):
         assert_frame_equal(eta0, eta1)        
 
     def test_surface_elevation_phasing(self):
-        S = wave.resource.bretschneider_spectrum(self.f,self.Tp,self.Hs)
+        S = wave.resource.jonswap_spectrum(self.f,self.Tp,self.Hs)
         eta0 = wave.resource.surface_elevation(S, self.t)        
         sig = inspect.signature(wave.resource.surface_elevation)
         seednum = sig.parameters['seed'].default
@@ -83,8 +86,8 @@ def test_surface_elevation_phasing(self):
 
 
     def test_surface_elevation_phases_np_and_pd(self):
-        S0 = wave.resource.bretschneider_spectrum(self.f,self.Tp,self.Hs)
-        S1 = wave.resource.bretschneider_spectrum(self.f,self.Tp,self.Hs*1.1)
+        S0 = wave.resource.jonswap_spectrum(self.f,self.Tp,self.Hs)
+        S1 = wave.resource.jonswap_spectrum(self.f,self.Tp,self.Hs*1.1)
         S = pd.concat([S0, S1], axis=1)
 
         phases_np = np.random.rand(S.shape[0], S.shape[1]) * 2 * np.pi
@@ -96,8 +99,8 @@ def test_surface_elevation_phases_np_and_pd(self):
         assert_frame_equal(eta_np, eta_pd)
 
     def test_surface_elevation_frequency_bins_np_and_pd(self):
-        S0 = wave.resource.bretschneider_spectrum(self.f,self.Tp,self.Hs)
-        S1 = wave.resource.bretschneider_spectrum(self.f,self.Tp,self.Hs*1.1)
+        S0 = wave.resource.jonswap_spectrum(self.f,self.Tp,self.Hs)
+        S1 = wave.resource.jonswap_spectrum(self.f,self.Tp,self.Hs*1.1)
         S = pd.concat([S0, S1], axis=1)
 
         eta0 = wave.resource.surface_elevation(S, self.t)
@@ -145,23 +148,12 @@ def test_surface_elevation_rmse(self):
 
         self.assertLess(rmse_sum, 0.02)
 
-    def test_jonswap_spectrum(self):
-        S = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)
-        Hm0 = wave.resource.significant_wave_height(S).iloc[0,0]
-        Tp0 = wave.resource.peak_period(S).iloc[0,0]
-        
-        errorHm0 = np.abs(self.Tp - Tp0)/self.Tp
-        errorTp0 = np.abs(self.Hs - Hm0)/self.Hs
-        
-        self.assertLess(errorHm0, 0.01)
-        self.assertLess(errorTp0, 0.01)
-    
     def test_plot_spectrum(self):            
         filename = abspath(join(testdir, 'wave_plot_spectrum.png'))
         if isfile(filename):
             os.remove(filename)
 
-        S = wave.resource.pierson_moskowitz_spectrum(self.f,self.Tp)
+        S = wave.resource.pierson_moskowitz_spectrum(self.f,self.Tp,self.Hs)
 
         plt.figure()
         wave.graphics.plot_spectrum(S)
@@ -351,7 +343,7 @@ def test_wave_celerity(self):
 
     def test_energy_flux_deep(self):
         # Dependent on mhkit.resource.BS spectrum
-        S = wave.resource.bretschneider_spectrum(self.f,self.Tp,self.Hs)
+        S = wave.resource.jonswap_spectrum(self.f,self.Tp,self.Hs)
         Te = wave.resource.energy_period(S)
         Hm0 = wave.resource.significant_wave_height(S)
         rho=1025

mhkit/wave/resource.py

@@ -63,16 +63,18 @@ def elevation_spectrum(eta, sample_rate, nnft, window='hann', detrend=True, nove
     return S
 
 
-def pierson_moskowitz_spectrum(f, Tp):
+def pierson_moskowitz_spectrum(f, Tp, Hs):
     """
-    Calculates Pierson-Moskowitz Spectrum from Tucker and Pitt (2001) 
+    Calculates Pierson-Moskowitz Spectrum from IEC TS 62600-2 ED2 Annex C.2 (2019)
     
     Parameters
     ------------
     f: numpy array
         Frequency [Hz]
     Tp: float/int
         Peak period [s]  
+    Hs: float/int
+        Significant wave height [m]         
     
     Returns
     ---------
@@ -86,59 +88,22 @@ def pierson_moskowitz_spectrum(f, Tp):
         pass
     assert isinstance(f, np.ndarray), 'f must be of type np.ndarray'
     assert isinstance(Tp, (int,float)), 'Tp must be of type int or float'
+    assert isinstance(Hs, (int,float)), 'Hs must be of type int or float'
 
     f.sort()
-    g = 9.81   
-    B_PM = (5/4)*(1/Tp)**(4)
-    A_PM = 0.0081*g**2*(2*np.pi)**(-4)
+    B_PM = (5/4)*(1/Tp)**4
+    A_PM = B_PM*(Hs/2)**2
     Sf  = A_PM*f**(-5)*np.exp(-B_PM*f**(-4)) 
 
     col_name = 'Pierson-Moskowitz ('+str(Tp)+'s)'
     S = pd.DataFrame(Sf, index=f, columns=[col_name])    
 
     return S
 
-def bretschneider_spectrum(f, Tp, Hs):
-    """
-    Calculates Bretschneider Sprectrum from Tucker and Pitt (2001)
-    
-    Parameters
-    ------------
-    f: numpy array
-        Frequency [Hz]
-    Tp: float/int
-        Peak period [s]
-    Hs: float/int
-        Significant wave height [m]        
-    
-    Returns
-    ---------    
-    S: pandas DataFrame
-        Spectral density [m^2/Hz] indexed by frequency [Hz]
 
+def jonswap_spectrum(f, Tp, Hs, gamma=None):
     """
-    try:
-        f = np.array(f)
-    except:
-        pass
-    assert isinstance(f, np.ndarray), 'f must be of type np.ndarray'
-    assert isinstance(Tp, (int,float)), 'Tp must be of type int or float'
-    assert isinstance(Hs, (int,float)), 'Hs must be of type int or float'
-
-    f.sort()
-    B_BS = (1.057/Tp)**4
-    A_BS = B_BS*(Hs/2)**2
-    Sf = A_BS*f**(-5)*np.exp(-B_BS*f**(-4))        
-
-    col_name = 'Bretschneider ('+str(Hs)+'m,'+str(Tp)+'s)'
-    S = pd.DataFrame(Sf, index=f, columns=[col_name])    
-    
-    return S
-
-
-def jonswap_spectrum(f, Tp, Hs, gamma=3.3):
-    """
-    Calculates JONSWAP spectrum from Hasselmann et al (1973)
+    Calculates JONSWAP Spectrum from IEC TS 62600-2 ED2 Annex C.2 (2019)
     
     Parameters
     ------------
@@ -164,11 +129,23 @@ def jonswap_spectrum(f, Tp, Hs, gamma=3.3):
     assert isinstance(f, np.ndarray), 'f must be of type np.ndarray'
     assert isinstance(Tp, (int,float)), 'Tp must be of type int or float'
     assert isinstance(Hs, (int,float)), 'Hs must be of type int or float'
-    assert isinstance(gamma, (int,float)), 'gamma must be of type int or float'
+    assert isinstance(gamma, (int,float, type(None))), \
+        'gamma must be of type int or float'
 
     f.sort()    
-    g = 9.81 
-    
+    B_PM = (5/4)*(1/Tp)**4
+    A_PM = B_PM*(Hs/2)**2
+    S_f  = A_PM*f**(-5)*np.exp(-B_PM*f**(-4)) 
+    
+    if not gamma:
+        TpsqrtHs = Tp/np.sqrt(Hs);
+        if TpsqrtHs <= 3.6:
+            gamma = 5;
+        elif TpsqrtHs > 5:
+            gamma = 1;
+        else:
+            gamma = np.exp(5.75 - 1.15*TpsqrtHs);
+        
     # Cutoff frequencies for gamma function
     siga = 0.07 
     sigb = 0.09 
@@ -179,9 +156,8 @@ def jonswap_spectrum(f, Tp, Hs, gamma=3.3):
     Gf = np.zeros(f.shape)
     Gf[lind] = gamma**np.exp(-(f[lind]-fp)**2/(2*siga**2*fp**2))
     Gf[hind] = gamma**np.exp(-(f[hind]-fp)**2/(2*sigb**2*fp**2))
-    S_temp = g**2*(2*np.pi)**(-4)*f**(-5)*np.exp(-(5/4)*(f/fp)**(-4))
-    alpha_JS = Hs**(2)/16/np.trapz(S_temp*Gf,f)
-    Sf = alpha_JS*S_temp*Gf # Wave Spectrum [m^2-s] 
+    C = 1- 0.287*np.log(gamma)
+    Sf = C*S_f*Gf
 
     col_name = 'JONSWAP ('+str(Hs)+'m,'+str(Tp)+'s)'
     S = pd.DataFrame(Sf, index=f, columns=[col_name])