1- #pragma once
2- #include < immintrin.h>
3- #include " floatdef.h"
4- #include < cmath>
5- #include < string>
6-
7- struct alignas (32 ) Particle {
8- real x, y, z;
9- real vx, vy, vz;
10- real m;
11- };
12-
13- inline void Gravity (Particle& a, Particle& b, real dt)
14- {
15- constexpr real G = real (1.0 );
16- real dx = b.x - a.x ;
17- real dy = b.y - a.y ;
18- real dz = b.z - a.z ;
19-
20- real distSq0 = dx*dx + dy*dy + dz*dz;
21- real dist0 = std::sqrt (distSq0);
22-
23- // adaptive softening = 1% of separation
24- real eps = dist0 * real (0.01 );
25-
26- real distSq = distSq0 + eps*eps;
27- real dist = std::sqrt (distSq);
28-
29- real invDist = real (1 ) / dist;
30- real invDist3 = invDist * invDist * invDist;
31-
32- real ax = G * b.m * dx * invDist3;
33- real ay = G * b.m * dy * invDist3;
34- real az = G * b.m * dz * invDist3;
35-
36- real bx = -G * a.m * dx * invDist3;
37- real by = -G * a.m * dy * invDist3;
38- real bz = -G * a.m * dz * invDist3;
39-
40-
41-
42- a.vx += ax * dt;
43- a.vy += ay * dt;
44- a.vz += az * dt;
45-
46- b.vx += bx * dt;
47- b.vy += by * dt;
48- b.vz += bz * dt;
49- }
50-
51-
1+ #pragma once
2+ #include < cmath>
3+ #include " floatdef.h"
4+
5+ struct alignas (32 ) Particle {
6+ real x, y, z;
7+ real vx, vy, vz;
8+ real m;
9+ };
10+
11+ inline void Gravity (Particle& a, Particle& b, real dt)
12+ {
13+ constexpr real G = real (1.0 );
14+ constexpr real eps2 = real (1e-4 );
15+
16+ // Relative position
17+ real dx = b.x - a.x ;
18+ real dy = b.y - a.y ;
19+ real dz = b.z - a.z ;
20+
21+ // Softened squared distance
22+ real r2 = dx*dx + dy*dy + dz*dz + eps2;
23+
24+ // Inverse distance and inverse distance cubed
25+ real invR = std::sqrt (real (1.0 ) / r2);
26+ real invR3 = invR * invR * invR;
27+
28+ // Accelerations (pairwise, symmetric)
29+ real ax = G * b.m * dx * invR3;
30+ real ay = G * b.m * dy * invR3;
31+ real az = G * b.m * dz * invR3;
32+
33+ real bx = -G * a.m * dx * invR3;
34+ real by = -G * a.m * dy * invR3;
35+ real bz = -G * a.m * dz * invR3;
36+
37+ a.vx += ax * dt;
38+ a.vy += ay * dt;
39+ a.vz += az * dt;
40+
41+ b.vx += bx * dt;
42+ b.vy += by * dt;
43+ b.vz += bz * dt;
44+ }
0 commit comments