clear;
close all;
clc;

% Create the figure and UI components
fig = figure('Position', [100, 100, 800, 800], 'Name', 'Truss Analysis', 'NumberTitle', 'off');
    
uicontrol(fig, 'Style', 'text', 'Position', [120, 750, 150, 20], 'String', 'Input File Path:', 'HorizontalAlignment', 'left');
uicontrol(fig, 'Style', 'edit', 'Position', [200, 750, 250, 20], 'Tag', 'input_file_path');
    
uicontrol(fig, 'Style', 'text', 'Position', [480, 750, 150, 20], 'String', 'Visual Scale:', 'HorizontalAlignment', 'left');
uicontrol(fig, 'Style', 'edit', 'Position', [550, 750, 100, 20], 'Tag', 'scale');
    
uicontrol(fig, 'Style', 'pushbutton', 'Position', [675, 750, 100, 20], 'String', 'Analyze', 'Callback', {@analyze, fig});
    
uicontrol(fig, 'Style', 'text', 'Position', [20, 150, 150, 20], 'String', 'Nodal Displacements:', 'HorizontalAlignment', 'left');
uicontrol(fig, 'Style', 'listbox', 'Position', [20, 20, 250, 120], 'Tag', 'nodal_displacements');
    
uicontrol(fig, 'Style', 'text', 'Position', [280, 150, 150, 20], 'String', 'Reaction Forces:', 'HorizontalAlignment', 'left');
uicontrol(fig, 'Style', 'listbox', 'Position', [280, 20, 250, 120], 'Tag', 'reaction_forces');
    
uicontrol(fig, 'Style', 'text', 'Position', [540, 150, 150, 20], 'String', 'Element Stresses:', 'HorizontalAlignment', 'left');
uicontrol(fig, 'Style', 'listbox', 'Position', [540, 20, 250, 120], 'Tag', 'element_stresses');

% Add an axes for the plot
ax = axes(fig, 'Units', 'pixels', 'Position', [60, 220, 700, 480], 'Tag', 'plot_axes');

function analyze(~, ~, fig)
    input_file_path = get(findobj(fig, 'Tag', 'input_file_path'), 'String');
    scale = str2double(get(findobj(fig, 'Tag', 'scale'), 'String'));
        
    % Run the truss analysis
    [nodes, deformed_nodes, fdisplacement, totF, stress, el] = truss_analysis(input_file_path, scale);
        
    % Update the GUI components with the results
    set(findobj(fig, 'Tag', 'nodal_displacements'), 'String', num2str(fdisplacement));
    set(findobj(fig, 'Tag', 'reaction_forces'), 'String', num2str(totF));
    set(findobj(fig, 'Tag', 'element_stresses'), 'String', num2str(stress));
        
    % Show the truss deformation plot
    axes(findobj(fig, 'Tag', 'plot_axes'));
    cla;
    hold on;

    num_elements = size(el, 1);
    for i = 1:num_elements
        node1 = el(i, 1);
        node2 = el(i, 2);
        plot([nodes(node1, 1), nodes(node2, 1)], [nodes(node1, 2), nodes(node2, 2)], 'k-o', 'LineWidth', 1.5);
    end

    % Plot deformed truss
    for i = 1:num_elements
        node1 = el(i, 1);
        node2 = el(i, 2);
        plot([deformed_nodes(node1, 1), deformed_nodes(node2, 1)], [deformed_nodes(node1, 2), deformed_nodes(node2, 2)], 'r--o', 'LineWidth', 1.5);
    end

    % Determine the minimum and maximum coordinates
    min_x = min(min(nodes(:, 1)), min(deformed_nodes(:, 1)));
    max_x = max(max(nodes(:, 1)), max(deformed_nodes(:, 1)));
    min_y = min(min(nodes(:, 2)), min(deformed_nodes(:, 2)));
    max_y = max(max(nodes(:, 2)), max(deformed_nodes(:, 2)));

    % Set the axes limits with some padding
    padding = 0.1;
    xlim([min_x - padding, max_x + padding]);
    ylim([min_y - padding, max_y + padding]);

    % Label nodes
    for i = 1:n
        text(nodes(i, 1) + 0.05, nodes(i, 2) + 0.05, num2str(i), 'Color', 'k', 'FontSize', 8);
        text(deformed_nodes(i, 1) + 0.05, deformed_nodes(i, 2) + 0.05, num2str(i), 'Color', 'r', 'FontSize', 8);
    end

    % Label lines
    for i = 1:num_elements
        node1 = el(i, 1);
        node2 = el(i, 2);
        midpoint_x = (nodes(node1, 1) + nodes(node2, 1)) / 2;
        midpoint_y = (nodes(node1, 2) + nodes(node2, 2)) / 2;
        text(midpoint_x, midpoint_y, num2str(i), 'Color', 'b', 'FontSize', 8);
    end


    legend('Original', 'Deformed');
    xlabel('X (m)');
    ylabel('Y (m)');
    title('Truss Deformation');
end

function [nodes, deformed_nodes, fdisplacement, totF, stress, el] = truss_analysis(input_file_path, scale)
    % Read input file
    fileID = fopen(input_file_path, 'r');
    
    if fileID == -1
        error('Error: Cannot open input file.');
    end

    % Read node coordinates
    nodes = [];
    line = fgets(fileID);
    while ~feof(fileID)
        line = fgets(fileID);
        if strncmp(line, '- ', 2)
            break;
        end
        nodes = [nodes; sscanf(line, '%f %f')'];
    end
    n = size(nodes, 1);

    % Read elements
    el = [];
    while ~feof(fileID)
        line = fgets(fileID);
        if strncmp(line, '- ', 2)
            break;
        end
        el = [el; sscanf(line, '%d %d')'];
    end
    m = size(el, 1);

    % Read element areas
    A = [];
    while ~feof(fileID)
        line = fgets(fileID);
        if strncmp(line, '- ', 2)
            break;
        end
        A = [A; sscanf(line, '%e')'];
    end

    % Read Young's modulus for each element
    E = [];
    while ~feof(fileID)
        line = fgets(fileID);
        if strncmp(line, '- ', 2)
            break;
        end
        E = [E; sscanf(line, '%e')'];
    end

    % Read displacement boundary conditions
    displacement = [];
    while ~feof(fileID)
        line = fgets(fileID);
        if strncmp(line, '- ', 2)
            break;
        end
        displacement = [displacement; sscanf(line, '%d')'];
    end
    displacement = displacement';

    % Read external forces
    forces = [];
    while ~feof(fileID)
        line = fgets(fileID);
        if strncmp(line, '- ', 2)
            break;
        end
        forces = [forces; sscanf(line, '%e')'];
    end
    forces = forces';

    % Close the input file
    fclose(fileID);

    % Assemble the global stiffness matrix and elements
    GK = zeros(2 * n, 2 * n);
    elements = zeros(m, 4);

    for i = 1:m
        node1 = el(i, 1);
        node2 = el(i, 2);
        x1 = nodes(node1, 1);
        y1 = nodes(node1, 2);
        x2 = nodes(node2, 1);
        y2 = nodes(node2, 2);

        L = sqrt((x2 - x1)^2 + (y2 - y1)^2);
        c = (x2 - x1) / L;
        s = (y2 - y1) / L;

        k = E(i) * A(i) / L;
        KL = k * [c^2, c*s, -c^2, -c*s;
                  c*s, s^2, -c*s, -s^2;
                  -c^2, -c*s, c^2, c*s;
                  -c*s, -s^2, c*s, s^2];

    index = [2 * node1 - 1, 2 * node1, 2 * node2 - 1, 2 * node2];
    GK(index, index) = GK(index, index) + KL;

    elements(i, :) = [x1, y1, x2, y2];
    end

% Calculate the external forces
loads = forces';

% Apply the displacement boundary conditions
fixed_nodes = find(displacement == 0);
free_nodes = find(displacement == 1);

GK_full = GK;

GK(fixed_nodes, :) = [];
GK(:, fixed_nodes) = [];

% Modify the loads vector according to the free nodes
loads = loads(free_nodes);

% Solve the system of equations
fdisplacement = linsolve(GK,loads);

% Reassemble fdisplacement to match the original size
full_fdisplacement = zeros(n * 2, 1);
full_fdisplacement(free_nodes) = fdisplacement;
fdisplacement = full_fdisplacement;

% Calculate the reaction forces
totF = GK_full*fdisplacement;

% Calculate stresses
stress = zeros(m, 1);
for i = 1:m
    node1 = el(i, 1);
    node2 = el(i, 2);
    L = norm(nodes(node1, :) - nodes(node2, :));
    cos_theta = (nodes(node2, 1) - nodes(node1, 1)) / L;
    sin_theta = (nodes(node2, 2) - nodes(node1, 2)) / L;

    local_node1_disp_x = fdisplacement(2 * node1 - 1);
    local_node1_disp_y = fdisplacement(2 * node1);
    local_node2_disp_x = fdisplacement(2 * node2 - 1);
    local_node2_disp_y = fdisplacement(2 * node2);

    stress(i) = E(i) * ((local_node2_disp_x - local_node1_disp_x) * cos_theta + (local_node2_disp_y - local_node1_disp_y) * sin_theta) / L;
end

% Calculate deformed nodes
deformed_nodes = nodes + reshape(fdisplacement, 2, n).' * scale;
end