//////////////////////////////////////////////////////////////////////////////////////////////////////
// Bind UI parameters to uniforms used in the various programs
// The metadata after the // on each line containing a uniform declaration is a JSON object which is used to
// generate a uniform variable for the shader, which is "bound" to (i.e. driven by) a UI slider or color picker
// (depending on whether "default" is a number or array).
//////////////////////////////////////////////////////////////////////////////////////////////////////

// "Physics"
uniform float gravity;          // {"name":"gravity",          "min":0.0, "max":1.0,   "step":0.01, "default":0.05}
uniform float buoyancy;         // {"name":"buoyancy",         "min":0.0, "max":0.5,   "step":0.01, "default":0.5}
uniform float radiationLoss;    // {"name":"radiationLoss",    "min":0.9, "max":1.0,   "step":0.01,  "default":0.99}

// Blast geometry
uniform float blast_height;     // {"name":"blast_height",     "min":0.1, "max":0.9,   "step":0.01, "default":0.25}
uniform float blast_radius;     // {"name":"blast_radius",     "min":0.0, "max":0.1,   "step":0.01, "default":0.1}
uniform float blast_velocity;   // {"name":"blast_velocity",   "min":0.0, "max":100.0, "step":0.1,  "default":50.0}
uniform float blast_heat_flux;  // {"name":"blast_heat_flux",  "min":0.0, "max":100.0, "step":1.0,  "default":100.0}

// Dust
uniform float dust_inflow_rate; // {"name":"dust_inflow_rate", "min":0.0, "max":10.0,   "step":0.01, "default":1.0}
uniform vec3  dust_absorption;  // {"name":"dust_absorption",  "default":[0.5,0.5,0.5], "scale":1.0}
uniform vec3  dust_scattering;  // {"name":"dust_scattering",  "default":[0.5,0.5,0.5], "scale":1.0}

// Rendering
uniform float TtoKelvin;        // {"name":"TtoKelvin",        "min":0.0, "max":300.0,  "step":0.01, "default":10.0}

float Tambient;

/******************************************************/
/*                 mandatory function                 */
/******************************************************/
void init()
{
    // Any global constants defined here are available in all functions
    Tambient = 1.0;
}

//////////////////////////////////////////////////////////////////////////////////////////////////////
// Specify the initial conditions for the simulation,
// i.e. populate all the relevant fields (velocity, temperature, debris density/albedo) at time 0.0
//////////////////////////////////////////////////////////////////////////////////////////////////////

/******************************************************/
/*                 mandatory function                 */
/******************************************************/
void initial_conditions(in vec3 wsP,               // world space center of current voxel
                        in vec3 L, in float dL,    // world-space extents of grid, and voxel-size
                        inout vec3 v,              // initial velocity
                        inout vec4 T,              // initial temperature
                        inout vec3 medium,         // initial per-channel medium density (extinction)
                        inout vec3 mediumAlbedo)   // initial per-channel medium albedo
{
    v = vec3(0.0);
    T = vec4(Tambient);
    medium = vec3(0.0);
    mediumAlbedo = vec3(0.0);
}

//////////////////////////////////////////////////////////////////////////////////////////////////////
// Update the velocity, temperature via either:
//  - specification of volumetric inflow/outflow rate due to sources/sinks (vInflow, Tinflow)
//  - modification in-place, i.e. Dirichlet boundary conditions (v, T)
// Also specify the injected medium density inflow rate, and its scattering albedo.
//////////////////////////////////////////////////////////////////////////////////////////////////////

/******************************************************/
/*                 mandatory function                 */
/******************************************************/

void inject(in vec3 wsP,                 // world space center of current voxel
            in float time,               // time (i.e. frame count times timestep value)
            in vec3 L, in float dL,      // world-space extents of grid, and voxel-size
            inout vec3 v,                // modify velocity in-place (defaults to no change)
            inout vec3 vInflow,          // velocity inflow rate (defaults to zero)
            inout vec4 T,                // modify temperature in-place (defaults to no change)
            inout vec4 Tinflow,          // temperature inflow rate (defaults to zero)
            inout vec3 mediumInflow,     // medium density inflow rate (defaults to zero)
            inout vec3 mediumAlbedo)     // medium albedo
{
    vec3 blast_center = vec3(0.5*L.x, blast_height*L.y, 0.5*L.z);
    vec3 dir = wsP - blast_center;
    float r = length(dir);
    dir /= r;
    float rt = r/(blast_radius*L.y);
    if (rt <= 1.0 && time<400.0)
    {
        // Within blast radius: inject velocity and temperature
        float radial_falloff = max(0.0, 1.0 - rt*rt*(3.0 - 2.0*rt));
        vInflow = dir * blast_velocity * radial_falloff;
        Tinflow.r = blast_heat_flux * radial_falloff;

        // Also inject absorbing/scattering "dust"
        vec3 dust_extinction = dust_absorption + dust_scattering;
        mediumInflow = dust_extinction * dust_inflow_rate * radial_falloff;
        mediumAlbedo = dust_scattering / dust_extinction;
    }
    else
    {
        // Apply thermal relaxation due to "radiation loss"
        T.r *= radiationLoss;
    }
}

//////////////////////////////////////////////////////////////////////////////////////////////////////
// Apply any external forces to the fluid
//////////////////////////////////////////////////////////////////////////////////////////////////////

/******************************************************/
/*                 mandatory function                 */
/******************************************************/

vec3 externalForces(in vec3 wsP,                       // world space center of current voxel
                    in float time,                     // time
                    in vec3 L, in float dL,            // world-space extents of grid, and voxel-size
                    in vec3 v, in float P, in vec4 T,  // velocity, pressure, temperature at current voxel
                    in vec3 medium)                    // medium density at current voxel
{
    // Boussinesq approximation (a la Fedkiw & Stam)
    float densityAvg = (medium.r + medium.g + medium.b)/3.0;
    float buoyancy_force = -densityAvg*gravity + buoyancy*(T.r - Tambient);
    return vec3(0.0, buoyancy_force, 0.0);
}

//////////////////////////////////////////////////////////////////////////////////////////////////////
// Specify regions which contain impenetrable, static collider material
//////////////////////////////////////////////////////////////////////////////////////////////////////

float sdSphere(vec3 X, in vec3 C, float r) { return length(X-C) - r; }

/******************************************************/
/*                 mandatory function                 */
/******************************************************/

float collisionSDF(in vec3 wsP,            // world space center of current voxel
                   in float time,          // time
                   in vec3 L, in float dL) // world-space extents of grid, and voxel-size
{
    // Return SDF of the collider surface.
    // (where the interior with SDF < 0.0 is a solid obstacle)
    return sdSphere(wsP, vec3(L.x/2.0, L.y/3.0, L.z/2.0), 0.5*L.x*collider_radius);
}

//////////////////////////////////////////////////////////////////////////////////////////////////////
// Specify the fluid emission field and phase function
//////////////////////////////////////////////////////////////////////////////////////////////////////

// Approximate map from temperature in Kelvin to blackbody emission color.
// Valid from 1000 to 40000 K (and additionally 0 for pure full white)
vec3 colorTemperatureToRGB(const in float temperature)
{
    // (implementation omitted here)
}

/******************************************************/
/*                 mandatory functions                */
/******************************************************/

// Specify how the temperature is mapped to the local emission radiance
vec3 temperatureToEmission(in vec4 T)
{
    vec3 emission = colorTemperatureToRGB(T.r * TtoKelvin) * pow(T.r/100.0, 4.0);
    return emission;
}

// Optionally remap the medium density (extinction) and albedo
void mediumRemap(inout vec3 medium,
                 inout vec3 mediumAlbedo)
{}

// Specify phase function of medium
float phaseFunction(float mu,         // cosine of angle between incident and scattered ray
                    float anisotropy) // anisotropy coefficient
{
    const float pi = 3.141592653589793;
    float g = anisotropy;
    float gSqr = g*g;
    return (1.0/(4.0*pi)) * (1.0 - gSqr) / pow(1.0 - 2.0*g*mu + gSqr, 1.5);
}