Stable Fluid for C# （）

FluidMath.cs

/**
 *  Jos Stam's Stable Fluids (for C#)
 **/

public static class FluidMath
{
    public static void VelStep(int N, ref float[] u, ref float[] v, ref float[] u0, ref float[] v0, float visc, float dt)
    {
        AddSource(N, ref u, ref u0, dt);
        AddSource(N, ref v, ref v0, dt);
        Swap(ref u0, ref u);
        Diffuse(N, 1, ref u, ref u0, visc, dt);
        Swap(ref v0, ref v);
        Diffuse(N, 2, ref v, ref v0, visc, dt);
        Project(N, ref u, ref v, ref u0, ref v0);
        Swap(ref u0, ref u);
        Swap(ref v0, ref v);
        Advect(N, 1, ref u, ref u0, ref u0, ref v0, dt);
        Advect(N, 2, ref v, ref v0, ref u0, ref v0, dt);
        Project(N, ref u, ref v, ref u0, ref v0);
    }

    public static void DensStep(int N, ref float[] x, ref float[] x0, ref float[] u, ref float[] v, float diff, float dt)
    {
        AddSource(N, ref x, ref x0, dt);
        Swap(ref x0, ref x);
        Diffuse(N, 0, ref x, ref x0, diff, dt);
        Swap(ref x0, ref x);
        Advect(N, 0, ref x, ref x0, ref u, ref v, dt);
    }
    
    public static int IX(int N, int i, int j)
    {
        return i + (N + 2) * j;
    }

    private static void Swap(ref float[] x, ref float[] x0)
    {
        float[] tmp = x0;
        x0 = x;
        x = tmp;
    }

    private static void AddSource(int N, ref float[] x, ref float[] s, float dt)
    {
        int size = (N + 2) * (N + 2);
        for (int i = 0; i < size; i++) x[i] += dt * s[i];
    }

    private static void SetBnd(int N, int b, ref float[] x)
    {
        for (int i = 1; i <= N; i++)
        {
            x[IX(N, 0    , i)] = b == 1 ? -x[IX(N, 1, i)] : x[IX(N, 1, i)];
            x[IX(N, N + 1, i)] = b == 1 ? -x[IX(N, N, i)] : x[IX(N, N, i)];
            x[IX(N, i    , 0)] = b == 2 ? -x[IX(N, i, 1)] : x[IX(N, i, 1)];
            x[IX(N, i, N + 1)] = b == 2 ? -x[IX(N, i, N)] : x[IX(N, i, N)];
        }

        x[IX(N, 0        , 0)] = 0.5f * (x[IX(N, 1, 0)]     + x[IX(N, 0, 1)]);
        x[IX(N, 0    , N + 1)] = 0.5f * (x[IX(N, 1, N + 1)] + x[IX(N, 0, N)]);
        x[IX(N, N + 1    , 0)] = 0.5f * (x[IX(N, N, 0)]     + x[IX(N, N + 1, 1)]);
        x[IX(N, N + 1, N + 1)] = 0.5f * (x[IX(N, N, N + 1)] + x[IX(N, N + 1, N)]);
    }

    private static void LinSolve(int N, int b, ref float[] x, ref float[] x0, float a, float c)
    {
        for (int k = 0; k < 20; k++)
        {
            for (int i = 1; i <= N; i++)
            {
                for (int j = 1; j <= N; j++)
                {
                    x[IX(N, i, j)] = (x0[IX(N, i, j)] + a * (x[IX(N, i - 1, j)] + x[IX(N, i + 1, j)] + x[IX(N, i, j - 1)] + x[IX(N, i, j + 1)])) / c;
                }
            }

            SetBnd(N, b, ref x);
        }
    }

    private static void Diffuse(int N, int b, ref float[] x, ref float[] x0, float diff, float dt)
    {
        float a = dt * diff * N * N;
        LinSolve(N, b, ref x, ref x0, a, 1 + 4 * a);
    }

    private static void Advect(int N, int b, ref  float[] d, ref float[] d0, ref float[] u, ref float[] v, float dt)
    {
        int i0, j0, i1, j1;
        float x, y, s0, t0, s1, t1, dt0;

        dt0 = dt * N;
        for (int i = 1; i <= N; i++)
        {
            for (int j = 1; j <= N; j++)
            {
                x = i - dt0 * u[IX(N, i, j)]; y = j - dt0 * v[IX(N, i, j)];
                if (x < 0.5f) x = 0.5f; if (x > N + 0.5f) x = N + 0.5f; i0 = (int)x; i1 = i0 + 1;
                if (y < 0.5f) y = 0.5f; if (y > N + 0.5f) y = N + 0.5f; j0 = (int)y; j1 = j0 + 1;
                s1 = x - i0; s0 = 1 - s1; t1 = y - j0; t0 = 1 - t1;
                d[IX(N, i, j)] = s0 * (t0 * d0[IX(N, i0, j0)] + t1 * d0[IX(N, i0, j1)]) + s1 * (t0 * d0[IX(N, i1, j0)] + t1 * d0[IX(N, i1, j1)]);
            }
        }

        SetBnd(N, b, ref d);
    }

    private static void Project(int N, ref float[] u, ref float[] v, ref float[] p, ref float[] div)
    {
        int i, j;

        for (i = 1; i <= N; i++)
        {
            for (j = 1; j <= N; j++)
            {
                div[IX(N, i, j)] = -0.5f * (u[IX(N, i + 1, j)] - u[IX(N, i - 1, j)] + v[IX(N, i, j + 1)] - v[IX(N, i, j - 1)]) / N;
                p[IX(N, i, j)] = 0;
            }
        }

        SetBnd(N, 0, ref div);
        SetBnd(N, 0, ref p);
        LinSolve(N, 0, ref p, ref div, 1, 4);

        for (i = 1; i <= N; i++)
        {
            for (j = 1; j <= N; j++)
            {
                u[IX(N, i, j)] -= 0.5f * N * (p[IX(N, i + 1, j)] - p[IX(N, i - 1, j)]);
                v[IX(N, i, j)] -= 0.5f * N * (p[IX(N, i, j + 1)] - p[IX(N, i, j - 1)]);
            }
        }

        SetBnd(N, 1, ref u);
        SetBnd(N, 2, ref v);
    }
}