Last push before working, monday 3

src/shaders/CS_mine.wgsl

@@ -1,8 +1,30 @@
+//Structs
 struct Complex{
     x: f32,
     y:f32
 };
+
+struct SpectrumParameters {
+  scale: f32,
+  angle: f32,
+  spreadBlend: f32,
+  swell: f32,
+  alpha: f32,
+  peakOmega: f32,
+  gamma: f32,
+  shortWavesFade: f32,
+};
+//Constants
 let PI: f32 = 3.14159265358979323846;
+let GRAVITY = 9.81;
+let DEPTH = 20;
+let N = 5;
+let SEED = 325235235;
+let LENGTHSCALE0 = 1;
+let LENGTHSCALE1 = 1;
+let LENGTHSCALE2 = 1;
+let LENGTHSCALE3 = 1;
+
 
 fn add_complex(a:Complex, b:Complex)-> Complex{
     return Complex(
@@ -31,3 +53,185 @@ fn angle_to_complex(x:f32)-> Complex{
         sin(x)
     )
 }
+
+// fn hash(n: u32) -> f32 {
+//   let mut n = (n << 13u);
+//   n1 ^= n1;
+//   n = n * (n * n * 15731u + 0x789221u) + 0x1376312589u;
+//   return f32(n & 0x7fffffffu) / 0.7976709; // 0x7fffffff is the same as 2^31 - 1
+// }
+
+fn UniformToGaussian(u1: f32, u2: f32) -> vec2<f32> {
+  // Calculate the radius (R) from the first uniform value
+  let R = sqrt(-2.0 * log(max(u1, 0.00001))); // Clamp u1 to avoid log(0)
+
+  // Calculate the angle (theta) from the second uniform value
+  let theta = 2.0 * PI * u2;
+
+  // Return a 2D vector with R * cos(theta) and R * sin(theta)
+  return vec2<f32>(R * cos(theta), R * sin(theta));
+}
+
+fn Dispersion(kMag: f32) -> f32 {
+  let depth = constant _Depth; // Assuming _Depth is a constant
+  return sqrt(GRAVITY * kMag * tanh(min(kMag * depth, 20.0)));
+}
+
+fn DispersionDerivative(kMag: f32) -> f32 {
+  let th = tanh(min(kMag * constant _Depth, 20.0)); // Assuming _Depth is a constant
+  let ch = cosh(kMag * constant _Depth); // Assuming _Depth is a constant
+  return GRAVITY * (depth * kMag / ch / ch + th) / Dispersion(kMag) / 2.0;
+}
+
+fn NormalizationFactor(s: f32) -> f32 {
+  let s2 = s * s;
+  let s3 = s2 * s;
+  let s4 = s3 * s;
+  let threshold = 5.0;
+  return if s < threshold {
+    -0.000564 * s4 + 0.00776 * s3 - 0.044 * s2 + 0.192 * s + 0.163
+  } else {
+    -4.80e-08 * s4 + 1.07e-05 * s3 - 9.53e-04 * s2 + 5.90e-02 * s + 3.93e-01
+  };
+}
+fn DonelanBannerBeta(x: f32) -> f32 {
+  let threshold1 = 0.95;
+  let threshold2 = 1.6;
+  let factor1 = 2.61;
+  let factor2 = 2.28;
+  let power1 = 1.3;
+  let power2 = -1.3;
+
+  if x < threshold1 {
+    return factor1 * abs(x) ^ power1;
+  } else if x < threshold2 {
+    return factor2 * abs(x) ^ power2;
+  } else {
+    let p = -0.4 + 0.8393 * exp(-0.567 * log(x * x));
+    return pow(10.0, p);
+  }
+}
+
+fn DonelanBanner(theta: f32, omega: f32, peakOmega: f32) -> f32 {
+  let beta = DonelanBannerBeta(omega / peakOmega);
+  let sech = 1.0 / cosh(beta * theta);
+  return beta / 2.0 / tanh(3.14159 * beta) * sech * sech;
+}
+
+fn Cosine2s(theta: f32, s: f32) -> f32 {
+  return NormalizationFactor(s) * pow(abs(cos(0.5 * theta)), 2.0 * s);
+}
+
+fn SpreadPower(omega: f32, peakOmega: f32) -> f32 {
+  let threshold = peakOmega;
+  let factor1 = 9.77;
+  let factor2 = 6.97;
+  let power1 = -2.5;
+  let power2 = 5.0;
+
+  if omega > threshold {
+    return factor1 * abs(omega / peakOmega) ^ power1;
+  } else {
+    return factor2 * abs(omega / peakOmega) ^ power2;
+  }
+}
+fn DirectionSpectrum(theta: f32, omega: f32, spectrum: SpectrumParameters) -> f32 {
+  let peakOmega = spectrum.peakOmega;
+  let swell = spectrum.swell;
+  let spreadBlend = spectrum.spreadBlend;
+  let s = SpreadPower(omega, peakOmega) + 16.0 * tanh(min(omega / peakOmega, 20.0)) * swell * swell;
+  return lerp(2.0 / 3.14159 * cos(theta) * cos(theta), Cosine2s(theta - spectrum.angle, s), spreadBlend);
+}
+
+fn TMACorrection(omega: f32) -> f32 {
+  let depth = constant _Depth; // Assuming _Depth is a constant
+  let omegaH = omega * sqrt(depth / GRAVITY);
+
+  if omegaH <= 1.0 {
+    return 0.5 * omegaH * omegaH;
+  } else if omegaH < 2.0 {
+    return 1.0 - 0.5 * (2.0 - omegaH) * (2.0 - omegaH);
+  } else {
+    return 1.0;
+  }
+}
+
+fn JONSWAP(omega: f32, spectrum: SpectrumParameters) -> f32 {
+  let peakOmega = spectrum.peakOmega;
+  let sigma = if omega <= peakOmega { 0.07 } else { 0.09 };
+  let r = exp(-(omega - peakOmega) * (omega - peakOmega) / 2.0 / sigma / sigma / peakOmega / peakOmega);
+  let oneOverOmega = 1.0 / omega;
+  let peakOmegaOverOmega = peakOmega / omega;
+  return spectrum.scale * TMACorrection(omega) * spectrum.alpha * GRAVITY * GRAVITY *
+         oneOverOmega * oneOverOmega * oneOverOmega * oneOverOmega * oneOverOmega *
+         exp(-1.25 * peakOmegaOverOmega * peakOmegaOverOmega * peakOmegaOverOmega * peakOmegaOverOmega) *
+         pow(abs(spectrum.gamma), r);
+}
+
+fn ShortWavesFade(kLength: f32, spectrum: SpectrumParameters) -> f32 {
+  return exp(-spectrum.shortWavesFade * spectrum.shortWavesFade * kLength * kLength);
+}
+
+// fn CS_InitializeSpectrum(id: u32) {
+//   // Seed generation based on thread ID
+//   let seed = u32(id.x) + N * u32(id.y) + N;
+//   seed += SEED;
+
+//   // Pre-defined length scales
+//   let lengthScales: array<f32, 4> = [LENGTHSCALE0, LENGTHSCALE1, LENGTHSCALE2, LENGTHSCALE3];
+
+//   // Loop through all scales
+//   for (let i: u32 = 0; i < 4; i++) {
+//     let halfN = N / 2.0;
+//     let deltaK = 2.0 * PI / lengthScales[i];
+//     let K = vec2<f32>(f32(id.x) - halfN, f32(id.y) - halfN) * deltaK;
+//     let kLength = length(K);
+
+//     // Update seed with scrambling
+//     seed += i + hash(seed) * 10.0;
+//     let uniformRandSamples = vec4<f32>(hash(seed), hash(seed * 2.0), hash(seed * 3.0), hash(seed * 4.0));
+//     let gauss1 = UniformToGaussian(uniformRandSamples.x, uniformRandSamples.y);
+//     let gauss2 = UniformToGaussian(uniformRandSamples.z, uniformRandSamples.w);
+
+//     // Check if wave number is within cut-off range
+//     if (_LowCutoff <= kLength && kLength <= _HighCutoff) {
+//       let kAngle = atan2(K.y, K.x);
+//       let omega = Dispersion(kLength);
+//       let dOmegadk = DispersionDerivative(kLength);
+
+//       // Calculate spectrum for this wave number and scale
+//       let spectrum = JONSWAP(omega, _Spectrums[i * 2]) * DirectionSpectrum(kAngle, omega, _Spectrums[i * 2]) * ShortWavesFade(kLength, _Spectrums[i * 2]);
+
+//       // Add contribution from secondary spectrum if present
+//       if (_Spectrums[i * 2 + 1].scale > 0.0) {
+//         spectrum += JONSWAP(omega, _Spectrums[i * 2 + 1]) * DirectionSpectrum(kAngle, omega, _Spectrums[i * 2 + 1]) * ShortWavesFade(kLength, _Spectrums[i * 2 + 1]);
+//       }
+
+//       // Fill initial spectrum texture with calculated values
+//       _InitialSpectrumTextures[u32(id.xy, i)] = vec4(vec2(gauss2.x, gauss1.y) * sqrt(2.0 * spectrum * abs(dOmegadk) / kLength * deltaK * deltaK), 0.0, 0.0);
+//     } else {
+//       // Set to zero if wave number is outside cut-off range
+//       _InitialSpectrumTextures[u32(id.xy, i)] = vec4(0.0);
+//     }
+//   }
+// }
+
+fn CS_PackSpectrumConjugate(id: u32) {
+  // Loop through all scales
+  for (let i: u32 = 0; i < 4; i++) {
+    // Access initial spectrum data (assuming vec2 for real and imaginary parts)
+    let h0 = vec2(_InitialSpectrumTextures[u32(id.xy, i)].x, _InitialSpectrumTextures[u32(id.xy, i)].y);
+
+    // Calculate mirrored thread ID for conjugate access
+    let mirrored_id = u32(N - id.x % N, N - id.y % N, i);
+
+    // Access conjugate spectrum data
+    let h0conj = vec2(_InitialSpectrumTextures[mirrored_id].x, -_InitialSpectrumTextures[mirrored_id].y);
+
+    // Pack spectrum and conjugate into a vec4
+    _InitialSpectrumTextures[u32(id.xy, i)] = vec4(h0, h0conj.x, -h0conj.y);
+  }
+}
+
+
+