Coordinate Conversion Equation Solution Issue

This program aims to represent the points of the new coordinate system by t, where the rotation angle theta and the translation length are given

the equation definitely has a solution, but I don't know why it cannot be solved.

theta_val = 30;

x_0_val = 10;

y_0_val = 20;

syms theta x_0 y_0 x_1 y_1 y__1 x__1 t x_2

equation1 = x_1 == x_0 + sin(theta + atan(x__1/y__1)) * sqrt(x__1^2 + y__1^2);

equation2 = y_1 == y_0 + cos(theta + atan(x__1/y__1)) * sqrt(x__1^2 + y__1^2);

equation1_sub = subs(equation1, [theta, x_0, y_0], [theta_val, x_0_val, y_0_val]);

equation2_sub = subs(equation2, [theta, x_0, y_0], [theta_val, x_0_val, y_0_val]);

ship_x = 2*t^2 + 2;

ship_y = 2*t + 3;

equation1_transformed = subs(equation1_sub, [x_1, y_1], [ship_x, ship_y]);

equation2_transformed = subs(equation2_sub, [x_1, y_1], [ship_x, ship_y]);

solution = solve([equation1_transformed, equation2_transformed], [x__1, y__1]);

Answer 1

[Torsten]

Here is a partial solution, but note that equation1_transformed and equation2_transformed are not satisfied for all values of t by the two solutions that the Symbolic Toolbox returns for the derived equations eqn1 and eqn2.

theta_val = 30*pi/180;
x_0_val = 10;
y_0_val = 20;
syms theta x_0 y_0 x_1 y_1 y__1 x__1 t x_2
equation1 = x_1 == x_0 + sin(theta + atan(x__1/y__1)) * sqrt(x__1^2 + y__1^2);
equation2 = y_1 == y_0 + cos(theta + atan(x__1/y__1)) * sqrt(x__1^2 + y__1^2);
equation1_sub = subs(equation1, [theta, x_0, y_0], [theta_val, x_0_val, y_0_val]);
equation2_sub = subs(equation2, [theta, x_0, y_0], [theta_val, x_0_val, y_0_val]);
ship_x = 2*t^2 + 2;
ship_y = 2*t + 3;
equation1_transformed = subs(equation1_sub, [x_1, y_1], [ship_x, ship_y])

equation1_transformed = 

equation2_transformed = subs(equation2_sub, [x_1, y_1], [ship_x, ship_y])

equation2_transformed = 

%solution = solve([equation1_transformed, equation2_transformed], [x__1, y__1]);
eqn1 = (2*t^2+2-10)^2+(2*t+3-20)^2 == x__1^2+y__1^2

eqn1 = 

eqn2 = tan(atan((2*t^2+2-10)/(2*t+3-20)) - pi/6) == x__1/y__1

eqn2 = 

solution = solve([eqn1,eqn2],[x__1,y__1])

solution = struct with fields:
    x__1: [2×1 sym]
    y__1: [2×1 sym]

solution.x__1

ans = 

solution.y__1

ans = 

simplify(subs([eqn1,eqn2],[x__1,y__1],[solution.x__1(1),solution.y__1(1)]))

ans = 

simplify(subs([eqn1,eqn2],[x__1,y__1],[solution.x__1(2),solution.y__1(2)]))

ans = 

simplify(subs([equation1_transformed,equation2_transformed],[x__1,y__1],[solution.x__1(1),solution.y__1(1)]))

ans = 

simplify(subs([equation1_transformed,equation2_transformed],[x__1,y__1],[solution.x__1(2),solution.y__1(2)]))

ans =