sciforge_lib/maths/vector/
fields.rs

1use super::types::{Particle, Vec3};
2use crate::constants::{G, K_COULOMB};
3
4pub fn gravitational_force(p1: &Particle, p2: &Particle) -> Vec3 {
5    let r = &p2.position - &p1.position;
6    let dist = r.magnitude();
7    if dist < 1e-30 {
8        return Vec3::zero();
9    }
10    r.normalize().scale(G * p1.mass * p2.mass / (dist * dist))
11}
12
13pub fn coulomb_force(p1: &Particle, p2: &Particle) -> Vec3 {
14    let r = &p2.position - &p1.position;
15    let dist = r.magnitude();
16    if dist < 1e-30 {
17        return Vec3::zero();
18    }
19    r.normalize()
20        .scale(-K_COULOMB * p1.charge * p2.charge / (dist * dist))
21}
22
23pub fn lorentz_force(charge: f64, velocity: &Vec3, e_field: &Vec3, b_field: &Vec3) -> Vec3 {
24    &e_field.scale(charge) + &velocity.cross(b_field).scale(charge)
25}
26
27pub fn spring_force(k: f64, rest_length: f64, p1: &Vec3, p2: &Vec3) -> Vec3 {
28    let r = p2 - p1;
29    let dist = r.magnitude();
30    if dist < 1e-30 {
31        return Vec3::zero();
32    }
33    r.normalize().scale(k * (dist - rest_length))
34}
35
36pub fn drag_force(cd: f64, rho: f64, area: f64, velocity: &Vec3) -> Vec3 {
37    let speed = velocity.magnitude();
38    if speed < 1e-30 {
39        return Vec3::zero();
40    }
41    velocity
42        .normalize()
43        .scale(-0.5 * cd * rho * area * speed * speed)
44}
45
46pub fn dipole_field(moment: &Vec3, position: &Vec3) -> Vec3 {
47    let r = position.magnitude();
48    if r < 1e-30 {
49        return Vec3::zero();
50    }
51    let r2 = r * r;
52    let r5 = r2 * r2 * r;
53    let mdotr = moment.dot(position);
54    let term1 = position.scale(3.0 * mdotr / r5);
55    let term2 = moment.scale(1.0 / (r2 * r));
56    &term1 - &term2
57}
58
59pub fn gradient_scalar_field(
60    f: impl Fn(f64, f64, f64) -> f64,
61    x: f64,
62    y: f64,
63    z: f64,
64    h: f64,
65) -> Vec3 {
66    let dx = (f(x + h, y, z) - f(x - h, y, z)) / (2.0 * h);
67    let dy = (f(x, y + h, z) - f(x, y - h, z)) / (2.0 * h);
68    let dz = (f(x, y, z + h) - f(x, y, z - h)) / (2.0 * h);
69    Vec3::new(dx, dy, dz)
70}
71
72pub fn divergence(
73    fx: impl Fn(f64, f64, f64) -> f64,
74    fy: impl Fn(f64, f64, f64) -> f64,
75    fz: impl Fn(f64, f64, f64) -> f64,
76    x: f64,
77    y: f64,
78    z: f64,
79    h: f64,
80) -> f64 {
81    let dfdx = (fx(x + h, y, z) - fx(x - h, y, z)) / (2.0 * h);
82    let dfdy = (fy(x, y + h, z) - fy(x, y - h, z)) / (2.0 * h);
83    let dfdz = (fz(x, y, z + h) - fz(x, y, z - h)) / (2.0 * h);
84    dfdx + dfdy + dfdz
85}
86
87pub fn curl(
88    fx: impl Fn(f64, f64, f64) -> f64,
89    fy: impl Fn(f64, f64, f64) -> f64,
90    fz: impl Fn(f64, f64, f64) -> f64,
91    x: f64,
92    y: f64,
93    z: f64,
94    h: f64,
95) -> Vec3 {
96    let dfz_dy = (fz(x, y + h, z) - fz(x, y - h, z)) / (2.0 * h);
97    let dfy_dz = (fy(x, y, z + h) - fy(x, y, z - h)) / (2.0 * h);
98    let dfx_dz = (fx(x, y, z + h) - fx(x, y, z - h)) / (2.0 * h);
99    let dfz_dx = (fz(x + h, y, z) - fz(x - h, y, z)) / (2.0 * h);
100    let dfy_dx = (fy(x + h, y, z) - fy(x - h, y, z)) / (2.0 * h);
101    let dfx_dy = (fx(x, y + h, z) - fx(x, y - h, z)) / (2.0 * h);
102    Vec3::new(dfz_dy - dfy_dz, dfx_dz - dfz_dx, dfy_dx - dfx_dy)
103}
104
105pub fn laplacian_scalar(f: impl Fn(f64, f64, f64) -> f64, x: f64, y: f64, z: f64, h: f64) -> f64 {
106    let fxyz = f(x, y, z);
107    let d2x = (f(x + h, y, z) - 2.0 * fxyz + f(x - h, y, z)) / (h * h);
108    let d2y = (f(x, y + h, z) - 2.0 * fxyz + f(x, y - h, z)) / (h * h);
109    let d2z = (f(x, y, z + h) - 2.0 * fxyz + f(x, y, z - h)) / (h * h);
110    d2x + d2y + d2z
111}
112
113pub fn potential_energy_field(
114    force: impl Fn(&Vec3) -> Vec3,
115    path_start: &Vec3,
116    path_end: &Vec3,
117    steps: usize,
118) -> f64 {
119    let mut total = 0.0;
120    for i in 0..steps {
121        let t0 = i as f64 / steps as f64;
122        let t1 = (i + 1) as f64 / steps as f64;
123        let p0 = lerp_vec3(path_start, path_end, t0);
124        let p1 = lerp_vec3(path_start, path_end, t1);
125        let f_mid = force(&lerp_vec3(path_start, path_end, 0.5 * (t0 + t1)));
126        let dr = &p1 - &p0;
127        total += f_mid.dot(&dr);
128    }
129    -total
130}
131
132fn lerp_vec3(a: &Vec3, b: &Vec3, t: f64) -> Vec3 {
133    Vec3::new(
134        a.x + t * (b.x - a.x),
135        a.y + t * (b.y - a.y),
136        a.z + t * (b.z - a.z),
137    )
138}
139
140pub fn streamline(
141    vx: impl Fn(f64, f64, f64) -> f64,
142    vy: impl Fn(f64, f64, f64) -> f64,
143    vz: impl Fn(f64, f64, f64) -> f64,
144    start: &Vec3,
145    dt: f64,
146    steps: usize,
147) -> Vec<Vec3> {
148    let mut points = Vec::with_capacity(steps + 1);
149    let mut p = start.clone();
150    points.push(p.clone());
151    for _ in 0..steps {
152        let v = Vec3::new(vx(p.x, p.y, p.z), vy(p.x, p.y, p.z), vz(p.x, p.y, p.z));
153        p = &p + &v.scale(dt);
154        points.push(p.clone());
155    }
156    points
157}
sciforge_lib/maths/vector/fields.rs

sciforge_lib/maths/vector/
fields.rs
