GerstnerWavesSumOf4.hlsl

/**
 * \brief Returns positional offset for a given point as a result of summing 4 gerstner waves 
 * \param positionWS point in world space
 * \param wavelengths wavelength of each of the 4 waves
 * \param amplitudes amplitudes of each of the 4 waves
 * \param directions direction of each of the 4 waves (each row = one direction). MUST BE NORMALIZED!
 * \param speed global speed multiplier. Individual wave speed depends on Wavelength
 * \param steepness Gerstner wave 'Steepness' parameter. Modulates the horizontal offset of points
 * \param normal returns the normal of given point.
 * \return positional offset of the given point
 */
float3 GerstnerWavesSumOf4(
    float3 positionWS,
    float4 wavelengths,
    float4 amplitudes,
    float4x2 directions,
    float speed,
    float steepness,
    out float3 normal)
{
    /*
     * these variable names are from the GPUGems piece.
     * w = frequency
     * phi = speed as phase-constant
     * Q = adjusted steepness
     *
     * Deviating from the GPUGems chapter becuase we have one
     * global speed multiplier and each wave moves at a speed proportional to
     * the squareroot of it's wavelength:
     *
     * phi = wave_speed * 2 / wave_length
     * wave_speed = global_speed / sqrt(wavelength)
     * phi = globalspeed * 2 * sqrt(wavelength) / wavelength = global_speed * 2 / sqrt(wavelength)
     */
    float4 w = 2 / wavelengths;
    float4 phi = speed * 2 * rsqrt(wavelengths);
    float4 phi_t = phi * -_Time.y;
    float4 Q = steepness / (w * max(.001, amplitudes) * 4);
    
    /*
     * Just grouping the argument that goes into the sine and cosine
     */
    float4 offsetInput = w * mul(directions, positionWS.xz) + phi_t;

    float4 s = sin(offsetInput);
    float4 c = cos(offsetInput);

    float3 offset = 0;
    offset.xz = mul(Q * amplitudes * c, directions);
    offset.y = dot(amplitudes * s, float4(1, 1, 1, 1));

    /*
     * NORMALS
     * 
     * Using notation from the GPUGems chapter
     * Dx = the X direction of each wave
     * Dy = the Y direction of each wave
     * P = the offset world position of each point
     */
    float2 P = positionWS.xz + offset.xz;
    float4 Dx = directions._m00_m10_m20_m30;
    float4 Dz = directions._m01_m11_m21_m31;
    float4 WA = w * amplitudes;
    float4 normalsInput = w * mul(directions, P) + phi_t;
    float4 S = sin(normalsInput);
    float4 C = cos(normalsInput);
    
    float3 bitangent = mul(
        float3x4(
            Q * Dx * Dx * WA * S,
            Dx * WA * C,
            Q * Dx * Dz * WA * S), float4(1, 1, 1, 1));
    // Tangent
    float3 tangent = mul(
        float3x4(
            Q * Dx * Dz * WA * S,
            Dz * WA * C,
            Q * Dz * Dz * WA * S), float4(1, 1, 1, 1));

    /*
     * The above matrix multiply takes care of the summation,
     * But these terms are outside of the summation:
     */
    bitangent.x = 1 - bitangent.x;
    bitangent.z = -bitangent.z;
    tangent.x = -tangent.x;
    tangent.z = 1 - tangent.z;


    /*
     * These need to be normalized because
     * the tangent space can be 'stretched' by the horizontal motion,
     * which means the tangent vectors are no longer length 1
     * But then it's better to just normalize the resulting normal
     * (because cross product will just multiply magnitude of result)
     */
    normal = normalize(cross(tangent, bitangent));

    return offset;
}

We don't really need amplitudes term in Q and WA, do we? This term will cancel out during the computation, (except line 71 and 77).