Untitled

import sympy as sp

def find_rotation_solutions(a, b, c):
    # Define the unknown matrix elements
    x, y, w, z, p, q = sp.symbols('x y w z p q', real=True)

    # Construct the matrix using the known and unknown values
    R = sp.Matrix([
        [a, x, y],
        [p, q, b],
        [w, c, z]
    ])

    # Create the system of equations based on orthogonality and determinant condition
    equations = []

    # Each row and each column should be a unit vector
    for i in range(3):
        equations.append(sum(element**2 for element in R.row(i)) - 1)
        equations.append(sum(element**2 for element in R.col(i)) - 1)

    # Dot product between any two distinct rows and columns should be zero
    for i in range(3):
        for j in range(i+1, 3):
            equations.append(R.row(i).dot(R.row(j)))
            equations.append(R.col(i).dot(R.col(j)))

    # Cross products for rows
    equations.append(R.row(0).cross(R.row(1)) - R.row(2))
    equations.append(R.row(1).cross(R.row(2)) - R.row(0))
    equations.append(R.row(2).cross(R.row(0)) - R.row(1))
    
    # Cross products for columns
    equations.append(R.col(0).cross(R.col(1)) - R.col(2))
    equations.append(R.col(1).cross(R.col(2)) - R.col(0))
    equations.append(R.col(2).cross(R.col(0)) - R.col(1))

    # The determinant should be 1
    equations.append(R.det() - 1)

    # Solve the system of equations
    solutions = sp.solve(equations)

    return solutions

a, b, c = 0.4835, 0.1690, 0.3256
solutions = find_rotation_solutions(a, b, c)

for idx, solution in enumerate(solutions, 1):
    print(f"Solution {idx}:")
    print(solution)
    print("---------")

print(f"Total solutions: {len(solutions)}")