|
| 1 | +type Number = f64; |
| 2 | + |
| 3 | +const G: Number = 6.67e-11; |
| 4 | +const TIMESTEP: Number = 0.25; |
| 5 | + |
| 6 | +struct Position { |
| 7 | + x: Number, |
| 8 | + y: Number, |
| 9 | + z: Number, |
| 10 | +} |
| 11 | +struct Velocity { |
| 12 | + dx: Number, |
| 13 | + dy: Number, |
| 14 | + dz: Number, |
| 15 | +} |
| 16 | +struct Force { |
| 17 | + fx: Number, |
| 18 | + fy: Number, |
| 19 | + fz: Number, |
| 20 | +} |
| 21 | +struct Acceleration { |
| 22 | + ax: Number, |
| 23 | + ay: Number, |
| 24 | + az: Number, |
| 25 | +} |
| 26 | + |
| 27 | +/// The main structure of this program. |
| 28 | +/// The state for a body is split across |
| 29 | +/// a position vector, |
| 30 | +/// a velocities vector, |
| 31 | +/// and a mass vector. |
| 32 | +/// |
| 33 | +/// Originally, this was implemented as an array of structures—now |
| 34 | +/// it's a structure of arrays. |
| 35 | +/// This was both for testing optimizations and for minute practice with EC(S). |
| 36 | +pub struct BodyStates { |
| 37 | + poss: Vec<Position>, |
| 38 | + vels: Vec<Velocity>, |
| 39 | + masses: Vec<Number>, |
| 40 | +} |
| 41 | + |
| 42 | +fn dist_squared(dx: Number, dy: Number, dz: Number) -> Number { |
| 43 | + (dx * dx) + (dy * dy) + (dz * dz) |
| 44 | +} |
| 45 | + |
| 46 | +fn force_d(mass1: Number, mass2: Number, distance_squared: Number) -> Number { |
| 47 | + (G * mass1 * mass2) / distance_squared |
| 48 | +} |
| 49 | + |
| 50 | +/// Given the position and mass of one body, |
| 51 | +/// calculate the force acting on it from all of the other bodies. |
| 52 | +fn forces_for_body<'a, I>(p: &Position, m: Number, reference: I) -> Force |
| 53 | +where |
| 54 | + I: IntoIterator<Item = (&'a Position, &'a Number)>, |
| 55 | +{ |
| 56 | + reference |
| 57 | + .into_iter() |
| 58 | + .map(|(ref otherpos, &othermass)| { |
| 59 | + let dx = p.x - otherpos.x; |
| 60 | + let dy = p.y - otherpos.y; |
| 61 | + let dz = p.z - otherpos.z; |
| 62 | + |
| 63 | + let d = dist_squared(dx, dy, dz); |
| 64 | + let f = force_d(m, othermass, d); |
| 65 | + |
| 66 | + Force { |
| 67 | + fx: (f * dx) / d, |
| 68 | + fy: (f * dy) / d, |
| 69 | + fz: (f * dz) / d, |
| 70 | + } |
| 71 | + }) |
| 72 | + .reduce(|acc: Force, f: Force| Force { |
| 73 | + fx: acc.fx + f.fx, |
| 74 | + fy: acc.fy + f.fy, |
| 75 | + fz: acc.fz + f.fz, |
| 76 | + }) |
| 77 | + .unwrap() |
| 78 | +} |
| 79 | + |
| 80 | +/// Calculate the accelerations for these bodies—used to update the bodies' velocities. |
| 81 | +fn accelerations(bs: &BodyStates) -> Vec<Acceleration> { |
| 82 | + bs.poss |
| 83 | + .iter() |
| 84 | + .zip(bs.masses.iter()) |
| 85 | + .map(|(ref p, &m)| { |
| 86 | + let reference = bs.poss.iter().zip(bs.masses.iter()); |
| 87 | + |
| 88 | + // calculate forces over all other things |
| 89 | + let Force { fx, fy, fz } = forces_for_body(p, m, reference); |
| 90 | + Acceleration { |
| 91 | + ax: fx / m, |
| 92 | + ay: fy / m, |
| 93 | + az: fz / m, |
| 94 | + } |
| 95 | + }) |
| 96 | + .collect() |
| 97 | +} |
| 98 | + |
| 99 | +/// Returns a new `Position` from a `Position` moving at a certain `Velocity`. |
| 100 | +fn move_position(p: &Position, v: &Velocity) -> Position { |
| 101 | + Position { |
| 102 | + x: p.x + v.dx * TIMESTEP, |
| 103 | + y: p.y + v.dy * TIMESTEP, |
| 104 | + z: p.z + v.dz * TIMESTEP, |
| 105 | + } |
| 106 | +} |
| 107 | + |
| 108 | +/// Returns a new `Velocity` from a `Velocity` accelerating at a certain `Acceleration`. |
| 109 | +fn update_velocity(v: &Velocity, a: &Acceleration) -> Velocity { |
| 110 | + Velocity { |
| 111 | + dx: v.dx + a.ax * TIMESTEP, |
| 112 | + dy: v.dy + a.ay * TIMESTEP, |
| 113 | + dz: v.dz + a.az * TIMESTEP, |
| 114 | + } |
| 115 | +} |
| 116 | + |
| 117 | +/// Computes the next `BodyStates`. |
| 118 | +pub fn compute_forces(bs: BodyStates) -> BodyStates { |
| 119 | + let accs = accelerations(&bs); |
| 120 | + BodyStates { |
| 121 | + poss: bs |
| 122 | + .poss |
| 123 | + .iter() |
| 124 | + .zip(bs.vels.iter()) |
| 125 | + .map(|(p, v)| move_position(p, v)) |
| 126 | + .collect(), |
| 127 | + vels: bs |
| 128 | + .vels |
| 129 | + .iter() |
| 130 | + .zip(accs.iter()) |
| 131 | + .map(|(v, a)| update_velocity(v, a)) |
| 132 | + .collect(), |
| 133 | + masses: bs.masses, |
| 134 | + } |
| 135 | +} |
| 136 | + |
| 137 | +/// Simple function to create a lot of bodies. |
| 138 | +/// Thank you, Larkins, for letting me use these umbers. |
| 139 | +pub fn init(count: usize) -> BodyStates { |
| 140 | + let range: Vec<Number> = (0..count).map(|i| i as Number).collect(); |
| 141 | + let ret = BodyStates { |
| 142 | + poss: range |
| 143 | + .iter() |
| 144 | + .map(|i| Position { |
| 145 | + x: 100. * (*i * 0.1), |
| 146 | + y: 200. * (*i * 0.1), |
| 147 | + z: 300. * (*i * 0.1), |
| 148 | + }) |
| 149 | + .collect(), |
| 150 | + vels: range |
| 151 | + .iter() |
| 152 | + .map(|i| Velocity { |
| 153 | + dx: 400. + *i, |
| 154 | + dy: 500. + *i, |
| 155 | + dz: 600. + *i, |
| 156 | + }) |
| 157 | + .collect(), |
| 158 | + masses: range.iter().map(|i| 10e6 * (*i + 100.2)).collect(), |
| 159 | + }; |
| 160 | + ret |
| 161 | +} |
0 commit comments