% Step 1: Load the simulation and plant data
%load('transfer_fcn.mat');       % Replace 'transfer_fcn.mat' with your simulation data file
%load('rc_circuit.mat');         % Replace 'rc_circuit.mat' with your orginal plant data file
%load('adj_rc_circuit.mat');     % Replace 'adj_rc_circuit.mat' with your adjusted plant data file

% Step 2: Define the lookup table resolution
rc_circuit_size = size(rc_circuit.signals.values);
num_points = rc_circuit_size(1);  % Adjust the number of points as needed

% Step 3: Interpolate the data to create the lookup table
lookup_time = linspace(min(rc_circuit.time), max(rc_circuit.time), num_points);
lookup_simulation = interp1(transfer_fcn.time, transfer_fcn.signals.values, lookup_time, 'linear');

% Step 4: Validate the lookup table
figure;
plot(rc_circuit.time, rc_circuit.signals.values, 'b', lookup_time, lookup_simulation, 'r--');
xlabel('Time');
ylabel('Response');
legend('Actual RC Circuit Output', 'Discrete Transfer Function Simulation');
title('Comparison between Actual RC Circuit output and Transfer Function Simulation');

% Step 5: Calculate the error between RC circuit data and the matched simulation
error = rc_circuit.signals.values - transfer_fcn.signals.values;

% Step 6: Find peak values and their indices in both the RC circuit data and the simulation
[~, plant_peak_indices] = findpeaks(rc_circuit.signals.values);
[~, sim_peak_indices] = findpeaks(lookup_simulation);

% Step 7: Calculate the error at the peak values
peak_errors = error(intersect(plant_peak_indices, sim_peak_indices));

% Step 8: Calculate the mean absolute error at peak values
mean_absolute_error_at_peak = mean(abs(peak_errors));

% Step 9: Display the size of the peak_errors variable and the mean absolute error at peak values
disp(['Size of peak_errors: ', num2str(size(peak_errors))]);

% Calculate the normalized mean absolute error at peak values relative to the maximum output of the RC circuit
normalized_mean_absolute_error = mean_absolute_error_at_peak / max(rc_circuit.signals.values);
disp(['Normalized mean absolute error at peak values: ', num2str(normalized_mean_absolute_error)]);

% Step 10: Calculate the mean peak error as a percentage
mean_peak_error_percentage = normalized_mean_absolute_error * 100;

% Step 12: Display the mean peak error as a percentage and calibration percentage
disp(['Mean peak error as a percentage: ', num2str(mean_peak_error_percentage), '%']);
