% Name: Rembrant Oyangoren Albeos
% Year: 2026
% Singly Reinforced Rectangular Beam Calculator

clc; clear;
disp('Singly Reinforced Beam Calculator');

% Inputs
fc = input('Concrete strength f''c (ksi): ');
fy = input('Steel yield strength fy (ksi): ');
b = input('Beam width b (in): ');
h = input('Beam depth h (in): ');
cover = input('Clear cover to centroid of bottom bars (in): ');
n_bars = input('Total number of bars: ');
a_bar = input('Area per bar (in^2): ');

% As calculation
As = n_bars * a_bar;
fprintf('Total As: %.3f in^2\n', As);

% Handle layers
layers = input('Double layer? (y/n): ', 's');
if lower(layers) == 'y'
    n_bot = input('Number of bars in bottom layer: ');
    n_top = input('Number of bars in top layer: ');
    spacing = input('Spacing between layers (in): ');
    
    % Centroid and depths
    g = (n_bot * a_bar * cover + n_top * a_bar * (cover + spacing)) / As;
    d = h - g;
    dt = h - cover;
    fprintf('Centroid g: %.3f in\n', g);
    fprintf('Effective d: %.3f in\n', d);
else
    d = h - cover;
    dt = d;
end

% Constants
Es = 29000;
Ecu = 0.003;

% Beta 1
if fc <= 4
    beta1 = 0.85;
elseif fc > 4 && fc <= 8
    beta1 = 0.85 - 0.05 * (fc - 4);
else
    beta1 = 0.65;
end

% Initial assumption: steel yields
T = As * fy;
a = T / (0.85 * fc * b);
c = a / beta1;

% Verify yield
eps_y = fy / Es;
eps_s = Ecu * (d - c) / c;

if eps_s < eps_y
    fprintf('Steel did not yield. Solving for exact c...\n');
    A_q = 0.85 * fc * b * beta1;
    B_q = As * Es * Ecu;
    C_q = -As * Es * Ecu * d;
    c = (-B_q + sqrt(B_q^2 - 4 * A_q * C_q)) / (2 * A_q);
    a = beta1 * c;
    eps_s = Ecu * (d - c) / c;
    T = As * (Es * eps_s);
end

% Nominal Moment
Mn = T * (d - a/2) / 12; % in kip-ft
fprintf('a: %.3f in\n', a);
fprintf('c: %.3f in\n', c);
fprintf('Mn: %.2f kip-ft\n', Mn);

% Minimum Area Check
As_min = max(3 * sqrt(fc * 1000), 200) / (fy * 1000) * b * d;
fprintf('As_min: %.3f in^2\n', As_min);
if As >= As_min
    disp('As >= As_min: OK');
else
    disp('As < As_min: FAILED');
end

% Reduction Factor
eps_t = Ecu * (dt - c) / c;
if eps_t >= 0.005
    phi = 0.90;
elseif eps_t <= 0.002
    phi = 0.65;
else
    phi = 0.65 + (eps_t - 0.002) * (250 / 3);
end
fprintf('phi: %.3f\n', phi);
fprintf('Design Strength (phi*Mn): %.2f kip-ft\n', phi * Mn);