cornerstone project

% Variables

r_out = input("Please type in the value of outer radius (m): ");
r_in = input("Please type in the value of inner radius (m): ");
F = input("Please type in the value of force (N): ");
d1 = input("Please type in the value of distance 1 (m): ");
d2 = input("Please type in the value of distance 2 (m): ");
len = input("Please type in the value of length (m): ");
% r_out = 0.015; r_in = 0.01;
% F = 900; len = 0.2;
% d2 = 0.25; d1 = 0.15;

% Main code
I = (pi / 4) * (r_out^4 - r_in^4);
J = (pi / 2) * (r_out^4 - r_in^4);

M_y = F * (d2 - d1);
M_x = 2 * len * F;

tau_xy = (M_x * r_out) / J;
sigma_x = -(M_y * r_out) / I;
sigma_y = 0;

% principal stresses
sigma_2 = ((2 * F * r_out) / (pi * (r_out^4 - r_in^4))) * ((d1 - d2) + sqrt(((d1 - d2)^2) + (4 * len^2)));
fprintf("Sigma 2: %.2f MPa\n", sigma_2/10^6)
sigma_1 = ((2 * F * r_out) / (pi * (r_out^4 - r_in^4))) * ((d1 - d2) - sqrt(((d1 - d2)^2) + (4 * len^2)));
fprintf("Sigma 1: %f MPa\n", sigma_1/10^6)

theta = (180 / pi) * 0.5 * atan((tau_xy) / (sigma_x/2));
fprintf("The principal angle is : %.2f, %.2f", theta, theta + 180);

angles = 0: 0.01: (2 * pi);
sigma_x_prime = (0.5 * (sigma_x + sigma_y)) + (0.5 * cos(2 * angles) * (sigma_x - sigma_y)) + tau_xy * sin(2 * angles);
sigma_y_prime = (0.5 * (sigma_x + sigma_y)) - (0.5 * cos(2 * angles) * (sigma_x - sigma_y)) - tau_xy * sin(2 * angles);
angles = angles * (180/ pi);
plot(angles, (sigma_x_prime / 10^6), angles, (sigma_y_prime / 10^6));