title: Matlab/CUDA: ocean wave simulation
body: <p>I've studied <a href="http://graphics.ucsd.edu/courses/rendering/2005/jdewall/tessendorf.pdf" rel="noreferrer">"Simulating Ocean Water"</a> article by Jerry Tessendorf and tried to program the Statistical Wave Model but I didn't get correct result and I don't understand why. </p>
  
  <p>In my program I tried only to create a wave height field at time <code>t = 0</code> without any further changes in time. After execution of my program I got not what I was expecting:<img src="https://i.stack.imgur.com/PSW7r.png" alt="enter image description here"></p>
  
  <p>Here's my source code:</p>
  
  <pre><code>clear all; close all; clc;
  rng(11); % setting seed for random numbers
  
  meshSize = 64; % field size
  windDir = [1, 0]; % ||windDir|| = 1
  patchSize = 64;
  A = 1e+4;
  g = 9.81; % gravitational constant
  windSpeed = 1e+2;
  
  x1 = linspace(-10, 10, meshSize+1); x = x1(1:meshSize);
  y1 = linspace(-10, 10, meshSize+1); y = y1(1:meshSize);
  [X,Y] = meshgrid(x, y);
  
  H0 = zeros(size(X)); % height field at time t = 0
  
  for i = 1:meshSize
      for j = 1:meshSize
          kx = 2.0 * pi / patchSize * (-meshSize / 2.0 + x(i)); % = 2*pi*n / Lx
          ky = 2.0 * pi / patchSize * (-meshSize / 2.0 + y(j)); % = 2*pi*m / Ly
          P = phillips(kx, ky, windDir, windSpeed, A, g); % phillips spectrum
          H0(i,j) = 1/sqrt(2) * (randn(1) + 1i * randn(1)) * sqrt(P);
      end
  end
  
  H0 = H0 + conj(H0);
  
  surf(X,Y,abs(ifft(H0)));
  axis([-10 10 -10 10 -10 10]);
  </code></pre>
  
  <p>And the <code>phillips</code> function:</p>
  
  <pre><code>function P = phillips(kx, ky, windDir, windSpeed, A, g)
      k_sq = kx^2 + ky^2;
      L = windSpeed^2 / g;
      k = [kx, ky] / sqrt(k_sq);
      wk = k(1) * windDir(1) + k(2) * windDir(2);
      P = A / k_sq^2 * exp(-1.0 / (k_sq * L^2)) * wk^2;
  end
  </code></pre>
  
  <p>Is there any matlab ocean simulation source code which could help me to understand my mistakes? Fast google search didn't get any results.</p>
  
  <hr>
  
  <p>Here's a "correct" result I got from "CUDA FFT Ocean Simulation". I didn't achieve this behavior in Matlab yet but I've ploted "surf" in matlab using data from  "CUDA FFT Ocean Simulation". Here's what it looks like:
  <img src="https://i.stack.imgur.com/TBrth.png" alt="enter image description here">
  <img src="https://i.stack.imgur.com/LbbQF.png" alt="enter image description here"></p>
  
  <hr>
  
  <p>I've made an experiment and got an interesting result:</p>
  
  <p>I've taken generated <code>h0</code> from "CUDA FFT Ocean Simulation". So I have to do ifft to transform from frequency domain to spatial domain to plot the graph. I've done it for the same <code>h0</code> using matlab <code>ifft</code> and using <code>cufftExecC2C</code> from CUDA library. Here's the result:</p>
  
  <p>CUDA ifft:
  <img src="https://i.stack.imgur.com/DXWXo.png" alt="enter image description here"></p>
  
  <p>Matlab ifft:
  <img src="https://i.stack.imgur.com/xXiVk.png" alt="enter image description here"></p>
  
  <p>Either I don't understand some aspects of realization of <code>cufftExecC2C</code> or <code>cufftExecC2C</code> and matlab ifft are different algorithms with different results.</p>
  
  <hr>
  
  <p>By the way parameters for generating such surface are:</p>
  
  <ul>
  <li><p>meshSize = 32</p></li>
  <li><p>A = 1e-7</p></li>
  <li><p>patchSize = 80</p></li>
  <li><p>windSpeed = 10</p></li>
  </ul>