Blender  V3.3
ocean_spectrum.c
Go to the documentation of this file.
1 /* SPDX-License-Identifier: GPL-2.0-or-later */
2 
7 #include "BKE_ocean.h"
8 #include "BLI_math.h"
9 #include "ocean_intern.h"
10 
11 #ifdef WITH_OCEANSIM
12 
13 /* -------------------------------------------------------------------- */
17 /*
18  * Original code from EncinoWaves project Copyright (c) 2015 Christopher Jon Horvath
19  * Modifications made to work within blender.
20  *
21  * Licensed under the Apache License, Version 2.0 (the "License");
22  * you may not use this file except in compliance with the License.
23  * You may obtain a copy of the License at
24  *
25  * http://www.apache.org/licenses/LICENSE-2.0
26  */
27 
32 static float alpha_beta_spectrum(const float alpha,
33  const float beta,
34  const float gamma,
35  const float omega,
36  const float peakomega)
37 {
38  return (alpha * sqrt(gamma) / pow(omega, 5.0)) * exp(-beta * pow(peakomega / omega, 4.0));
39 }
40 
41 static float peak_sharpen(const float omega, const float m_peakomega, const float m_gamma)
42 {
43  const float peak_sharpening_sigma = (omega < m_peakomega) ? 0.07 : 0.09;
44  const float peak_sharpening = pow(
45  m_gamma, exp(-sqrt((omega - m_peakomega) / (peak_sharpening_sigma * m_peakomega)) / 2.0));
46 
47  return peak_sharpening;
48 }
49 
53 static float ocean_spectrum_wind_and_damp(const Ocean *oc,
54  const float kx,
55  const float kz,
56  const float val)
57 {
58  const float k2 = kx * kx + kz * kz;
59  const float k_mag_inv = 1.0f / k2;
60  const float k_dot_w = (kx * k_mag_inv * oc->_wx) + (kz * k_mag_inv * oc->_wz);
61 
62  /* Bias towards wind dir. */
63  float newval = val * pow(fabs(k_dot_w), oc->_wind_alignment);
64 
65  /* Eliminate wavelengths smaller than cutoff. */
66  // val *= exp(-k2 * m_cutoff);
67 
68  /* Reduce reflected waves. */
69  if (k_dot_w < 0.0f) {
70  if (oc->_wind_alignment > 0.0) {
71  newval *= oc->_damp_reflections;
72  }
73  }
74 
75  return newval;
76 }
77 
78 static float jonswap(const Ocean *oc, const float k2)
79 {
80  /* Get our basic JONSWAP value from #alpha_beta_spectrum. */
81  const float k_mag = sqrt(k2);
82 
83  const float m_omega = GRAVITY * k_mag * tanh(k_mag * oc->_depth);
84  const float omega = sqrt(m_omega);
85 
86  const float m_fetch = oc->_fetch_jonswap;
87 
88  /* Strictly, this should be a random value from a Gaussian (mean 3.3, variance 0.67),
89  * clamped 1.0 to 6.0. */
90  float m_gamma = oc->_sharpen_peak_jonswap;
91  if (m_gamma < 1.0) {
92  m_gamma = 1.00;
93  }
94  if (m_gamma > 6.0) {
95  m_gamma = 6.0;
96  }
97 
98  const float m_windspeed = oc->_V;
99 
100  const float m_dimensionlessFetch = fabs(GRAVITY * m_fetch / sqrt(m_windspeed));
101  const float m_alpha = 0.076 * pow(m_dimensionlessFetch, -0.22);
102 
103  const float m_tau = M_PI * 2;
104  const float m_peakomega = m_tau * 3.5 * fabs(GRAVITY / oc->_V) *
105  pow(m_dimensionlessFetch, -0.33);
106 
107  const float beta = 1.25f;
108 
109  float val = alpha_beta_spectrum(m_alpha, beta, GRAVITY, omega, m_peakomega);
110 
111  /* Peak sharpening. */
112  val *= peak_sharpen(m_omega, m_peakomega, m_gamma);
113 
114  return val;
115 }
116 
117 float BLI_ocean_spectrum_piersonmoskowitz(const Ocean *oc, const float kx, const float kz)
118 {
119  const float k2 = kx * kx + kz * kz;
120 
121  if (k2 == 0.0f) {
122  /* No DC component. */
123  return 0.0f;
124  }
125 
126  /* Get Pierson-Moskowitz value from #alpha_beta_spectrum. */
127  const float peak_omega_PM = 0.87f * GRAVITY / oc->_V;
128 
129  const float k_mag = sqrt(k2);
130  const float m_omega = GRAVITY * k_mag * tanh(k_mag * oc->_depth);
131 
132  const float omega = sqrt(m_omega);
133  const float alpha = 0.0081f;
134  const float beta = 1.291f;
135 
136  float val = alpha_beta_spectrum(alpha, beta, GRAVITY, omega, peak_omega_PM);
137 
138  val = ocean_spectrum_wind_and_damp(oc, kx, kz, val);
139 
140  return val;
141 }
142 
143 float BLI_ocean_spectrum_texelmarsenarsloe(const Ocean *oc, const float kx, const float kz)
144 {
145  const float k2 = kx * kx + kz * kz;
146 
147  if (k2 == 0.0f) {
148  /* No DC component. */
149  return 0.0f;
150  }
151 
152  float val = jonswap(oc, k2);
153 
154  val = ocean_spectrum_wind_and_damp(oc, kx, kz, val);
155 
156  /* TMA modifications to JONSWAP. */
157  const float m_depth = oc->_depth;
158  const float gain = sqrt(m_depth / GRAVITY);
159 
160  const float k_mag = sqrt(k2);
161 
162  const float m_omega = GRAVITY * k_mag * tanh(k_mag * oc->_depth);
163  const float omega = sqrt(m_omega);
164 
165  const float kitaigorodskiiDepth_wh = omega * gain;
166  const float kitaigorodskiiDepth = 0.5 + (0.5 * tanh(1.8 * (kitaigorodskiiDepth_wh - 1.125)));
167 
168  val *= kitaigorodskiiDepth;
169 
170  val = ocean_spectrum_wind_and_damp(oc, kx, kz, val);
171 
172  return val;
173 }
174 
175 float BLI_ocean_spectrum_jonswap(const Ocean *oc, const float kx, const float kz)
176 {
177  const float k2 = kx * kx + kz * kz;
178 
179  if (k2 == 0.0f) {
180  /* No DC component. */
181  return 0.0f;
182  }
183 
184  float val = jonswap(oc, k2);
185 
186  val = ocean_spectrum_wind_and_damp(oc, kx, kz, val);
187 
188  return val;
189 }
190 
193 #endif /* WITH_OCEANSIM */
float BLI_ocean_spectrum_piersonmoskowitz(const struct Ocean *oc, float kx, float kz)
float BLI_ocean_spectrum_jonswap(const struct Ocean *oc, float kx, float kz)
float BLI_ocean_spectrum_texelmarsenarsloe(const struct Ocean *oc, float kx, float kz)
sqrt(x)+1/max(0
#define M_PI
Definition: BLI_math_base.h:20
btVector3 m_depth
ccl_device_inline float2 fabs(const float2 &a)
Definition: math_float2.h:222
ccl_device_inline float3 exp(float3 v)
Definition: math_float3.h:392
ccl_device_inline float3 pow(float3 v, float e)
Definition: math_float3.h:533
INLINE Rall1d< T, V, S > tanh(const Rall1d< T, V, S > &arg)
Definition: rall1d.h:422
ccl_device_inline float beta(float x, float y)
Definition: util/math.h:775