Finally, physically-accurate simulation of one-fourth-of-a-cylinder

src/l3/r4/r4.c

@@ -6,6 +6,15 @@
 
 #include "linux/input-event-codes.h"
 
+
+typedef mat3 quad_form3_t;
+
+float quad_form3_mul_vec(quad_form3_t Q, vec3 M) {
+    return M.x * M.x * Q.x.x * Q.x.x + M.y * M.y * Q.y.y * Q.y.y + M.z * M.z * Q.z.z * Q.z.z +
+        2 * (M.x * M.y * Q.x.y + M.x * M.z * Q.x.z + M.y * M.x * Q.y.z);
+}
+
+
 typedef struct{
     vec3 pos;
     vec3 color;
@@ -16,7 +25,7 @@ typedef struct{
 
 typedef struct {
     float mass;  /* In kg, Not zero */
-    float moment_of_inertia; /* In kg*m^2 */
+    quad_form3_t inertia_moment; /* Qadratic form, yields kg*m^2 */
     /* Center of mass relative to "center of mesh" */
     vec3 mass_center;
 } RigidBodyConfig;
@@ -126,7 +135,10 @@ void RigidBodyState_when_shot(RigidBodyState* self, vec3 imp, float m2, vec3 v){
     self->speed = vec3_add_vec3(self->speed, linear_speed_gain);
 
 
-    vec3 angular_speed_gain = vec3_mul_scal(vec3_cross(v, IO), m2 / self->p.moment_of_inertia);
+    vec3 www = vec3_mul_scal(vec3_cross(v, IO), m2);
+    vec3 w = vec3_normalize(www);
+    float I = quad_form3_mul_vec(self->p.inertia_moment, w);
+    vec3 angular_speed_gain = vec3_div_by_scal(www, I);
     self->angular_speed = vec3_add_vec3(self->angular_speed, angular_speed_gain);
 }
 
@@ -325,8 +337,15 @@ void run_app(){
             .specular_texture_path = vcstr("./src/l3/textures/log_10_2_6_specular.png")
         });
     AliceGenericMeshHand_resize_instance_arr(st.alice, &st.ROA_mesh->el, 1);
+    const float gamma_l_c = 4.f / 3 / M_PIf;
     st.ROA_state = (RigidBodyState){
-        .p.mass = 50.f, .p.moment_of_inertia = 5000.f, .p.mass_center = (vec3){5, 0.77f, -0.77f},
+        .p.mass = 1.f * 1/4 * M_PIf * 2*2 * 10,
+        .p.inertia_moment = (mat3){
+            .x = { 16.926308911016392, 1.5411267097545367e-13, -9.491580156388568e-15},
+            .y = { 1.5411267097545367e-13,  259.0981544556158,  2.539003832491438},
+            .z = {-9.491580156388568e-15, 2.539003832491438, 259.09815445561975},
+        },
+        .p.mass_center = (vec3){5, 2 * gamma_l_c, -2 * gamma_l_c},
         .pos = (vec3){11.f, 3, 4}, .speed = {0},
         .rot = mat3_E, .angular_speed = (vec3){0}};