R = input('Enter radii [R1, R2, R3]: ');
theta = input('Enter angles [theta1, theta2, theta3] in degrees: ');
stopping_criterion = input('Enter stopping criterion for Newton-Raphson iterations: ');
disp('Choose a solution method:');
disp('1. Gauss-Elimination with partial pivoting');
disp('3. Naïve Gauss-Jordan');
method = input('Enter the method number: ');
if method < 1 || method > 3
error('Invalid method number');
f = @(params) [params(1)/(1 + params(2) * sind(theta(1) + params(3))) - R(1);
params(1)/(1 + params(2) * sind(theta(2) + params(3))) - R(2);
params(1)/(1 + params(2) * sind(theta(3) + params(3))) - R(3)];
relative_errors_C = zeros(100, 1);
relative_errors_e = zeros(100, 1);
relative_errors_alpha = zeros(100, 1);
J = numerical_jacobian(f, [C, e, alpha]);
error('The Jacobian matrix is not invertible. The Newton-Raphson method cannot be used.');
update = gauss_elimination_pivoting(J, f([C, e, alpha]));
update = cramers_rule_solver(J, f([C, e, alpha]));
update = naive_gauss_jordan(J, f([C, e, alpha]));
if norm(update) < stopping_criterion
alpha = alpha - update(3);
relative_errors_C(iteration) = abs(update(1) / C);
relative_errors_e(iteration) = abs(update(2) / e);
relative_errors_alpha(iteration) = abs(update(3) / alpha);
fprintf('Iteration %d: Relative Error in C = %f, e = %f, alpha = %f\n', iteration, relative_errors_C(iteration), relative_errors_e(iteration), relative_errors_alpha(iteration));
plot(1:iteration, relative_errors_C(1:iteration), 'o-');
ylabel('Relative Error in C');
plot(1:iteration, relative_errors_e(1:iteration), 'o-');
ylabel('Relative Error in e');
plot(1:iteration, relative_errors_alpha(1:iteration), 'o-');
ylabel('Relative Error in alpha');
fprintf('Final values:\nC = %f\n e = %f\n alpha = %f\n', C, e, alpha);
theta_values = linspace(0, 360, 1000);
R_values = C./(1 + e * sind(theta_values + alpha));
polarplot(deg2rad(theta_values), R_values);
polarplot(deg2rad(theta), R, 'ro');
title('Satellite Orbit');
legend('Orbit', 'Data Points');
function J = numerical_jacobian(f, params)
params_plus(i) = params_plus(i) + delta;
params_minus(i) = params_minus(i) - delta;
J(:, i) = (f(params_plus) - f(params_minus)) / (2 * delta);