// ============================================================================
// EROSION DIRECTION WRANGLE – Determines sand sliding direction based on slope,
// and user-defined thresholds. This is the first step in the 
// erosion-sedimentation cycle.
// ============================================================================

float grid_resolution = chf("grid_resolution"); // Size of each grid cell
float slope_threshold = chf("slope_threshold"); // Minimum slope required for sand to slide
float scale_factor = chf("scale_factor"); // Scale factor to avoid extreme results

// --- Compute the effective slope threshold for sliding ---
float sliding_threshold = grid_resolution * slope_threshold;

// --- Cache current cell height and grid position ---
int x = @ix;
int y = @iy;

/*
   Define 4 cardinal directions (left, right, up, down).
   These are used to evaluate slope and determine possible sliding directions.
*/

vector dirs[] = array(
    set(-1, 0, 0),  // links
    set( 1, 0, 0),  // rechts
    set( 0, -1, 0), // unten
    set( 0,  1, 0)  // oben
);

vector slide_dirs[] = array();

/*
   Loop through all neighboring cells to evaluate slope.
   If the slope toward a neighbor exceeds the adjusted threshold,
   that direction is considered a valid candidate for sand sliding.
*/

foreach (vector dir; dirs) {
    int nx = x + int(dir.x);
    int ny = y + int(dir.y);
    vector neighbor_pos = set(nx, ny, 0); // Grid position

    float h_neighbor = volumeindex(0, "height", neighbor_pos); // Neighbor height
    int unstable = (@height - h_neighbor) > sliding_threshold;

    if (unstable) {
        append(slide_dirs, dir); // Save direction
    }
}

// --- If there are valid sliding directions, choose one.

int direction_code = -1; // Stable condition

if (len(slide_dirs) > 0) {
    float seed_x = @Time + v@P.x * @height;
    float seed_y = @Time + v@P.y * @height;
    float rand = noise(set(seed_x, seed_y, 0));

    int index = int(fit01(rand, 0, len(slide_dirs) - 0.001));
    vector chosen_dir = slide_dirs[index];

    if (chosen_dir == set(-1, 0, 0)) direction_code = 0;
    else if (chosen_dir == set(1, 0, 0)) direction_code = 1;
    else if (chosen_dir == set(0, -1, 0)) direction_code = 2;
    else if (chosen_dir == set(0, 1, 0)) direction_code = 3;
}

@delta = float(direction_code); // Direction as integer layer
@mask = direction_code != -1 ? 1.0 : 0.0; // Mark instable cells

//--- Save important values as detail attributes for further use.

setdetailattrib(geoself(), "sliding_threshold", sliding_threshold, "set");
setdetailattrib(geoself(), "scale_factor", scale_factor, "set");

// ============================================================================
// EROSION APPLICATION – Applies sand movement based on previously computed d
// irections. This is the second step in the erosion-sedimentation cycle.
// ============================================================================

float sliding_threshold = detail(0, "sliding_threshold") * detail(0, "scale_factor");

/*
   Define 4 cardinal directions (left, right, up, down).
   These are used to evaluate slope and determine possible sliding directions.
*/

vector dirs[] = array(
    set(-1, 0, 0), // Left
    set( 1, 0, 0), // Right
    set( 0, -1, 0), // Bottom
    set( 0, 1, 0)  // Top
);

float height_increment = 0;
int x = @ix;
int y = @iy;

int opposite[] = array(1, 0, 3, 2); // Flip directions: left -> right, bottom -> top for correct results

// --- Read the delta layer values and apply them to the currently process cell ---

for (int i = 0; i < len(dirs); i++) {
    int nx = x + int(dirs[i].x);
    int ny = y + int(dirs[i].y);
    float neighbor_delta = volumeindex(0, "delta", set(nx, ny, 0));

    if (int(neighbor_delta) == opposite[i]) {
        height_increment += sliding_threshold;
    }
}

if (@delta != -1) {
    height_increment -= sliding_threshold;
}

// --- Apply accumulated height change ---

@height += height_increment;