#version 430

// Runge-Kutta 4th order - Stage 1 compute shader
// Computes k1 = h * f(t, y)

layout(local_size_x = 256, local_size_y = 1, local_size_z = 1) in;

layout(std430, binding = 0) buffer YBuffer {
    float y[];
};

layout(std430, binding = 1) buffer K1Buffer {
    float k1[];
};

uniform float t;
uniform float h;
uniform int n;

// Example ODE function: dy/dt = f(t, y)
// This would be replaced with actual problem-specific function
vec4 ode_function(float time, vec4 state) {
    // Example: harmonic oscillator
    // y[0] = position, y[1] = velocity
    // dy[0]/dt = y[1]
    // dy[1]/dt = -y[0]
    return vec4(state.y, -state.x, 0.0, 0.0);
}

void main() {
    uint index = gl_GlobalInvocationID.x;
    
    if (index >= n) return;
    
    // Load current state
    vec4 current_state = vec4(
        y[index * 4 + 0],
        y[index * 4 + 1], 
        y[index * 4 + 2],
        y[index * 4 + 3]
    );
    
    // Compute derivative at current state
    vec4 derivative = ode_function(t, current_state);
    
    // Store k1 = h * f(t, y)
    k1[index * 4 + 0] = h * derivative.x;
    k1[index * 4 + 1] = h * derivative.y;
    k1[index * 4 + 2] = h * derivative.z;
    k1[index * 4 + 3] = h * derivative.w;
}