#ifndef __GLITTER_INC
#define __GLITTER_INC

#define PI 3.1415926535897932384626433832795028841971

float gaussian(float3 x, float3 mu, float sigma) {
  float3 x_minus_mu = x - mu;
  return exp(-dot(x_minus_mu, x_minus_mu) / (2 * sigma * sigma)) / (sigma * sqrt(2 * PI));
}

float det(float2x2 m) {
  return m[0][0] * m[1][1] - m[0][1] * m[1][0];
}

float dirac(float3 x) {
  // Ternary operator short-circuits in slang >:(
  return any(abs(x) < 1e-4) ? 1e4 : 0;
}

// Maps a point on the unit hemisphere to an inscribed circle.
float2 hemisphere_to_inscribed_circle(float3 p) {
  // TODO
  return 0;
}

float3 rand3_hash(float3 p)
{
    // Improved Murmurhash3 by Squirrel Eiserloh (GDC 2017)
    p = float3(dot(p, float3(127.1, 311.7, 74.7)),
               dot(p, float3(269.5, 183.3, 246.1)),
               dot(p, float3(113.5, 271.9, 124.6)));
    return -1.0 + 2.0 * frac(sin(p) * 43758.5453123);
}

public float D_Kemppinen(
  float3 worldPos,
  float3 H_tan,  // halfway vector in tangent space
  float amount,
  float roughness)
{
  float2x2 J_T__omega_h = 0;  // TODO
  float det_J_T__omega_h = det(J_T__omega_h);
  
  float sum = 0;
  // x is the current world space position.
  float3 x = worldPos;
  // omega_h is the halfway vector *in TBN coordinates.*
  float3 omega_h = H_tan;
  float3 T_omega_h = float3(
      hemisphere_to_inscribed_circle(omega_h),
      0);
  uint K_p = 1;  // TODO
  for (uint i = 0; i < K_p; ++i) {
    // x_i is the position of the facet.
    float3 x_i = worldPos;  // TODO randomize position
    // mu_i is the orientation of the facet.
    float3 mu_i = abs(rand3_hash(x_i));

    float a_star = roughness;

    float cur_gauss = gaussian(
        T_omega_h,
        mu_i,
        a_star * sqrt(det_J_T__omega_h));

    sum += cur_gauss * dirac(x_i - x);
  }

  float scale = det_J_T__omega_h / K_p;
  return sum * scale;
}

#endif  // __GLITTER_INC