Puede haber una manera, sigue siendo una manera de programa de ordenador, pero puede ser más rápido.
El método consistiría en dejar variar la dirección del eje en pequeños pasos 0-1 para $a$ , $b$ , $c$ en un bucle anidado y calcular el momento angular total de un eje. Para ello tendríamos que encontrar el vector de la línea roja para cada partícula $Q_i$
Si el COM es (0,0,0) y la ecuación del eje es $A = \left[\begin{array}{ccc}0\\0\\0\end{array}\right]+\left[\begin{array}{ccc}a\\b\\c\end{array}\right]t$
donde $t$ es un parámetro, entonces el valor de $t$ para que la línea roja sea perpendicular al eje se encuentra de la siguiente manera:
$$\left[\begin{array}{ccc}a\\b\\c\end{array}\right].\left[\begin{array}{ccc}x_i-at\\y_i-bt\\z_i-ct\end{array}\right]=0$$
que conduce al valor de $t_i$ para cada partícula
$$t_i = \frac{ax_i+by_i+cz_i}{a^2+b^2+c^2}$$
que encontraría el derecho $t_i$ entonces la línea roja tiene vector $$r_i = \left[\begin{array}{ccc}x_i-at_i\\y_i-bt_i\\z_i-ct_i\end{array}\right]$$
y el momento angular de la partícula $Q_i$ es entonces $mv_i \times r_i$ (producto cruzado), luego suma para todas las partículas usando la adición de vectores y encuentra la magnitud del resultado.
Guarda el valor y repite en el bucle anidado, y guarda los valores que superen al mejor anterior.
Digamos que el mejor $\left[\begin{array}{ccc}a\\b\\c\end{array}\right]$ fue $\left[\begin{array}{ccc}0.8\\0.4\\0.2\end{array}\right]$ a continuación, repita con un tamaño de paso más pequeño entre $\left[\begin{array}{ccc}0.7-0.9\\0.3 - 0.5\\0.1-0.3\end{array}\right]$
a continuación, reducir el tamaño del paso, repetir, etc ... tal vez podría hacerse automática. El programa se centraría entonces en la dirección del mejor eje.