#version 330 core
#include "noise.glsl"
#include "planet-height.glsl"
#include "camera-ray.glsl"
#include "cloud-common.glsl"
in vec2 vUv;
uniform sampler2D uScene;
uniform sampler2D uDepth;
uniform float uFar;
uniform float uSeaRadius;
uniform float uCameraHeight;
uniform vec3 uAxisX;
uniform vec3 uAxisY;
uniform vec3 uAxisZ;
uniform vec3 uWaveVectors[32];
uniform float uWavePhases[32];
uniform float uWaveAmplitudes[32];
uniform vec3 uDetailOrigins[3];
uniform vec3 uWaveSettings; // maximum displacement, choppiness, pixel angular size
uniform vec3 uSurfaceDetail; // ripple strength, geometry range (metres), whitecap threshold
uniform vec3 uWaterOptics; // roughness, reflection, refraction
uniform vec3 uDeepColor;
uniform vec3 uShallowColor;
uniform vec3 uAbsorption;
uniform float uFoamStrength;
uniform int uLightCount;
uniform vec3 uLightDirections[2];
uniform vec3 uLightColors[2];
uniform float uLightFlux[2];
uniform float uTime;
layout(location=0) out vec4 fragColor;
layout(location=1) out vec4 reflectionNormal;
layout(location=2) out vec4 reflectionFallback;
const float PI=3.14159265359;
vec3 toLocal(vec3 v) { return vec3(dot(v,uAxisX),dot(v,uAxisY),dot(v,uAxisZ)); }
vec3 toWorld(vec3 v) { return uAxisX*v.x+uAxisY*v.y+uAxisZ*v.z; }
float sceneDistance(vec2 uv,vec3 ray) {
    return (exp2(texture(uDepth,uv).r*log2(uFar+1.0))-1.0)/max(.001,dot(ray,uForward));
}
// Stable roots avoid subtracting two planet-sized floats when the camera is metres above water.
vec2 oceanSphere(vec3 ro,vec3 rd,float expansion) {
    float b=dot(ro+vec3(0,uSeaRadius,0),rd);
    float c=dot(ro,ro)+2.0*uSeaRadius*(ro.y-expansion)-expansion*expansion;
    float h=b*b-c;
    if(h<0.0)return vec2(1e25,-1e25);
    float q=-b-(b<0.0?-1.0:1.0)*sqrt(h);
    float a=q, d=c/(abs(q)<1e-12?1e-12:q);
    return vec2(min(a,d),max(a,d));
}
float radialHeight(vec3 p) {
    return (dot(p,p)+2.0*uSeaRadius*p.y)/(length(p+vec3(0,uSeaRadius,0))+uSeaRadius);
}
float wavePhase(int i,vec3 p,out vec3 phaseGradient) {
    int partner=(i/4)*4+(i+1)%4;
    float crossing=dot(p,uWaveVectors[partner])+uWavePhases[partner];
    // Crossing waves bend and break up the crests without a tiled height texture.
    phaseGradient=uWaveVectors[i]+uWaveVectors[partner]*(.55*cos(crossing));
    return dot(p,uWaveVectors[i])+uWavePhases[i]+.55*sin(crossing);
}
float waveHeight(vec3 p,float footprint) {
    float height=0.0;
    for(int i=0;i<16;i++) {
        float width=length(uWaveVectors[i])*footprint;
        float weight=exp(-.5*width*width);
        vec3 phaseGradient;
        float base=wavePhase(i,p,phaseGradient);
        float phase=base;
        for(int j=0;j<2;j++)phase=base-uWaveSettings.y*.48*cos(phase);
        height+=sin(phase)*uWaveAmplitudes[i]*weight;
    }
    return height;
}
vec3 waveNormal(vec3 p,vec3 dx,vec3 dy,float footprint,out float compression,out float variance) {
    vec3 gradient=vec3(0); compression=0.0;
    // Waves below the smallest normal band still contribute microscopic surface roughness.
    variance=.004*uSurfaceDetail.x*uSurfaceDetail.x;
    vec3 up=normalize(p+vec3(0,uSeaRadius,0));
    for(int i=0;i<32;i++) {
        // Integrate each direction over the actual pixel footprint, including grazing views.
        vec3 phaseGradient;
        float base=wavePhase(i,p,phaseGradient);
        float width=length(vec2(dot(phaseGradient,dx),dot(phaseGradient,dy)));
        float weight=exp(-.5*width*width);
        float phase=base;
        float q=uWaveSettings.y*.48;
        for(int j=0;j<3;j++)phase=base-q*cos(phase);
        float a=uWaveAmplitudes[i]*weight;
        gradient+=phaseGradient*(a*cos(phase)/max(.15,1.0-q*sin(phase)));
        vec3 tangent=uWaveVectors[i]-up*dot(up,uWaveVectors[i]);
        float slope=uWaveAmplitudes[i]*length(tangent);
        // Unresolved waves become BRDF roughness instead of disappearing into a flat mirror.
        variance+=.5*slope*slope*(1.0-weight*weight);
        compression+=sin(phase)*slope*weight*uWaveSettings.y;
    }
    gradient-=up*dot(gradient,up);
    // Capillary ripples use a planet-anchored noise field and disappear below the pixel footprint.
    vec3 bodyPoint=mat3(uWorldToBody)*toWorld(p);
    vec3 bodyX=mat3(uWorldToBody)*uAxisX,bodyZ=mat3(uWorldToBody)*uAxisZ;
    for(int detail=0;detail<3;detail++) {
        float frequency=detail==0?.045:(detail==1?.35:1.4);
        float width=footprint*frequency*2.0;
        float weight=exp(-.5*width*width);
        // Modulo origins on the CPU preserve centimetre detail without planet-sized floats.
        vec3 drift=vec3(uTime*.018,0,uTime*.011)*sqrt(frequency);
        vec3 q=bodyPoint*frequency+uDetailOrigins[detail]+drift;
        float gx=(noise3(q+bodyX*.04)-noise3(q-bodyX*.04))/.08;
        float gz=(noise3(q+bodyZ*.04)-noise3(q-bodyZ*.04))/.08;
        float strength=(detail==0?.055:(detail==1?.14:.08))*uSurfaceDetail.x;
        gradient+=vec3(gx,0,gz)*weight*strength;
        variance+=strength*strength*.35*(1.0-weight*weight);
    }
    return normalize(up-gradient);
}
vec3 skyReflection(vec3 direction,vec3 up) {
    float daylight=0.0;
    for(int l=0;l<2;l++){if(l>=uLightCount)break;daylight+=sqrt(uLightFlux[l])*smoothstep(-.08,.16,dot(up,uLightDirections[l]));}
    float elevation=max(0.0,dot(direction,up));
    vec3 sky=mix(vec3(.56,.73,.87),vec3(.12,.32,.64),sqrt(elevation));
    for(int l=0;l<2;l++) {
        if(l>=uLightCount)break;
        sky+=uLightColors[l]*pow(max(0.0,dot(direction,uLightDirections[l])),24.0)*.22*uLightFlux[l];
    }
    return sky*(.006+daylight*1.35);
}
void main() {
    // Misses discard, preserving earlier oceans and their reflection metadata.
    reflectionNormal=vec4(0); reflectionFallback=vec4(0);
    vec3 scene=texture(uScene,vUv).rgb;
    float originalDepth=texture(uDepth,vUv).r;
    gl_FragDepth=originalDepth;
    vec3 worldRay=cameraRay(vUv), rd=toLocal(worldRay), ro=vec3(0,uCameraHeight,0);
    float opaque=sceneDistance(vUv,worldRay);
    float displacement=uWaveSettings.x;
    bool underwater=uCameraHeight<displacement && uCameraHeight<waveHeight(ro,.02);
    vec2 envelope=oceanSphere(ro,rd,displacement);
    vec2 base=oceanSphere(ro,rd,0.0);
    float t=base.x>0.0?base.x:base.y;
    bool hit=envelope.y>0.0 && t>0.0;
    float geometryRange=uSurfaceDetail.y;
    if(hit && abs(uCameraHeight)<geometryRange && t<geometryRange) {
        vec2 inner=oceanSphere(ro,rd,-displacement);
        float low=max(.015,envelope.x),high=inner.x>0.0?inner.x:envelope.y;
        if(underwater) {low=max(.015,inner.y);high=envelope.y;}
        float reference=radialHeight(ro+rd*low)-waveHeight(ro+rd*low,max(.03,low*uWaveSettings.z));
        float endValue=radialHeight(ro+rd*high)-waveHeight(ro+rd*high,max(.03,high*uWaveSettings.z));
        if(reference*endValue<0.0 && high>low) {
            for(int i=0;i<18;i++) {
                float mid=(low+high)*.5;
                vec3 p=ro+rd*mid;
                float height=radialHeight(p)-waveHeight(p,max(.03,mid*uWaveSettings.z));
                if(height*reference>0.0)low=mid;else high=mid;
            }
            float solved=(low+high)*.5;
            float geometry=1.0-smoothstep(geometryRange*.35,geometryRange,t);
            geometry*=1.0-smoothstep(geometryRange*.5,geometryRange,abs(uCameraHeight));
            t=mix(t,solved,geometry);
        }
    }
    vec3 p=ro+rd*max(0.0,t);
    // Evaluate derivatives before any divergent discard at the shoreline/horizon.
    vec3 dx=dFdx(p),dy=dFdy(p);
    vec3 radial=normalize(uAxisY+toWorld(p)/uSeaRadius);
    float waterMask=uHasOceanMask!=0?oceanMaskValue(bodyDirection(radial)):1.0;
    hit=hit && waterMask>.01 && t<opaque+.15;
    if(!hit) {
        float maskHere=uHasOceanMask!=0?oceanMaskValue(bodyDirection(uAxisY)):1.0;
        if(underwater && maskHere>.5) {
            vec3 transmission=exp(-uAbsorption*min(opaque,400.0));
            scene=scene*transmission+uShallowColor*(1.0-transmission)*.5;
            fragColor=vec4(scene,1);return;
        }
        discard;
    }
    float footprint=max(.02,max(length(dx),length(dy)));
    float compression,variance;
    vec3 n=waveNormal(p,dx,dy,footprint,compression,variance), up=normalize(p+vec3(0,uSeaRadius,0));
    float alpha=clamp(sqrt(pow(uWaterOptics.x,4.0)+variance*.7),.001,1.0);
    float roughness=clamp(sqrt(alpha),.04,1.0);
    vec3 view=-rd;
    if(dot(n,view)<0.0)n=-n;
    float ndv=max(.02,dot(n,view));
    // Average Fresnel over unresolved slopes, rather than applying the flat-surface grazing spike.
    float environmentNdv=min(1.0,sqrt(ndv*ndv+2.0*variance));
    float fresnel=.0204+(max(.0204,1.0-roughness*.65)-.0204)*pow(1.0-environmentNdv,5.0);
    vec3 reflectionDirection=reflect(rd,n);
    vec3 reflected=skyReflection(normalize(mix(reflectionDirection,up,roughness*.3)),up);
    reflected=cloudReflection(toWorld(p-ro),toWorld(reflectionDirection),reflected,roughness*roughness);
    vec2 refractUv=clamp(vUv+vec2(dot(toWorld(n),uRight),dot(toWorld(n),uUp))*uWaterOptics.z/(1.0+t*.012),vec2(.001),vec2(.999));
    float floorDistance=sceneDistance(refractUv,cameraRay(refractUv));
    float thickness=max(.0,min(400.0,floorDistance-t));
    vec3 transmission=exp(-uAbsorption*thickness);
    vec3 scatter=mix(uShallowColor,uDeepColor,1.0-exp(-thickness*.075));
    float skyDaylight=0.0;
    for(int l=0;l<2;l++){if(l>=uLightCount)break;skyDaylight+=sqrt(uLightFlux[l])*smoothstep(-.08,.16,dot(up,uLightDirections[l]));}
    // Diffuse skylight keeps deep water blue even when a cloud blocks the direct sun.
    vec3 lighting=vec3(.004)+vec3(.12,.2,.3)*skyDaylight, highlights=vec3(0);
    for(int l=0;l<2;l++) {
        if(l>=uLightCount)break;
        vec3 light=uLightDirections[l];
        float sunHeight=dot(up,light), visible=smoothstep(-.02,.03,sunHeight);
        visible*=cloudShadow(toWorld(p-ro),toWorld(light));
        float ndl=max(.001,dot(n,light));
        lighting+=uLightColors[l]*uLightFlux[l]*visible*(.25+.75*max(0.0,sunHeight));
        vec3 halfVector=normalize(view+light);
        float ndh=max(0.0,dot(n,halfVector));
        float brdfAlpha=max(.025,alpha);
        float a2=brdfAlpha*brdfAlpha, den=ndh*ndh*(a2-1.0)+1.0;
        float distribution=a2/(PI*den*den);
        float k=(brdfAlpha+1.0)*(brdfAlpha+1.0)/8.0;
        float geometry=ndv/(ndv*(1.0-k)+k)*ndl/(ndl*(1.0-k)+k);
        float f=.0204+.9796*pow(1.0-max(0.0,dot(view,halfVector)),5.0);
        highlights+=uLightColors[l]*uLightFlux[l]*visible*min(12.0,distribution*geometry*f/(4.0*ndv));
        scatter+=uShallowColor*pow(max(0.0,dot(view,-light)),4.0)*max(0.0,compression)*.28*visible;
    }
    vec3 refracted=texture(uScene,refractUv).rgb*transmission+scatter*(1.0-transmission)*lighting;
    vec3 color=mix(refracted,reflected*uWaterOptics.y,fresnel)+highlights;
    float shore=(1.0-smoothstep(.2,3.5,thickness))*(.5+.5*sin(uTime*1.3+t*.03));
    vec3 bodyDetail=mat3(uWorldToBody)*toWorld(p)*.045+uDetailOrigins[0];
    float breakup=noise3(bodyDetail+vec3(uTime*.025,0,uTime*.016))*.5+.5;
    float breaking=smoothstep(uSurfaceDetail.z,uSurfaceDetail.z+.12,compression);
    float foam=clamp((breaking*smoothstep(.25,.75,breakup)+shore*.7)*uFoamStrength,0.0,.9);
    color=mix(color,vec3(.72,.86,.88)*lighting,foam);
    if(underwater) {
        vec3 through=exp(-uAbsorption*min(t,400.0));
        color=color*through+uShallowColor*(1.0-through)*.5;
    }
    fragColor=vec4(max(color,vec3(0)),1);
    float reflectance=underwater?0.0:fresnel*uWaterOptics.y*(1.0-foam);
    reflectionNormal=vec4(toWorld(n)*reflectance,1);
    reflectionFallback=vec4(reflected*reflectance,roughness);
    gl_FragDepth=clamp(log2(1.0+t*max(.001,dot(worldRay,uForward)))/log2(uFar+1.0),0.0,1.0);
}