%% ========================================================================
%  PURDUE BAJA RACING -- STEERING & BRAKES
%  Brake Rotor Thermal Analysis
%  ------------------------------------------------------------------------
%  Implements the DRB "Thermal Analysis" flowchart:
%   1) Total braking energy from vehicle speed/mass (or telemetry)
%   2) Energy split between rotor and pad (thermal effusivity)
%   3) Single-stop rotor temperature rise (conservative, adiabatic)
%   4) Endurance simulation: repeated stops + convective cooling
%      (worst case: rotor treated as stationary, i.e. no forced convection)
%   5) Thermal stress / expansion check -> failure temperature
%   6) Rotor sizing: minimum mass for thermal capacity, pocket removal sweep
%
%  Every section marked "EDIT ME" needs YOUR numbers (CAD mass, dyno data,
%  material specs). Everything else is generic and shouldn't need changes.
% ========================================================================

clear; clc; close all;

%% ------------------------------------------------------------------
%  VARIABLE / ABBREVIATION GLOSSARY
%  ------------------------------------------------------------------
%  STRUCTS (group related variables):
%    veh     = vehicle              endur = endurance (event)
%    rotor   = brake rotor          pad   = brake pad
%    mu      = friction coefficient env   = environment
%
%  VEHICLE / EVENT FIELDS:
%    v_top, v_bottom     = velocity at top/bottom of a brake event (m/s)
%    duration_s          = duration, seconds
%    brake_events_hr     = number of brake events per hour
%    brake_time_s        = duration of one braking event (s)
%
%  MATERIAL PROPERTY FIELDS (standard engineering symbols):
%    rho (rho)  = density                      k     = thermal conductivity
%    cp         = specific heat capacity        E     = elastic (Young's) modulus
%    alpha      = coefficient of thermal expansion
%    nu         = Poisson's ratio                yield = yield strength
%    OD / ID    = outer / inner diameter
%    A_contact  = pad-rotor contact area (interface)
%    A_surface  = total exposed surface area (for convection)
%
%  ENERGY / HEAT VARIABLES:
%    Q               = heat energy, Joules (standard thermo symbol)
%    Q_total_single  = total kinetic energy dissipated per brake stop
%    Q_rotor_single  = rotor's share of that energy
%    Q_in_profile    = time-varying heat input array (endurance sim)
%    Q_loss          = convective heat lost to ambient air
%    e_rotor, e_pad  = thermal effusivity of rotor / pad material
%    gamma_d, gamma_p = fraction of heat into disc (rotor) vs pad
%                       (subscript d = disc/rotor, p = pad -- matches DRB slide notation)
%
%  TEMPERATURE VARIABLES:
%    dT           = delta-T, a temperature CHANGE (not absolute temp)
%    T_sweep_C    = array of temperatures swept through, in Celsius
%    T_rotor      = rotor temperature over time, in Celsius
%    T_rotor_C    = same array (kept for naming consistency with other _C vars)
%    T_peak_C     = peak rotor temperature reached
%    T_fail_C     = estimated rotor failure temperature
%    T_ambient_C  = ambient air temperature, Celsius
%
%  SIMULATION / TIME VARIABLES:
%    dt          = time step size (s)
%    t_vec       = time vector (array of time points)
%    N           = number of time steps
%    n_events    = number of simulated brake events
%    event_times = randomly generated timestamps for each event
%    idx_peak, idx_start, idx_end = array indices (position in a list)
%    h_conv      = convective heat transfer coefficient, W/(m^2*K)
%
%  STRESS / SIZING VARIABLES:
%    grad_frac      = gradient fraction (how much of bulk dT becomes
%                     a radial temperature gradient across the rotor)
%    sigma_thermal  = thermal stress (standard stress symbol, sigma)
%    FOS            = factor of safety
%    m_min          = minimum rotor mass (for thermal capacity)
%    mass_removed_pct = percent of rotor mass removed via pocketing
%  ------------------------------------------------------------------

%% ------------------------------------------------------------------
%  SECTION 1: VEHICLE + EVENT INPUTS                          [EDIT ME]
%  ------------------------------------------------------------------
veh.mass          = 177.8;    % kg, from Overall DCR Falcon Baseline: actual measured "Weight: 392 lbs comp-ready" (measured, not goal target)
veh.v_top         = 15.83;    % m/s, from Overall DCR Falcon Baseline: actual measured "Top Speed: 35.4 2WD, 35.1 4WD" -- used 35.4 MPH (2WD)
veh.v_bottom      = 2;        % m/s, speed after braking event (not full stop) -- still no documented value, EDIT
veh.wheel_radius  = 0.279;    % m, tire rolling radius (~11 in)

% Endurance event assumptions (used in Section 4 if you don't import
% real telemetry -- replace this block with actual dyno/data-log import)
endur.duration_s      = 4*3600;   % 4 hr endurance, seconds
endur.brake_events_hr = 60;       % derived, not guessed: Overall DCR "Purdue AVG Lap Time 3:51.4" (231 s)
                                   % -> 3600/231 = 15.6 laps/hr; assuming ~4 brake
                                   % zones/lap (engineering estimate, not documented) -> ~60 events/hr.
                                   % EDIT once we have real endurance brake-pressure event counts.
endur.brake_time_s    = 1.5;      % avg duration of each braking event (still undocumented, EDIT)

%% ------------------------------------------------------------------
%  SECTION 2: ROTOR MATERIAL + GEOMETRY                        [EDIT ME]
%  ------------------------------------------------------------------
rotor.material   = 'Steel 4130 (annealed)';
rotor.rho        = 7850;      % kg/m^3
rotor.k          = 42.7;      % W/(m*K) thermal conductivity
rotor.cp         = 477;       % J/(kg*K) specific heat
rotor.E          = 205e9;     % Pa, elastic modulus
rotor.alpha      = 12.2e-6;   % 1/K, coefficient of thermal expansion
rotor.nu         = 0.29;      % Poisson's ratio
rotor.yield      = 460e6;     % Pa, yield strength (adjust for temper)
rotor.mass       = 0.451;      % kg, CURRENT rotor mass from CAD (EDIT)
rotor.OD         = 0.1778;     % m, outer swept diameter (EDIT)
rotor.ID         = 0.0782;     % m, inner swept diameter (EDIT)
rotor.thickness  = 0.002286;     % m, rotor thickness (EDIT)
rotor.A_contact  = pi/4*(rotor.OD^2 - rotor.ID^2) * 0.15; % m^2, approx
                                % swept pad contact area (EDIT: from pad geo)
rotor.A_surface  = 2*(pi/4*(rotor.OD^2-rotor.ID^2)) * 1.3; % both faces
                                % +30% fudge for vents/pockets, EDIT if known

%% ------------------------------------------------------------------
%  SECTION 3: PAD MATERIAL + GEOMETRY            (CAD)              [EDIT ME]
%  ------------------------------------------------------------------
pad.material  = 'Sintered metallic (GMP-type)';
pad.rho       = 5000;         % kg/m^3   (typical sintered pad, EDIT if known)
pad.k         = 2.5;          % W/(m*K)  (composite pads: ~0.5-3, sintered higher)
pad.cp        = 800;          % J/(kg*K)
pad.A_contact = rotor.A_contact; % m^2, assume same interface area
mu.pad        = 0.35;         % coefficient of friction (from your DRJ testing, EDIT)

%% ------------------------------------------------------------------
%  SECTION 4: ENVIRONMENT
%  ------------------------------------------------------------------
env.T_ambient_C = 25;  % C, ambient reference temperature (Kelvin removed --
                        % thermal deltaT is identical in K or C, so tracking
                        % Celsius directly through the simulation removes an
                        % unnecessary conversion step and possible confusion)

%% ------------------------------------------------------------------
%  SECTION 5: BRAKE DISTRIBUTION ACROSS ROTORS                [EDIT ME]
%  ------------------------------------------------------------------
%  IMPORTANT PHYSICS FIX: the vehicle's total KE is shared across
%  multiple rotors (front-left, front-right, rear), not absorbed by
%  a single rotor. Standard Baja layout: 2 front rotors (one per
%  wheel) + 1 rear inboard rotor (shared across the rear axle via
%  the differential) -- EDIT n_rear_rotors if your car uses dual
%  rear rotors instead.
%
%  Bias source: S&B design poster, "Actual Value" for brake pressure
%  bias = 65R-35F. (DCR goal range was 63R-37F; poster shows the
%  achieved value.)
brake.bias_rear      = 0.65;  % rear fraction of total brake force, from design poster
brake.bias_front     = 1 - brake.bias_rear;
brake.n_front_rotors = 2;     % EDIT if not one rotor per front wheel
brake.n_rear_rotors  = 1;     % EDIT if dual rear rotors instead of single inboard

% Per-rotor energy fraction of total vehicle KE, for front vs rear:
%
% Torque correction: the 65/35 figure is a HYDRAULIC PRESSURE bias --
% it sets the ratio of clamp FORCE front vs rear, not the ratio of
% braking TORQUE or dissipated energy directly. Torque = F_clamp * mu *
% r_eff, so if the front and rear rotors have different effective
% friction radii, equal-ratio clamp force does not produce the same
% ratio of braking work. This section weights the bias fraction by
% effective pad radius to get the actual energy split, rather than
% assuming force bias = energy bias.
brake.r_eff_rear  = (rotor.OD/2 + rotor.ID/2)/2; % m, mean friction radius, computed from measured rear rotor geometry
brake.r_eff_front = brake.r_eff_rear; % m, EDIT once front rotor OD/ID is known --
                                       % placeholder equal to rear radius, so the
                                       % torque correction has no effect until a
                                       % real front rotor geometry is entered

torque_weight_front = brake.bias_front * brake.r_eff_front;
torque_weight_rear  = brake.bias_rear  * brake.r_eff_rear;
torque_weight_total = torque_weight_front + torque_weight_rear;

brake.energy_frac_front_axle = torque_weight_front / torque_weight_total;
brake.energy_frac_rear_axle  = torque_weight_rear  / torque_weight_total;

brake.frac_per_front_rotor = brake.energy_frac_front_axle / brake.n_front_rotors;
brake.frac_per_rear_rotor  = brake.energy_frac_rear_axle  / brake.n_rear_rotors;

fprintf('=== SECTION 5: Brake Distribution ===\n');
fprintf('Front: %.0f%% bias / %d rotors -> %.1f%% of vehicle KE per front rotor\n', ...
    brake.bias_front*100, brake.n_front_rotors, brake.frac_per_front_rotor*100);
fprintf('Rear:  %.0f%% bias / %d rotor(s) -> %.1f%% of vehicle KE per rear rotor\n\n', ...
    brake.bias_rear*100, brake.n_rear_rotors, brake.frac_per_rear_rotor*100);

% This script models the REAR rotor, since it sees the larger
% per-rotor energy share (single rotor, higher bias) -- the worst
% case of the two. Swap to brake.frac_per_front_rotor to check the
% front rotor instead.
brake.frac_this_rotor = brake.frac_per_rear_rotor; % EDIT to check front rotor instead

%% ========================================================================
%  STAGE 1: TOTAL BRAKING ENERGY (single event, from vehicle KE)
%  ========================================================================
Q_total_single = 0.5*veh.mass*(veh.v_top^2 - veh.v_bottom^2); % Joules, WHOLE VEHICLE

% Energy reaching THIS rotor (rear, worst case) before the rotor/pad split:
Q_rotor_event = brake.frac_this_rotor * Q_total_single; % Joules, one rotor's share

fprintf('=== STAGE 1: Single-Stop Braking Energy ===\n');
fprintf('Vehicle KE dissipated per brake event: %.1f J\n', Q_total_single);
fprintf('Energy reaching this rotor (rear, %.1f%% share): %.1f J\n\n', ...
    brake.frac_this_rotor*100, Q_rotor_event);

%% ========================================================================
%  STAGE 2: ENERGY SPLIT VIA THERMAL EFFUSIVITY
%  ========================================================================
%  gamma_d = (e_d*A_d) / (e_d*A_d + e_p*A_p)
%  e = sqrt(rho * k * cp)   [thermal effusivity, W*s^0.5/(m^2*K)]

e_rotor = sqrt(rotor.rho * rotor.k * rotor.cp);
e_pad   = sqrt(pad.rho   * pad.k   * pad.cp);

gamma_d = (e_rotor*rotor.A_contact) / ...
          (e_rotor*rotor.A_contact + e_pad*pad.A_contact);
gamma_p = 1 - gamma_d;

fprintf('=== STAGE 2: Effusivity-Based Energy Split ===\n');
fprintf('Rotor effusivity e_d = %.1f  |  Pad effusivity e_p = %.1f\n', e_rotor, e_pad);
fprintf('Fraction of heat into ROTOR: %.1f%%\n', gamma_d*100);
fprintf('Fraction of heat into PAD:   %.1f%%\n\n', gamma_p*100);

Q_rotor_single = gamma_d * Q_rotor_event; % J into THIS rotor's material (rear, worst case)

%% ------------------------------------------------------------------
%  Sensitivity: how does the split change with rotor temperature?
%  (per your "see sensitivity of energy split with temperature" box)
%  ------------------------------------------------------------------
T_sweep_C = linspace(25, 500, 50);
gamma_d_sweep = zeros(size(T_sweep_C));

for i = 1:length(T_sweep_C)
    % Simple linear derating of steel k and cp with temperature.
    % Replace with real material-datasheet curves if you have them.
    k_T  = rotor.k  * (1 - 0.00035*(T_sweep_C(i)-25));   % k drops w/ temp
    cp_T = rotor.cp * (1 + 0.00040*(T_sweep_C(i)-25));   % cp rises w/ temp
    e_rotor_T = sqrt(rotor.rho * k_T * cp_T);
    gamma_d_sweep(i) = (e_rotor_T*rotor.A_contact) / ...
        (e_rotor_T*rotor.A_contact + e_pad*pad.A_contact);
end

%% ========================================================================
%  STAGE 3: SINGLE-STOP ROTOR TEMPERATURE RISE (conservative/adiabatic)
%  ========================================================================
%  Per the DRB note: "assume all the other energy went into the brake
%  pad (conservative)" -- i.e. this is a floor-level check assuming the
%  rotor sees ONLY its effusivity-split share and no time to cool.

dT_single = Q_rotor_single / (rotor.mass * rotor.cp);   % K rise, one event

fprintf('=== STAGE 3: Single-Stop Rotor Temp Rise (adiabatic) ===\n');
fprintf('Energy into rotor this stop: %.1f J\n', Q_rotor_single);
fprintf('Rotor temperature rise:      %.2f K\n\n', dT_single);

%% ========================================================================
%  STAGE 4: ENDURANCE SIMULATION (repeated stops + convective cooling)
%  ========================================================================
%  Lumped capacitance model:  m*cp*dT/dt = Q_in(t) - h*A*(T - T_amb)
%  Worst case: rotor treated as STATIONARY (no forced convection from
%  wheel rotation) -- matches your DRB assumption.
%
%  Replace the synthetic brake-event generator below with your real
%  brake-pressure / wheel-speed data log when available (see
%  IMPORT_TELEMETRY_TEMPLATE at bottom of this file).

dt = 0.1; % s, simulation time step
t_vec = 0:dt:endur.duration_s;
N = length(t_vec);
T_rotor = zeros(1,N);
T_rotor(1) = env.T_ambient_C;  % C, simulation now tracks Celsius directly

% --- Build synthetic brake-event energy input timeline -----------------
n_events = round(endur.brake_events_hr * endur.duration_s/3600);
event_times = sort(rand(1,n_events)*endur.duration_s);
Q_in_profile = zeros(1,N);
for ev = 1:n_events
    idx_start = find(t_vec >= event_times(ev), 1, 'first');
    idx_end   = find(t_vec >= event_times(ev)+endur.brake_time_s, 1, 'first');
    if isempty(idx_end), idx_end = N; end
    if isempty(idx_start), continue; end
    n_steps = max(idx_end - idx_start, 1);
    % Deliver this event's rotor-side energy as a power pulse over brake_time_s
    Q_in_profile(idx_start:idx_end) = Q_in_profile(idx_start:idx_end) + ...
        Q_rotor_single / (n_steps*dt);   % W
end

% --- Convective cooling coefficient ---------------------------------
% NOTE: the DRB slides specify "assume rotor is stationary (worst
% case)" -- that assumption is correctly kept fully conservative for
% the single-event checks in Stage 3 (adiabatic, zero cooling) and
% Stage 5 (failure temp). For THIS endurance time-march, using zero
% forced convection for the full 4-hour race is unrealistic: the car
% is actually moving and getting airflow between brake zones, so a
% "parked in still air the whole race" cooling rate makes the sim
% diverge into non-physical territory. Using a light-forced-convection
% estimate instead (still well below full ram-air cooling at speed).
h_conv = 50; % W/(m^2*K), light forced convection while driving between
             % brake events -- EDIT if you have real data (e.g. from
             % your planned dyno airflow correlation)

% --- Time march -----------------------------------------------------
for i = 2:N
    Q_loss = h_conv * rotor.A_surface * (T_rotor(i-1) - env.T_ambient_C); % W
    dTdt = (Q_in_profile(i-1) - Q_loss) / (rotor.mass * rotor.cp);
    T_rotor(i) = T_rotor(i-1) + dTdt*dt;
end

T_rotor_C = T_rotor;  % already Celsius, no conversion needed
[T_peak_C, idx_peak] = max(T_rotor_C);

fprintf('=== STAGE 4: Endurance Simulation ===\n');
fprintf('Simulated %d brake events over %.1f hr\n', n_events, endur.duration_s/3600);
fprintf('Peak rotor temperature: %.1f C at t = %.0f s\n', T_peak_C, t_vec(idx_peak));
fprintf('Final rotor temperature: %.1f C\n\n', T_rotor_C(end));

%% ========================================================================
%  STAGE 5: THERMAL STRESS / EXPANSION -- FAILURE TEMPERATURE CHECK
%  ========================================================================
%  Real rotor cracking is driven by the RADIAL temperature gradient
%  (hot outer band vs cooler hub), not the bulk temp rise alone.
%  grad_frac = fraction of bulk deltaT assumed to exist as a radial
%  gradient across the rotor -- conservative default 0.6 (EDIT if you
%  have thermal-camera data from mock endurance to refine this).

grad_frac = 0.6; % EDIT ME once you have thermocouple/IR data

dT_bulk_from_ambient = T_peak_C - env.T_ambient_C;
dT_gradient = grad_frac * dT_bulk_from_ambient;

% Constrained thermal (biaxial) stress approximation for a disc with a
% local hot band relative to a cooler, more rigid hub region:
%   sigma_thermal = E*alpha*deltaT / (1-nu)
sigma_thermal = rotor.E * rotor.alpha * dT_gradient / (1 - rotor.nu); % Pa

FOS_thermal = rotor.yield / sigma_thermal;

% Solve for the bulk temperature rise that would drop FOS to 1.0 (failure)
dT_gradient_fail = rotor.yield*(1-rotor.nu) / (rotor.E*rotor.alpha);
dT_bulk_fail = dT_gradient_fail / grad_frac;
T_fail_C = env.T_ambient_C + dT_bulk_fail;

fprintf('=== STAGE 5: Thermal Stress Check ===\n');
fprintf('Peak bulk deltaT from ambient: %.1f K\n', dT_bulk_from_ambient);
fprintf('Assumed radial gradient (%.0f%% of bulk): %.1f K\n', grad_frac*100, dT_gradient);
fprintf('Resulting thermal stress: %.1f MPa\n', sigma_thermal/1e6);
fprintf('Thermal FOS at current peak temp: %.2f\n', FOS_thermal);
fprintf('Estimated rotor bulk temp at FAILURE (FOS=1): %.0f C\n\n', T_fail_C);

% --- Structural energy capacity --------------------------------------
% Context: this is not a heat-soak or endurance-duration limit -- it is
% the maximum SINGLE-EVENT energy the rotor can structurally absorb
% before the resulting thermal gradient stress reaches yield (FOS=1).
% It is the direct energy-domain equivalent of the failure temperature
% above, expressed as a hard energy ceiling instead of a temperature,
% so it can be compared directly against Stage 1's per-event input
% (Q_rotor_single) or against a torque/energy-based duty cycle.
Q_structural_capacity = rotor.mass * rotor.cp * dT_bulk_fail; % J

fprintf('Rotor structural energy capacity (single event, FOS=1): %.0f J\n', Q_structural_capacity);
fprintf('Current single-event energy input: %.0f J (%.1f%% of structural capacity)\n\n', ...
    Q_rotor_single, Q_rotor_single/Q_structural_capacity*100);

%% ========================================================================
%  STAGE 6: ROTOR SIZING -- MIN MASS + POCKET REMOVAL SWEEP
%  ========================================================================
%  Given a max allowable bulk temp rise (from a target FOS, e.g. 1.2 per
%  your spec sheet), back-solve the minimum rotor mass for thermal
%  capacity, then sweep how much mass you can pocket out before you
%  blow that budget.

target_FOS = 1.2; % EDIT ME -- pull from the requirements doc (spec 4.1 uses
                   % 1.2 for hydraulic pressure; reuse or set your own thermal FOS)

dT_gradient_allow = (rotor.yield/target_FOS) * (1-rotor.nu) / (rotor.E*rotor.alpha);
dT_bulk_allow = dT_gradient_allow / grad_frac;

% Minimum rotor mass so that the worst-case single-event dT stays under budget
m_min = Q_rotor_single / (rotor.cp * dT_bulk_allow);

fprintf('=== STAGE 6: Rotor Thermal Sizing ===\n');
fprintf('Allowable bulk deltaT for FOS=%.1f: %.1f K\n', target_FOS, dT_bulk_allow);
fprintf('Minimum rotor mass for thermal capacity: %.3f kg\n', m_min);
fprintf('Current CAD rotor mass: %.3f kg  (margin: %.1f%%)\n\n', ...
    rotor.mass, (rotor.mass-m_min)/m_min*100);

% --- Reference rotor cross-check --------------------------------------
% Compares the computed minimum thermal mass against a known/benchmark
% rotor (e.g. a prior-year rotor mass, or an off-the-shelf reference)
% as a sanity check on whether the sizing result is realistic relative
% to hardware that has already run in competition, rather than trusting
% the calculator output in isolation.
reference_rotor_mass = 0.45; % kg, EDIT -- set to a known/prior rotor mass for comparison

fprintf('=== Reference Rotor Cross-Check ===\n');
if m_min > reference_rotor_mass
    fprintf(['WARNING: computed minimum thermal mass (%.3f kg) exceeds the ' ...
        'reference rotor mass (%.3f kg) -- the reference rotor may be ' ...
        'thermally undersized for this duty cycle.\n\n'], m_min, reference_rotor_mass);
else
    fprintf('Reference rotor (%.3f kg) clears the minimum thermal mass (%.3f kg) by %.1f%% margin.\n\n', ...
        reference_rotor_mass, m_min, (reference_rotor_mass-m_min)/m_min*100);
end

% --- Pocket removal sweep: how much can you cut before violating m_min?
mass_removed_pct = 0:1:40;   % sweep 0-40% mass removed via pockets
mass_candidates = rotor.mass * (1 - mass_removed_pct/100);

% Gradient-mass coupling: removing material (pocketing) reduces the
% cross-section available to conduct heat away from the hot band toward
% the hub, so the fraction of bulk deltaT that shows up as a radial
% gradient is modeled as increasing as mass is removed, rather than
% held constant at the single-point grad_frac used above. Linear
% inverse-mass scaling, capped at 1.0 (fully localized gradient) --
% EDIT this model once IR/thermocouple data exists across multiple
% pocket configurations to replace the assumed scaling.
grad_frac_candidates = min(1, grad_frac .* (rotor.mass ./ mass_candidates));

dT_candidates = Q_rotor_single ./ (mass_candidates * rotor.cp);
dT_gradient_candidates = grad_frac_candidates .* dT_candidates;
sigma_candidates = rotor.E * rotor.alpha .* dT_gradient_candidates / (1 - rotor.nu);
FOS_candidates = rotor.yield ./ sigma_candidates;

max_removable_idx = find(mass_candidates >= m_min, 1, 'last');
fprintf('Max mass you can remove via pockets and stay above m_min: %.0f%%\n\n', ...
    mass_removed_pct(max_removable_idx));

%% ========================================================================
%  PLOTS
%  ========================================================================
figure('Name','Baja Brake Thermal Analysis','Position',[100 100 1100 750]);

subplot(2,2,1);
plot(T_sweep_C, gamma_d_sweep*100, 'LineWidth', 2);
xlabel('Rotor Temperature (C)'); ylabel('% Energy into Rotor');
title('Energy Split Sensitivity vs Temperature'); grid on;

subplot(2,2,2);
plot(t_vec/60, T_rotor_C, 'LineWidth', 1.2);
xlabel('Time (min)'); ylabel('Rotor Temp (C)');
title('Endurance Rotor Temperature (Stationary/Worst-Case Cooling)'); grid on;
yline(T_fail_C, 'r--', 'Est. Failure Temp');

subplot(2,2,3);
plot(mass_removed_pct, mass_candidates, 'LineWidth', 2); hold on;
yline(m_min, 'r--', 'Min Mass for Thermal Cap.');
xlabel('% Mass Removed (pockets)'); ylabel('Rotor Mass (kg)');
title('Pocket Removal vs Thermal Mass Limit'); grid on; legend('Candidate mass','Location','best');

subplot(2,2,4);
plot(mass_removed_pct, FOS_candidates, 'LineWidth', 2); hold on;
yline(target_FOS, 'r--', sprintf('Target FOS = %.1f', target_FOS));
xlabel('% Mass Removed (pockets)'); ylabel('Thermal FOS');
title('Thermal FOS vs Pocket Removal'); grid on;

sgtitle('Purdue Baja Racing -- Brake Rotor Thermal Analysis');

%% ========================================================================
%  FIGURE 2: MASS-GRADIENT BALANCE
%  ========================================================================
%  Shows the coupling added in Stage 6: as pocketing removes mass, the
%  radial temperature gradient grows, working against the FOS margin
%  from two directions at once (less thermal mass to absorb energy, AND
%  a steeper gradient per degree of bulk temperature rise).
figure('Name','Mass-Gradient Balance','Position',[150 150 900 700]);

subplot(2,1,1);
yyaxis left
plot(mass_removed_pct, dT_candidates, 'LineWidth', 2);
ylabel('Bulk deltaT (K)');
yyaxis right
plot(mass_removed_pct, grad_frac_candidates*100, 'LineWidth', 2);
ylabel('Assumed Radial Gradient (% of bulk)');
xlabel('% Mass Removed (pockets)');
title('Bulk Temperature Rise vs Radial Gradient Assumption'); grid on;

subplot(2,1,2);
plot(mass_removed_pct, mass_candidates, 'LineWidth', 2); hold on;
yline(m_min, 'r--', 'Min Mass for Thermal Cap.');
yline(reference_rotor_mass, 'b--', 'Reference Rotor Mass');
xlabel('% Mass Removed (pockets)'); ylabel('Rotor Mass (kg)');
title('Candidate Mass vs Thermal Minimum and Reference Rotor'); grid on;
legend('Candidate mass','Location','best');

%% ========================================================================
%  IMPORT_TELEMETRY_TEMPLATE
%  ------------------------------------------------------------------------
%  When we have real dyno / brake-pressure data logs, replace the
%  synthetic event generator in STAGE 4 with something like this:
%
%   data = readtable('endurance_log.csv');   % columns: time_s, speed_mps,
%                                             % brake_pressure_psi, wheel_rpm
%   t_vec = data.time_s;
%   dt = mean(diff(t_vec));
%
%   % Convert brake line pressure -> clamp force -> torque -> power:
%   A_piston   = pi/4*(piston_dia_m)^2 * n_pistons;   % m^2, EDIT
%   F_clamp    = data.brake_pressure_psi*6894.76 .* A_piston; % N
%   Torque_br  = mu.pad * F_clamp * effective_pad_radius_m;   % N*m, EDIT
%   omega      = data.wheel_rpm * 2*pi/60;                    % rad/s
%   Power_in   = Torque_br .* omega;                          % W (total, both pads)
%   Q_in_profile = gamma_d * Power_in;                        % W into rotor
%
%  Then feed Q_in_profile directly into the time-march loop in Stage 4.
% ========================================================================