MassSpringDamperSimulator.m

function MassSpringDamperSimulator()
% Mass Spring Damper Simulator (1-DOF)
% Interactive simulation with base excitation
% Replicates the web-based simulator functionality

    % Create figure
    fig = figure('Name', 'Mass Spring Damper Simulator (1-DOF)', ...
                 'NumberTitle', 'off', ...
                 'Position', [100, 100, 1400, 800], ...
                 'Color', [0.95, 0.95, 0.98]);
    
    % Main title
    uicontrol('Style', 'text', ...
              'String', 'Mass Spring Damper Simulator', ...
              'FontSize', 24, ...
              'FontWeight', 'bold', ...
              'Units', 'normalized', ...
              'Position', [0.25, 0.93, 0.5, 0.05], ...
              'BackgroundColor', [0.95, 0.95, 0.98], ...
              'HorizontalAlignment', 'center');
    
    uicontrol('Style', 'text', ...
              'String', 'Interactive 1-DOF system with base excitation. Tweak parameters and watch the animation update in real time.', ...
              'FontSize', 10, ...
              'Units', 'normalized', ...
              'Position', [0.15, 0.89, 0.7, 0.03], ...
              'BackgroundColor', [0.95, 0.95, 0.98], ...
              'HorizontalAlignment', 'center');
    
    % Parameters panel
    paramPanel = uipanel('Title', 'Parameters', ...
                        'FontSize', 12, ...
                        'FontWeight', 'bold', ...
                        'BackgroundColor', [1, 1, 1], ...
                        'Units', 'normalized', ...
                        'Position', [0.02, 0.05, 0.22, 0.82]);
    
    % Visualization panel
    visPanel = uipanel('Title', 'Mass spring damper - 1 DOF', ...
                      'FontSize', 12, ...
                      'FontWeight', 'bold', ...
                      'BackgroundColor', [1, 1, 1], ...
                      'Units', 'normalized', ...
                      'Position', [0.26, 0.05, 0.72, 0.82]);
    
    % Initialize parameters
    params = struct();
    params.mass = 1000.0;           % kg
    params.springK = 48000.0;       % N/m
    params.dampingC = 4000.0;       % N*s/m
    params.amplitudeA = 0.1;        % m
    params.frequencyF = 1.0;        % Hz
    params.L0 = 0.7;                % m
    params.blockH = 0.4;            % m
    params.blockA = 0.8;            % m
    params.timeTF = 15.0;           % s
    params.frameRate = 30;          % fps (reduced from 60 for stability)
    params.isPaused = false;
    params.time = 0;
    params.position = 0;
    params.velocity = 0;
    
    % Create sliders and labels
    yPos = 0.88;
    deltaY = 0.11;
    
    % Mass slider
    createSlider(paramPanel, 'Mass [kg]', 100, 5000, params.mass, ...
                 [0.05, yPos, 0.9, 0.05], @(val) updateParam('mass', val));
    yPos = yPos - deltaY;
    
    % Spring k slider
    createSlider(paramPanel, 'Spring k [N/m]', 1000, 100000, params.springK, ...
                 [0.05, yPos, 0.9, 0.05], @(val) updateParam('springK', val));
    yPos = yPos - deltaY;
    
    % Damping c slider
    createSlider(paramPanel, 'Damping c [N*s/m]', 0, 10000, params.dampingC, ...
                 [0.05, yPos, 0.9, 0.05], @(val) updateParam('dampingC', val));
    yPos = yPos - deltaY;
    
    % Amplitude A slider
    createSlider(paramPanel, 'Amplitude A [m]', 0, 0.5, params.amplitudeA, ...
                 [0.05, yPos, 0.9, 0.05], @(val) updateParam('amplitudeA', val));
    yPos = yPos - deltaY;
    
    % Frequency f slider
    createSlider(paramPanel, 'Frequency f [Hz]', 0.1, 5, params.frequencyF, ...
                 [0.05, yPos, 0.9, 0.05], @(val) updateParam('frequencyF', val));
    yPos = yPos - deltaY;
    
    % L0 and Block h sliders (side by side)
    createSlider(paramPanel, 'L0 [m]', 0.3, 1.5, params.L0, ...
                 [0.05, yPos, 0.42, 0.05], @(val) updateParam('L0', val));
    createSlider(paramPanel, 'Block h [m]', 0.2, 0.8, params.blockH, ...
                 [0.53, yPos, 0.42, 0.05], @(val) updateParam('blockH', val));
    yPos = yPos - deltaY;
    
    % Block a and Time TF sliders (side by side)
    createSlider(paramPanel, 'Block a [m]', 0.3, 1.5, params.blockA, ...
                 [0.05, yPos, 0.42, 0.05], @(val) updateParam('blockA', val));
    createSlider(paramPanel, 'Time TF [s]', 5, 30, params.timeTF, ...
                 [0.53, yPos, 0.42, 0.05], @(val) updateParam('timeTF', val));
    yPos = yPos - deltaY;
    
    % Frame rate control
    uicontrol('Parent', paramPanel, 'Style', 'text', ...
              'String', 'Frame rate [fps]', ...
              'FontSize', 10, ...
              'Units', 'normalized', ...
              'Position', [0.05, yPos+0.015, 0.6, 0.03], ...
              'BackgroundColor', [1, 1, 1], ...
              'HorizontalAlignment', 'left');
    
    uicontrol('Parent', paramPanel, 'Style', 'edit', ...
              'String', num2str(params.frameRate), ...
              'FontSize', 10, ...
              'Units', 'normalized', ...
              'Position', [0.7, yPos+0.015, 0.25, 0.03], ...
              'Callback', @(src, ~) updateFrameRate(str2double(src.String)));
    
    % Pause and Reset buttons
    pauseBtn = uicontrol('Parent', paramPanel, 'Style', 'pushbutton', ...
              'String', 'Pause', ...
              'FontSize', 11, ...
              'Units', 'normalized', ...
              'Position', [0.05, 0.02, 0.42, 0.06], ...
              'Callback', @togglePause);
    
    uicontrol('Parent', paramPanel, 'Style', 'pushbutton', ...
              'String', 'Reset', ...
              'FontSize', 11, ...
              'Units', 'normalized', ...
              'Position', [0.53, 0.02, 0.42, 0.06], ...
              'Callback', @resetSimulation);
    
    % Create axes for visualization
    ax = axes('Parent', visPanel, ...
              'Units', 'normalized', ...
              'Position', [0.1, 0.1, 0.85, 0.85]);
    hold(ax, 'on');
    grid(ax, 'on');
    axis(ax, 'equal');
    xlabel(ax, 'x [m]', 'FontSize', 11);
    ylabel(ax, 'y [m]', 'FontSize', 11);
    
    % Timer variable
    tmr = [];
    
    % Initialize simulation
    resetSimulation();
    
    % Start animation timer with safer settings
    tmr = timer('ExecutionMode', 'fixedRate', ...
                'Period', max(0.001, 1/params.frameRate), ...  % Ensure minimum 1ms
                'BusyMode', 'drop', ...
                'TimerFcn', @updateAnimation);
    
    try
        start(tmr);
        disp('Animation started successfully');
    catch ME
        disp(['Timer error: ' ME.message]);
    end
    
    % Cleanup on close
    fig.CloseRequestFcn = @cleanupFigure;
    
    function cleanupFigure(~,~)
        if ~isempty(tmr) && isvalid(tmr)
            stop(tmr);
            delete(tmr);
        end
        delete(fig);
    end
    
    function createSlider(parent, label, minVal, maxVal, initVal, pos, callback)
        % Label with value on same line
        labelStr = sprintf('%s: %.3f', label, initVal);
        valDisp = uicontrol('Parent', parent, 'Style', 'text', ...
                           'String', labelStr, ...
                           'FontSize', 9, ...
                           'Units', 'normalized', ...
                           'Position', [pos(1), pos(2)+0.055, pos(3), 0.025], ...
                           'BackgroundColor', [1, 1, 1], ...
                           'HorizontalAlignment', 'left');
        
        % Slider
        uicontrol('Parent', parent, 'Style', 'slider', ...
                  'Min', minVal, 'Max', maxVal, ...
                  'Value', initVal, ...
                  'Units', 'normalized', ...
                  'Position', pos, ...
                  'Callback', @(src, ~) sliderCallback(src, valDisp, callback, label));
    end
    
    function sliderCallback(src, valDisp, callback, label)
        val = get(src, 'Value');
        labelStr = sprintf('%s: %.3f', label, val);
        set(valDisp, 'String', labelStr);
        callback(val);
    end
    
    function updateParam(paramName, value)
        params.(paramName) = value;
    end
    
    function updateFrameRate(value)
        if value > 0 && value <= 60
            params.frameRate = value;
            if ~isempty(tmr) && isvalid(tmr)
                stop(tmr);
                tmr.Period = max(0.001, 1/value);
                start(tmr);
            end
        end
    end
    
    function togglePause(~, ~)
        params.isPaused = ~params.isPaused;
        if params.isPaused
            set(pauseBtn, 'String', 'Resume');
        else
            set(pauseBtn, 'String', 'Pause');
        end
    end
    
    function resetSimulation(~, ~)
        params.time = 0;
        params.position = 0;
        params.velocity = 0;
        params.isPaused = false;
        set(pauseBtn, 'String', 'Pause');
        
        % Draw initial state
        drawSystem();
    end
    
    function drawSystem()
        % Clear axes
        cla(ax);
        hold(ax, 'on');
        grid(ax, 'on');
        
        % Calculate positions
        omega = 2 * pi * params.frequencyF;
        baseDisp = params.amplitudeA * sin(omega * params.time);
        
        yBase = baseDisp;
        yMass = yBase + params.L0 + params.position;
        
        % Set axis limits dynamically
        xRange = 2.5;
        yMin = min(-0.5, yBase - 0.5);
        yMax = max(params.L0 + params.blockH + 1.0, yMass + params.blockH/2 + 0.5);
        axis(ax, [-xRange/2, xRange/2, yMin, yMax]);
        
        % Draw ground (fixed reference)
        groundY = -0.3;
        groundWidth = 2.0;
        plot(ax, [-groundWidth/2, groundWidth/2], [groundY, groundY], ...
             'k-', 'LineWidth', 3);
        
        % Draw hatching for ground
        for i = -10:10
            x = i * 0.1;
            if abs(x) <= groundWidth/2
                plot(ax, [x, x-0.05], [groundY, groundY-0.1], 'k-', 'LineWidth', 1);
            end
        end
        
        % Draw base platform
        baseWidth = params.blockA;
        baseHeight = 0.15;
        rectangle('Parent', ax, ...
                  'Position', [-baseWidth/2, yBase-baseHeight/2, baseWidth, baseHeight], ...
                  'FaceColor', [0.7, 0.7, 0.7], ...
                  'EdgeColor', 'k', ...
                  'LineWidth', 2);
        
        % Draw spring (zigzag pattern)
        springCoils = 12;
        springWidth = 0.15;
        springX = -0.3;
        numPoints = springCoils * 2 + 1;
        springY = linspace(yBase + baseHeight/2, yMass - params.blockH/2, numPoints);
        springXCoords = zeros(1, numPoints);
        springXCoords(1) = springX;
        springXCoords(end) = springX;
        for i = 2:numPoints-1
            if mod(i, 2) == 0
                springXCoords(i) = springX + springWidth;
            else
                springXCoords(i) = springX - springWidth;
            end
        end
        plot(ax, springXCoords, springY, 'k-', 'LineWidth', 2);
        
        % Draw damper
        damperX = 0.3;
        totalDamperLen = yMass - params.blockH/2 - (yBase + baseHeight/2);
        damperPistonLen = totalDamperLen * 0.7;
        damperCylLen = totalDamperLen * 0.3;
        
        % Damper piston (line)
        plot(ax, [damperX, damperX], [yBase + baseHeight/2, yBase + baseHeight/2 + damperPistonLen], ...
             'k-', 'LineWidth', 3);
        
        % Damper cylinder
        cylWidth = 0.15;
        rectangle('Parent', ax, ...
                  'Position', [damperX-cylWidth/2, yBase+baseHeight/2+damperPistonLen, cylWidth, damperCylLen], ...
                  'FaceColor', [0.6, 0.6, 0.6], ...
                  'EdgeColor', 'k', ...
                  'LineWidth', 2);
        
        % Damper connection to mass
        plot(ax, [damperX, damperX], ...
             [yBase+baseHeight/2+damperPistonLen+damperCylLen, yMass-params.blockH/2], ...
             'k-', 'LineWidth', 3);
        
        % Draw mass (block)
        rectangle('Parent', ax, ...
                  'Position', [-params.blockA/2, yMass-params.blockH/2, ...
                               params.blockA, params.blockH], ...
                  'FaceColor', [0.4, 0.3, 0.8], ...
                  'EdgeColor', 'k', ...
                  'LineWidth', 2);
        
        % Draw guide lines (dashed)
        xLim = get(ax, 'XLim');
        yEquilibrium = params.L0;
        
        plot(ax, xLim, [yEquilibrium, yEquilibrium], 'k--', 'LineWidth', 1);
        plot(ax, xLim, [yBase, yBase], 'k--', 'LineWidth', 1);
        plot(ax, xLim, [yMass, yMass], 'k--', 'LineWidth', 1);
        
        drawnow;
    end
    
    function updateAnimation(~, ~)
        if params.isPaused
            return;
        end
        
        try
            % Update time
            dt = 1/params.frameRate;
            params.time = params.time + dt;
            
            % Reset if time exceeds simulation time
            if params.time >= params.timeTF
                params.time = 0;
                params.position = 0;
                params.velocity = 0;
            end
            
            % Calculate base excitation
            omega = 2 * pi * params.frequencyF;
            baseDisp = params.amplitudeA * sin(omega * params.time);
            baseVel = params.amplitudeA * omega * cos(omega * params.time);
            
            % Calculate relative displacement
            relDisp = params.position - baseDisp;
            relVel = params.velocity - baseVel;
            
            % Forces
            springForce = -params.springK * relDisp;
            dampingForce = -params.dampingC * relVel;
            
            % Acceleration
            acceleration = (springForce + dampingForce) / params.mass;
            
            % Update velocity and position (Euler integration)
            params.velocity = params.velocity + acceleration * dt;
            params.position = params.position + params.velocity * dt;
            
            % Redraw system
            drawSystem();
        catch ME
            disp(['Animation error: ' ME.message]);
        end
    end
end