1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
|
use std::io::Write;
use std::{fs::File, path::PathBuf};
use cgmath::{Angle, ElementWise, InnerSpace, Matrix3, MetricSpace, Rad, Vector3, Zero};
type Vector = Vector3<f32>;
#[derive(Debug, Clone, Copy)]
struct Particle {
position: Vector,
position_old: Vector,
// velocity: Vector,
acceleration: Vector,
mass: f32,
fixed: bool,
}
impl Particle {
pub fn new(position: Vector, mass: f32, fixed: bool) -> Self {
Self {
position,
position_old: position,
acceleration: Vector::zero(),
mass,
fixed,
}
}
fn apply_force(&mut self, force: Vector) {
self.acceleration += force / self.mass;
}
fn tick(&mut self, dt: f32) {
if !self.fixed {
// Verlet method
let velocity = self.position - self.position_old;
self.position_old = self.position;
self.position += velocity + self.acceleration * dt * dt;
}
self.acceleration = Vector::zero();
}
}
#[derive(Debug, Clone, Copy)]
struct Spring {
rest_length: f32,
stiffness: f32,
particle_ids: [usize; 2],
}
#[derive(Debug, Clone)]
struct Simulation {
particles: Vec<Particle>,
springs: Vec<Spring>,
}
impl Simulation {
fn tick(&mut self, dt: f32) {
for spring in self.springs.iter() {
let [a, b] = spring.particle_ids;
let dist_vec = self.particles[b].position - self.particles[a].position;
let force_magnitude = spring.stiffness * (dist_vec.magnitude() - spring.rest_length);
let normalized = dist_vec.normalize();
self.particles[a].apply_force(normalized.mul_element_wise(force_magnitude));
self.particles[b].apply_force(normalized.mul_element_wise(force_magnitude * -1.0));
}
for p in &mut self.particles {
p.tick(dt);
}
}
}
struct PPM {
prefix: PathBuf,
width: usize,
height: usize,
// pixels_per_unit: usize,
}
impl PPM {
fn render_particles(&self, particles: &Vec<Particle>) -> String {
let mut s = format!("P3\n{} {}\n255\n", self.width, self.height);
let white = "255 255 255 ";
let black = "0 0 0 ";
for pixel_row in 0..self.height {
for pixel_col in 0..self.width {
let point = Vector::new(pixel_col as f32, 0.0, (pixel_row as f32) * -1.0) + Vector::new(self.width as f32 / -2.0, 0.0, self.height as f32 / 2.0);
let color = match particles.iter().any(|p| p.position.distance(point) <= p.mass / 2.0) {
true => black,
false => white,
};
s += color;
}
}
s
}
fn save_frame(&self, particles: &Vec<Particle>, time: f32) {
let file_name = format!("frame-{:08.3}", time);
let path = self.prefix.join(file_name);
let mut file = File::create(path).unwrap();
let data = self.render_particles(particles);
file.write(data.as_bytes()).unwrap();
}
}
fn main() {
let rotation_matrix = Matrix3::from_axis_angle(Vector::unit_y(), -Rad::full_turn() / 12.0);
let mut simulation = Simulation {
particles: vec![
Particle::new(Vector::zero(), 5.0, true),
Particle::new(rotation_matrix * Vector::unit_x() * 40.0, 30.0, false),
Particle::new(rotation_matrix * Vector::unit_x() * 40.0 + rotation_matrix * rotation_matrix * Vector::unit_x() * 50.0, 15.0, false)
],
springs: vec![
Spring {
particle_ids: [0, 1],
rest_length: 40.0,
stiffness: 400.0,
},
Spring {
particle_ids: [1, 2],
rest_length: 50.0,
stiffness: 5.0,
},
],
};
let ppm = PPM {
prefix: PathBuf::from("./out"),
width: 150,
height: 800,
};
let gravity = Vector::new(0.0, 0.0, -9.8);
let time_step = 0.01;
let mut tick = 1;
let mut time = 0.0;
while time < 130.0 {
{
let p = &mut simulation.particles[0];
p.apply_force(gravity * -p.mass);
}
let mut all_particles = Vec::new();
for p in &mut simulation.particles {
p.apply_force(gravity * p.mass);
}
simulation.tick(time_step);
all_particles.append(&mut simulation.particles.clone());
if tick % 10 == 0 {
ppm.save_frame(&all_particles, time);
}
time += time_step;
tick +=1;
}
}
|
