#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 <cerrno>
#include <cctype>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <filesystem>
#include <fstream>
#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.46;
    constexpr double SpherePushMetersPerSecondSquared=14.0;
    constexpr double Pi=3.14159265358979323846;

    struct FOptions
    {
        std::string Device="auto";
        std::filesystem::path Output="Saved/Linux/l162/living-plane";
        int Size=128;
        int Steps=360;
        int CaptureEvery=1;
        bool Breathing=false;
    };

    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,rg32f) readonly uniform image2D uSource;
layout(binding=1,rg32f) 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 uForcePosition;
uniform float uForceSigma;
uniform float uForceStrength;
uniform int uBreathing;
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);
    float h=imageLoad(uSource,p).r;
    float v=imageLoad(uSource,p).g;
    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);
    vec2 offset=world-uForcePosition;
    float sigma2=max(uForceSigma*uForceSigma,1.0e-8);
    float interaction=-uForceStrength*exp(-dot(offset,offset)/(2.0*sigma2));
    float breathing=0.0;
    if(uBreathing!=0)
    {
        float spatial=sin(world.x*0.72)*sin(world.y*0.72);
        breathing=2.4*spatial*sin(uTime*1.45);
    }
    float acceleration=uWaveSpeed*uWaveSpeed*laplacian-uFoundationStiffness*h
        +interaction+breathing;
    float nextVelocity=(v+acceleration*uDt)*exp(-uDamping*uDt);
    float nextHeight=h+nextVelocity*uDt;
    imageStore(uTarget,p,vec4(nextHeight,nextVelocity,0.0,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 mat4 uViewProjection;
out vec3 vWorld;
out vec3 vNormal;
out float vHeight;
flat out int vSphere;
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));
        height=texelFetch(uField,texel,0).r;
        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=(texelFetch(uField,right,0).r-texelFetch(uField,left,0).r)/dx;
        float dhdy=(texelFetch(uField,up,0).r-texelFetch(uField,down,0).r)/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;
    vSphere=uRenderSphere;
    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 vSphere;
uniform vec3 uCameraPosition;
out vec4 outColor;
void main()
{
    vec3 n=normalize(vNormal);
    vec3 l=normalize(vec3(-0.45,-0.55,0.72));
    vec3 v=normalize(uCameraPosition-vWorld);
    float diffuse=max(dot(n,l),0.0);
    float spec=pow(max(dot(reflect(-l,n),v),0.0),vSphere==0?58.0:38.0);
    if(vSphere!=0)
    {
        vec3 base=vec3(0.98,0.43,0.13);
        float light=0.24+0.84*diffuse;
        outColor=vec4(base*light+vec3(0.65,0.43,0.22)*spec,1.0);
    }
    else
    {
        float crest=clamp(0.50+vHeight*0.8,0.0,1.0);
        vec3 deep=vec3(0.035,0.22,0.34);
        vec3 crestColor=vec3(0.13,0.62,0.70);
        vec3 base=mix(deep,crestColor,crest);
        float light=0.28+0.88*diffuse;
        outColor=vec4(base*light+vec3(0.31,0.62,0.68)*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; };

    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)
    {
        const FVec3 Target={0.0,0.0,0.0};
        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=100.0;
        const double F=1.0/std::tan(42.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{8.8,-12.2,9.4};

        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,256);
            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);

            ViewProjection=MakeViewProjection(Width,Height,Camera);
            glBindVertexArray(0);
            return glGetError()==GL_NO_ERROR;
        }

        bool Render(GLuint FieldTexture,const std::filesystem::path& Output,
                    double SphereX,double SphereY,double SphereZ)
        {
            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);
            glUniform3f(glGetUniformLocation(Program,"uSphereCenter"),
                static_cast<float>(SphereX),static_cast<float>(SphereY),static_cast<float>(SphereZ));
            glUniform1f(glGetUniformLocation(Program,"uSphereRadius"),
                static_cast<float>(SphereRadiusMeters));
            glBindVertexArray(SphereVao);
            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 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] [--size=N] [--steps=N] "
                    "[--out=DIR] [--capture-every=N] [--breathing=on|off]\n"
                    "  --size is the square height-field and PNG resolution.\n"
                    "  --capture-every=0 writes only timing.json. Default is 1.\n";
                std::exit(0);
            }
            else if(Take("--device",Value))Options.Device=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("--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;}
            }
            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.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;}
        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)];
    }

    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,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());
        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"
            <<"  \"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"
            <<"  \"frames_written\": "<<Frames<<",\n"
            <<"  \"first_step_equation_audit\": \""<<(FirstStepAudit?"pass":"fail")<<"\",\n"
            <<"  \"field_dynamics_audit\": \""<<(DynamicsAudit?"pass":"fail")<<"\",\n"
            <<"  \"timing_interpretation\": \""
            <<(Context.ExecutionDevice=="gpu"?"GPU GL_TIME_ELAPSED per compute dispatch"
                :"software-renderer GL_TIME_ELAPSED values; not GPU performance timings")<<"\",\n"
            <<"  \"field_min_m\": "<<HeightMin<<",\n"
            <<"  \"field_max_m\": "<<HeightMax<<",\n"
            <<"  \"field_mean_m\": "<<HeightMean<<",\n"
            <<"  \"field_rms_m\": "<<HeightRms<<",\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*2,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_RG32F,Size,Size,0,GL_RG,GL_FLOAT,Zeroes.data());
        }

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

        const GLuint QueryCount=static_cast<GLuint>(Options.Steps);
        std::vector<GLuint> Queries(QueryCount);
        glGenQueries(static_cast<GLsizei>(QueryCount),Queries.data());
        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 LocationForcePosition=glGetUniformLocation(ComputeProgram,"uForcePosition");
        const GLint LocationForceSigma=glGetUniformLocation(ComputeProgram,"uForceSigma");
        const GLint LocationForceStrength=glGetUniformLocation(ComputeProgram,"uForceStrength");
        const GLint LocationBreathing=glGetUniformLocation(ComputeProgram,"uBreathing");
        if(LocationGrid<0||LocationDt<0||LocationDx<0||LocationPlane<0||LocationWaveSpeed<0
            ||LocationDamping<0||LocationFoundation<0||LocationTime<0
            ||LocationForcePosition<0||LocationForceSigma<0
            ||LocationForceStrength<0||LocationBreathing<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);

        int Current=0;
        const auto Dispatch=[&](int Target,double Time,double SphereX,double SphereY,
                                double ForceStrength,GLuint Query)
        {
            glUseProgram(ComputeProgram);
            glUniform1f(LocationTime,static_cast<float>(Time));
            glUniform2f(LocationForcePosition,static_cast<float>(SphereX),static_cast<float>(SphereY));
            glUniform1f(LocationForceStrength,static_cast<float>(ForceStrength));
            glBindImageTexture(0,Fields[static_cast<std::size_t>(Current)],0,GL_FALSE,0,
                GL_READ_ONLY,GL_RG32F);
            glBindImageTexture(1,Fields[static_cast<std::size_t>(Target)],0,GL_FALSE,0,
                GL_WRITE_ONLY,GL_RG32F);
            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 with zero forcing before timing.
        for(int Warm=0;Warm<8;++Warm)
            Dispatch(1-Current,-8.0*Dt,0.0,0.0,0.0,0);
        glFinish();
        for(const GLuint Field:Fields)
        {
            glBindTexture(GL_TEXTURE_2D,Field);
            glTexSubImage2D(GL_TEXTURE_2D,0,0,0,Size,Size,GL_RG,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 SphereX,double SphereY,double SphereZ)->bool
        {
            if(!RenderPtr)return true;
            char Filename[32]{};
            std::snprintf(Filename,sizeof(Filename),"frame_%06d.png",FrameIndex++);
            if(!Render.Render(Fields[static_cast<std::size_t>(Current)],
                Options.Output/Filename,SphereX,SphereY,SphereZ))return false;
            ++FramesWritten;
            return true;
        };

        const int MoveSteps=std::max(2,static_cast<int>(std::lround(Options.Steps*0.58)));
        const auto SphereAt=[&](int Step)
        {
            const double MoveT=std::clamp(static_cast<double>(Step)/std::max(1,MoveSteps-1),0.0,1.0);
            const double Smooth=MoveT*MoveT*(3.0-2.0*MoveT);
            const double X=-4.8+9.6*Smooth;
            const bool Rest=Step>=MoveSteps;
            const double RestTime=Rest?static_cast<double>(Step-MoveSteps)*Dt:0.0;
            const double Bob=Rest?std::sin(RestTime*2.0*Pi*1.25):0.0;
            const double Z=0.36+0.08*Bob;
            const double Force=SpherePushMetersPerSecondSquared*(1.0-0.32*Bob);
            return std::array<double,4>{X,0.0,Z,Force};
        };

        const auto Initial=SphereAt(0);
        if(!SaveFrame(Initial[0],Initial[1],Initial[2]))
        {Render.Destroy();glDeleteQueries(static_cast<GLsizei>(Queries.size()),Queries.data());
            glDeleteTextures(2,Fields.data());glDeleteProgram(ComputeProgram);return false;}

        std::vector<float> FirstField;
        bool FirstStepAudit=false;
        for(int Step=0;Step<Options.Steps;++Step)
        {
            const auto Sphere=SphereAt(Step);
            const int Target=1-Current;
            Dispatch(Target,static_cast<double>(Step)*Dt,Sphere[0],Sphere[1],Sphere[3],
                Queries[static_cast<std::size_t>(Step)]);
            if(Step==0)
            {
                glFinish();
                FirstField.resize(static_cast<std::size_t>(Size)*Size*2);
                glBindTexture(GL_TEXTURE_2D,Fields[static_cast<std::size_t>(Current)]);
                glGetTexImage(GL_TEXTURE_2D,0,GL_RG,GL_FLOAT,FirstField.data());
                const int X=std::clamp(static_cast<int>(std::floor(
                    (Sphere[0]+PlaneLengthMeters*0.5)/Dx)),0,Size-1);
                const int Y=std::clamp(static_cast<int>(std::floor(
                    (Sphere[1]+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;
                const double Distance2=(CellX-Sphere[0])*(CellX-Sphere[0])
                    +(CellY-Sphere[1])*(CellY-Sphere[1]);
                const double Force=-Sphere[3]*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)*2;
                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))
            {
                const auto FrameSphere=SphereAt(Step);
                if(!SaveFrame(FrameSphere[0],FrameSphere[1],FrameSphere[2]))
                {Render.Destroy();glDeleteQueries(static_cast<GLsizei>(Queries.size()),Queries.data());
                    glDeleteTextures(2,Fields.data());glDeleteProgram(ComputeProgram);return false;}
            }
        }
        glFinish();

        std::vector<GLuint64> Times(static_cast<std::size_t>(Options.Steps));
        for(int Step=0;Step<Options.Steps;++Step)
            glGetQueryObjectui64v(Queries[static_cast<std::size_t>(Step)],GL_QUERY_RESULT,
                &Times[static_cast<std::size_t>(Step)]);
        std::vector<float> FinalField(static_cast<std::size_t>(Size)*Size*2);
        glBindTexture(GL_TEXTURE_2D,Fields[static_cast<std::size_t>(Current)]);
        glGetTexImage(GL_TEXTURE_2D,0,GL_RG,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;
        for(std::size_t I=0;I<FinalField.size();I+=2)
        {
            Finite=Finite&&std::isfinite(FinalField[I])&&std::isfinite(FinalField[I+1]);
            Heights.push_back(FinalField[I]);
            Velocities.push_back(FinalField[I+1]);
        }
        const auto [HeightMinIt,HeightMaxIt]=std::minmax_element(Heights.begin(),Heights.end());
        const bool DynamicsAudit=Finite&&*HeightMinIt<-1.0e-7&&*HeightMaxIt-*HeightMinIt<2.0;
        const bool Good=FirstStepAudit&&DynamicsAudit&&glGetError()==GL_NO_ERROR;
        if(!DynamicsAudit)
            std::cerr<<"[LivingPlane] field audit failed: finite="<<Finite<<" min="
                <<*HeightMinIt<<" max="<<*HeightMaxIt<<"\n";
        if(!WriteTimingJson(Options,Context,Renderer,Vendor,Version,Dt,Times,Heights,Velocities,
            FramesWritten,FirstStepAudit,DynamicsAudit))
        {Render.Destroy();glDeleteQueries(static_cast<GLsizei>(Queries.size()),Queries.data());
            glDeleteTextures(2,Fields.data());glDeleteProgram(ComputeProgram);return false;}

        double MeanNs=0.0;
        for(const GLuint64 Ns:Times)MeanNs+=static_cast<double>(Ns);
        MeanNs/=Times.size();
        std::cout<<"[LivingPlane] size="<<Size<<" steps="<<Options.Steps<<" dt_ms="
            <<Dt*1000.0<<" GL_TIME_ELAPSED_mean_ms="<<MeanNs*1.0e-6
            <<" median_ms="<<Quantile([&]{std::vector<double> V;for(auto N:Times)V.push_back(N*1.0e-6);return V;}(),0.5)
            <<" field_min_m="<<*HeightMinIt<<" field_max_m="<<*HeightMaxIt
            <<" frames="<<FramesWritten<<" equation_audit="<<(FirstStepAudit?"PASS":"FAIL")
            <<" status="<<(Good?"PASS":"FAIL")<<" out="<<Options.Output.string()<<"\n";

        Render.Destroy();
        glDeleteQueries(static_cast<GLsizei>(Queries.size()),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;
}