How to take input values from users symbolically while running a code?

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);

%%in this code,i want to take theta1 and theta2 values from the user before from doing parametric substitution(t) given below and want to show eq in terms of theta1 and theta2 and phi

eq=norm(coupler)-106;

disp(eq);

%%after doing this ,put values of theta 1 and theta2 in

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

%%and continue further as below

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(1));

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);

%%at last i want to display value of angle for both theta1 and theta2 equal to 10 degrees

%%please help someone ,i am getting no real solution if i am taking theta1 and theta2 from user as 10 degrees but no error if i am explicitly defining like in this code

Answer 1

[Athanasios Paraskevopoulos]

I think this is what you asked for

% Function to generate rotation matrix
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);

% Predefined values for theta1 and theta2
theta1_deg = 10; % Replace with desired value
theta2_deg = 10; % Replace with desired value
theta1 = deg2rad(theta1_deg); % Convert to radians
theta2 = deg2rad(theta2_deg); % Convert to radians

% Initial data and vectors
o2 = [0, 0, 0];
ain = [26, 0, 0];
input_axis = [0, 1, 0]; % Axis of rotation for ain (y-axis)
a_rotated = rot_matrix(input_axis, theta1) * (ain' - o2') + o2';
a_final = a_rotated';

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];
output_axis = [0, 1, 0]; % Axis of rotation for bin (y-axis)
b1_rotated = rot_matrix(output_axis, theta2) * (bin' - o4') + o4';
b1_final = b1_rotated';

% Calculate the coupler
coupler = c_final_rotated' - b1_final;
coupler = subs(coupler, conj(phi), phi);

% Show the equation in terms of theta1, theta2, and phi
eq = norm(coupler) - 106;
disp('Equation in terms of theta1, theta2, and phi:');

Equation in terms of theta1, theta2, and phi:

disp(eq);

% Continue further calculations with substituted values of theta1 and theta2
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]);
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));
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);

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

Angles in degrees before adjustment:

disp(angles_deg);

phi = double(angles_rad(1));
c1_position = double(rot_matrix(naxis, phi) * (c_rotated - a_rotated) + a_rotated);
p = (c1_position' - a_final)';
angle = acosd(p(2) / norm(p));

disp('Final angle:');

Final angle:

disp(angle);

   88.7282

[Aman]

@Athanasios Paraskevopoulos sir thank you for your answer,but i wanted to disp(eq) in terms of variables theta1 and theta2 and phi .and after displaying eq in the term of these,i want to solve for phi by doing parametric substitution by asking user to input values of theta1 and theta2 after displaying eq .

[Athanasios Paraskevopoulos]

`@Aman I think you mean something like the following code

% Function to generate rotation matrix
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);

% Initial data and vectors
o2 = [0, 0, 0];
ain = [26, 0, 0];
input_axis = [0, 1, 0]; % Axis of rotation for ain (y-axis)

% Symbolic variables for theta1, theta2, and phi
syms theta1 theta2 phi;

% Calculate rotated vectors symbolically
a_rotated = rot_matrix(input_axis, theta1) * (ain' - o2') + o2';
a_final = a_rotated';

cin = [122.95, -20, 0];
c_rotated = rot_matrix(input_axis, theta1) * (cin' - o2') + o2';
naxis = [sin(theta1), 0, cos(theta1)];
c_final_rotated = rot_matrix(naxis, phi) * (c_rotated - a_rotated) + a_rotated;

bin = [29.5, 30, 0];
o4 = [13.5, 30, 0];
output_axis = [0, 1, 0]; % Axis of rotation for bin (y-axis)
b1_rotated = rot_matrix(output_axis, theta2) * (bin' - o4') + o4';
b1_final = b1_rotated';

% Calculate the coupler symbolically
coupler = c_final_rotated' - b1_final;
coupler = subs(coupler, conj(phi), phi);

% Show the equation in terms of theta1, theta2, and phi
eq = norm(coupler) - 106;
disp('Equation in terms of theta1, theta2, and phi:');

Equation in terms of theta1, theta2, and phi:

disp(eq);

% Predefined values for theta1 and theta2
theta1_deg = 10; % Replace with desired value
theta2_deg = 10; % Replace with desired value
theta1_val = deg2rad(theta1_deg); % Convert to radians
theta2_val = deg2rad(theta2_deg); % Convert to radians

% Substitute the values of theta1 and theta2
coupler_subs = subs(coupler, {theta1, theta2}, {theta1_val, theta2_val});

% Continue further calculations with substituted values of theta1 and theta2
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_subs, [cos(phi), sin(phi)], [cos_phi, sin_phi]);
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));
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);

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

Angles in degrees before adjustment:

disp(angles_deg);

if ~isempty(angles_rad)
    phi = double(angles_rad(1));
    c1_position = double(rot_matrix(naxis, phi) * (c_rotated - a_rotated) + a_rotated);
    p = (c1_position' - a_final)';
    angle = acosd(p(2) / norm(p));
    disp('Final angle:');
    disp(angle);
else
    disp('No real solution found for phi.');
end

No real solution found for phi.

[Aman]

@Athanasios Paraskevopoulos sir can you answer please now why the code provided by you before is giving the real solutions but the later one code doesn't for same values of theta1 and theta2 ?

[Athanasios Paraskevopoulos]

@Amani I just did what you asked for. You wanted disp(eq) in terms of variables theta1 and theta2 and phi and then solve for phi by doing parametric substitution. The first one does not give you that

[Aman]

@Athanasios Paraskevopoulos thank you sir,but why i am not getting real solutions in 2nd code instead i am giving the same input values,if you can help,it will be nice.

Thanks!!