// Define the custom parameters

float base_angle_deg = chf("wind_angle");
float variation_pos_deg = chf("angle_variation_pos");
float variation_neg_deg = chf("angle_variation_neg");
float turbulence_strength = chf("turbulence_strength");
float wind_strength = chf("wind_strength");
float wind_shear_factor = chf("wind_shear_factor");

// Calculate the wind direction with variance
// This part converts the angle's degrees into radians.
// The offset calculates a random value from time, position and an arbitrary vector. The result is then remapped to values between variation_neg_deg and variation_pos_deg.
// The get final wind direction, the values are added.

float base_angle_rad = radians(base_angle_deg + 180.0);
float angle_offset = fit01(rand(@Time + dot(@P, set(12.345, 67.89, 0))), -variation_neg_deg, variation_pos_deg);
float final_angle = base_angle_rad + radians(angle_offset);

// Main wind vector with turbulence
// The sine and cosine values of final_angle are assembled to a vector.
// A curlnoise function adds small-scale turbulence to the wind.
// The three vectors wind_dir_base, wind_turbulence, turbulence_strength are normalized to get a direction instead of a magnitude.

vector wind_dir_base = set(cos(final_angle), 0.0, sin(final_angle));
vector wind_turbulence = curlnoise(@P * 0.1 + @Time);
vector wind_dir = normalize(wind_dir_base + wind_turbulence * turbulence_strength);

// Wind shear
// The wind shear has a 90° clockwise rotation in the XZ plane.
// The wind shear strength is the product of wind_strength and wind_shear_factor.

vector cross_wind_dir = normalize(set(-wind_dir.z, 0.0, wind_dir.x));
float cross_strength  = wind_strength * wind_shear_factor;

// target_pos determines the target position of the current voxel's sand

vector target_pos = @P + wind_dir * wind_strength + cross_wind_dir * cross_strength;

// Probability factor
// To get a more random appearance, the script calculates a transport probability for each voxel.

float prob_factor = fit01(rand(@Time + dot(@P, wind_dir)), 0.5, 1.0);
float transport_prob = wind_strength * prob_factor;

// This part samples the terrain's height values at the target position and the current voxel's position.

float target_height = volumesample(0, "height", target_pos);
float current_height = volumesample(0, "height", @P);

// Set the terrain's height at the current voxel
// The script checks if a voxel's sand will be moving based on the previous samples and calculations.
// If the sand moves, the script determines if the sand is removed or accumulated.
// If sand is removed, the amount of sand is directly converted into the terrain's actual height.

if (current_height > target_height) {
    if (rand(@P + @Time) < transport_prob) {
        float sand_amount = current_height - target_height;

        if (sand_amount > 0.0) {
            current_height -= sand_amount;
            target_height  += sand_amount;

            @height = current_height;
        }
    }
}