How to simulate code in symbolic form?

o2 = [0, 0, 0]; % Origin for ain

ain = [26, 0, 0]; % Initial vector for ain

input_axis = [0, 1, 0]; % Axis of rotation for ain (y-axis)

theta1 = deg2rad(10); % Angle of rotation for ain in radians

% Rotation matrix function

rot_matrix = @(axis, theta) cos(theta) * eye(3) + ...

sin(theta) * [0, -axis(3), axis(2); axis(3), 0, -axis(1); -axis(2), axis(1), 0] + ...

(1 - cos(theta)) * (axis' * axis);

% Compute the rotated vector for ain

a_rotated = rot_matrix(input_axis, theta1) * (ain' - o2') + o2';

a_final = a_rotated';

disp(norm(a_final));

cin=[122.95, -20, 0];

c_rotated= rot_matrix(input_axis, theta1) * (cin' - o2') + o2';

naxis=[sin(theta1),0,cos(theta1)];

syms phi;

c_final_rotated=rot_matrix(naxis,phi)*(c_rotated-a_rotated)+a_rotated;

bin = [29.5, 30, 0];

o4 = [13.5, 30, 0]; % Origin for bin

output_axis = [0, 1, 0]; % Axis of rotation for bin (y-axis)

theta2 = deg2rad(10); % Angle of rotation for bin in radians

% Compute the rotated vector for bin

b1_rotated = rot_matrix(output_axis, theta2) * (bin' - o4') + o4';

b1_final = b1_rotated';

disp(norm(b1_final-o4));

coupler = c_final_rotated'-b1_final;

coupler = subs(coupler, conj(phi), phi);

syms t;

cos_phi = (1 - t^2) / (1 + t^2);

sin_phi = 2 * t / (1 + t^2);

% Substitute parametric forms into coupler components

coupler_parametric = subs(coupler, [cos(phi), sin(phi)], [cos_phi, sin_phi]);

% Display the parametric coupler

disp('Parametric form of coupler:');

disp(coupler_parametric);

syms targetvalue % it might be 3.5 ...

normsq = expand(sum(coupler_parametric.^2) - targetvalue^2);

normpoly = simplify(normsq*(t^2+1)^2);

vpa(expand(normpoly),4);

tsolve = solve(normpoly,t,'maxdegree',4,'returnconditions',true);

h=vpa(subs(tsolve.t,targetvalue, 106));

%disp(h);

real_solutions = h(imag(h) == 0);

disp('Real roots:');

disp(real_solutions);

angles_rad = 2 * atan(real_solutions);

angles_deg = rad2deg(angles_rad);

% Display angles in degrees

disp('Angles in degrees before adjustment:');

disp(angles_deg);

phi=double(angles_rad(2));

c1_position = double(rot_matrix(naxis,phi) * (c_rotated - a_rotated) + a_rotated);

p=(c1_position'-a_final)';

%q=(c1_position'-b1_final)';

angle=acosd(p(2)/norm(p));

disp(angle);

%%I want to run this code for symbolic variables including o2,ain,o4,bin,input_axis,output_axis,naxis,theta1,theta2.

can someone tell how can i do all above calculations in symbolic form

Answer 1

[Aquatris]

Here is one way, though it cannot solve it at one point I think so throws an error when asked to display the solution

% assume real to prevent complicating things
syms o2 ain input_axis bin o4 output_axis [1 3] real
syms theta1 theta2 real

% Rotation matrix function
rot_matrix = @(axis, theta) cos(theta) * eye(3) + ...
    sin(theta) * [0, -axis(3), axis(2); axis(3), 0, -axis(1); -axis(2), axis(1), 0] + ...
    (1 - cos(theta)) * (axis' * axis);
% Compute the rotated vector for ain
a_rotated = rot_matrix(input_axis, theta1) * (ain' - o2') + o2';
a_final = a_rotated';
disp(norm(a_final));

cin=[122.95, -20, 0];
c_rotated= rot_matrix(input_axis, theta1) * (cin' - o2') + o2';
naxis=[sin(theta1),0,cos(theta1)];
syms phi;
c_final_rotated=rot_matrix(naxis,phi)*(c_rotated-a_rotated)+a_rotated;
% Compute the rotated vector for bin
b1_rotated = rot_matrix(output_axis, theta2) * (bin' - o4') + o4';
b1_final = b1_rotated';
disp(norm(b1_final-o4));

coupler = c_final_rotated'-b1_final;
coupler = subs(coupler, conj(phi), phi);
syms t;
cos_phi = (1 - t^2) / (1 + t^2);
sin_phi = 2 * t / (1 + t^2);
% Substitute parametric forms into coupler components
coupler_parametric = subs(coupler, [cos(phi), sin(phi)], [cos_phi, sin_phi]);
% Display the parametric coupler
disp('Parametric form of coupler:');

Parametric form of coupler:

disp(coupler_parametric);

syms targetvalue % it might be 3.5 ...
normsq = expand(sum(coupler_parametric.^2) - targetvalue^2);
normpoly = simplify(normsq*(t^2+1)^2);
vpa(expand(normpoly),4);
tsolve = solve(normpoly,t,'maxdegree',4,'returnconditions',true);
h=vpa(subs(tsolve.t,targetvalue, 106));
%disp(h);
real_solutions = h(imag(h) == 0);
disp('Real roots:');

Real roots:

disp(real_solutions);
angles_rad = 2 * atan(real_solutions);
angles_deg = rad2deg(angles_rad);
% Display angles in degrees
disp('Angles in degrees before adjustment:');

Angles in degrees before adjustment:

disp(angles_deg);
phi=double(angles_rad(2));

Index exceeds the number of array elements. Index must not exceed 0.

Error in indexing (line 962)
R_tilde = builtin('subsref',L_tilde,Idx);

c1_position = double(rot_matrix(naxis,phi) * (c_rotated - a_rotated) + a_rotated);
p=(c1_position'-a_final)';
%q=(c1_position'-b1_final)';
angle=acosd(p(2)/norm(p));
disp(angle);

[Aman]

@Aquatris Thank you very much sir, if i give these symbolic variables a numerical value after "disp(coupler_parametric);" ,then can i get solution?

Please Take these values to substitute in symbolic variables

o2 = [0, 0, 0]; % Origin for ain

ain = [26, 0, 0]; % Initial vector for ain

input_axis = [0, 1, 0]; % Axis of rotation for ain (y-axis)

bin = [29.5, 30, 0];

o4 = [13.5, 30, 0]; % Origin for bin

output_axis = [0, 1, 0]; % Axis of rotation for bin (y-axis)

[Aquatris]

It depends on a couple of things, one of which is if there exist a solution. you can plugin values using subs() function and try.