Elastisk kollision (2D)

Kollisioner af objekter foretages nemmest i CM systemet...

% Collide in CM
clf; clc; clear;

s1 = struct(); % Første kugle
s1.r = [-10, -5];
s1.radius = 2;
s1.m = s1.radius * 10;
s1.p = s1.m .* [8, -2];

s2 = struct(); % Anden kugle
s2.r = [10, 2];
s2.radius = 2;
s2.m = s2.radius * 10; % Lad massen være skaleret med radius
s2.p = s1.m .* [-5, 4];


dt = 0.1; % tidsskridt
boks_sider = 20; % (m)

t = 0; % start tid (s)
tf = 100; % slut tid (s)
 
while t < tf
    t = t + dt; % Opdater tid
    
    % Opdater position
    s1.r = s1.r + dt.* s1.p./s1.m;    
    s2.r = s2.r + dt.* s2.p./s2.m;
    
    % Tegn cirkler
    [x,y,z] = cylinder(s1.radius,200);
    plot(s1.r(1) + x(1,:),s1.r(2) + y(1,:), 'k')
    hold on
    
    [x,y,z] = cylinder(s2.radius,200);
    plot(s2.r(1) + x(1,:),s2.r(2) + y(1,:), 'r')
    
 
    % Tegn boks
    plot([-boks_sider -boks_sider], [-boks_sider boks_sider], 'k') % venstre side 
    plot([boks_sider boks_sider], [-boks_sider boks_sider], 'k') % højde side
    plot([-boks_sider, boks_sider], [boks_sider boks_sider], 'k') % top
    plot([-boks_sider, boks_sider], [-boks_sider -boks_sider], 'k') % top


    hold off
    axis equal
    axis([-21, 21, -21, 21])
    
    dr = norm(s1.r - s2.r);
    if dr <= (s1.radius + s2.radius) % Kollision!
        
                 hold on
         % Plot hastighedsvektor
         s1.p_unit = s1.p ./norm(s1.p) * 2;
         s2.p_unit = s2.p ./norm(s2.p) * 2;
         pu = plot([s1.r(1), s1.r(1) + s1.p_unit(1)], [s1.r(2), s1.r(2) + s1.p_unit(2)],'g','LineWidth',1.2)        
         plot([s2.r(1), s2.r(1) + s2.p_unit(1)], [s2.r(2), s2.r(2) + s2.p_unit(2)],'g','LineWidth',1.2)

         % Plot linje mellem de to OPDATEREDE punkter
         pos_cm_p = plot([s1.r(1), s2.r(1)], [s1.r(2), s2.r(2)],'.-b')
                 
         % Plot centre of mass point
         r_cm = (s1.m * s1.r + s2.m * s2.r)./(s1.m + s2.m)
         cm_plot = plot(r_cm(1), r_cm(2),'*k')
          
        % Take a step back, we had a collision..
        % TODO: Vi bør beregne det eksakte kollisionstidspunkt..
        

        % Calculate the velocity
        s1.v = s1.p./s1.m;
        s2.v = s2.p./s2.m;

        % Center of mass frame
        M = s1.m + s2.m;
        v_cm = ( (s1.m .* s1.v) + (s2.m .* s2.v) ) ./ M;
        
        s1.v = -(s1.v - v_cm) + v_cm;
        s2.v = -(s2.v - v_cm) + v_cm;
        
        % Back to momentum!

        s1.p = s1.v .* s1.m;
        s2.p = s2.v .* s2.m;
        
        
        %% Move along, we are done now!
        
        s1.p_unit = s1.p ./norm(s1.p) * 2;
        s2.p_unit = s2.p ./norm(s2.p) * 2;
        plt_mom_after = plot([s1.r(1), s1.r(1) + s1.p_unit(1)], [s1.r(2), s1.r(2) + s1.p_unit(2)],'y','LineWidth',1.2)        
        plot([s2.r(1), s2.r(1) + s2.p_unit(1)], [s2.r(2), s2.r(2) + s2.p_unit(2)],'y','LineWidth',1.2)


         % Plot linje mellem de to punkter et-tidsskridt tilbage
         plt_pos_after_coll = plot([s1.r(1), s2.r(1)], [s1.r(2), s2.r(2)],'.--y')

         pause(0.1)
    end
    
    % Kollision med væg S1
    if s1.r(2)+s1.radius >= boks_sider
        s1.p = [s1.p(1) -s1.p(2)];
    end
     
    if s1.r(1)+s1.radius >= boks_sider
        s1.p = [-s1.p(1) s1.p(2)];
    end
     
    if s1.r(2)-s1.radius <= -boks_sider
        s1.p = [s1.p(1) -s1.p(2)];
    end
     
    if s1.r(1)-s1.radius <= -boks_sider
        s1.p = [-s1.p(1) s1.p(2)];
    end
    
    % Kollision med væg S2
    if s2.r(2)+s2.radius >= boks_sider
        s2.p = [s2.p(1) -s2.p(2)];
    end

    if s2.r(1)+s2.radius >= boks_sider
        s2.p = [-s2.p(1) s2.p(2)];
    end

    if s2.r(2)-s2.radius <= -boks_sider
        s2.p = [s2.p(1) -s2.p(2)];
    end

    if s2.r(1)-s2.radius <= -boks_sider
        s2.p = [-s2.p(1) s2.p(2)];
    end
 
    pause(0.01)
    
end
Download kildekode