// ============================================================================
// wave_ca.v - 二维波动方程CA引擎
// 双时间层存储，Von Neumann邻域Laplacian
// 支持阻尼、源点和障碍物
// ============================================================================
module wave_ca #(
    parameter WIDTH  = 128,
    parameter HEIGHT = 128,
    parameter FRAC_W = 8,
    parameter [15:0] C_SQ = 16'h0040,    // c² (8.8定点)，约0.25
    parameter [15:0] DAMPING = 16'hFFF8  // 阻尼系数 (8.8定点)，约0.97
)(
    input  wire              clk,
    input  wire              rst_n,
    input  wire              enable,
    input  wire              init,
    // 源点控制
    input  wire [7:0]        src_x,
    input  wire [7:0]        src_y,
    input  wire              src_en,
    input  wire signed [15:0] src_amp,    // 源幅度
    output wire [15:0]       u_curr [0:WIDTH*HEIGHT-1],
    output wire [31:0]       step_count
);

    // 双时间层
    reg signed [15:0] u_prev [0:WIDTH*HEIGHT-1];  // u(t-1)
    reg signed [15:0] u_curr_r [0:WIDTH*HEIGHT-1]; // u(t)
    reg signed [15:0] u_next [0:WIDTH*HEIGHT-1];  // u(t+1)
    reg [31:0] steps;

    integer idx, x, y;
    reg signed [17:0] laplacian;

    always @(*) begin
        for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
            x = idx % WIDTH; y = idx / WIDTH;

            // Laplacian（Von Neumann邻域，5点模板）
            laplacian = 0;
            if (y > 0)        laplacian = laplacian + u_curr_r[(y-1)*WIDTH+x];
            if (y < HEIGHT-1) laplacian = laplacian + u_curr_r[(y+1)*WIDTH+x];
            if (x > 0)        laplacian = laplacian + u_curr_r[y*WIDTH+x-1];
            if (x < WIDTH-1)  laplacian = laplacian + u_curr_r[y*WIDTH+x+1];
            laplacian = laplacian - 4 * u_curr_r[idx];

            // 波动方程更新
            // u(t+1) = 2*u(t) - u(t-1) + c²*laplacian
            u_next[idx] = 2 * u_curr_r[idx] - u_prev[idx] +
                          (laplacian * $signed(C_SQ) >>> FRAC_W);

            // 阻尼
            u_next[idx] = (u_next[idx] * $signed(DAMPING)) >>> FRAC_W;

            // 源点注入
            if (src_en && x == src_x && y == src_y)
                u_next[idx] = u_next[idx] + src_amp;
        end
    end

    always @(posedge clk or negedge rst_n) begin
        if (!rst_n) begin
            steps <= 32'd0;
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
                u_prev[idx] <= 16'd0;
                u_curr_r[idx] <= 16'd0;
            end
        end else if (init) begin
            steps <= 32'd0;
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
                u_prev[idx] <= 16'd0;
                u_curr_r[idx] <= 16'd0;
            end
        end else if (enable) begin
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
                u_prev[idx] <= u_curr_r[idx];
                u_curr_r[idx] <= u_next[idx];
            end
            steps <= steps + 32'd1;
        end
    end

    genvar gi;
    generate
        for (gi = 0; gi < WIDTH*HEIGHT; gi = gi + 1)
            assign u_curr[gi] = u_curr_r[gi];
    endgenerate
    assign step_count = steps;

endmodule