function F = formapolinomioinverso(nombrefuncion, K, grado, X)
%formapolinomioinverso
%   F = formapolinomioinverso(funcion)
%
%   Construye una aproximación en forma de recetas en tiempo discreto de
%   unas recetas en tiempo continuo, definidas en un script o con el
%   comando inline; un ejemplo de este ultimo caso seria  
%      F = inline('[-280+575/k,-12; -120,-8+20/k; -280+400/k, -12]')
%
%   F = formapolinomioinverso(funcion, K, grado)
%      funcion   nombre de la función con las recetas en tiempo continuo
%      K         factores sobre los que se calcula la aproximación
%      grado     número de pasos temporales para el tiempo discreto
%      F         aproximación con recetas en tiempo discreto
%
%   Si K es un escalar la aproximación se hace en el intervalo [K/2, K*3/2]

if nargin < 2, K = []; end
if nargin < 3, grado = []; end
if nargin < 4, X = []; end

if isempty(K), K = 1; end                                %si K no está definido es 1
if size(K(:),1) == 1, K = (K/2:K/5:K*3/2); end           %si K es un escalar la aproximación se efectua en [K/2, K*3/2] con 6 puntos
puntos = size(K, 2);                                     %determina el número de factores para los que se hace la aproximación
if isempty(grado), grado = fix((puntos+3)/2); end        %define el grado del polinomio inverso por defecto
if grado > puntos, grado = puntos; end                   %el grado debe ser menor o igual al numero de puntos analizados
                                                         %
if isempty(X)                                            %si el problema está escrito con una matriz de flujos actualizados
    F0 = feval(nombrefuncion, K(1));                         %calcula la magnitud de los flujos actualizados para K(1)
else                                                     %si el problema está escrito con funciones de producción actualizadas
    F0 = feval(nombrefuncion, K(1), X);                      %calcula la magnitud de los flujos actualizados para K(1) y X
end                                                      %
[filas, columnas] = size(F0);                            %establece el tamaño del problema
Fk = zeros(filas, columnas, puntos);                     %crea la variable para almacenar los resultados
Fk(:, :, 1) = F0;                                        %almacena la magnitud de los flujos actualizados para K(1)
for c1 = 2:puntos                                        %para el resto de K analizadas
   if isempty(X)                                             %si el problema está escrito con una matriz de flujos actualizados
       Fk(:,:,c1) = feval(nombrefuncion, K(c1));                         %calcula la magnitud de los flujos actualizados para K(1)
   else                                                      %si el problema está escrito con funciones de producción actualizadas
       Fk(:,:,c1) = feval(nombrefuncion, K(c1), X);                      %calcula la magnitud de los flujos actualizados para K(1) y X
   end                                                       %
end                                                      %
x = 1 ./ K;                                              %define x como 1 / K analizadas
y = zeros(1, puntos);                                    %define y
for c1 = 1:filas                                         %para cada fila
    for c2 = 1:columnas                                      %para cada columna
        y(:) = Fk(c1, c2, :);                                     %extrae y
        F(c1, c2, :) = ajusta(x, y, puntos, grado);               %calcula los coeficientes del polinomio inverso	
    end                                                      %
end                                                      %






function A = ajusta(x, y, fil, col)                      %función para el cálculo de los coeficientes de un polinomio inverso
M = zeros(fil, col);                                     %crea la matriz del sistema lineal
for c1 = 1:fil                                           %para cada fila
    for c2 = 1:col                                           %para cada columna
        M(c1, col-c2+1) = x(c1) .^ (col-c2);                     %establece la matriz
    end                                                      %
end                                                      %
A = (M \ y')';                                           %resuelve el sistema lineal para establecer los coeficientes del polinomio
