% Programa para que procesa los modelos de simulación
% Versión 01. San Salvador de Jujuy, 12/01/2025. Enrique Tarifa.
disp('Resolvedor v01, 2025');
disp('');
disp('Resolviendo el modelo...');

function v = vector(leyenda)
  % Devuelve el vector columna correspondiente a la variable leyenda.
  global LX LY tpts X Y
  indicex = find(strcmp(LX, leyenda)); % Índice del elemento
  if length(indicex) == 1
    v = X(:,indicex);
  else
    indicey = find(strcmp(LY, leyenda)); % Índice del elemento
    if length(indicey) == 1
      v = Y(:,indicey);
    else
      disp(['Error: Variable "' leyenda '" no encontrada.']);
	  error('Código de error: %d - Descripción del error', 1); % Detener con un mensaje
    endif
  endif
endfunction % vector
%---------------------------------------

function graficar(LV, titulo, rotulox, rotuloy, limitesy)
  % Crea una figura
  % LV: Arreglo de celdas fila que contiene los textos para las leyendas de las variables a graficar.
  % titulo: Título de la figura.
  % rotulox: Rótulo para la abscisa.
  % rotuloy: Rótulo para la ordenada.
  % limitesy: Vector fila con el límite inferior y el superior para la ordenada. Es opcional.
  global tpts

  colores = ['r' 'g' 'b' 'c' 'm' 'y' 'k'];

  figure;

  % Variables
  hold on; % Mantiene la figura para superponer la siguiente gráfica.
  for i = 1:length(LV)
    plot(tpts, vector(LV{i}), colores(mod(i-1,length(LV)) + 1), 'LineWidth', 2); % Línea con espesor 2
  endfor

  % Título del gráfico
  title(titulo);

  % Configurar los ejes
  xlabel(rotulox); % Título del eje x
  ylabel(rotuloy); % Título del eje y

  % Verificar y asignar valores predeterminados
  if nargin == 5
    ylim(limitesy); % Rango del eje y
  endif

  % Mostrar la cuadrícula
  grid on;

  % Añadir la leyenda
  legend(LV, 'Location', 'northeast'); % Leyenda en la esquina superior derecha

endfunction % graficar
%---------------------------------------

function exportar(archivo, LV)
  % Exporta todos los resultados o las variables seleccionadas
  % archivo: Nombre del archivo al que se exportará. Debe incluir la extensión.
  % LV: Lista de var  iables a exportar. Es opcional. Si no se incluye, se exportan todas.

  global LX LY tpts X Y
  
  sep = ';'; % Separador de columnas
  
  % Preparación de las variables a exportar
  if nargin == 1 % no se envió LV, se exportan todas las variables
    LV = [LX LY];
  endif

  A = tpts;
  for i = 1:length(LV)
    A = [A vector(LV{i})];
  endfor

  % Abrir el archivo para escritura
  fid = fopen(archivo, 'w');

  % Escribir la primera fila (encabezados)
  fprintf(fid, '%s%s', 'tpts', sep) % Escribe el rótulo del tiempo
  fprintf(fid, '%s\n', strjoin(LV, sep)) % Escribir encabezados separados por sep

  % Escribir los datos de la matriz X
  for i = 1:size(A, 1) % Iterar sobre las filas de A
      fprintf(fid, '%s\n', strjoin(arrayfun(@num2str, A(i, :), 'UniformOutput', false), sep)); % Escribir elementos de la fila separados por sep
  endfor

  % Cerrar el archivo
  fclose(fid);

  disp(['Archivo exportado como "' archivo '" en el directorio de trabajo.']);
end % exportar
%---------------------------------------

function [tpts X Y] = simulacion(tfin,dt,Xini)
  % Realiza la simulación.

  % Resolución
  nts = ceil(tfin/dt + 1); % redondea por exceso
  tpts = linspace(0, tfin, nts)';
  [tpts X] = ode15s(@ODEs,tpts,Xini');

  % Cálculo de las variables dependientes
  for i = 1:size(tpts,1)
    Y(i,:) = AEs(tpts(i),X(i,:)');
  endfor

endfunction % simulacion

%=============== Simulación =================
global LX LY tpts X Y

% Inicialización
[tfin dt Xini LX LY] = inicializacion;

% Resolución
[tpts X Y] = simulacion(tfin,dt,Xini);

% Análisis
analizar(LX,LY,tpts,X,Y);

disp('');
disp('Simulación finalizada.');