fn next_velocity()

in rayon/rayon-demo/src/nbody/nbody.rs [195:331]

• 99 lines of code
• 15 McCabe index (conditional complexity)

fn next_velocity(time: usize, prev: &Body, bodies: &[Body]) -> (Vector3<f64>, Vector3<f64>) {
    let time = time as f64;
    let center = Point3 {
        x: (time / 22.0).cos() * -4200.0,
        y: (time / 14.0).sin() * 9200.0,
        z: (time / 27.0).sin() * 6000.0,
    };

    // pull to center
    let max_distance = 3400.0;
    let pull_strength = 0.042;

    // zones
    let zone = 400.0;
    let repel = 100.0;
    let align = 300.0;
    let attract = 100.0;

    let (speed_limit, attract_power);
    if time < 500.0 {
        speed_limit = 2000.0;
        attract_power = 100.9;
    } else {
        speed_limit = 0.2;
        attract_power = 20.9;
    }

    let zone_sqrd = 3.0 * (zone * zone);

    let mut acc = Vector3::zero();
    let mut acc2 = Vector3::zero();

    let dir_to_center = center - prev.position;
    let dist_to_center = dir_to_center.magnitude();

    // orient to center
    if dist_to_center > max_distance {
        let velc = if time < 200.0 {
            0.2
        } else {
            (dist_to_center - max_distance) * pull_strength
        };

        let diff = (dir_to_center / dist_to_center) * velc;
        acc += diff;
    }

    // TODO -- parallelize this loop too? would be easy enough, it's just a reduction
    // Just have to factor things so that `tick_seq` doesn't do it in parallel.
    let zero: Vector3<f64> = Vector3::zero();
    let (diff, diff2) = bodies
        .iter()
        .fold((zero, zero), |(mut diff, mut diff2), body| {
            let r = body.position - prev.position;

            // make sure we are not testing the particle against its own position
            let are_same = r == Vector3::zero();

            let dist_sqrd = r.magnitude2();

            if dist_sqrd < zone_sqrd && !are_same {
                let length = dist_sqrd.sqrt();
                let percent = dist_sqrd / zone_sqrd;

                if dist_sqrd < repel {
                    let f = (repel / percent - 1.0) * 0.025;
                    let normal = (r / length) * f;
                    diff += normal;
                    diff2 += normal;
                } else if dist_sqrd < align {
                    let thresh_delta = align - repel;
                    let adjusted_percent = (percent - repel) / thresh_delta;
                    let q = (0.5 - (adjusted_percent * PI * 2.0).cos() * 0.5 + 0.5) * 100.9;

                    // normalize vel2 and multiply by factor
                    let vel2_length = body.velocity2.magnitude();
                    let vel2 = (body.velocity2 / vel2_length) * q;

                    // normalize own velocity
                    let vel_length = prev.velocity.magnitude();
                    let vel = (prev.velocity / vel_length) * q;

                    diff += vel2;
                    diff2 += vel;
                }

                if dist_sqrd > attract {
                    // attract
                    let thresh_delta2 = 1.0 - attract;
                    let adjusted_percent2 = (percent - attract) / thresh_delta2;
                    let c =
                        (1.0 - ((adjusted_percent2 * PI * 2.0).cos() * 0.5 + 0.5)) * attract_power;

                    // normalize the distance vector
                    let d = (r / length) * c;

                    diff += d;
                    diff2 -= d;
                }
            }

            (diff, diff2)
        });

    acc += diff;
    acc2 += diff2;

    // Speed limits
    if time > 500.0 {
        let acc_squared = acc.magnitude2();
        if acc_squared > speed_limit {
            acc *= 0.015;
        }

        let acc_squared2 = acc2.magnitude2();
        if acc_squared2 > speed_limit {
            acc2 *= 0.015;
        }
    }

    let mut new = prev.velocity + acc;
    let mut new2 = prev.velocity2 + acc2;

    if time < 500.0 {
        let acs = new2.magnitude2();
        if acs > speed_limit {
            new2 *= 0.15;
        }

        let acs2 = new.magnitude2();
        if acs2 > speed_limit {
            new *= 0.15;
        }
    }

    (new, new2)
}