%% test_heatup_rate.m — measure the actual heatup-rate trajectory from ctrl_heatup
clear; addpath('controllers'); addpath('../reachability');
plant = pke_params();
pred  = load_predicates(plant);
x0 = [1e-3; (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * 1e-3; ...
      pred.constants.T_standby; pred.constants.T_standby; pred.constants.T_standby];
ref = struct('T_start', pred.constants.T_standby, ...
             'T_target', plant.T_c0, ...
             'ramp_rate', 28/3600);
[t, X, U] = pke_solver(plant, @(t) 0, @ctrl_heatup, ref, [0 5400], x0);

a_f = 0.4587; a_c = -0.9587; a_cold = 0.5000;
dTcdt = a_f*X(:,8) + a_c*X(:,9) + a_cold*X(:,10);     % C/s

fprintf('\n=== Heatup rate statistics ===\n');
fprintf('  max dT_c/dt:  %.4f C/s  = %6.1f C/hr\n', max(dTcdt), max(dTcdt)*3600);
fprintf('  min dT_c/dt:  %.4f C/s  = %6.1f C/hr\n', min(dTcdt), min(dTcdt)*3600);
fprintf('  rate limit:   +/- 0.01389 C/s = +/- 50.0 C/hr\n');
fprintf('  violates upper? %s\n', string(any(dTcdt > 0.01389)));
fprintf('  violates lower? %s\n', string(any(dTcdt < -0.01389)));