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

// --- Retrieve simulation parameters from the UI ---
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"); // Multiplier controlling erosion intensity
float sedimentation_scale = chf("sedimentation_scale"); // Controls how much sediment accumulates in low areas
float wind_angle_deg = chf("wind_direction"); // Wind direction in degrees (0° = +X, 90° = +Y)
float wind_bias_strength = chf("wind_bias_strength"); // Strength of wind influence on sand movement

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

/* 
   Convert wind angle from degrees to radians and generate a unit vector
   pointing in the wind direction. This vector will be used to bias sand
   movement toward the wind.
*/
float wind_angle_rad = radians(wind_angle_deg);
vector base_wind_dir = set(cos(wind_angle_rad), sin(wind_angle_rad), 0);

/*
   Define the 8 possible neighbor directions:
   - 4 cardinal directions (left, right, up, down)
   - 4 diagonal directions (top-left, top-right, bottom-left, bottom-right)
   These are used to evaluate slope and determine possible sliding directions.
*/
vector dirs[] = array(
    set(-1,  0, 0), set( 1,  0, 0),
    set( 0, -1, 0), set( 0,  1, 0),
    set(-1, -1, 0), set( 1, -1, 0),
    set(-1,  1, 0), set( 1,  1, 0)
);

// --- Initialize array to store valid sliding directions ---
int slide_indices[];

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

/*
   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.
*/
for (int i = 0; i < len(dirs); i++) {
    vector dir = dirs[i];
    int nx = x + int(dir.x);
    int ny = y + int(dir.y);
    float h_neighbor = volumeindex(0, "height", set(nx, ny, 0));

    float slope = h_self - h_neighbor;
    float distance = length(dir); // 1.0 for cardinal, ~1.414 for diagonal
    float adjusted_threshold = sliding_threshold * distance;

    if (slope > adjusted_threshold) {
        append(slide_indices, i); // Store index of valid direction
    }
}

/*
   If there are valid sliding directions, choose one.
   The selection is biased toward the wind direction using a dot product.
   Directions aligned with the wind receive higher weights.
*/
int direction_code = -1;
if (len(slide_indices) > 0) {
    float weights[] = {};

    for (int i = 0; i < len(slide_indices); i++) {
        vector dir = dirs[slide_indices[i]];
        float bias = dot(normalize(dir), base_wind_dir); // Alignment with wind
        float weight = 1 + bias * wind_bias_strength;    // Wind bias scaling
        append(weights, weight);
    }

    /*
       Randomly select one direction from the weighted list.
       This introduces a stochastic approach while still favoring wind-aligned movement.
    */
    int pick = sample_discrete(weights, rand(@ix * 17 + @iy * 43 + @Time));
    direction_code = slide_indices[pick];
}

// --- Store the chosen direction as a float attribute for use in the next wrangle ---
@delta = float(direction_code);

/*
   On the first frame, store global parameters as detail attributes.
   These will be accessed in Wrangle 2 to ensure consistent values across voxels.
*/
if (@Frame == 1) {
    setdetailattrib(0, "sliding_threshold", sliding_threshold, "set");
    setdetailattrib(0, "scale_factor", scale_factor, "set");
    setdetailattrib(0, "sedimentation_scale", sedimentation_scale, "set");
}

// ============================================================================
// EROSION APPLICATION & SEDIMENTATION WRANGLE – Applies sand movement based on
// previously computed directions, and adds sediment in low-lying regions.
// This is the second step in the erosion-sedimentation cycle.
// ============================================================================

// --- Retrieve global parameters stored in Wrangle W1 ---
float sliding_threshold     = detail(0, "sliding_threshold") * detail(0, "scale_factor");
float sedimentation_scale   = detail(0, "sedimentation_scale");

/*
   Define the 8 neighboring directions:
   - 4 cardinal directions (left, right, up, down)
   - 4 diagonal directions (top-left, top-right, bottom-left, bottom-right)
   These are used to evaluate incoming and outgoing sand flow.
*/
vector dirs[] = array(
    set(-1,  0, 0), set( 1,  0, 0),
    set( 0, -1, 0), set( 0,  1, 0),
    set(-1, -1, 0), set( 1, -1, 0),
    set(-1,  1, 0), set( 1,  1, 0)
);

// --- Cache current voxel position ---
int x = @ix;
int y = @iy;

// --- Initialize height change accumulator ---
float height_increment = 0;

// --- Retrieve own sliding direction from previous wrangle ---
int my_delta = int(@delta);
vector my_dir = my_delta >= 0 ? dirs[my_delta] : {0, 0, 0};

/*
   Loop through all neighbors to check if any of them are sliding sand
   toward this cell. If a neighbor's sliding direction points toward us,
   we receive sand from them.
*/
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 (neighbor_delta >= 0) {
        vector neighbor_dir = dirs[int(neighbor_delta)];

        /*
           If the neighbor's direction is opposite to the current direction,
           it means sand is flowing into this cell.
        */
        if (length(neighbor_dir + dirs[i]) < 0.01) {
            height_increment += sliding_threshold;
        }
    }
}

/*
   If this cell is actively sliding sand away (i.e. has a valid direction),
   subtract sand from its height.
*/
if (my_delta != -1) {
    height_increment -= sliding_threshold;
}

/*
   If this cell is not sliding (i.e. stable), check if it lies in a depression.
   If the average height of its neighbors is higher, we simulate sedimentation
   by adding sand proportionally.
*/
if (@delta == -1) {
    float h_self = @height;
    float h_sum = 0;
    int count = 0;

    for (int i = 0; i < len(dirs); i++) {
        int nx = x + int(dirs[i].x);
        int ny = y + int(dirs[i].y);
        float h_neighbor = volumeindex(0, "height", set(nx, ny, 0));
        h_sum += h_neighbor;
        count++;
    }

    float h_avg = h_sum / count;
    float sediment_amount = clamp(h_avg - h_self, 0, sliding_threshold);

    // Add sediment scaled by user-defined factor
    height_increment += sediment_amount * sedimentation_scale;
}

// --- Apply accumulated height change ---
@height += height_increment;