vec3 GGX1(
    vec3 N,
    vec3 V,
    vec3 L,
    float roughness,
    vec3 F0)
{
    // half‐vector + dot products
    vec3 H = normalize(V + L);
    float NdotV = clamp(dot(N, V), 0.0, 1.0);
    float NdotL = clamp(dot(N, L), 0.0, 1.0);
    float NdotH = clamp(dot(N, H), 0.0, 1.0);
    float VdotH = clamp(dot(V, H), 0.0, 1.0);

    // sun angular radius: blend between dry & wet limits, then clamp
    float minSun = 0.005; // ~0.57°
    float maxSun = 0.200; // ~11.5°
    float sunAngRad = mix(minSun, maxSun, clamp(rainStrength, 0.0, 1.0));
    sunAngRad = clamp(sunAngRad, minSun, maxSun);

    // combined variance = roughness² + sunAngRad², with a floor
    float rou = clamp(roughness, 0.0, 1.0);
    float rou2 = rou * rou + sunAngRad * sunAngRad;
    rou2 = max(rou2, 1e-6);
    float a = sqrt(rou2);
    float a2 = a * a;

    // GGX NDF
    float denomD = NdotH * NdotH * (a2 - 1.0) + 1.0;
    denomD = max(denomD, 1e-6);
    float D = a2 / (PI * denomD * denomD);

    // Smith‐Schlick geometry
    float k = (a + 1.0) * (a + 1.0) * 0.125;
    float Gv = NdotV / (NdotV * (1.0 - k) + k);
    float Gl = NdotL / (NdotL * (1.0 - k) + k);
    float G = Gv * Gl;

    // Schlick Fresnel
    vec3 F = mix(F0, vec3(1.0), pow(1.0 - VdotH, 5.0));

    // keep denominators safe
    float denom = max(4.0 * NdotV * NdotL, 1e-6);

    // final specular, modulated by N·L
    return D * G * F * (NdotL / denom);
}

#define HALF_ANGLE 0.05 // radians
#define TILT_ANGLE 1.5  // radians

// Helper: smoothstep-interpolated 5-point curve
float sampleCurve(float r, float v0, float v1, float v2, float v3, float v4)
{
    if (r < 0.25)
        return mix(v0, v1, smoothstep(0.0, 0.25, r));
    else if (r < 0.5)
        return mix(v1, v2, smoothstep(0.25, 0.5, r));
    else if (r < 0.75)
        return mix(v2, v3, smoothstep(0.5, 0.75, r));
    else
        return mix(v3, v4, smoothstep(0.75, 1.0, r));
}

// Generalized L^p norm (Minkowski) for flexible diamond-to-circle shaping
float minkowskiNorm(vec2 v, float p)
{
    return pow(pow(abs(v.x), p) + pow(abs(v.y), p), 1.0 / p);
}

#define HALF_ANGLE 0.05 // radians
#define TILT_ANGLE 1.55 // radians
vec3 GGX3(
    vec3 normal,
    vec3 viewDir,
    vec3 lightDir,
    float roughness,
    vec3 baseReflectance)
{
    // Safety normalization
    normal = normalize(normal);
    viewDir = normalize(viewDir);
    lightDir = normalize(lightDir);

    vec3 halfVector = normalize(lightDir + viewDir);
    float NdotH = max(dot(normal, halfVector), 0.0),
          NdotL = max(dot(normal, lightDir), 0.0),
          NdotV = max(dot(normal, viewDir), 0.0),
          VdotH = max(dot(viewDir, halfVector), 0.0);

    // Reflection vector
    vec3 R = reflect(-viewDir, normal);

    // Construct lobe-aligned basis
    vec3 up = abs(lightDir.y) < 0.99 ? vec3(0, 1, 0) : vec3(1, 0, 0);
    vec3 R0 = cross(up, lightDir);
    if (length(R0) < 1e-4)
        R0 = vec3(1, 0, 0); // Fallback axis
    R0 = normalize(R0);

    float ca = cos(TILT_ANGLE), sa = sin(TILT_ANGLE);
    vec3 X = normalize(R0 * ca + cross(lightDir, R0) * sa);
    vec3 Y = normalize(cross(lightDir, X));

    // Project R into X/Y basis
    float r_dot_center = max(dot(R, lightDir), 1e-4); // Avoid divide-by-zero
    float r_dot_x = clamp(dot(R, X) / r_dot_center, -10.0, 10.0);
    float r_dot_y = clamp(dot(R, Y) / r_dot_center, -10.0, 10.0);

    // GGX-style roughness shaping

    vec3 fresnel = mix(baseReflectance, vec3(1.0), pow(1.0 - VdotH, 5.0));
    float adjustedRoughness = roughness + 1e-2;
    float alpha = adjustedRoughness * adjustedRoughness;

    float NdotH2 = NdotH * NdotH;

    float distributionDenom = NdotH2 * (alpha - 1.0) + 1.0;
    float distributionTerm = alpha / (PI * distributionDenom * distributionDenom);
    float k = (adjustedRoughness + 1.0) * (adjustedRoughness + 1.0) * 0.125;
    float geomVisibilityV = NdotV / (NdotV * (1.0 - k) + k);
    float geomVisibilityL = NdotL / (NdotL * (1.0 - k) + k);
    float denominator = max(4.0 * NdotV * NdotL, 1e-7);

    // GGX-style roughness shaping
    float baseExtent = tan(HALF_ANGLE);

    float height = baseExtent * mix(3.0, 128.0, alpha);
    float width = baseExtent * mix(2.0, 64.0, alpha);
    float feather = mix(1.0, 5.0, alpha);
    float energyScale = 1.0 / mix(1.0, 10.0, alpha);

    float normPower = mix(1.0, 1.5, alpha);
    float shapeVal = clamp(minkowskiNorm(vec2(r_dot_x / height, r_dot_y / width), normPower), 0.0, 1.0);

    float fade = clamp(1.0 - smoothstep(1.0 - feather, 1.0, shapeVal), 0.0, 1.0);

    // Fresnel (Schlick approximation)
    float G_V = NdotV / (NdotV * (1.0 - alpha) + alpha);
    float G_L = NdotL / (NdotL * (1.0 - alpha) + alpha);
    float G = G_V * G_L;

    //   return max(distributionTerm * geomVisibilityV * geomVisibilityL * fresnel * (NdotL / denominator), 0.0);

    // Final shaded specular output
    return clamp(fade * geomVisibilityV * geomVisibilityL * fresnel * (NdotL / denominator), 0, 32);
}
vec3 GGX1(vec3 normal, vec3 viewDir, vec3 lightDir, float roughness, vec3 baseReflectance); // Declare if needed

vec3 GGX1_Disc(
    vec3 normal,
    vec3 viewDir,
    vec3 lightDir,
    float roughness,
    vec3 baseReflectance,
    float thetaOffset)
{
    float r2 = roughness * roughness;
    float theta2 = thetaOffset * thetaOffset;
    float boostedRoughness = sqrt(r2 + theta2);
    return GGX1(normal, viewDir, lightDir, boostedRoughness, baseReflectance);
}
vec2 sampleCMJ(int i, int N, float noise)
{
    int sx = i % int(sqrt(float(N)));
    int sy = i / int(sqrt(float(N)));
    float jx = fract(sin(float(i) * 12.9898 + noise * 78.233) * 43758.5453);
    float jy = fract(sin(float(i) * 93.9898 + noise * 67.345) * 12345.6789);
    return (vec2(sx, sy) + vec2(jx, jy)) / sqrt(float(N));
}

vec3 GGX1_MonteCarlo(
    vec3 normal,
    vec3 viewDir,
    vec3 lightDir,
    float roughness,
    vec3 baseReflectance)
{
    vec3 X = normalize(cross(lightDir, abs(lightDir.y) < 0.99 ? vec3(0, 1, 0) : vec3(1, 0, 0)));
    vec3 Y = normalize(cross(lightDir, X));

    float extent = tan(HALF_ANGLE);
    vec3 offsetX = X * extent;
    vec3 offsetY = Y * extent;

    float noise = blueNoise(gl_FragCoord.xy);

    // Shape-GGX blend based on roughness
    float shapeBlend = clamp((0.75 - roughness) / (0.75 - 0.1), 0.0, 1.0);
    const int maxSamples = 16;
    int sampleCount = int(mix(1.0, float(maxSamples), shapeBlend) + 0.5);

    // GGX1 directly (used either as fallback or for blending)
    vec3 ggxCenter = GGX1(normal, viewDir, lightDir, roughness, baseReflectance);

    // Early out: roughness high enough, just use standard GGX
    if (sampleCount <= 1 || shapeBlend <= 0.01)
        return ggxCenter;

    vec3 result = vec3(0.0);

    for (int i = 0; i < maxSamples; ++i)
    {
        if (i >= sampleCount)
            break;

        vec2 jitteredUV = sampleCMJ(i, sampleCount, noise);
        float sx = (jitteredUV.x - 0.5) * 2.0;
        float sy = (jitteredUV.y - 0.5) * 2.0;

        // Rotate 45° to align square samples diagonally
        float d = 0.70710678;
        float sx_rot = d * (sx - sy);
        float sy_rot = d * (sx + sy);

        vec3 sampleDir = normalize(lightDir + sx_rot * offsetX + sy_rot * offsetY);

        float cosTheta = clamp(dot(lightDir, sampleDir), 0.0, 1.0);
        float thetaOffset = sqrt(2.0 * (1.0 - cosTheta));

        float falloff = pow(1.0 - clamp(length(vec2(sx_rot, sy_rot)), 0.0, 1.0), 2.0);
        float thetaBoost = thetaOffset * falloff;

        result += GGX1_Disc(normal, viewDir, sampleDir, roughness, baseReflectance, thetaBoost);
    }

    vec3 shaped = result / float(sampleCount);
    return mix(ggxCenter, shaped, shapeBlend);
}

vec3 GGX2(
    vec3 N,
    vec3 V,
    vec3 sunDir,
    float roughness,
    vec3 F0)
{
#ifdef GGX_SQUARE
    vec3 outval = GGX1_MonteCarlo(
        N,
        V,
        sunDir,
        roughness,
        F0);
#else
    vec3 outval = GGX1(
        N,
        V,
        sunDir,
        roughness,
        F0);
#endif

    return outval;
}