# 海洋知识数学笔记2

private void GenerateMesh()
{
int indiceCount = 0;
int halfResolution = resolution / 2;
for (int i = 0; i < resolution; i++)
{
float horizontalPosition = (i - halfResolution) * unitWidth;
for (int j = 0; j < resolution; j++)
{
int currentIdx = i * (resolution) + j;
float verticalPosition = (j - halfResolution) * unitWidth;
vertices[currentIdx] = new Vector3(horizontalPosition + (resolution % 2 == 0? unitWidth / 2f : 0f), 0f, verticalPosition + (resolution % 2 == 0? unitWidth / 2f : 0f));
normals[currentIdx] = new Vector3(0f, 1f, 0f);
verttilde[currentIdx] = htilde0(i, j);
htiled0就是
complex cOcean::hTilde_0(int n_prime, int m_prime) {
complex r = gaussianRandomVariable();
return r * sqrt(phillips(n_prime, m_prime) / 2.0f);
}

Vector2 temp = htilde0(resolution - i, resolution - j);
vertConj[currentIdx] = new Vector2(temp.x, -temp.y);
uvs[currentIdx] = new Vector2(i * 1.0f / (resolution - 1), j * 1.0f / (resolution - 1));
if (j == resolution - 1)
continue;
if (i != resolution - 1)
{
indices[indiceCount++] = currentIdx;
indices[indiceCount++] = currentIdx + 1;
indices[indiceCount++] = currentIdx + resolution;
}
if (i != 0)
{
indices[indiceCount++] = currentIdx;
indices[indiceCount++] = currentIdx - resolution + 1;
indices[indiceCount++] = currentIdx + 1;
}
}
}
mesh.vertices = vertices;
mesh.SetIndices(indices, MeshTopology.Triangles, 0);
mesh.normals = normals;
mesh.uv = uvs;
filter.mesh = mesh;
}

private Vector3 Displacement(Vector2 x, float t, out Vector3 nor)
{
Vector2 h = new Vector2(0f, 0f);
Vector2 d = new Vector2(0f, 0f);
Vector3 n = Vector3.zero;
Vector2 c, htilde_c, k;
float kx, kz, k_length, kDotX;
for (int i = 0; i < resolution; i++)
{
kx = 2 * PI * (i - resolution / 2.0f) / length;
for (int j = 0; j < resolution; j++)
{
kz = 2 * PI * (j - resolution / 2.0f) / length;
k = new Vector2(kx, kz);
k_length = k.magnitude;
kDotX = Vector2.Dot(k, x);
c = new Vector2(Mathf.Cos(kDotX), Mathf.Sin(kDotX));
Vector2 temp = htilde(t, i, j);
hTilde是根据n m和时间t，然后获得顶点的htiled0,     以及他们的共轭。
complex cOcean::hTilde(float t, int n_prime, int m_prime) {
int index = m_prime * Nplus1 + n_prime;

complex htilde0(vertices[index].a,  vertices[index].b);
complex htilde0mkconj(vertices[index]._a, vertices[index]._b);

float omegat = dispersion(n_prime, m_prime) * t;
这个就是根据公式去计算，这里w定义不同，代码里是开根号gk
float cOcean::dispersion(int n_prime, int m_prime) {
float w_0 = 2.0f * M_PI / 200.0f;
float kx = M_PI * (2 * n_prime - N) / length;
float kz = M_PI * (2 * m_prime - N) / length;
return floor(sqrt(g * sqrt(kx * kx + kz * kz)) / w_0) * w_0;
}
float cos_ = cos(omegat);
float sin_ = sin(omegat);

complex c0(cos_,  sin_);
complex c1(cos_, -sin_);

complex res = htilde0 * c0 + htilde0mkconj * c1;

return htilde0 * c0 + htilde0mkconj*c1;
}
htilde_c = new Vector2(temp.x * c.x - temp.y * c.y, temp.x * c.y + temp.y * c.x);
h += htilde_c;
n += new Vector3(-kx * htilde_c.y, 0f, -kz * htilde_c.y);
if (k_length < EPSILON)
continue;
d += new Vector2(kx / k_length * htilde_c.y, -kz / k_length * htilde_c.y);
}
}
nor = Vector3.Normalize(Vector3.up - n);
return new Vector3(d.x, h.x, d.y);
}

现在还有一个猜想，就是算出来的实部有用，虚部直接舍弃。

https://www.keithlantz.net/2011/10/ocean-simulation-part-one-using-the-discrete-fourier-transform/

https://www.keithlantz.net/2011/11/ocean-simulation-part-two-using-the-fast-fourier-transform/

http://www.cnblogs.com/wantnon/p/7062364.html

https://www.cnblogs.com/wantnon/p/6985141.html

  public void GenerateWaveData(int componentsPerOctave, ref float[] wavelengths, ref float[] anglesDeg, ref float[] phases)
{
int totalComponents = NUM_OCTAVES * componentsPerOctave;

if (wavelengths == null || wavelengths.Length != totalComponents) wavelengths = new float[totalComponents];
if (anglesDeg == null || anglesDeg.Length != totalComponents) anglesDeg = new float[totalComponents];
if (phases == null || phases.Length != totalComponents) phases = new float[totalComponents];

float minWavelength = Mathf.Pow(2f, SMALLEST_WL_POW_2);
float invComponentsPerOctave = 1f / componentsPerOctave;

for (int octave = 0; octave < NUM_OCTAVES; octave++)
{
for (int i = 0; i < componentsPerOctave; i++)
{
int index = octave * componentsPerOctave + i;

// stratified random sampling - should give a better range of wavelengths, and also means i can generated the
// wavelengths in sorted order!
分层采样，这里用了最小的波浪长度，每一组件都不一样
float minWavelengthi = minWavelength + invComponentsPerOctave * minWavelength * i;
最大波长就是取i = 最大值，如果两倍的最小值要更小，那么就取两倍的最小值
float maxWavelengthi = Mathf.Min(minWavelengthi + invComponentsPerOctave * minWavelength, 2f * minWavelength);
wavelengths[index] = Mathf.Lerp(minWavelengthi, maxWavelengthi, Random.value);

float rnd;

rnd = (i + Random.value) * invComponentsPerOctave;
anglesDeg[index] = (2f * rnd - 1f) * _waveDirectionVariance;

rnd = (i + Random.value) * invComponentsPerOctave;
phases[index] = 2f * Mathf.PI * rnd;
}

minWavelength *= 2f;
}
}

 public float GetAmplitude(float wavelength, float componentsPerOctave)
{
if (wavelength <= 0.001f)
{
Debug.LogError("Wavelength must be >= 0f");
return 0f;
}

float wl_pow2 = Mathf.Log(wavelength) / Mathf.Log(2f);
wl_pow2 = Mathf.Clamp(wl_pow2, SMALLEST_WL_POW_2, SMALLEST_WL_POW_2 + NUM_OCTAVES - 1f);

int index = (int)(wl_pow2 - SMALLEST_WL_POW_2);

if (index >= _powerLog.Length)
{
Debug.LogError("Out of bounds index");
return 0f;
}

if (_powerDisabled[index])
{
return 0f;
}

// The amplitude calculation follows this nice paper from Frechot:
// https://hal.archives-ouvertes.fr/file/index/docid/307938/filename/frechot_realistic_simulation_of_ocean_surface_using_wave_spectra.pdf
float wl_lo = Mathf.Pow(2f, Mathf.Floor(wl_pow2));
float k_lo = 2f * Mathf.PI / wl_lo;
float omega_lo = k_lo * ComputeWaveSpeed(wl_lo);
float wl_hi = 2f * wl_lo;
float k_hi = 2f * Mathf.PI / wl_hi;
float omega_hi = k_hi * ComputeWaveSpeed(wl_hi);

float domega = (omega_lo - omega_hi) / componentsPerOctave;

float a_2 = 2f * Mathf.Pow(10f, _powerLog[index]) * domega;
var a = Mathf.Sqrt(a_2);
return a;
}

float ComputeWaveSpeed(float wavelength)
{
// wave speed of deep sea ocean waves: https://en.wikipedia.org/wiki/Wind_wave
// https://en.wikipedia.org/wiki/Dispersion_(water_waves)#Wave_propagation_and_dispersion
float g = 9.81f;
float k = 2f * Mathf.PI / wavelength;
//float h = max(depth, 0.01);
//float cp = sqrt(abs(tanh_clamped(h * k)) * g / k);
float cp = Mathf.Sqrt(g / k);
return cp;
}

10-24
01-28 532

09-04 1万+
04-01 291
12-17 17万+
05-03 1018
02-13 7296
04-16 2732
09-21 1246
04-06 5万+
08-31 412