#ifndef TSUNAMI_CORE_H
#define TSUNAMI_CORE_H

#include <stdint.h>
#include <stddef.h>

#define TSUNAMI_MAGIC "TSUNAMI_SIM_V1"
#define TSUNAMI_MAGIC_SIZE 16
#define SCRATCH_LAND_VALUE 100.0f
#define DEFAULT_DT_SECONDS 15 /* v125: use a true 15-second model tick instead of 30s + 2 internal substeps */
#define WAVE_INTERNAL_SUBSTEPS 1 /* v125: one numerical wave update per 15-second model tick */
#define WAVE_FACE_DEPTH_FLUX_ENABLED 1 /* v137: use face-depth flux form for bathymetry transitions */
#define WAVE_FACE_DEPTH_MODE_GEOMETRIC 1 /* 1=sqrt(d0*d1), 0=arithmetic mean */
#define WAVE_OUTPUT_SHOALING_ENABLED 0 /* solver/display field: no extra depth^-1/4 multiplier; avoid double-counting shoaling already produced by the variable-depth solver. */
#define STATION_SHORELINE_EQUIV_SHOALING_ENABLED 1 /* stations/JESIM remain 1 m shoreline-equivalent gauges, because finite grid cells may sample 100 m+ water. */
#define OCEAN_MAX_HEIGHT_ARRIVAL_SHORELINE_EQUIV_ENABLED 1 /* Max-height colors use actual grid-cell eta, but arrival-time contours keep the 1 m shoreline-equivalent threshold basis. */
#define DEFAULT_STEPS 5760 /* keep the default model horizon close to the old 24h span */
#define DEFAULT_COAST_GUARD_CELLS 0
#define OUTER_ABSORB_CELLS 80
#define ABSORB_MARGIN_CELLS OUTER_ABSORB_CELLS
#define ABSORB_MAX_STRENGTH 0.46f
#define OUTER_ABSORB_DEPTH_M 6000.0f
#define EAST_VISIBLE_GATE_CELLS 1
#define EAST_GATE_MIN_DEPTH_M 4000.0f
#define EAST_GATE_DEPTH_SCAN_CELLS 120
#define STATION_OBSERVATION_DEPTH_M 1.0f
#define STATION_OFFSHORE_SAMPLE_RADIUS_CELLS 2
#define GLOBAL_HEIGHT_DECAY_KEEP_6H 0.775f /* global decay: 22.5% amplitude loss per 6h */
#define GLOBAL_HEIGHT_DECAY_PERIOD_S 21600.0f
#define COAST_REFLECTION_DEFAULT_ENABLED 1
#define COAST_REFLECTION_DEFAULT_ENERGY_LOSS 0.55f
#define COAST_REFLECTION_DISTRIBUTE_STEPS 24.0f
#define FORCE_WATER_MIN_DEPTH_M 50.0f
#define FORCE_WATER_DEPTH_SCAN_CELLS 96
#define NARROW_SEA_DEFAULT_ENABLED 0
#define NARROW_SEA_DEFAULT_ENERGY_LOSS 0.0f
#define SOURCE_RAMP_DEFAULT_SECONDS 120.0f
#define OCEAN_CELL_ARRIVAL_THRESHOLD_M 0.05f  // arrival contour threshold on 1m-depth shoreline-equivalent height

typedef struct {
    uint32_t width;
    uint32_t height;
    uint32_t station_count;
    float *depth_m;
    float *effective_depth_m; // depth actually used by the solver/display after force-water depth repair
    uint8_t *base_land;     // DEP-derived base mask before shoreline flood-fill correction
    uint8_t *land;          // corrected land mask used for drawing; station cells are water
    uint8_t *solid;         // land plus optional coastline guard; used by the solver
    int8_t  *user_override; // -1 = none, 0 = force water, 1 = force land
    int32_t *station_index;
    int32_t *station_sample_index;  // nearest ocean cell used for station height sampling
    int32_t *station_id;
    uint8_t *station_enabled; // 0 = hidden/ignored observation point; keeps original list indices aligned
    uint8_t *edge_absorb_gate; // far-east edge cells used only as an open boundary into the hidden absorber
    float   *edge_absorb_depth_m; // smoothed/extended depth for the east open-boundary edge
} GridData;

typedef struct {
    uint32_t width;
    uint32_t height;
    uint32_t padded_height;      // visible display buffer height: height + 2
    uint32_t outer_pad;          // hidden absorbing ocean cells around the visible map
    uint32_t solver_width;       // width + 2 * outer_pad
    uint32_t solver_height;      // height + 2 * outer_pad
    size_t cells;                // visible cells only
    size_t padded_cells;         // visible display cells with top/bottom guard rows
    size_t solver_cells;         // hidden solver domain cells
    float *height_norm;          // solver_width * solver_height
    float *height_prev;          // solver_width * solver_height
    float *velocity;             // solver_width * solver_height
    float *display_height_m;     // width * (height + 2), visible map only
    float *display_max_m;        // width * (height + 2), max positive wave rise seen at each visible ocean cell
    int32_t *display_arrival_s;   // width * (height + 2), first arrival time at each visible ocean cell
    float *source_height_m;
    int32_t *source_delay_s;
    float *source_applied_m;     // cumulative source height already injected by smooth rupture ramp
    float *station_max_m;       // maximum positive water-level rise; negative drawdown is still used for arrival detection/output
    int32_t *station_arrival_s;
} Simulation;

typedef struct {
    int grid_x;
    int grid_y;
    float height_m;
    float slip_m;   // per-cell fault slip used for M0/Mw; may differ between painted cells
    float stamp_dip_deg;       // dip setting at the moment this source node was painted/loaded
    float moment_area_factor;  // per-node M0 area correction frozen at paint/load time
    int fault_id;              // 0=legacy/free brush, >0 finite-fault stamp group id
    int fault_patch_id;        // 0=not a stamp patch, >0 stamp sequence id
    float stamp_strike_deg;    // strike setting frozen for finite-fault stamp nodes
    float stamp_depth_km;      // top/reference depth frozen for finite-fault stamp nodes
    float stamp_width_km;      // true down-dip width frozen for stamp nodes
    float stamp_length_km;     // along-strike patch length frozen for stamp nodes
    int delay_s;
} SourceSpec;

typedef struct {
    SourceSpec *items;
    size_t count;
    size_t capacity;
} SourceList;

void tsunami_fatal(const char *message);
void tsunami_fatal_errno(const char *message);
void *tsunami_xcalloc(size_t count, size_t size);

GridData load_grid_data(const char *path);
void free_grid_data(GridData *g);
void rebuild_land_mask_from_coast_stations(GridData *g);
void rebuild_solid_mask(GridData *g, int coast_guard_cells);
void clear_land_overrides(GridData *g);
void set_land_override(GridData *g, int grid_x, int grid_y, int mode);
size_t count_land_overrides(const GridData *g);
int save_land_override_csv(const char *path, const GridData *g);
int load_land_override_csv(GridData *g, const char *path);
Simulation create_simulation(const GridData *g);
void free_simulation(Simulation *s);
void reset_simulation(const GridData *g, Simulation *s);
void simulate_one_step(const GridData *g, Simulation *s, int time_s);
void set_source_ramp_seconds(float seconds);
float source_ramp_seconds(void);
void update_station_stats(const GridData *g, Simulation *s, int time_s);
void write_station_file(const char *path, const GridData *g, const Simulation *s);
void write_forecast_zone_file(const char *path, const GridData *g, const Simulation *s);
const char *forecast_zone_name_from_id(int region_id);
const char *tsunami_height_color_bin(float height_m);
const char *tsunami_height_color_hex(float height_m);
float shoaling_factor(float depth_m);
float output_shoaling_factor(float depth_m);
float visible_effective_depth_m(const GridData *g, int visible_idx);
float station_shoreline_equiv_factor(float source_depth_m);
float station_observed_height_m(const GridData *g, const Simulation *s, uint32_t station_i);
float global_height_decay_per_step(void);
void set_coastal_reflection_absorption(int enabled, float energy_loss);
int coastal_reflection_absorption_enabled(void);
float coastal_reflection_absorption_energy_loss(void);
float coastal_reflection_keep_per_step(void);
void set_narrow_sea_damping(int enabled, float energy_loss);
int narrow_sea_damping_enabled(void);
float narrow_sea_damping_energy_loss(void);
float narrow_sea_keep_per_step(void);
int station_is_enabled(const GridData *g, uint32_t station_i);
uint32_t count_enabled_stations(const GridData *g);

void source_list_append(SourceList *list, int grid_x, int grid_y, float height_m, int delay_s);
void source_list_clear(SourceList *list);
void source_list_free(SourceList *list);
int parse_source_csv_line(const char *text, int *x, int *y, float *height_m, int *delay_s);
int load_sources_csv(SourceList *list, const char *path);
void apply_sources_to_simulation(const GridData *g, Simulation *s, const SourceList *sources);
void add_source_grid(const GridData *g, Simulation *s, int grid_x, int grid_y, float height_m, int delay_s);

static inline size_t tsunami_index(const GridData *g, int x, int y) {
    return (size_t)y * (size_t)g->width + (size_t)x;
}

#endif
