create_model_fnum.asv

function [sys, mats] = create_model_fnum(side_length, Fo, Co, By)
    
% This function creates the dynamical model
% It's much faster than the livescript version, which used symbolic
% variables. This ver

% Fo - Fourier number
% Co - Temperature-dependent heat loss to environment
% By - determines the effect of the inputs 'u', unit [K * W^-1]

    disp('Creating model')
    disp(cd)
    
    filename = "models" + filesep + "heatPlateN_" + string(side_length)+...
        "_Fo" + replace(string(Fo),'.','-') + "_Co" + ...
        replace(string(Co),'.','-') + "_By" + replace(string(By),'.','-') + ".mat";
    if isfile(filename)
        loaded_file = load(filename);
        sys = loaded_file.sys;
        mats = loaded_file.mats;
        disp('Model already exists, loading model from ' + filename + '.')
        return
    end
    
    % index (0,0) is assumed to be bottom left, x moves right, y moves up
    % tested only for side length 3, 99% sure it'll break for any other size
    
    n_elements = side_length ^ 2;

%     diag_main = spdiag(ones(n_elements,1)*(-4)) + spdiag(ones(n_elements,1) * Fo^(-1));
    A = spdiag(n_elements,0)*(-4) + spdiag(n_elements,0)/Fo;

    %Assuming all points are interior points for now
%     diag_sub1 = spdiag(ones(n_elements-1,1), -1);
%     diag_super1 = spdiag(ones(n_elements-1,1), 1);
%     diag_sub2 = spdiag(ones(n_elements-side_length,1), -side_length);
%     diag_super2 = spdiag(ones(n_elements-side_length,1), side_length);
    diag_sub1 = spdiag(n_elements, -1);
    diag_super1 = spdiag(n_elements, 1);
    diag_sub2 = spdiag(n_elements, -side_length);
    diag_super2 = spdiag(n_elements, side_length);
    A = A + diag_sub1 + diag_super1 + diag_sub2 + diag_super2;

    %Indices of edge nodes (excluding corner nodes) and corner nodes
    % Automatically finding edges for NxN plate 
    %   sides must be equal, could probably be generalized
    edge_l = 1:side_length:n_elements;
    edge_r = side_length:side_length:n_elements;
    edge_u = (n_elements-side_length+1):1:n_elements;
    edge_d = 1:1:side_length;

    A = remove_connections(A, edge_l, -1);
    A = remove_connections(A, edge_r, 1);
    A = remove_connections(A, edge_u, side_length-1);
    A = remove_connections(A, edge_d, -(side_length-1));
    A = A * Fo;
    
    B = spdiag(n_elements, 0) * By;
    
    C = spdiag(n_elements, 0) * Co;

    sys = @(Tk, u) A*Tk - C*Tk + B*u;
    
    mats.A = A-C;
    mats.B = B;
    
    save(filename, "sys", "mats");
    
    disp('Model created')
end



function A_new = remove_connections(A, idxs, index_dist)
     % A -> heat equation dynamics matrix, assumed as A*Fo^-1
     % index_dist -> row-wise index distance from diagonal
     [h, w] = size(A);

     for idx = idxs
         A(idx, idx) = A(idx, idx) + 1;  % add 1 to diagonal
         col_idx_to_remove = idx + index_dist;
         if (col_idx_to_remove<1) || (col_idx_to_remove>w)
             continue
         else
             A(idx, col_idx_to_remove) = 0;           
         end
     end
     A_new = sparse(A);
end