clc; clear all; close all;

pkg load control signal

%% Configuración de Octave

set(0, "defaultAxesColorOrder", [
        1.0000 0.0000 0.0000;  % red
        0.0000 0.0000 1.0000;  % blue
        1.0000 0.4980 0.0549;  % orange
        0.1725 0.6275 0.1725;  % green
        0.5804 0.4039 0.7412;  % purple
        0.5490 0.3373 0.2941;  % brown
        0.8902 0.4667 0.7608;  % pink
        0.4980 0.4980 0.4980;  % gray
]);

set(0, "DefaultFigureColor", [1 1 1]);
set(0, "DefaultAxesColor", [1 1 1]);
set(0, "DefaultAxesXColor", [0 0 0]);
set(0, "DefaultAxesYColor", [0 0 0]);
set(0, "DefaultTextColor", [0 0 0]);
set(0, "DefaultAxesGridColor", [0 0 0]);
set(0, 'DefaultAxesGridAlpha', 1.00);
set(0, "DefaultAxesMinorGridColor", [0 0 0]);
set(0, 'DefaultAxesMinorGridAlpha', 0.20);
set(0, "DefaultLineLinewidth", 3.00);
set(0, "DefaultAxesFontSize", 16);
set(0, "DefaultTextFontSize", 16);
set(0, "DefaultAxesLineWidth", 1.00);
set(0, "DefaultAxesXGrid", "on");
set(0, "DefaultAxesYGrid", "on");
set(0, "DefaultAxesZGrid", "on");
set(0, "DefaultAxesXMinorGrid", "on");
set(0, "DefaultAxesYMinorGrid", "on");
set(0, "DefaultAxesXMinorTick", "on");
set(0, "DefaultAxesYMinorTick", "on");
set(0, 'DefaultAxesGridAlpha', 0.50);
set(0, "defaultAxesFontName", "Nimbus Sans");
set(0, "defaultaxesminorgridlinestyle", ":");

set(0, "defaultfigurepaperunits", "inches");
set(0, "defaultfigurepapersize", [10 4]);
set(0, "defaultfigurepaperposition", [0 0 10 4]);

% No abrir ventanas para ninguna figura
set(0, "defaultfigurevisible", "off");

%% Funcion para leer archivss

function [steps_data, step_names] = ltspice_step_file(filename)
    % IMPORTAR_LTSPICE_STEPS Lee archivos de LTspice exportados 
    % con directivas .STEP
    %
    % Uso: [data, nombres] = ltspice_step_file("archivo.csv")
    
    % Abrir el archivo en modo lectura
    fid = fopen(filename, 'r');

    if fid == -1
        error('No se pudo abrir el archivo: %s. Revisa la ruta.', filename);
    end

    steps_data = {};  % Celda donde guardaremos los datos de cada paso
    step_names = {};  % Nombres de los pasos (ej: "F=250")
    current_step = 0; % Contador de pasos detectados
    temp_data = [];

    % Leer línea por línea
    while !feof(fid)
        linea = fgetl(fid);
        
        % Si la línea está vacía o es la cabecera original, se ignora
        if isempty(linea) || startsWith(linea, 'time')
            continue;
        end
        
        % Detectar si la línea contiene información de un nuevo paso (Step)
        if !isempty(strfind(linea, 'Step Information:'))
            % Si ya veníamos acumulando datos de un paso anterior, los guardamos
            if current_step > 0 && !isempty(temp_data)
                steps_data{current_step} = temp_data;
                temp_data = [];
            end
            
            current_step = current_step + 1;
            
            % Extraer el identificador del parámetro usando expresiones regulares
            tokens = regexp(linea, 'Step Information:\s*([^\(]+)', 'tokens');
            if !isempty(tokens)
                step_names{current_step} = strtrim(tokens{1}{1});
            else
                step_names{current_step} = sprintf('Step %d', current_step);
            end
            
        else
            % Si no es un header de Step, es una línea con números (Tiempo y Voltaje)
            valores = sscanf(linea, '%f %f');
            if length(valores) == 2
                temp_data = [temp_data; valores'];
            end
        end
    end

    % No olvidar guardar el último paso leído al llegar al final del archivo
    if current_step > 0 && !isempty(temp_data)
        steps_data{current_step} = temp_data;
    end

    fclose(fid);
    printf('Se importaron con éxito %d pasos desde: %s\n', ...
        length(steps_data), filename);
end

% f = logspace(0, 5, 1000);
% w = 2*pi*f;

%% Transferencia ideal

H1 = tf([1 250], [1 500])
H2 = tf([4000^2], [1 4000^2/9000 4000^2])
H = 12*H1*H2

%% Transferencia normalizada

H1_norm = tf([1 252.52], [1 505.05])
H2_norm = tf([4014.80^2], [1 4014.80/2.264 4014.80^2])
H_norm = 12*H1_norm*H2_norm

%% Diagrama de Bode

% Datos simulaciones LTSpice

filename = "data/sim_bode.csv";
fid = fopen(filename, "r");
fgetl(fid); % Saltear encabezado
data = textscan(fid, "%f %s", "Delimiter", "\t");
fclose(fid);
freq = data{1};
w = 2*pi*freq;
response = data{2};
N = numel(freq);
mag = zeros(N, 1);
phase = zeros(N, 1);

for k = 1:N
    tokens = regexp(response{k}, ...
        '\(([-+0-9.eE]+)dB,([-+0-9.eE]+)°\)', ...
        'tokens');

    mag_circ(k)   = str2double(tokens{1}{1});
    pha_circ(k) = str2double(tokens{1}{2});
end

% Transferencia ideal
[mag_tf, pha_tf] = bode(H, w);

% La funcion `bode` devuelve el eje `y` en escala lineal
mag_tf = squeeze(mag_tf);
pha_tf = squeeze(pha_tf);
mag_tf = 20*log10(mag_tf);

% Transferencia normalizada
[mag_tf_norm, pha_tf_norm] = bode(H_norm, w);

% La funcion `bode` devuelve el eje `y` en escala lineal
mag_tf_norm = squeeze(mag_tf_norm);
pha_tf_norm = squeeze(pha_tf_norm);
mag_tf_norm = 20*log10(mag_tf_norm);

% Magnitud

ax1 = subplot(2, 1, 1);
set(gcf, "paperunits", "inches");
set(gcf, "papersize", [10 5]);
set(gcf, "paperposition", [0 0 10 5]);

grid on;
hold on;

y1 = semilogx(w, mag_tf, "LineWidth", 3.50);
y2 = semilogx(w, mag_tf_norm, "LineWidth", 2.75);
y3 = semilogx(w, mag_circ, "LineWidth", 2.00);

ylim([-35 35])
yticks([-20 0 20]);
xlim([1 1e5])
ylabel("MAGNITUD [dB]");
set(ax1, "XTickLabel", []); % Ocultar x-tick del gráfico superior

% Fase

ax2 = subplot(2, 1, 2)
grid on;
hold on;

semilogx(w, pha_tf, "LineWidth", 3.00);
semilogx(w, pha_tf_norm, "LineWidth", 2.75);
semilogx(w, pha_circ, "LineWidth", 2.00);

ylim([-225 45])
yticks([-180 -90 0]);
xlim([1 1e5])
xlabel("FRECUENCIA  [rad/s]");
ylabel("FASE  [grados]");

% Ejes invisibles para `legend`
ax0 = axes("Position", [0 0 1 1], ...
           "Visible", "off", ...
           "Units", "normalized");

lgd = legend(ax0, [y1 y2 y3], ...
             {"H", "H normalizada", "Circuito"});

set(lgd, "position", [0.16 0.54 0.10 0.05]);
set(lgd, "numcolumns", 1);

% Reducir espacio entre sub-figuras
set(ax1, "Position", [0.13 0.58 0.80 0.38]);
set(ax2, "Position", [0.13 0.13 0.80 0.38]);

print("plot/raw/sim_bode.png", "-dpng", "-r500");

%% Respuesta al impulso (solo para las transferencias)

t = linspace(0, 0.010, 10000);
[h1, t1] = impulse(H, t);
[h2, t2] = impulse(H_norm, t);

figure();
grid on;
hold on;

y1 = plot(t1*1e3, h1, 'linewidth', 3.00, 'displayname', 'G1');
y2 = plot(t2*1e3, h2, 'linewidth', 2.25, 'displayname', 'G2');

xlabel("TIEMPO  t  [ms]");
ylabel("AMPLITUD");
xlim([0 10])
xticks(0:1:16);
ylim([-35000 35000])

legend("location", "southeast");

legend(gca, [y1 y2], {"H", "H normalizada"});

print('plot/raw/impulse.png', '-dpng', '-r500');

%% Respuesta al escalon

% Datos simulacion LTSpice
filename = "data/sim_step.csv";
fid = fopen(filename, 'r');
fgetl(fid); % Saltear encabezado
data = fscanf(fid, '%f %f', [2, Inf]); % Leer dos columnas
fclose(fid);
data = data'; 

t0  = data(:, 1);
u0 = data(:, 2);

% Transferencias
[u1, t1] = step(H, t0);
[u2, t2] = step(H_norm, t0);

% Plot
figure();
grid on;
hold on;

y1 = plot(t1*1e3, u1, 'linewidth', 3.50);
y2 = plot(t2*1e3, u2, 'linewidth', 2.75);
y0 = plot(t0*1e3, u0, 'linewidth', 2.00);

xlabel("TIEMPO  t  [ms]");
ylabel("AMPLITUD");
legend('location', 'northeast');

legend(gca, [y0 y1 y2], ...
    {"Circuito", 'H_{1}', 'H_{1} normalizada'});

xlim([0 10])
xticks(0:1:16);
ylim([0 20])

print('plot/raw/step.png', '-dpng', '-r500');

%% Respuesta a señales senoidales

set(0, "defaultAxesColorOrder", [
        0.1725 0.6275 0.1725;  % green
        0.0000 0.0000 1.0000;  % blue
        1.0000 0.4980 0.0549;  % orange
        1.0000 0.0000 0.0000;  % red
]);

% Abrir el archivo en modo lectura
file_path = "data/sim_sin.csv";
[data, steps_name] = ltspice_step_file(file_path);

freqs = [250, 636.6, 10e3]; % Frequencias [Hz]

for i = 1:length(data)
    f = freqs(i);

    figure();
    grid on;
    hold on;

    xlabel("TIEMPO  t  [ms]");
    ylabel("AMPLITUD");

    switch i
        case 1
            sin_signal_lgd = sprintf("Señal senoidal 250 Hz");
            save_filename = sprintf("plot/raw/sim_sine_resp_250.png");
            % t = 0:1e-6:500e-3;
            xlim([0 50])
            xticks(0:5:50);
            ylim([-4 4])
            yticks([-3 -1.5 0 1.5 3]);
        case 2
            sin_signal_lgd = sprintf("Señal senoidal 637 Hz");
            save_filename = sprintf("plot/raw/sim_sine_resp_637.png", f);
            % t = 0:1e-6:10e-3;
            xlim([0 10])
            xticks(0:1:10);
            
            ylim([-4 4])
            yticks([-3 -1.5 0 1.5 3]);
        case 3
            sin_signal_lgd = sprintf("Señal senoidal 10 kHz");
            save_filename = sprintf("plot/raw/sim_sine_resp_10k.png");
            % t = 0:1e-9:100e-6;
            xlim([0 10])
            % xticks(0:0.1:0.5);
            xticks(0:1:10);
            ylim([-2.5 2.5])
            yticks([-2 -1 0 1 2]);
        otherwise
    end

    step_data = data{i};
    t = step_data(:, 1);
    y_circ = step_data(:, 2);

    u = .1*sin(2*pi*f*t);
    y      = lsim(H, u, t);
    y_norm = lsim(H_norm, u, t);

    y0 = plot(t*1e3, u, 'linewidth', 3.00);
    y1 = plot(t*1e3, y, 'linewidth', 3.75);
    y2 = plot(t*1e3, y_norm, 'linewidth', 3.00);
    y3 = plot(t*1e3, y_circ, 'linewidth', 2.25);

    legend("location", "southeast");

    legend(gca, [y0 y1 y2 y3], ...
        {sin_signal_lgd,...
            "H", "H normalizada", "Circuito"});

    print(save_filename, "-dpng", "-r500");
end

%% Respuesta a señales cuadradas

file_path = "data/sim_squ.csv";
[data, steps_name] = ltspice_step_file(file_path);

freqs = [10, 250, 2.5e3]; % Frequencias [Hz]

for i = 1:length(data)
    f = freqs(i);

    figure();
    grid on;
    hold on;

    xlabel("TIEMPO  t  [ms]");
    ylabel("AMPLITUD");
    % xticks(0:1:16);

    switch i
        case 1
            sig_lgd = sprintf("Señal cuadrada 10 Hz");
            save_file_name = sprintf("plot/raw/sim_square_resp_10.png");
            % t = 0:1e-6:500e-3;
            xlim([0 300])
            % xticks(0:1:10);
            ylim([-3 3])
            % yticks([-12 -6 0 6 12]);
        case 2
            sig_lgd = sprintf("Señal cuadrada 250 Hz");
            save_file_name = sprintf("plot/raw/sim_square_resp_250.png");
            % t = 0:1e-6:10e-3;
            xlim([0 25])
            xticks(0:5:25);
            ylim([-6 6])
        case 3
            sig_lgd = sprintf("Señal cuadrada 2.5 kHz");
            save_file_name = sprintf("plot/raw/sim_square_resp_2k5.png");
            % t = 0:1e-9:10e-6;
            xlim([0 5])
            xticks(0:1:5);
            ylim([-5 5])
            yticks([-4 -2 0 2 4]);
        otherwise
    end

    step_data = data{i};
    t = step_data(:, 1);
    y_circ = step_data(:, 2);

    u = .1*square(2*pi*f*t);
    y      = lsim(H, u, t);
    y_norm = lsim(H_norm, u, t);

    y0 = plot(t*1e3, u, 'linewidth', 3.00);
    y1 = plot(t*1e3, y, 'linewidth', 3.50);
    y2 = plot(t*1e3, y_norm, 'linewidth', 2.75);
    y3 = plot(t*1e3, y_circ, 'linewidth', 2.00);

    % legend("location", "southeast");

    legend(gca, [y0 y1 y2 y3], ...
        {sig_lgd, "H", "H normalizada", "Circuito"});

    print(save_file_name, "-dpng", "-r500");
end

% legend('location', 'southeast');
% xlim([0 10])
% ylim([0 20])

% print('tf_step.png', '-dpng', '-r500');