#define GL_GLEXT_PROTOTYPES 1

#include <EGL/egl.h>
#include <EGL/eglext.h>
#include <GL/glcorearb.h>
#include <png.h>

#include <algorithm>
#include <array>
#include <bit>
#include <cerrno>
#include <chrono>
#include <cctype>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <set>
#include <span>
#include <string>
#include <string_view>
#include <unistd.h>
#include <utility>
#include <vector>

namespace
{
    constexpr double PlaneLengthMeters=12.0;
    constexpr double WaveSpeedMetersPerSecond=4.5;
    constexpr double DampingPerSecond=1.35;
    constexpr double FoundationStiffnessPerSecondSquared=28.0;
    constexpr double SphereRadiusMeters=0.43;
    constexpr double SphereSigmaMeters=0.58;
    constexpr double SpherePushMetersPerSecondSquared=36.0;
    constexpr GLuint TimerQueryRingSize=64;
    constexpr double Pi=3.14159265358979323846;

    struct FBodyState
    {
        double X=0.0,Y=0.0,Z=0.0,Force=0.0;
    };

    struct FOptions
    {
        std::string Device="auto";
        std::filesystem::path Output="Saved/Linux/l164/living-plane";
        std::string Scene="wake";
        int Size=128;
        int Steps=360;
        int CaptureEvery=1;
        bool TimerQueries=true;
        bool Breathing=false;
        bool BreathingConfigured=false;
        bool Imprint=false;
        bool ImprintConfigured=false;
        double BreathingStrength=4.8;
        double ImprintCaptureMetersPerSecond=0.15;
        double ImprintRecoverySeconds=6.0;
        double ImprintMaxDepthMeters=0.34;
    };

    struct FDeviceCandidate
    {
        EGLDeviceEXT Handle=EGL_NO_DEVICE_EXT;
        std::string Vendor;
        std::string Extensions;
        std::string RenderNode;
        int Score=0;
    };

    struct FContext
    {
        EGLDisplay Display=EGL_NO_DISPLAY;
        EGLSurface Surface=EGL_NO_SURFACE;
        EGLContext Context=EGL_NO_CONTEXT;
        EGLDeviceEXT Device=EGL_NO_DEVICE_EXT;
        std::string DeviceLabel="surfaceless Mesa";
        std::string RenderNode;
        std::string ExecutionDevice;
        std::string Error;

        ~FContext() { Reset(); }

        void Reset()
        {
            if(Display!=EGL_NO_DISPLAY)
            {
                eglMakeCurrent(Display,EGL_NO_SURFACE,EGL_NO_SURFACE,EGL_NO_CONTEXT);
                if(Context!=EGL_NO_CONTEXT)eglDestroyContext(Display,Context);
                if(Surface!=EGL_NO_SURFACE)eglDestroySurface(Display,Surface);
                eglTerminate(Display);
            }
            Display=EGL_NO_DISPLAY;
            Surface=EGL_NO_SURFACE;
            Context=EGL_NO_CONTEXT;
        }

        bool TryDisplay(EGLDisplay Candidate,const std::string& Label)
        {
            Reset();
            if(Candidate==EGL_NO_DISPLAY)
            {Error="eglGetPlatformDisplay returned EGL_NO_DISPLAY";return false;}
            Display=Candidate;
            EGLint Major=0,Minor=0;
            if(eglInitialize(Display,&Major,&Minor)!=EGL_TRUE)
            {Error="eglInitialize failed: "+EglError();Reset();return false;}
            if(eglBindAPI(EGL_OPENGL_API)!=EGL_TRUE)
            {Error="eglBindAPI(EGL_OPENGL_API) failed: "+EglError();Reset();return false;}

            const EGLint ConfigAttributes[]={EGL_SURFACE_TYPE,EGL_PBUFFER_BIT,
                EGL_RENDERABLE_TYPE,EGL_OPENGL_BIT,EGL_RED_SIZE,8,EGL_GREEN_SIZE,8,
                EGL_BLUE_SIZE,8,EGL_ALPHA_SIZE,8,EGL_DEPTH_SIZE,24,EGL_NONE};
            EGLConfig Config{};
            EGLint ConfigCount=0;
            if(eglChooseConfig(Display,ConfigAttributes,&Config,1,&ConfigCount)!=EGL_TRUE
                ||ConfigCount==0)
            {Error="no desktop OpenGL pbuffer config: "+EglError();Reset();return false;}

            const EGLint SurfaceAttributes[]={EGL_WIDTH,16,EGL_HEIGHT,16,EGL_NONE};
            Surface=eglCreatePbufferSurface(Display,Config,SurfaceAttributes);
            const EGLint ContextAttributes[]={EGL_CONTEXT_MAJOR_VERSION_KHR,4,
                EGL_CONTEXT_MINOR_VERSION_KHR,3,EGL_CONTEXT_OPENGL_PROFILE_MASK_KHR,
                EGL_CONTEXT_OPENGL_CORE_PROFILE_BIT_KHR,EGL_NONE};
            Context=eglCreateContext(Display,Config,EGL_NO_CONTEXT,ContextAttributes);
            if(Surface==EGL_NO_SURFACE||Context==EGL_NO_CONTEXT
                ||eglMakeCurrent(Display,Surface,Surface,Context)!=EGL_TRUE)
            {Error="OpenGL 4.3 context creation failed: "+EglError();Reset();return false;}

            GLint GlMajor=0,GlMinor=0;
            glGetIntegerv(GL_MAJOR_VERSION,&GlMajor);
            glGetIntegerv(GL_MINOR_VERSION,&GlMinor);
            if(GlMajor<4||(GlMajor==4&&GlMinor<3))
            {Error="driver returned OpenGL "+std::to_string(GlMajor)+"."+
                std::to_string(GlMinor)+", need 4.3";Reset();return false;}

            const auto* Renderer=reinterpret_cast<const char*>(glGetString(GL_RENDERER));
            if(!Renderer)
            {Error="glGetString(GL_RENDERER) returned null";Reset();return false;}
            if(IsSoftwareRenderer(Renderer))ExecutionDevice="cpu";
            else ExecutionDevice="gpu";
            DeviceLabel=Label;
            return true;
        }

        static std::string EglError()
        {
            char Text[24]{};
            std::snprintf(Text,sizeof(Text),"0x%04x",eglGetError());
            return Text;
        }

        static std::string Lower(std::string_view Input)
        {
            std::string Result(Input);
            std::transform(Result.begin(),Result.end(),Result.begin(),[](unsigned char C)
            {return static_cast<char>(std::tolower(C));});
            return Result;
        }

        static bool IsSoftwareRenderer(std::string_view Renderer)
        {
            const std::string LowerRenderer=Lower(Renderer);
            return LowerRenderer.find("llvmpipe")!=std::string::npos
                ||LowerRenderer.find("softpipe")!=std::string::npos
                ||LowerRenderer.find("software rasterizer")!=std::string::npos;
        }
    };

    std::string JsonEscape(std::string_view Text)
    {
        std::string Result;
        for(const unsigned char C:Text)
        {
            switch(C)
            {
                case '"':Result+="\\\"";break;
                case '\\':Result+="\\\\";break;
                case '\n':Result+="\\n";break;
                case '\r':Result+="\\r";break;
                case '\t':Result+="\\t";break;
                default:if(C<0x20)Result+=' ';else Result+=static_cast<char>(C);break;
            }
        }
        return Result;
    }

    bool HasExtension(std::string_view Extensions,std::string_view Wanted)
    {
        std::size_t At=0;
        while(At<Extensions.size())
        {
            while(At<Extensions.size()&&Extensions[At]==' ')++At;
            const std::size_t End=Extensions.find(' ',At);
            const std::string_view Part=Extensions.substr(At,
                End==std::string_view::npos?Extensions.size()-At:End-At);
            if(Part==Wanted)return true;
            if(End==std::string_view::npos)break;
            At=End+1;
        }
        return false;
    }

    std::vector<FDeviceCandidate> FindHardwareDevices(std::string& Error)
    {
        using TQueryDevices=PFNEGLQUERYDEVICESEXTPROC;
        using TQueryDeviceString=PFNEGLQUERYDEVICESTRINGEXTPROC;
        using TGetPlatformDisplay=PFNEGLGETPLATFORMDISPLAYEXTPROC;
        auto QueryDevices=reinterpret_cast<TQueryDevices>(eglGetProcAddress("eglQueryDevicesEXT"));
        auto QueryDeviceString=reinterpret_cast<TQueryDeviceString>(
            eglGetProcAddress("eglQueryDeviceStringEXT"));
        auto GetPlatformDisplay=reinterpret_cast<TGetPlatformDisplay>(
            eglGetProcAddress("eglGetPlatformDisplayEXT"));
        const auto* ClientExtensions=reinterpret_cast<const char*>(
            eglQueryString(EGL_NO_DISPLAY,EGL_EXTENSIONS));
        if(!QueryDevices||!QueryDeviceString||!GetPlatformDisplay||!ClientExtensions
            ||!HasExtension(ClientExtensions,"EGL_EXT_device_enumeration")
            ||!HasExtension(ClientExtensions,"EGL_EXT_platform_device"))
        {Error="EGL device enumeration/platform-device extensions are unavailable";return {};}

        EGLint Count=0;
        if(QueryDevices(0,nullptr,&Count)!=EGL_TRUE||Count<=0)
        {Error="eglQueryDevicesEXT found no EGL devices";return {};}
        std::vector<EGLDeviceEXT> Devices(static_cast<std::size_t>(Count));
        if(QueryDevices(Count,Devices.data(),&Count)!=EGL_TRUE)
        {Error="eglQueryDevicesEXT failed while listing devices";return {};}

        std::vector<FDeviceCandidate> Candidates;
        for(EGLint I=0;I<Count;++I)
        {
            FDeviceCandidate Candidate;
            Candidate.Handle=Devices[static_cast<std::size_t>(I)];
            const auto* Vendor=QueryDeviceString(Candidate.Handle,EGL_VENDOR);
            const auto* Extensions=QueryDeviceString(Candidate.Handle,EGL_EXTENSIONS);
            if(Vendor)Candidate.Vendor=reinterpret_cast<const char*>(Vendor);
            if(Extensions)Candidate.Extensions=reinterpret_cast<const char*>(Extensions);
            if(!HasExtension(Candidate.Extensions,"EGL_EXT_device_drm_render_node"))continue;
            const auto* Node=QueryDeviceString(Candidate.Handle,EGL_DRM_RENDER_NODE_FILE_EXT);
            if(Node)Candidate.RenderNode=reinterpret_cast<const char*>(Node);
            if(Candidate.RenderNode.empty()
                ||access(Candidate.RenderNode.c_str(),R_OK|W_OK)!=0)continue;

            const std::string VendorLower=FContext::Lower(Candidate.Vendor);
            const bool IsIntel=VendorLower.find("intel")!=std::string::npos
                ||Candidate.RenderNode.find("renderD128")!=std::string::npos;
            // This host's integrated Alder Lake device is exposed by Mesa as Intel
            // on renderD128. EGL does not expose Vulkan's physical-device-type enum.
            Candidate.Score=IsIntel?100:50;
            Candidates.push_back(std::move(Candidate));
        }
        std::stable_sort(Candidates.begin(),Candidates.end(),[](const auto& A,const auto& B)
        {return A.Score>B.Score;});
        if(Candidates.empty()&&Error.empty())
            Error="no accessible EGL DRM render node was found";
        return Candidates;
    }

    bool CreateContext(const std::string& Request,FContext& Out)
    {
        if(Request=="cpu")
        {
            setenv("LIBGL_ALWAYS_SOFTWARE","1",1);
            const auto GetPlatformDisplay=reinterpret_cast<PFNEGLGETPLATFORMDISPLAYEXTPROC>(
                eglGetProcAddress("eglGetPlatformDisplayEXT"));
#ifdef EGL_PLATFORM_SURFACELESS_MESA
            EGLDisplay Display=EGL_NO_DISPLAY;
            if(GetPlatformDisplay)Display=GetPlatformDisplay(EGL_PLATFORM_SURFACELESS_MESA,
                EGL_DEFAULT_DISPLAY,nullptr);
            if(Display==EGL_NO_DISPLAY)Display=eglGetDisplay(EGL_DEFAULT_DISPLAY);
            if(!Out.TryDisplay(Display,"Mesa surfaceless (LIBGL_ALWAYS_SOFTWARE=1)"))
            {std::cerr<<"[LivingPlane] CPU context failed: "<<Out.Error<<"\n";return false;}
#else
            const EGLDisplay Display=eglGetDisplay(EGL_DEFAULT_DISPLAY);
            if(!Out.TryDisplay(Display,"EGL default display (LIBGL_ALWAYS_SOFTWARE=1"))
            {std::cerr<<"[LivingPlane] CPU context failed: "<<Out.Error<<"\n";return false;}
#endif
            if(Out.ExecutionDevice!="cpu")
            {std::cerr<<"[LivingPlane] --device=cpu selected a non-software renderer\n";return false;}
            return true;
        }

        unsetenv("LIBGL_ALWAYS_SOFTWARE");
        std::string EnumerationError;
        const auto Candidates=FindHardwareDevices(EnumerationError);
        if(Candidates.empty())
        {
            if(Request=="gpu")
            {std::cerr<<"[LivingPlane] --device=gpu unavailable: "<<EnumerationError
                <<" (or no accessible DRM render node)\n";return false;}
        }
        else
        {
            auto GetPlatformDisplay=reinterpret_cast<PFNEGLGETPLATFORMDISPLAYEXTPROC>(
                eglGetProcAddress("eglGetPlatformDisplayEXT"));
            for(const auto& Candidate:Candidates)
            {
                // Automatic selection is intentionally limited to this host's
                // integrated Intel target; a discrete-only host falls back to CPU.
                if(Request=="auto"&&Candidate.Score<100)continue;
                EGLDisplay Display=GetPlatformDisplay(EGL_PLATFORM_DEVICE_EXT,
                    Candidate.Handle,nullptr);
                const std::string Label=Candidate.Vendor+" "+Candidate.RenderNode;
                if(Out.TryDisplay(Display,Label))
                {
                    Out.Device=Candidate.Handle;
                    Out.RenderNode=Candidate.RenderNode;
                    if(Out.ExecutionDevice=="gpu")return true;
                    Out.Reset();
                    Out.Error="hardware EGL device produced a software renderer";
                }
            }
            if(Request=="gpu")
            {std::cerr<<"[LivingPlane] --device=gpu context failed: "<<Out.Error<<"\n";return false;}
        }

        setenv("LIBGL_ALWAYS_SOFTWARE","1",1);
        auto GetPlatformDisplay=reinterpret_cast<PFNEGLGETPLATFORMDISPLAYEXTPROC>(
            eglGetProcAddress("eglGetPlatformDisplayEXT"));
        EGLDisplay SoftwareDisplay=EGL_NO_DISPLAY;
#ifdef EGL_PLATFORM_SURFACELESS_MESA
        if(GetPlatformDisplay)SoftwareDisplay=GetPlatformDisplay(EGL_PLATFORM_SURFACELESS_MESA,
            EGL_DEFAULT_DISPLAY,nullptr);
#endif
        if(SoftwareDisplay==EGL_NO_DISPLAY)SoftwareDisplay=eglGetDisplay(EGL_DEFAULT_DISPLAY);
        if(!Out.TryDisplay(SoftwareDisplay,"Mesa surfaceless (LIBGL_ALWAYS_SOFTWARE=1)"))
        {std::cerr<<"[LivingPlane] llvmpipe fallback failed: "<<Out.Error<<"\n";return false;}
        if(Out.ExecutionDevice!="cpu")
        {std::cerr<<"[LivingPlane] software fallback did not return a software renderer\n";return false;}
        return true;
    }

    GLuint CompileShader(GLenum Type,const char* Source)
    {
        const GLuint Shader=glCreateShader(Type);
        glShaderSource(Shader,1,&Source,nullptr);
        glCompileShader(Shader);
        GLint Success=GL_FALSE;
        glGetShaderiv(Shader,GL_COMPILE_STATUS,&Success);
        if(Success==GL_FALSE)
        {
            GLint Length=0;
            glGetShaderiv(Shader,GL_INFO_LOG_LENGTH,&Length);
            std::string Log(static_cast<std::size_t>(std::max(Length,1)),0);
            glGetShaderInfoLog(Shader,Length,nullptr,Log.data());
            std::cerr<<"[LivingPlane] shader compilation failed:\n"<<Log<<"\n";
            glDeleteShader(Shader);
            return 0;
        }
        return Shader;
    }

    GLuint LinkProgram(std::span<const GLuint> Shaders)
    {
        const GLuint Program=glCreateProgram();
        for(const GLuint Shader:Shaders)glAttachShader(Program,Shader);
        glLinkProgram(Program);
        GLint Success=GL_FALSE;
        glGetProgramiv(Program,GL_LINK_STATUS,&Success);
        if(Success==GL_FALSE)
        {
            GLint Length=0;
            glGetProgramiv(Program,GL_INFO_LOG_LENGTH,&Length);
            std::string Log(static_cast<std::size_t>(std::max(Length,1)),0);
            glGetProgramInfoLog(Program,Length,nullptr,Log.data());
            std::cerr<<"[LivingPlane] program link failed:\n"<<Log<<"\n";
            glDeleteProgram(Program);
            return 0;
        }
        return Program;
    }

    GLuint MakeComputeProgram()
    {
        static constexpr char Source[]=R"GLSL(#version 430 core
layout(local_size_x=16,local_size_y=16) in;
layout(binding=0,rgba32f) readonly uniform image2D uSource;
layout(binding=1,rgba32f) writeonly uniform image2D uTarget;
uniform int uGridSize;
uniform float uDt;
uniform float uDx;
uniform float uPlaneLength;
uniform float uWaveSpeed;
uniform float uDamping;
uniform float uFoundationStiffness;
uniform float uTime;
uniform vec2 uForcePositions[2];
uniform float uForceSigma;
uniform float uForceStrengths[2];
uniform int uBreathing;
uniform float uBreathingStrength;
uniform int uImprint;
uniform float uImprintCaptureRate;
uniform float uImprintRecoverySeconds;
uniform float uImprintMaxDepth;
void main()
{
    ivec2 p=ivec2(gl_GlobalInvocationID.xy);
    if(p.x>=uGridSize||p.y>=uGridSize)return;
    ivec2 lo=ivec2(0);
    ivec2 hi=ivec2(uGridSize-1);
    vec4 state=imageLoad(uSource,p);
    float h=state.r;
    float v=state.g;
    float imprint=state.b;
    float l=imageLoad(uSource,ivec2(max(p.x-1,lo.x),p.y)).r;
    float r=imageLoad(uSource,ivec2(min(p.x+1,hi.x),p.y)).r;
    float d=imageLoad(uSource,ivec2(p.x,max(p.y-1,lo.y))).r;
    float u=imageLoad(uSource,ivec2(p.x,min(p.y+1,hi.y))).r;
    float laplacian=(l+r+d+u-4.0*h)/(uDx*uDx);
    vec2 world=(vec2(p)+vec2(0.5))*uDx-vec2(uPlaneLength*0.5);
    float sigma2=max(uForceSigma*uForceSigma,1.0e-8);
    float interaction=0.0;
    float contact=0.0;
    for(int i=0;i<2;++i)
    {
        vec2 offset=world-uForcePositions[i];
        float weight=exp(-dot(offset,offset)/(2.0*sigma2));
        interaction-=uForceStrengths[i]*weight;
        if(uForceStrengths[i]>0.0)contact=max(contact,weight);
    }
    float breathing=0.0;
    if(uBreathing!=0)
    {
        float swell=sin(world.x*0.48-uTime*0.22)
            *cos(world.y*0.42+uTime*0.16);
        breathing=uBreathingStrength*swell*sin(uTime*0.48);
    }
    float nextImprint=imprint;
    if(uImprint!=0)
    {
        nextImprint=clamp(imprint*exp(-uDt/max(uImprintRecoverySeconds,1.0e-4))
            -contact*uImprintCaptureRate*uDt,-uImprintMaxDepth,0.0);
    }
    float foundationTarget=uImprint!=0?imprint:0.0;
    float acceleration=uWaveSpeed*uWaveSpeed*laplacian
        -uFoundationStiffness*(h-foundationTarget)
        +interaction+breathing;
    float nextVelocity=(v+acceleration*uDt)*exp(-uDamping*uDt);
    float nextHeight=h+nextVelocity*uDt;
    imageStore(uTarget,p,vec4(nextHeight,nextVelocity,nextImprint,0.0));
})GLSL";
        const GLuint Shader=CompileShader(GL_COMPUTE_SHADER,Source);
        if(!Shader)return 0;
        const std::array<GLuint,1> Shaders={Shader};
        const GLuint Program=LinkProgram(Shaders);
        glDeleteShader(Shader);
        return Program;
    }

    GLuint MakeRenderProgram()
    {
        static constexpr char VertexSource[]=R"GLSL(#version 430 core
layout(location=0) in vec3 aPosition;
layout(location=1) in vec2 aUv;
layout(location=2) in vec3 aNormal;
uniform sampler2D uField;
uniform int uGridSize;
uniform float uPlaneLength;
uniform int uRenderSphere;
uniform vec3 uSphereCenter;
uniform float uSphereRadius;
uniform int uSphereMaterial;
uniform mat4 uViewProjection;
out vec3 vWorld;
out vec3 vNormal;
out float vHeight;
flat out int vSphereMaterial;
float TerrainHeight(vec2 p)
{
    return 0.105*sin(p.x*0.62+p.y*0.13)*cos(p.y*0.57)
        +0.045*sin(p.x*1.04-p.y*0.78);
}
void main()
{
    vec3 world;
    vec3 normal;
    float height;
    if(uRenderSphere==0)
    {
        ivec2 texel=clamp(ivec2(round(aUv*float(uGridSize-1))),ivec2(0),
            ivec2(uGridSize-1));
        ivec2 left=ivec2(max(texel.x-1,0),texel.y);
        ivec2 right=ivec2(min(texel.x+1,uGridSize-1),texel.y);
        ivec2 down=ivec2(texel.x,max(texel.y-1,0));
        ivec2 up=ivec2(texel.x,min(texel.y+1,uGridSize-1));
        float sampleDx=uPlaneLength/float(uGridSize);
        vec2 leftWorld=(vec2(left)+vec2(0.5))*sampleDx-vec2(uPlaneLength*0.5);
        vec2 rightWorld=(vec2(right)+vec2(0.5))*sampleDx-vec2(uPlaneLength*0.5);
        vec2 downWorld=(vec2(down)+vec2(0.5))*sampleDx-vec2(uPlaneLength*0.5);
        vec2 upWorld=(vec2(up)+vec2(0.5))*sampleDx-vec2(uPlaneLength*0.5);
        float leftHeight=texelFetch(uField,left,0).r+TerrainHeight(leftWorld);
        float rightHeight=texelFetch(uField,right,0).r+TerrainHeight(rightWorld);
        float downHeight=texelFetch(uField,down,0).r+TerrainHeight(downWorld);
        float upHeight=texelFetch(uField,up,0).r+TerrainHeight(upWorld);
        height=texelFetch(uField,texel,0).r+TerrainHeight(
            (vec2(texel)+vec2(0.5))*sampleDx-vec2(uPlaneLength*0.5));
        float dx=max(float(right.x-left.x)*uPlaneLength/float(uGridSize),1.0e-6);
        float dy=max(float(up.y-down.y)*uPlaneLength/float(uGridSize),1.0e-6);
        float dhdx=(rightHeight-leftHeight)/dx;
        float dhdy=(upHeight-downHeight)/dy;
        world=vec3((aUv-vec2(0.5))*uPlaneLength,height);
        normal=normalize(vec3(-dhdx,-dhdy,1.0));
    }
    else
    {
        world=uSphereCenter+aPosition*uSphereRadius;
        normal=normalize(aNormal);
        height=world.z;
    }
    vWorld=world;
    vNormal=normal;
    vHeight=height;
    vSphereMaterial=uRenderSphere==0?0:uSphereMaterial;
    gl_Position=uViewProjection*vec4(world,1.0);
})GLSL";
        static constexpr char FragmentSource[]=R"GLSL(#version 430 core
in vec3 vWorld;
in vec3 vNormal;
in float vHeight;
flat in int vSphereMaterial;
uniform vec3 uCameraPosition;
out vec4 outColor;
void main()
{
    vec3 n=normalize(vNormal);
    vec3 l=normalize(vec3(-0.62,-0.54,0.32));
    vec3 v=normalize(uCameraPosition-vWorld);
    float diffuse=max(dot(n,l),0.0);
    float slope=1.0-clamp(n.z,0.0,1.0);
    float spec=pow(max(dot(reflect(-l,n),v),0.0),vSphereMaterial==0?72.0:38.0);
    if(vSphereMaterial!=0)
    {
        vec3 base=vSphereMaterial==1?vec3(0.96,0.34,0.10):vec3(0.18,0.70,0.92);
        float light=0.24+0.84*diffuse;
        outColor=vec4(base*light+vec3(0.65,0.43,0.22)*spec,1.0);
    }
    else
    {
        float low= smoothstep(-0.46,-0.08,vHeight);
        float high=smoothstep(0.02,0.32,vHeight);
        vec3 base=mix(vec3(0.025,0.075,0.13),vec3(0.055,0.30,0.31),low);
        base=mix(base,vec3(0.28,0.46,0.28),high*0.82);
        float contourPhase=fract((vHeight+0.60)/0.075);
        float contourDistance=min(contourPhase,1.0-contourPhase);
        float contour=smoothstep(0.025,0.16,contourDistance);
        base=mix(base*0.35,base,contour);
        base+=vec3(0.20,0.22,0.12)*smoothstep(0.08,0.36,slope);
        float light=0.34+0.98*diffuse;
        outColor=vec4(base*light+vec3(0.30,0.52,0.48)*spec,1.0);
    }
})GLSL";
        const GLuint VS=CompileShader(GL_VERTEX_SHADER,VertexSource);
        const GLuint FS=CompileShader(GL_FRAGMENT_SHADER,FragmentSource);
        if(!VS||!FS)
        {
            if(VS)glDeleteShader(VS);
            if(FS)glDeleteShader(FS);
            return 0;
        }
        const std::array<GLuint,2> Shaders={VS,FS};
        const GLuint Program=LinkProgram(Shaders);
        glDeleteShader(VS);
        glDeleteShader(FS);
        return Program;
    }

    struct FVec3 { double X=0.0,Y=0.0,Z=0.0; };

    double TerrainHeight(double X,double Y)
    {
        return 0.105*std::sin(X*0.62+Y*0.13)*std::cos(Y*0.57)
            +0.045*std::sin(X*1.04-Y*0.78);
    }

    double Smooth01(double T)
    {
        T=std::clamp(T,0.0,1.0);
        return T*T*(3.0-2.0*T);
    }

    std::array<FBodyState,2> GetBodiesAt(double Time)
    {
        constexpr double RollEnd=1.62;
        constexpr double RestEnd=2.55;
        constexpr double LiftEnd=2.95;
        constexpr double SecondaryStart=2.70;
        constexpr double SecondaryEnd=5.55;

        FBodyState Primary;
        Primary.Y=-0.15;
        if(Time<RollEnd)
        {
            Primary.X=-3.55+4.05*Smooth01(Time/RollEnd);
            Primary.Z=TerrainHeight(Primary.X,Primary.Y)+0.41;
            Primary.Force=SpherePushMetersPerSecondSquared;
        }
        else if(Time<RestEnd)
        {
            const double RestTime=Time-RollEnd;
            Primary.X=0.50;
            Primary.Z=TerrainHeight(Primary.X,Primary.Y)+0.41
                +0.065*std::sin(RestTime*2.0*Pi*1.15);
            Primary.Force=SpherePushMetersPerSecondSquared
                *(0.82+0.18*std::sin(RestTime*2.0*Pi*1.15));
        }
        else if(Time<LiftEnd)
        {
            const double Lift=Smooth01((Time-RestEnd)/(LiftEnd-RestEnd));
            Primary.X=0.50;
            Primary.Z=TerrainHeight(Primary.X,Primary.Y)+0.41+1.10*Lift;
        }
        else
        {
            Primary.X=0.50;
            Primary.Z=TerrainHeight(Primary.X,Primary.Y)+1.51;
        }

        FBodyState Secondary;
        Secondary.Y=1.05+0.14*std::sin(Time*0.72);
        const double Travel=Smooth01((Time-SecondaryStart)/(SecondaryEnd-SecondaryStart));
        Secondary.X=5.15-10.0*Travel;
        Secondary.Z=TerrainHeight(Secondary.X,Secondary.Y)+0.41;
        if(Time>=SecondaryStart&&Time<=SecondaryEnd)
            Secondary.Force=SpherePushMetersPerSecondSquared*0.78;
        return {Primary,Secondary};
    }

    FVec3 Sub(FVec3 A,FVec3 B){return {A.X-B.X,A.Y-B.Y,A.Z-B.Z};}
    FVec3 Cross(FVec3 A,FVec3 B)
    {return {A.Y*B.Z-A.Z*B.Y,A.Z*B.X-A.X*B.Z,A.X*B.Y-A.Y*B.X};}
    double Dot(FVec3 A,FVec3 B){return A.X*B.X+A.Y*B.Y+A.Z*B.Z;}
    FVec3 Normalize(FVec3 V)
    {
        const double Length=std::sqrt(Dot(V,V));
        return Length>1.0e-12?FVec3{V.X/Length,V.Y/Length,V.Z/Length}:FVec3{};
    }

    std::array<float,16> MakeViewProjection(int Width,int Height,FVec3 Camera,FVec3 Target)
    {
        const FVec3 Forward=Normalize(Sub(Target,Camera));
        const FVec3 Side=Normalize(Cross(Forward,{0.0,0.0,1.0}));
        const FVec3 Up=Cross(Side,Forward);
        std::array<float,16> View{};
        View[0]=static_cast<float>(Side.X);View[4]=static_cast<float>(Side.Y);
        View[8]=static_cast<float>(Side.Z);View[12]=static_cast<float>(-Dot(Side,Camera));
        View[1]=static_cast<float>(Up.X);View[5]=static_cast<float>(Up.Y);
        View[9]=static_cast<float>(Up.Z);View[13]=static_cast<float>(-Dot(Up,Camera));
        View[2]=static_cast<float>(-Forward.X);View[6]=static_cast<float>(-Forward.Y);
        View[10]=static_cast<float>(-Forward.Z);View[14]=static_cast<float>(Dot(Forward,Camera));
        View[15]=1.0f;

        const double Near=0.1,Far=60.0;
        const double F=1.0/std::tan(39.0*Pi/360.0);
        const double Aspect=static_cast<double>(Width)/static_cast<double>(Height);
        std::array<float,16> Projection{};
        Projection[0]=static_cast<float>(F/Aspect);
        Projection[5]=static_cast<float>(F);
        Projection[10]=static_cast<float>((Far+Near)/(Near-Far));
        Projection[11]=-1.0f;
        Projection[14]=static_cast<float>((2.0*Far*Near)/(Near-Far));

        std::array<float,16> Result{};
        for(int Column=0;Column<4;++Column)for(int Row=0;Row<4;++Row)
            for(int K=0;K<4;++K)
                Result[static_cast<std::size_t>(Column*4+Row)]+=
                    Projection[static_cast<std::size_t>(K*4+Row)]
                    *View[static_cast<std::size_t>(Column*4+K)];
        return Result;
    }

    struct FGridVertex { float U=0.0f,V=0.0f; };
    struct FSphereVertex { float Position[3]{},Normal[3]{}; };

    void MakeGrid(int Segments,std::vector<FGridVertex>& Vertices,
                  std::vector<std::uint32_t>& Indices)
    {
        const int Side=Segments+1;
        Vertices.reserve(static_cast<std::size_t>(Side)*Side);
        for(int Y=0;Y<=Segments;++Y)for(int X=0;X<=Segments;++X)
            Vertices.push_back({static_cast<float>(X)/Segments,static_cast<float>(Y)/Segments});
        Indices.reserve(static_cast<std::size_t>(Segments)*Segments*6);
        for(int Y=0;Y<Segments;++Y)for(int X=0;X<Segments;++X)
        {
            const std::uint32_t A=static_cast<std::uint32_t>(Y*Side+X);
            const std::uint32_t B=A+1;
            const std::uint32_t D=A+static_cast<std::uint32_t>(Side);
            const std::uint32_t C=D+1;
            Indices.insert(Indices.end(),{A,B,C,A,C,D});
        }
    }

    void MakeSphere(std::vector<FSphereVertex>& Vertices,
                    std::vector<std::uint32_t>& Indices)
    {
        constexpr int Rings=24,Segments=40;
        for(int Ring=0;Ring<=Rings;++Ring)
        {
            const double Theta=Pi*static_cast<double>(Ring)/Rings;
            for(int Segment=0;Segment<=Segments;++Segment)
            {
                const double Phi=2.0*Pi*static_cast<double>(Segment)/Segments;
                const float X=static_cast<float>(std::sin(Theta)*std::cos(Phi));
                const float Y=static_cast<float>(std::sin(Theta)*std::sin(Phi));
                const float Z=static_cast<float>(std::cos(Theta));
                Vertices.push_back({{X,Y,Z},{X,Y,Z}});
            }
        }
        const int Side=Segments+1;
        for(int Ring=0;Ring<Rings;++Ring)for(int Segment=0;Segment<Segments;++Segment)
        {
            const std::uint32_t A=static_cast<std::uint32_t>(Ring*Side+Segment);
            const std::uint32_t B=A+static_cast<std::uint32_t>(Side);
            const std::uint32_t C=B+1;
            const std::uint32_t D=A+1;
            Indices.insert(Indices.end(),{A,B,C,A,C,D});
        }
    }

    bool WritePng(const std::filesystem::path& Path,int Width,int Height,
                  const std::vector<std::uint8_t>& Pixels)
    {
        std::error_code Error;
        std::filesystem::create_directories(Path.parent_path(),Error);
        if(Error)return false;
        FILE* File=std::fopen(Path.c_str(),"wb");
        if(!File)return false;
        png_structp Png=png_create_write_struct(PNG_LIBPNG_VER_STRING,nullptr,nullptr,nullptr);
        png_infop Info=Png?png_create_info_struct(Png):nullptr;
        if(!Png||!Info)
        {if(Png)png_destroy_write_struct(&Png,nullptr);std::fclose(File);return false;}
        if(setjmp(png_jmpbuf(Png)))
        {png_destroy_write_struct(&Png,&Info);std::fclose(File);return false;}
        png_init_io(Png,File);
        png_set_IHDR(Png,Info,static_cast<png_uint_32>(Width),static_cast<png_uint_32>(Height),
            8,PNG_COLOR_TYPE_RGBA,PNG_INTERLACE_NONE,PNG_COMPRESSION_TYPE_DEFAULT,
            PNG_FILTER_TYPE_DEFAULT);
        png_write_info(Png,Info);
        const std::size_t Stride=static_cast<std::size_t>(Width)*4;
        for(int Y=Height-1;Y>=0;--Y)
            png_write_row(Png,const_cast<png_bytep>(Pixels.data()+static_cast<std::size_t>(Y)*Stride));
        png_write_end(Png,nullptr);
        png_destroy_write_struct(&Png,&Info);
        return std::fclose(File)==0;
    }

    struct FRenderResources
    {
        GLuint Program=0,Framebuffer=0,ColorTexture=0,DepthBuffer=0;
        GLuint GridVao=0,GridVbo=0,GridEbo=0,SphereVao=0,SphereVbo=0,SphereEbo=0;
        GLsizei GridIndexCount=0,SphereIndexCount=0;
        int Width=0,Height=0,GridSize=0;
        std::array<float,16> ViewProjection{};
        FVec3 Camera{};

        void Destroy()
        {
            if(SphereEbo)glDeleteBuffers(1,&SphereEbo);
            if(SphereVbo)glDeleteBuffers(1,&SphereVbo);
            if(SphereVao)glDeleteVertexArrays(1,&SphereVao);
            if(GridEbo)glDeleteBuffers(1,&GridEbo);
            if(GridVbo)glDeleteBuffers(1,&GridVbo);
            if(GridVao)glDeleteVertexArrays(1,&GridVao);
            if(DepthBuffer)glDeleteRenderbuffers(1,&DepthBuffer);
            if(ColorTexture)glDeleteTextures(1,&ColorTexture);
            if(Framebuffer)glDeleteFramebuffers(1,&Framebuffer);
            if(Program)glDeleteProgram(Program);
            *this=FRenderResources{};
        }

        bool Initialize(int Size,int InWidth,int InHeight)
        {
            Width=InWidth;Height=InHeight;GridSize=Size;
            Program=MakeRenderProgram();
            if(!Program)return false;
            glGenTextures(1,&ColorTexture);
            glBindTexture(GL_TEXTURE_2D,ColorTexture);
            glTexImage2D(GL_TEXTURE_2D,0,GL_RGBA8,Width,Height,0,GL_RGBA,GL_UNSIGNED_BYTE,nullptr);
            glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER,GL_LINEAR);
            glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER,GL_LINEAR);
            glGenRenderbuffers(1,&DepthBuffer);
            glBindRenderbuffer(GL_RENDERBUFFER,DepthBuffer);
            glRenderbufferStorage(GL_RENDERBUFFER,GL_DEPTH_COMPONENT24,Width,Height);
            glGenFramebuffers(1,&Framebuffer);
            glBindFramebuffer(GL_FRAMEBUFFER,Framebuffer);
            glFramebufferTexture2D(GL_FRAMEBUFFER,GL_COLOR_ATTACHMENT0,GL_TEXTURE_2D,ColorTexture,0);
            glFramebufferRenderbuffer(GL_FRAMEBUFFER,GL_DEPTH_ATTACHMENT,GL_RENDERBUFFER,DepthBuffer);
            if(glCheckFramebufferStatus(GL_FRAMEBUFFER)!=GL_FRAMEBUFFER_COMPLETE)
            {std::cerr<<"[LivingPlane] offscreen framebuffer is incomplete\n";return false;}

            const int Segments=std::min(Size,512);
            std::vector<FGridVertex> GridVertices;
            std::vector<std::uint32_t> GridIndices;
            MakeGrid(Segments,GridVertices,GridIndices);
            GridIndexCount=static_cast<GLsizei>(GridIndices.size());
            glGenVertexArrays(1,&GridVao);glBindVertexArray(GridVao);
            glGenBuffers(1,&GridVbo);glBindBuffer(GL_ARRAY_BUFFER,GridVbo);
            glBufferData(GL_ARRAY_BUFFER,static_cast<GLsizeiptr>(GridVertices.size()*sizeof(FGridVertex)),
                GridVertices.data(),GL_STATIC_DRAW);
            glEnableVertexAttribArray(1);
            glVertexAttribPointer(1,2,GL_FLOAT,GL_FALSE,sizeof(FGridVertex),nullptr);
            glGenBuffers(1,&GridEbo);glBindBuffer(GL_ELEMENT_ARRAY_BUFFER,GridEbo);
            glBufferData(GL_ELEMENT_ARRAY_BUFFER,
                static_cast<GLsizeiptr>(GridIndices.size()*sizeof(std::uint32_t)),
                GridIndices.data(),GL_STATIC_DRAW);

            std::vector<FSphereVertex> SphereVertices;
            std::vector<std::uint32_t> SphereIndices;
            MakeSphere(SphereVertices,SphereIndices);
            SphereIndexCount=static_cast<GLsizei>(SphereIndices.size());
            glGenVertexArrays(1,&SphereVao);glBindVertexArray(SphereVao);
            glGenBuffers(1,&SphereVbo);glBindBuffer(GL_ARRAY_BUFFER,SphereVbo);
            glBufferData(GL_ARRAY_BUFFER,
                static_cast<GLsizeiptr>(SphereVertices.size()*sizeof(FSphereVertex)),
                SphereVertices.data(),GL_STATIC_DRAW);
            glEnableVertexAttribArray(0);
            glVertexAttribPointer(0,3,GL_FLOAT,GL_FALSE,sizeof(FSphereVertex),
                reinterpret_cast<void*>(offsetof(FSphereVertex,Position)));
            glEnableVertexAttribArray(2);
            glVertexAttribPointer(2,3,GL_FLOAT,GL_FALSE,sizeof(FSphereVertex),
                reinterpret_cast<void*>(offsetof(FSphereVertex,Normal)));
            glGenBuffers(1,&SphereEbo);glBindBuffer(GL_ELEMENT_ARRAY_BUFFER,SphereEbo);
            glBufferData(GL_ELEMENT_ARRAY_BUFFER,
                static_cast<GLsizeiptr>(SphereIndices.size()*sizeof(std::uint32_t)),
                SphereIndices.data(),GL_STATIC_DRAW);

            glBindVertexArray(0);
            return glGetError()==GL_NO_ERROR;
        }

        bool Render(GLuint FieldTexture,const std::filesystem::path& Output,
                    const FBodyState& Primary,const FBodyState& Secondary)
        {
            const FVec3 Target={Primary.X,Primary.Y,0.0};
            Camera={Primary.X+2.65,Primary.Y-4.45,3.65};
            ViewProjection=MakeViewProjection(Width,Height,Camera,Target);
            glBindFramebuffer(GL_FRAMEBUFFER,Framebuffer);
            glViewport(0,0,Width,Height);
            glDisable(GL_CULL_FACE);
            glEnable(GL_DEPTH_TEST);
            glDepthFunc(GL_LEQUAL);
            glClearColor(0.018f,0.055f,0.095f,1.0f);
            glClearDepth(1.0);
            glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
            glUseProgram(Program);
            glUniform1i(glGetUniformLocation(Program,"uField"),0);
            glUniform1i(glGetUniformLocation(Program,"uGridSize"),GridSize);
            glUniform1f(glGetUniformLocation(Program,"uPlaneLength"),
                static_cast<float>(PlaneLengthMeters));
            glUniformMatrix4fv(glGetUniformLocation(Program,"uViewProjection"),1,GL_FALSE,
                ViewProjection.data());
            glUniform3f(glGetUniformLocation(Program,"uCameraPosition"),
                static_cast<float>(Camera.X),static_cast<float>(Camera.Y),static_cast<float>(Camera.Z));
            glActiveTexture(GL_TEXTURE0);
            glBindTexture(GL_TEXTURE_2D,FieldTexture);
            glUniform1i(glGetUniformLocation(Program,"uRenderSphere"),0);
            glBindVertexArray(GridVao);
            glDrawElements(GL_TRIANGLES,GridIndexCount,GL_UNSIGNED_INT,nullptr);
            glUniform1i(glGetUniformLocation(Program,"uRenderSphere"),1);
            glUniform1f(glGetUniformLocation(Program,"uSphereRadius"),
                static_cast<float>(SphereRadiusMeters));
            glBindVertexArray(SphereVao);
            glUniform1i(glGetUniformLocation(Program,"uSphereMaterial"),1);
            glUniform3f(glGetUniformLocation(Program,"uSphereCenter"),
                static_cast<float>(Primary.X),static_cast<float>(Primary.Y),
                static_cast<float>(Primary.Z));
            glDrawElements(GL_TRIANGLES,SphereIndexCount,GL_UNSIGNED_INT,nullptr);
            glUniform1i(glGetUniformLocation(Program,"uSphereMaterial"),2);
            glUniform3f(glGetUniformLocation(Program,"uSphereCenter"),
                static_cast<float>(Secondary.X),static_cast<float>(Secondary.Y),
                static_cast<float>(Secondary.Z));
            glDrawElements(GL_TRIANGLES,SphereIndexCount,GL_UNSIGNED_INT,nullptr);
            glFinish();
            std::vector<std::uint8_t> Pixels(static_cast<std::size_t>(Width)*Height*4);
            glPixelStorei(GL_PACK_ALIGNMENT,1);
            glReadPixels(0,0,Width,Height,GL_RGBA,GL_UNSIGNED_BYTE,Pixels.data());
            if(glGetError()!=GL_NO_ERROR||!WritePng(Output,Width,Height,Pixels))
            {std::cerr<<"[LivingPlane] frame render/write failed: "<<Output<<"\n";return false;}
            return true;
        }
    };

    bool ParseInt(std::string_view Text,int& Value)
    {
        if(Text.empty())return false;
        char* End=nullptr;
        errno=0;
        const std::string Copy(Text);
        const long Parsed=std::strtol(Copy.c_str(),&End,10);
        if(errno||!End||*End!='\0'||Parsed<std::numeric_limits<int>::min()
            ||Parsed>std::numeric_limits<int>::max())return false;
        Value=static_cast<int>(Parsed);
        return true;
    }

    bool ParseDouble(std::string_view Text,double& Value)
    {
        if(Text.empty())return false;
        char* End=nullptr;
        errno=0;
        const std::string Copy(Text);
        const double Parsed=std::strtod(Copy.c_str(),&End);
        if(errno||!End||*End!='\0'||!std::isfinite(Parsed))return false;
        Value=Parsed;
        return true;
    }

    bool ParseOptions(int Argc,char** Argv,FOptions& Options)
    {
        for(int I=1;I<Argc;++I)
        {
            const std::string_view Arg(Argv[I]);
            const auto Take=[&](std::string_view Key,std::string_view& Value)->bool
            {
                if(Arg.starts_with(Key)&&Arg.size()>Key.size()&&Arg[Key.size()]=='=')
                {Value=Arg.substr(Key.size()+1);return true;}
                return false;
            };
            std::string_view Value;
            if(Arg=="--help"||Arg=="-h")
            {
                std::cout<<"Usage: living-plane [--device=auto|gpu|cpu] [--scene=wake|alive] "
                    "[--size=N] [--steps=N] [--out=DIR] [--capture-every=N]\n"
                    "  [--queries=on|off] [--breathing=on|off] [--imprint=on|off]\n"
                    "  [--breathing-strength=MPS2] [--imprint-capture=MPS] "
                    "[--imprint-recovery=SECONDS] [--imprint-depth=M]\n"
                    "  --capture-every=0 writes only timing.json. --queries=off measures "
                    "batched dispatch wall time without per-step query overhead.\n";
                std::exit(0);
            }
            else if(Take("--device",Value))Options.Device=Value;
            else if(Take("--scene",Value))Options.Scene=Value;
            else if(Take("--size",Value))
            {if(!ParseInt(Value,Options.Size)){std::cerr<<"invalid --size\n";return false;}}
            else if(Take("--steps",Value))
            {if(!ParseInt(Value,Options.Steps)){std::cerr<<"invalid --steps\n";return false;}}
            else if(Take("--capture-every",Value))
            {if(!ParseInt(Value,Options.CaptureEvery)){std::cerr<<"invalid --capture-every\n";return false;}}
            else if(Take("--out",Value))Options.Output=Value;
            else if(Take("--queries",Value))
            {
                if(Value=="on")Options.TimerQueries=true;
                else if(Value=="off")Options.TimerQueries=false;
                else {std::cerr<<"--queries must be on or off\n";return false;}
            }
            else if(Take("--breathing",Value))
            {
                if(Value=="on")Options.Breathing=true;
                else if(Value=="off")Options.Breathing=false;
                else {std::cerr<<"--breathing must be on or off\n";return false;}
                Options.BreathingConfigured=true;
            }
            else if(Take("--imprint",Value))
            {
                if(Value=="on")Options.Imprint=true;
                else if(Value=="off")Options.Imprint=false;
                else {std::cerr<<"--imprint must be on or off\n";return false;}
                Options.ImprintConfigured=true;
            }
            else if(Take("--breathing-strength",Value))
            {
                if(!ParseDouble(Value,Options.BreathingStrength))
                {std::cerr<<"invalid --breathing-strength\n";return false;}
            }
            else if(Take("--imprint-capture",Value))
            {
                if(!ParseDouble(Value,Options.ImprintCaptureMetersPerSecond))
                {std::cerr<<"invalid --imprint-capture\n";return false;}
            }
            else if(Take("--imprint-recovery",Value))
            {
                if(!ParseDouble(Value,Options.ImprintRecoverySeconds))
                {std::cerr<<"invalid --imprint-recovery\n";return false;}
            }
            else if(Take("--imprint-depth",Value))
            {
                if(!ParseDouble(Value,Options.ImprintMaxDepthMeters))
                {std::cerr<<"invalid --imprint-depth\n";return false;}
            }
            else {std::cerr<<"unknown option: "<<Arg<<"\n";return false;}
        }
        if(Options.Device!="auto"&&Options.Device!="gpu"&&Options.Device!="cpu")
        {std::cerr<<"--device must be auto, gpu, or cpu\n";return false;}
        if(Options.Scene!="wake"&&Options.Scene!="alive")
        {std::cerr<<"--scene must be wake or alive\n";return false;}
        if(!Options.BreathingConfigured)Options.Breathing=Options.Scene=="alive";
        if(!Options.ImprintConfigured)Options.Imprint=Options.Scene=="alive";
        if(Options.Size<16||Options.Size>4096)
        {std::cerr<<"--size must be in [16,4096]\n";return false;}
        if(Options.Steps<1||Options.Steps>100000)
        {std::cerr<<"--steps must be in [1,100000]\n";return false;}
        if(Options.CaptureEvery<0||Options.CaptureEvery>100000)
        {std::cerr<<"--capture-every must be in [0,100000]\n";return false;}
        if(Options.BreathingStrength<0.0||Options.BreathingStrength>30.0
            ||Options.ImprintCaptureMetersPerSecond<0.0
            ||Options.ImprintCaptureMetersPerSecond>1.0
            ||Options.ImprintRecoverySeconds<0.1||Options.ImprintRecoverySeconds>60.0
            ||Options.ImprintMaxDepthMeters<0.0||Options.ImprintMaxDepthMeters>1.0)
        {std::cerr<<"alive-mode parameters are outside their supported ranges\n";return false;}
        return true;
    }

    double Quantile(std::vector<double> Values,double Q)
    {
        if(Values.empty())return 0.0;
        std::sort(Values.begin(),Values.end());
        const std::size_t I=static_cast<std::size_t>(std::ceil(Q*Values.size())-1.0);
        return Values[std::min(I,Values.size()-1)];
    }

    std::uint64_t FieldChecksum(const std::vector<float>& Field)
    {
        // Hash float bit patterns in canonical little-endian byte order so the
        // result does not depend on the host's native byte order.
        std::uint64_t Hash=14695981039346656037ull;
        for(const float Value:Field)
        {
            const std::uint32_t Bits=std::bit_cast<std::uint32_t>(Value);
            for(unsigned Byte=0;Byte<4;++Byte)
            {
                Hash^=static_cast<std::uint8_t>((Bits>>(Byte*8))&0xffu);
                Hash*=1099511628211ull;
            }
        }
        return Hash;
    }

    bool WriteTimingJson(const FOptions& Options,const FContext& Context,
        std::string_view Renderer,std::string_view Vendor,std::string_view Version,
        double Dt,const std::vector<GLuint64>& Times,const std::vector<double>& Heights,
        const std::vector<double>& Velocities,std::size_t Frames,int TimerQueryBits,
        int TimerQueryWaits,bool TimerResultsAvailable,double DispatchWallMsPerStep,
        std::uint64_t FieldChecksumValue,bool FirstStepAudit,bool DynamicsAudit)
    {
        std::error_code Error;
        std::filesystem::create_directories(Options.Output,Error);
        if(Error){std::cerr<<"cannot create output directory: "<<Error.message()<<"\n";return false;}
        std::ofstream Json(Options.Output/"timing.json",std::ios::trunc);
        if(!Json){std::cerr<<"cannot write timing.json\n";return false;}
        std::vector<double> Ms;
        Ms.reserve(Times.size());
        for(const GLuint64 Ns:Times)Ms.push_back(static_cast<double>(Ns)*1.0e-6);
        const double Mean=Ms.empty()?0.0:std::accumulate(Ms.begin(),Ms.end(),0.0)/Ms.size();
        const double Min=Ms.empty()?0.0:*std::min_element(Ms.begin(),Ms.end());
        const double Max=Ms.empty()?0.0:*std::max_element(Ms.begin(),Ms.end());
        const double HeightMin=*std::min_element(Heights.begin(),Heights.end());
        const double HeightMax=*std::max_element(Heights.begin(),Heights.end());
        const double HeightMean=std::accumulate(Heights.begin(),Heights.end(),0.0)/Heights.size();
        double HeightSquare=0.0;
        for(const double H:Heights)HeightSquare+=H*H;
        const double HeightRms=std::sqrt(HeightSquare/Heights.size());
        char FieldChecksumHex[17]{};
        std::snprintf(FieldChecksumHex,sizeof(FieldChecksumHex),"%016llx",
            static_cast<unsigned long long>(FieldChecksumValue));
        const double Dx=PlaneLengthMeters/Options.Size;
        const double SpectralStep=Dt*std::sqrt(8.0*WaveSpeedMetersPerSecond*
            WaveSpeedMetersPerSecond/(Dx*Dx)+FoundationStiffnessPerSecondSquared);
        const double VelocityMax=*std::max_element(Velocities.begin(),Velocities.end(),
            [](double A,double B){return std::abs(A)<std::abs(B);});

        Json<<"{\n  \"schema\": 1,\n"
            <<"  \"requested_device\": \""<<JsonEscape(Options.Device)<<"\",\n"
            <<"  \"execution_device\": \""<<Context.ExecutionDevice<<"\",\n"
            <<"  \"egl_device\": \""<<JsonEscape(Context.DeviceLabel)<<"\",\n"
            <<"  \"drm_render_node\": \""<<JsonEscape(Context.RenderNode)<<"\",\n"
            <<"  \"gl_vendor\": \""<<JsonEscape(Vendor)<<"\",\n"
            <<"  \"gl_renderer\": \""<<JsonEscape(Renderer)<<"\",\n"
            <<"  \"gl_version\": \""<<JsonEscape(Version)<<"\",\n"
            <<"  \"size\": "<<Options.Size<<",\n"
            <<"  \"steps\": "<<Options.Steps<<",\n"
            <<"  \"scene\": \""<<Options.Scene<<"\",\n"
            <<"  \"plane_length_m\": "<<PlaneLengthMeters<<",\n"
            <<"  \"dt_seconds\": "<<Dt<<",\n"
            <<"  \"spectral_step_ratio\": "<<SpectralStep<<",\n"
            <<"  \"wave_speed_m_per_s\": "<<WaveSpeedMetersPerSecond<<",\n"
            <<"  \"damping_per_s\": "<<DampingPerSecond<<",\n"
            <<"  \"foundation_stiffness_per_s2\": "<<FoundationStiffnessPerSecondSquared<<",\n"
            <<"  \"breathing\": "<<(Options.Breathing?"true":"false")<<",\n"
            <<"  \"breathing_strength_m_per_s2\": "<<Options.BreathingStrength<<",\n"
            <<"  \"imprint\": "<<(Options.Imprint?"true":"false")<<",\n"
            <<"  \"imprint_capture_m_per_s\": "<<Options.ImprintCaptureMetersPerSecond<<",\n"
            <<"  \"imprint_recovery_seconds\": "<<Options.ImprintRecoverySeconds<<",\n"
            <<"  \"imprint_max_depth_m\": "<<Options.ImprintMaxDepthMeters<<",\n"
            <<"  \"frames_written\": "<<Frames<<",\n"
            <<"  \"first_step_equation_audit\": \""<<(FirstStepAudit?"pass":"fail")<<"\",\n"
            <<"  \"field_dynamics_audit\": \""<<(DynamicsAudit?"pass":"fail")<<"\",\n"
            <<"  \"timer_queries_enabled\": "<<(Options.TimerQueries?"true":"false")<<",\n"
            <<"  \"timer_query_ring_size\": "<<(Options.TimerQueries
                ?std::min<std::size_t>(TimerQueryRingSize,Times.size()):0)<<",\n"
            <<"  \"timer_query_counter_bits\": "<<TimerQueryBits<<",\n"
            <<"  \"timer_query_results_available_after_completion\": ";
        if(Options.TimerQueries)Json<<(TimerResultsAvailable?"true":"false");
        else Json<<"null";
        Json<<",\n"
            <<"  \"timer_query_waits_before_reuse\": "<<TimerQueryWaits<<",\n"
            <<"  \"dispatch_wall_ms_per_step\": ";
        if(DispatchWallMsPerStep>=0.0)Json<<DispatchWallMsPerStep;
        else Json<<"null";
        Json<<",\n  \"dispatch_wall_window\": \""
            <<(Options.CaptureEvery==0
                ?"simulation loop plus final glFinish; no frame capture; first-step equation readback included"
                :"not suitable for step cost because frame capture is enabled")<<"\",\n"
            <<"  \"field_min_m\": "<<HeightMin<<",\n"
            <<"  \"field_max_m\": "<<HeightMax<<",\n"
            <<"  \"field_mean_m\": "<<HeightMean<<",\n"
            <<"  \"field_rms_m\": "<<HeightRms<<",\n"
            <<"  \"field_checksum_fnv1a64\": \""<<FieldChecksumHex<<"\",\n"
            <<"  \"max_abs_velocity_m_per_s\": "<<std::abs(VelocityMax)<<",\n"
            <<"  \"gl_time_elapsed_ms_per_step\": [";
        for(std::size_t I=0;I<Times.size();++I)
        {
            if(I)Json<<",";
            if(I%8==0)Json<<"\n    ";
            Json<<static_cast<double>(Times[I])*1.0e-6;
        }
        Json<<"\n  ],\n  \"gl_time_elapsed_summary_ms\": {\n"
            <<"    \"mean\": "<<Mean<<",\n"
            <<"    \"median\": "<<Quantile(Ms,0.50)<<",\n"
            <<"    \"p95\": "<<Quantile(Ms,0.95)<<",\n"
            <<"    \"min\": "<<Min<<",\n"
            <<"    \"max\": "<<Max<<"\n  }\n}\n";
        return static_cast<bool>(Json);
    }

    bool Run(const FOptions& Options)
    {
        FContext Context;
        if(!CreateContext(Options.Device,Context))return false;
        const auto* RendererPtr=reinterpret_cast<const char*>(glGetString(GL_RENDERER));
        const auto* VendorPtr=reinterpret_cast<const char*>(glGetString(GL_VENDOR));
        const auto* VersionPtr=reinterpret_cast<const char*>(glGetString(GL_VERSION));
        if(!RendererPtr||!VendorPtr||!VersionPtr)
        {std::cerr<<"[LivingPlane] OpenGL identity query failed\n";return false;}
        const std::string Renderer=RendererPtr,Vendor=VendorPtr,Version=VersionPtr;
        if(Options.Device=="gpu"&&Context.ExecutionDevice!="gpu")
        {std::cerr<<"[LivingPlane] refused software renderer for requested GPU run: "
            <<Renderer<<"\n";return false;}
        if(Options.Device=="cpu"&&Context.ExecutionDevice!="cpu")
        {std::cerr<<"[LivingPlane] refused non-software renderer for requested CPU run: "
            <<Renderer<<"\n";return false;}
        std::cout<<"[LivingPlane] device="<<Context.ExecutionDevice<<" renderer=\""
            <<Renderer<<"\" version=\""<<Version<<"\" egl_device=\""
            <<Context.DeviceLabel<<"\" render_node=\""<<Context.RenderNode<<"\"\n";

        const GLuint ComputeProgram=MakeComputeProgram();
        if(!ComputeProgram)return false;
        const int Size=Options.Size;
        const double Dx=PlaneLengthMeters/Size;
        const double Dt=1.30/std::sqrt(8.0*WaveSpeedMetersPerSecond*
            WaveSpeedMetersPerSecond/(Dx*Dx)+FoundationStiffnessPerSecondSquared);
        std::array<GLuint,2> Fields{};
        glGenTextures(2,Fields.data());
        const std::vector<float> Zeroes(static_cast<std::size_t>(Size)*Size*4,0.0f);
        for(const GLuint Field:Fields)
        {
            glBindTexture(GL_TEXTURE_2D,Field);
            glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER,GL_NEAREST);
            glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER,GL_NEAREST);
            glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_WRAP_S,GL_CLAMP_TO_EDGE);
            glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_WRAP_T,GL_CLAMP_TO_EDGE);
            glTexImage2D(GL_TEXTURE_2D,0,GL_RGBA32F,Size,Size,0,GL_RGBA,GL_FLOAT,Zeroes.data());
        }

        FRenderResources Render;
        if(Options.CaptureEvery>0&&!Render.Initialize(Size,Size,Size))
        {Render.Destroy();glDeleteTextures(2,Fields.data());glDeleteProgram(ComputeProgram);return false;}

        const int QueryCount=Options.TimerQueries
            ?std::min<int>(static_cast<int>(TimerQueryRingSize),Options.Steps):0;
        std::vector<GLuint> Queries(static_cast<std::size_t>(QueryCount));
        if(QueryCount>0)glGenQueries(QueryCount,Queries.data());
        GLint TimerQueryBits=0;
        if(Options.TimerQueries)glGetQueryiv(GL_TIME_ELAPSED,GL_QUERY_COUNTER_BITS,&TimerQueryBits);

        const GLint LocationGrid=glGetUniformLocation(ComputeProgram,"uGridSize");
        const GLint LocationDt=glGetUniformLocation(ComputeProgram,"uDt");
        const GLint LocationDx=glGetUniformLocation(ComputeProgram,"uDx");
        const GLint LocationPlane=glGetUniformLocation(ComputeProgram,"uPlaneLength");
        const GLint LocationWaveSpeed=glGetUniformLocation(ComputeProgram,"uWaveSpeed");
        const GLint LocationDamping=glGetUniformLocation(ComputeProgram,"uDamping");
        const GLint LocationFoundation=glGetUniformLocation(ComputeProgram,"uFoundationStiffness");
        const GLint LocationTime=glGetUniformLocation(ComputeProgram,"uTime");
        const GLint LocationForcePositions=glGetUniformLocation(ComputeProgram,"uForcePositions[0]");
        const GLint LocationForceSigma=glGetUniformLocation(ComputeProgram,"uForceSigma");
        const GLint LocationForceStrengths=glGetUniformLocation(ComputeProgram,"uForceStrengths[0]");
        const GLint LocationBreathing=glGetUniformLocation(ComputeProgram,"uBreathing");
        const GLint LocationBreathingStrength=glGetUniformLocation(ComputeProgram,"uBreathingStrength");
        const GLint LocationImprint=glGetUniformLocation(ComputeProgram,"uImprint");
        const GLint LocationImprintCapture=glGetUniformLocation(ComputeProgram,"uImprintCaptureRate");
        const GLint LocationImprintRecovery=glGetUniformLocation(ComputeProgram,"uImprintRecoverySeconds");
        const GLint LocationImprintDepth=glGetUniformLocation(ComputeProgram,"uImprintMaxDepth");
        if(LocationGrid<0||LocationDt<0||LocationDx<0||LocationPlane<0||LocationWaveSpeed<0
            ||LocationDamping<0||LocationFoundation<0||LocationTime<0
            ||LocationForcePositions<0||LocationForceSigma<0||LocationForceStrengths<0
            ||LocationBreathing<0||LocationBreathingStrength<0||LocationImprint<0
            ||LocationImprintCapture<0||LocationImprintRecovery<0||LocationImprintDepth<0)
        {std::cerr<<"[LivingPlane] missing compute shader uniforms\n";return false;}

        glUseProgram(ComputeProgram);
        glUniform1i(LocationGrid,Size);
        glUniform1f(LocationDt,static_cast<float>(Dt));
        glUniform1f(LocationDx,static_cast<float>(Dx));
        glUniform1f(LocationPlane,static_cast<float>(PlaneLengthMeters));
        glUniform1f(LocationWaveSpeed,static_cast<float>(WaveSpeedMetersPerSecond));
        glUniform1f(LocationDamping,static_cast<float>(DampingPerSecond));
        glUniform1f(LocationFoundation,static_cast<float>(FoundationStiffnessPerSecondSquared));
        glUniform1f(LocationForceSigma,static_cast<float>(SphereSigmaMeters));
        glUniform1i(LocationBreathing,Options.Breathing?1:0);
        glUniform1f(LocationBreathingStrength,static_cast<float>(Options.BreathingStrength));
        glUniform1i(LocationImprint,Options.Imprint?1:0);
        glUniform1f(LocationImprintCapture,
            static_cast<float>(Options.ImprintCaptureMetersPerSecond));
        glUniform1f(LocationImprintRecovery,static_cast<float>(Options.ImprintRecoverySeconds));
        glUniform1f(LocationImprintDepth,static_cast<float>(Options.ImprintMaxDepthMeters));

        int Current=0;
        const auto Dispatch=[&](int Target,double Time,
                                const std::array<FBodyState,2>& Bodies,GLuint Query)
        {
            glUseProgram(ComputeProgram);
            glUniform1f(LocationTime,static_cast<float>(Time));
            const std::array<float,4> Positions={
                static_cast<float>(Bodies[0].X),static_cast<float>(Bodies[0].Y),
                static_cast<float>(Bodies[1].X),static_cast<float>(Bodies[1].Y)};
            const std::array<float,2> Strengths={
                static_cast<float>(Bodies[0].Force),static_cast<float>(Bodies[1].Force)};
            glUniform2fv(LocationForcePositions,2,Positions.data());
            glUniform1fv(LocationForceStrengths,2,Strengths.data());
            glBindImageTexture(0,Fields[static_cast<std::size_t>(Current)],0,GL_FALSE,0,
                GL_READ_ONLY,GL_RGBA32F);
            glBindImageTexture(1,Fields[static_cast<std::size_t>(Target)],0,GL_FALSE,0,
                GL_WRITE_ONLY,GL_RGBA32F);
            if(Query)glBeginQuery(GL_TIME_ELAPSED,Query);
            glDispatchCompute(static_cast<GLuint>((Size+15)/16),
                static_cast<GLuint>((Size+15)/16),1);
            if(Query)glEndQuery(GL_TIME_ELAPSED);
            glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT|GL_TEXTURE_FETCH_BARRIER_BIT);
            Current=Target;
        };

        // Warm the runtime-compiled compute kernel before measurement.
        const std::array<FBodyState,2> NoBodies{};
        for(int Warm=0;Warm<8;++Warm)Dispatch(1-Current,-8.0*Dt,NoBodies,0);
        glFinish();
        for(const GLuint Field:Fields)
        {
            glBindTexture(GL_TEXTURE_2D,Field);
            glTexSubImage2D(GL_TEXTURE_2D,0,0,0,Size,Size,GL_RGBA,GL_FLOAT,Zeroes.data());
        }
        glMemoryBarrier(GL_TEXTURE_UPDATE_BARRIER_BIT|GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);
        Current=0;

        FRenderResources* RenderPtr=Options.CaptureEvery>0?&Render:nullptr;
        std::size_t FramesWritten=0;
        int FrameIndex=0;
        auto SaveFrame=[&](double Time)->bool
        {
            if(!RenderPtr)return true;
            char Filename[32]{};
            std::snprintf(Filename,sizeof(Filename),"frame_%06d.png",FrameIndex++);
            const auto Bodies=GetBodiesAt(Time);
            if(!Render.Render(Fields[static_cast<std::size_t>(Current)],
                Options.Output/Filename,Bodies[0],Bodies[1]))return false;
            ++FramesWritten;
            return true;
        };

        if(!SaveFrame(0.0))
        {Render.Destroy();if(QueryCount>0)glDeleteQueries(QueryCount,Queries.data());
            glDeleteTextures(2,Fields.data());glDeleteProgram(ComputeProgram);return false;}

        std::vector<GLuint64> Times(Options.TimerQueries
            ?static_cast<std::size_t>(Options.Steps):0);
        int TimerQueryWaits=0;
        bool TimerResultsAvailable=Options.TimerQueries;
        auto ReadTimerQuery=[&](int QueryStep)
        {
            const GLuint Query=Queries[static_cast<std::size_t>(QueryStep%QueryCount)];
            GLint Available=GL_FALSE;
            glGetQueryObjectiv(Query,GL_QUERY_RESULT_AVAILABLE,&Available);
            if(Available!=GL_TRUE)++TimerQueryWaits;
            glGetQueryObjectui64v(Query,GL_QUERY_RESULT,
                &Times[static_cast<std::size_t>(QueryStep)]);
            glGetQueryObjectiv(Query,GL_QUERY_RESULT_AVAILABLE,&Available);
            TimerResultsAvailable=TimerResultsAvailable&&Available==GL_TRUE;
        };

        std::vector<float> FirstField;
        bool FirstStepAudit=false;
        const auto WallStart=std::chrono::steady_clock::now();
        for(int Step=0;Step<Options.Steps;++Step)
        {
            if(QueryCount>0&&Step>=QueryCount)ReadTimerQuery(Step-QueryCount);
            const double Time=static_cast<double>(Step)*Dt;
            const auto Bodies=GetBodiesAt(Time);
            const int Target=1-Current;
            const GLuint Query=QueryCount>0
                ?Queries[static_cast<std::size_t>(Step%QueryCount)]:0;
            Dispatch(Target,Time,Bodies,Query);
            if(Step==0)
            {
                glFinish();
                FirstField.resize(static_cast<std::size_t>(Size)*Size*4);
                glBindTexture(GL_TEXTURE_2D,Fields[static_cast<std::size_t>(Current)]);
                // glGetTexImage is a texture-update access. The image/fetch
                // barriers used for simulation and rendering do not cover this
                // readback visibility path.
                glMemoryBarrier(GL_TEXTURE_UPDATE_BARRIER_BIT);
                glGetTexImage(GL_TEXTURE_2D,0,GL_RGBA,GL_FLOAT,FirstField.data());
                const int X=std::clamp(static_cast<int>(std::floor(
                    (Bodies[0].X+PlaneLengthMeters*0.5)/Dx)),0,Size-1);
                const int Y=std::clamp(static_cast<int>(std::floor(
                    (Bodies[0].Y+PlaneLengthMeters*0.5)/Dx)),0,Size-1);
                const double CellX=(X+0.5)*Dx-PlaneLengthMeters*0.5;
                const double CellY=(Y+0.5)*Dx-PlaneLengthMeters*0.5;
                double Force=0.0;
                for(const FBodyState& Body:Bodies)
                {
                    const double Distance2=(CellX-Body.X)*(CellX-Body.X)
                        +(CellY-Body.Y)*(CellY-Body.Y);
                    Force-=Body.Force*std::exp(-Distance2/
                        (2.0*SphereSigmaMeters*SphereSigmaMeters));
                }
                const double Expected=Force*Dt*Dt*std::exp(-DampingPerSecond*Dt);
                const std::size_t Index=(static_cast<std::size_t>(Y)*Size+X)*4;
                const double Actual=FirstField[Index];
                const double Tolerance=std::max(2.0e-7,std::abs(Expected)*2.5e-3);
                FirstStepAudit=std::isfinite(Actual)&&std::abs(Actual-Expected)<=Tolerance;
                if(!FirstStepAudit)
                    std::cerr<<"[LivingPlane] first-step equation audit failed: expected "
                        <<Expected<<" actual "<<Actual<<" tolerance "<<Tolerance<<"\n";
            }
            if(RenderPtr&&(Options.CaptureEvery==1
                ||((Step+1)%Options.CaptureEvery==0)||Step==Options.Steps-1))
            {
                if(!SaveFrame(static_cast<double>(Step+1)*Dt))
                {Render.Destroy();if(QueryCount>0)glDeleteQueries(QueryCount,Queries.data());
                    glDeleteTextures(2,Fields.data());glDeleteProgram(ComputeProgram);return false;}
            }
        }
        glFinish();
        const auto WallEnd=std::chrono::steady_clock::now();
        const double WallMs=std::chrono::duration<double,std::milli>(WallEnd-WallStart).count();
        const double DispatchWallMsPerStep=Options.CaptureEvery==0
            ?WallMs/Options.Steps:-1.0;

        if(QueryCount>0)
            for(int Step=std::max(0,Options.Steps-QueryCount);Step<Options.Steps;++Step)
                ReadTimerQuery(Step);
        const bool QueryAudit=!Options.TimerQueries||TimerResultsAvailable;

        std::vector<float> FinalField(static_cast<std::size_t>(Size)*Size*4);
        glBindTexture(GL_TEXTURE_2D,Fields[static_cast<std::size_t>(Current)]);
        glMemoryBarrier(GL_TEXTURE_UPDATE_BARRIER_BIT);
        glGetTexImage(GL_TEXTURE_2D,0,GL_RGBA,GL_FLOAT,FinalField.data());
        std::vector<double> Heights,Velocities;
        Heights.reserve(static_cast<std::size_t>(Size)*Size);
        Velocities.reserve(static_cast<std::size_t>(Size)*Size);
        bool Finite=true;
        double ImprintMin=0.0;
        for(std::size_t I=0;I<FinalField.size();I+=4)
        {
            Finite=Finite&&std::isfinite(FinalField[I])&&std::isfinite(FinalField[I+1])
                &&std::isfinite(FinalField[I+2]);
            Heights.push_back(FinalField[I]);
            Velocities.push_back(FinalField[I+1]);
            ImprintMin=std::min(ImprintMin,static_cast<double>(FinalField[I+2]));
        }
        const auto [HeightMinIt,HeightMaxIt]=std::minmax_element(Heights.begin(),Heights.end());
        const bool DynamicsAudit=Finite&&*HeightMinIt<-1.0e-7&&*HeightMaxIt-*HeightMinIt<2.0
            &&(!Options.Imprint||ImprintMin>=-Options.ImprintMaxDepthMeters-1.0e-5);
        const std::uint64_t FinalFieldChecksum=FieldChecksum(FinalField);
        const bool Good=FirstStepAudit&&DynamicsAudit&&QueryAudit&&glGetError()==GL_NO_ERROR;
        if(!DynamicsAudit)
            std::cerr<<"[LivingPlane] field audit failed: finite="<<Finite<<" min="
                <<*HeightMinIt<<" max="<<*HeightMaxIt<<" imprint_min="<<ImprintMin<<"\n";
        if(!QueryAudit)
            std::cerr<<"[LivingPlane] timer-query result availability audit failed after glFinish\n";
        if(!WriteTimingJson(Options,Context,Renderer,Vendor,Version,Dt,Times,Heights,Velocities,
            FramesWritten,TimerQueryBits,TimerQueryWaits,TimerResultsAvailable,
            DispatchWallMsPerStep,FinalFieldChecksum,FirstStepAudit,DynamicsAudit))
        {Render.Destroy();if(QueryCount>0)glDeleteQueries(QueryCount,Queries.data());
            glDeleteTextures(2,Fields.data());glDeleteProgram(ComputeProgram);return false;}

        const std::vector<double> QueryMs=[&]
        {
            std::vector<double> Result;
            Result.reserve(Times.size());
            for(const GLuint64 Ns:Times)Result.push_back(static_cast<double>(Ns)*1.0e-6);
            return Result;
        }();
        std::cout<<"[LivingPlane] scene="<<Options.Scene<<" size="<<Size
            <<" steps="<<Options.Steps<<" dt_ms="<<Dt*1000.0
            <<" dispatch_wall_ms_per_step=";
        if(DispatchWallMsPerStep>=0.0)std::cout<<DispatchWallMsPerStep;
        else std::cout<<"n/a(capture enabled)";
        if(Options.TimerQueries)
            std::cout<<" GL_TIME_ELAPSED_median_ms="<<Quantile(QueryMs,0.50)
                <<" timer_query_bits="<<TimerQueryBits<<" timer_query_waits="<<TimerQueryWaits;
        std::cout<<" field_min_m="<<*HeightMinIt<<" field_max_m="<<*HeightMaxIt
            <<" imprint_min_m="<<ImprintMin<<" frames="<<FramesWritten
            <<" field_checksum_fnv1a64="<<std::hex<<std::setw(16)<<std::setfill('0')
            <<FinalFieldChecksum<<std::dec<<std::setfill(' ')
            <<" equation_audit="<<(FirstStepAudit?"PASS":"FAIL")
            <<" timer_query_audit="<<(Options.TimerQueries
                ?(QueryAudit?"PASS":"FAIL"):"NOT_RUN")
            <<" status="<<(Good?"PASS":"FAIL")<<" out="<<Options.Output.string()<<"\n";

        Render.Destroy();
        if(QueryCount>0)glDeleteQueries(QueryCount,Queries.data());
        glDeleteTextures(2,Fields.data());
        glDeleteProgram(ComputeProgram);
        return Good;
    }
}

int main(int Argc,char** Argv)
{
    FOptions Options;
    if(!ParseOptions(Argc,Argv,Options))return 2;
    if(!Run(Options))return 1;
    return 0;
}