diff options
Diffstat (limited to 'v_windows/v/examples/pendulum_sim/sim.v')
-rw-r--r-- | v_windows/v/examples/pendulum_sim/sim.v | 366 |
1 files changed, 366 insertions, 0 deletions
diff --git a/v_windows/v/examples/pendulum_sim/sim.v b/v_windows/v/examples/pendulum_sim/sim.v new file mode 100644 index 0000000..f8ef11d --- /dev/null +++ b/v_windows/v/examples/pendulum_sim/sim.v @@ -0,0 +1,366 @@ +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +// sim.v * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +// created by: jordan bonecutter * * * * * * * * * * * * * * * * * * * +// jpbonecutter@gmail.com * * * * * * * * * * * * * * * * * * * * * * +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +// +// I wrote the pendulum simulator to learn V, I think it could be a +// good addition to the examples directory. +// Essentially, the pendulum sim runs a simulation of a pendulum with +// a metallic tip swinging over three magnets. +// I run this simulation with the initial position at each pixel in an +// image and color the pixel according to the magnet over which it +// finally rests. +// I used some fun features in V like coroutines, channels, +// struct embedding, mutability, methods, and the like. +import math +import os +import term +import runtime + +// customisable through setting VJOBS +const parallel_workers = runtime.nr_jobs() + +const width = 800 + +const height = 600 + +struct Vec3D { + x f64 + y f64 + z f64 +} + +fn (v Vec3D) add(v2 Vec3D) Vec3D { + return Vec3D{ + x: v.x + v2.x + y: v.y + v2.y + z: v.z + v2.z + } +} + +fn (v Vec3D) dot(v2 Vec3D) f64 { + return (v.x * v2.x) + (v.y * v2.y) + (v.z * v2.z) +} + +fn (v Vec3D) scale(scalar f64) Vec3D { + return Vec3D{ + x: v.x * scalar + y: v.y * scalar + z: v.z * scalar + } +} + +fn (v Vec3D) norm_squared() f64 { + return v.dot(v) +} + +fn (v Vec3D) norm() f64 { + return math.sqrt(v.norm_squared()) +} + +struct SimState { +mut: + position Vec3D + velocity Vec3D + accel Vec3D +} + +// magnets lie at [ +// math.cos(index * 2 * math.pi / 3) * magnet_spacing +// math.sin(index * 2 * math.pi / 3) * magnet_spacing +// -magnet_height +// ] +struct SimParams { + rope_length f64 + bearing_mass f64 + magnet_spacing f64 + magnet_height f64 + magnet_strength f64 + gravity f64 +} + +fn (params SimParams) get_rope_vector(state SimState) Vec3D { + rope_origin := Vec3D{ + x: 0 + y: 0 + z: params.rope_length + } + + return state.position.add(rope_origin.scale(-1)) +} + +fn (mut state SimState) satisfy_rope_constraint(params SimParams) { + mut rope_vector := params.get_rope_vector(state) + rope_vector = rope_vector.scale(params.rope_length / rope_vector.norm()) + state.position = Vec3D{ + x: 0 + y: 0 + z: params.rope_length + }.add(rope_vector) +} + +fn (params SimParams) get_grav_force(state SimState) Vec3D { + return Vec3D{ + x: 0 + y: 0 + z: -params.bearing_mass * params.gravity + } +} + +fn (params SimParams) get_magnet_position(theta f64) Vec3D { + return Vec3D{ + x: math.cos(theta) * params.magnet_spacing + y: math.sin(theta) * params.magnet_spacing + z: -params.magnet_height + } +} + +fn (params SimParams) get_magnet_force(theta f64, state SimState) Vec3D { + magnet_position := params.get_magnet_position(theta) + mut diff := magnet_position.add(state.position.scale(-1)) + distance_squared := diff.norm_squared() + diff = diff.scale(1.0 / math.sqrt(distance_squared)) + return diff.scale(params.magnet_strength / distance_squared) +} + +fn (params SimParams) get_magnet_dist(theta f64, state SimState) f64 { + return params.get_magnet_position(theta).add(state.position.scale(-1)).norm() +} + +fn (params SimParams) get_magnet1_force(state SimState) Vec3D { + return params.get_magnet_force(0.0 * math.pi / 3.0, state) +} + +fn (params SimParams) get_magnet2_force(state SimState) Vec3D { + return params.get_magnet_force(2.0 * math.pi / 3.0, state) +} + +fn (params SimParams) get_magnet3_force(state SimState) Vec3D { + return params.get_magnet_force(4.0 * math.pi / 3.0, state) +} + +fn (params SimParams) get_tension_force(state SimState, f_passive Vec3D) Vec3D { + rope_vector := params.get_rope_vector(state) + rope_vector_norm := rope_vector.scale(1.0 / rope_vector.norm()) + return rope_vector_norm.scale(-1.0 * rope_vector_norm.dot(f_passive)) +} + +fn (mut state SimState) increment(delta_t f64, params SimParams) { + // basically just add up all forces => + // get an accelleration => + // add to velocity => + // ensure rope constraint is satisfied + + // force due to gravity + f_gravity := params.get_grav_force(state) + + // force due to each magnet + f_magnet1 := params.get_magnet1_force(state) + + // force due to each magnet + f_magnet2 := params.get_magnet2_force(state) + + // force due to each magnet + f_magnet3 := params.get_magnet3_force(state) + + // passive forces + f_passive := f_gravity.add(f_magnet1.add(f_magnet2.add(f_magnet3))) + + // force due to tension of the rope + f_tension := params.get_tension_force(state, f_passive) + + // sum up all the fores + f_sum := f_tension.add(f_passive) + + // get the acceleration + accel := f_sum.scale(1.0 / params.bearing_mass) + state.accel = accel + + // update the velocity + state.velocity = state.velocity.add(accel.scale(delta_t)) + + // update the position + state.position = state.position.add(state.velocity.scale(delta_t)) + + // ensure the position satisfies rope constraint + state.satisfy_rope_constraint(params) +} + +fn (state SimState) done() bool { + return state.velocity.norm() < 0.05 && state.accel.norm() < 0.01 +} + +struct PPMWriter { +mut: + file os.File +} + +struct ImageSettings { + width int + height int +} + +struct Pixel { + r byte + g byte + b byte +} + +fn (mut writer PPMWriter) start_for_file(fname string, settings ImageSettings) { + writer.file = os.create(fname) or { panic("can't create file $fname") } + writer.file.writeln('P6 $settings.width $settings.height 255') or {} +} + +fn (mut writer PPMWriter) next_pixel(p Pixel) { + writer.file.write([p.r, p.g, p.b]) or {} +} + +fn (mut writer PPMWriter) finish() { + writer.file.close() +} + +fn sim_runner(mut state SimState, params SimParams) Pixel { + // do the simulation! + for _ in 0 .. 1000 { + state.increment(0.0005, params) + if state.done() { + println('done!') + break + } + } + + // find the closest magnet + m1_dist := params.get_magnet_dist(0, state) + m2_dist := params.get_magnet_dist(2.0 * math.pi / 3.0, state) + m3_dist := params.get_magnet_dist(4.0 * math.pi / 3.0, state) + + if m1_dist < m2_dist && m1_dist < m3_dist { + return Pixel{ + r: 255 + g: 0 + b: 0 + } + } else if m2_dist < m1_dist && m2_dist < m3_dist { + return Pixel{ + r: 0 + g: 255 + b: 0 + } + } else { + return Pixel{ + r: 0 + g: 0 + b: 255 + } + } +} + +struct SimResult { + id u64 + p Pixel +} + +struct SimRequest { + id u64 + params SimParams +mut: + initial SimState +} + +fn sim_worker(request_chan chan SimRequest, result_chan chan SimResult) { + // serve sim requests as they come in + for { + mut request := <-request_chan or { break } + + result_chan <- SimResult{ + id: request.id + p: sim_runner(mut request.initial, request.params) + } + } +} + +struct ValidPixel { + Pixel +mut: + valid bool +} + +fn image_worker(mut writer PPMWriter, result_chan chan SimResult, total_pixels u64) { + // as new pixels come in, write them to the image file + mut current_index := u64(0) + mut pixel_buf := []ValidPixel{len: int(total_pixels), init: ValidPixel{ + valid: false + }} + for { + result := <-result_chan or { break } + pixel_buf[result.id].Pixel = result.p + pixel_buf[result.id].valid = true + + for current_index < total_pixels && pixel_buf[current_index].valid { + writer.next_pixel(pixel_buf[current_index].Pixel) + current_index++ + } + + if current_index >= total_pixels { + break + } + } +} + +fn main() { + params := SimParams{ + rope_length: 0.25 + bearing_mass: 0.03 + magnet_spacing: 0.05 + magnet_height: 0.03 + magnet_strength: 10.0 + gravity: 4.9 + } + + mut writer := PPMWriter{} + writer.start_for_file('test.ppm', ImageSettings{ + width: width + height: height + }) + defer { + writer.finish() + } + + result_chan := chan SimResult{} + request_chan := chan SimRequest{} + + // start a worker on each core + for _ in 0 .. parallel_workers { + go sim_worker(request_chan, result_chan) + } + + go fn (request_chan chan SimRequest, params SimParams) { + mut index := u64(0) + println('') + for y in 0 .. height { + term.clear_previous_line() + println('Line: $y') + for x in 0 .. width { + // setup initial conditions + mut state := SimState{} + state.position = Vec3D{ + x: 0.1 * ((f64(x) - 0.5 * f64(width - 1)) / f64(width - 1)) + y: 0.1 * ((f64(y) - 0.5 * f64(height - 1)) / f64(height - 1)) + z: 0.0 + } + state.velocity = Vec3D{} + state.satisfy_rope_constraint(params) + request_chan <- SimRequest{ + id: index + initial: state + params: params + } + index++ + } + } + request_chan.close() + }(request_chan, params) + + image_worker(mut writer, result_chan, width * height) +} |