Various Symbolic calculation for Physical Consistency

physconst.py

from sympy import *
from sympy.matrices import *

x1, y1, x2, y2, x3, y3, x4, y4 = symbols('x1 y1 x2 y2 x3 y3 x4 y4');

r1 = x1*x1+y1*y1;
r2 = x2*x2+y2*y2;
r3 = x3*x3+y3*y3;
r4 = x4*x4+y4*y4;

#A = Matrix([1,1,1,1],[x1,x2,x3,x4],[y1,y2,y3,y4],[r1,r2,r3,r4]);

A = eye(4);

A[0,0] = A[0,1] = A[0,2] = A[0,3] = 1;

A[1,0] = x1;
A[1,1] = x2;
A[1,2] = x3;
A[1,3] = x4;

A[2,0] = y1;
A[2,1] = y2;
A[2,2] = y3;
A[2,3] = y4;

A[3,0] = r1;
A[3,1] = r2;
A[3,2] = r3;
A[3,3] = r4;

I = eye(3);

Ixx, Iyy, Izz = symbols('Ixx Iyy Izz')
r, p, y = symbols('r p y')

R = rot_axis1(r)*rot_axis2(p)*rot_axis3(y);

I[0,0] = Ixx
I[1,1] = Iyy
I[2,2] = Izz

rI = R*I*R.transpose();

Ixy, Ixz, Iyz = symbols('Ixy Ixz Iyz')
I[1,0] = I[0,1] = Ixy;
I[2,0] = I[0,2] = Ixz;
I[2,1] = I[1,2] = Iyz;

#cholesky decomposition
l00, l10, l11, l20, l21, l22 = symbols('l00 l10 l11 l20 l21 l22');
L = eye(3);
L[0,0] = l00;
L[1,0] = l10;
L[1,1] = l11;
L[2,0] = l20;
L[2,1] = l21;
L[2,2] = l22;

L*L.transpose();

#beyer 1987
Ib = I;
Ib[0,0] = I[0,0]-x;
Ib[1,1] = I[1,1]-x;
Ib[2,2] = I[2,2]-x;