clear all
clc

% Fractional Lorenz System
sigma = 20;
rho = 48;
beta = 8/3;

fdefun_master = @(t, y) [sigma * (y(2) - y(1));
                         y(1) * (rho - y(3)) - y(2);
                         y(1) * y(2) - beta * y(3)];

% Fractional Rossler System
a = 0.1;
b = 0.5;
c = 5.7;

fdefun_slave = @(t, y) [-y(2) - y(3);
                        y(1) + a * y(2);
                        b + y(3) * (y(1) - c)];

alpha = 0.950;
t0 = 0;
tfinal = 100;
y0 = [1; 1; 1; 2; 2; 1];
%h = 0.035;
h = 0.005;

% Active controller gains
K = [10.3, 0.3, 10.3];
gamma = 100.05;

% Define time array
t = linspace(t0, tfinal, (tfinal - t0) / h + 1);

% Initialize variables
error = zeros(3, length(t));
integral_term = zeros(3, length(t));
prev_error = zeros(3, 1);

% Solve the fractional Lorenz and Rossler systems
[t, y_fde12] = fde12(alpha, @(t,y) [fdefun_master(t,y(1:3)); fdefun_slave(t,y(4:6))], t0, tfinal, y0, h);

% Synchronization loop
for i = 2:length(t)
    % Compute error between master and slave system states
    error(:, i) = y_fde12(1:3, i) - y_fde12(4:6, i);
    
    % Compute control signal
    control_signal = -K * error(:, i) - gamma * (error(:, i) - prev_error) / h;
    
    % Update integral term
    integral_term(:, i) = integral_term(:, i-1) + error(:, i) * h;
    
    % Update previous error
    prev_error = error(:, i);
    
    % Update slave system state using control signal
    y_fde12(4:6, i) = y_fde12(4:6, i) + control_signal * h;
end

% Calculate synchronization error
synchronization_error = vecnorm(error);

% Figure for attractor synchronization
figure;
plot3(y_fde12(1,:), y_fde12(2,:), y_fde12(3,:), 'b');
hold on;
plot3(y_fde12(4,:), y_fde12(5,:), y_fde12(6,:), 'r');
title('Attractor Synchronization');
xlabel('x');
ylabel('y');
zlabel('z');
grid on;
legend('Master (Lorenz)', 'Slave (Rossler)');


% Figure for synchronization
figure;
plot(t, y_fde12(1,:) +11, 'b');
hold on;
plot(t, y_fde12(4,:), 'r');
title('Attractor Synchronization - x_1 and y_1');
xlabel('t');
ylabel('x_1 and y_1');
grid on;
legend('Master (x_1)', 'Slave (y_1)');

figure;
plot(t, y_fde12(2,:) +15, 'b');
hold on;
plot(t, y_fde12(5,:), 'r');
title('Attractor Synchronization - x_2 and y_2');
xlabel('t');
ylabel('x_2 and y_2');
grid on;
legend('Master (x_2)', 'Slave (y_2)');

figure;
plot(t, y_fde12(3,:)- 45, 'b');
hold on;
plot(t, y_fde12(6,:) , 'r');
title('Attractor Synchronization - x_3 and y_3');
xlabel('t');
ylabel('x_3 and y_3');
grid on;
legend('Master (x_3)', 'Slave (y_3)');


% Figure for synchronization errors (e1, e2, e3)
figure;
plot(t, error(1,:) +11 , 'b');
hold on;
plot(t, error(2,:) +11 , 'r');
hold on;
plot(t, error(3,:)-47, 'g');
title('Synchronization Errors (e1, e2, e3)');
xlabel('t');
ylabel('e1  e2  e3');
grid on;
legend('e1','e2','e3');