Partikel i en boks

% Partikel i en boks
%
% Opgaver:
%   - Gennemse kode (kommenter ting ud du ikke forstår, hvad sker der? Slå kommandoer op med help eller doc)
%   - Prøv at se om du forstår alle linjer
%   - Hvilken betydning har begyndelsesbetingelserne?
%   - Der er mulighed for at tilføre en konstant acceleration, hvad vil det
%   gøre for partiklens bevægelse, prøv både and indføre forskellige
%   værdier i x og y komponenterne.
%   - Simuleringen kører i maksimalt 40000 iterationer, er dette den
%   smarteste måde at bestemme længden af simuleringen? prøv at indføre en variabel t som for hver iteration bliver adderet med dt.
%   - Hvad har størrelsen af tidsskridtet af betydning for simuleringen?
%   - Beregn den kinetiske energi af partiklen undervejs, ændrer den sig?
%   - Randbetingelserne er blot simpel reflektion af hastighedsvektorens x
%   eller y komponent, kan du tilføje en dæmpningsfaktor, som gør at
%   partiklen taber energi når den rammer kanten?
%   - Der er en termineringsbetingelse i løkken som afbryder simuleringen
%   hvis hastigheden af partiklen når under 0.01, denne sætning kan
%   omskrives til at stoppe ved en specifik kinetisk energi, prøv dette...



clf; clear; clc;

dt = 0.001; % Tidsskridt (s)

x = [10 10]; % (x,y) Begydelsesposition (m)
v = [6 -8.5]; % (x,y) Begyndelseshastighed (m/s)
a = [0, 0]; % (x,y) Konstant acceleration (m/s^2)

plot(x(1),x(2),'or')
hold on


boks_sider = 20; % (m)

% Tegn boks
plot([0 0], [0 boks_sider], 'k') % venstre side 
plot([boks_sider boks_sider], [0 boks_sider], 'k') % højde side
plot([0, boks_sider], [boks_sider boks_sider], 'k') % top
plot([0, boks_sider], [0 0], 'k') % bund

axis([-2 boks_sider+2 -2 boks_sider+2])
axis('equal') 

for i=1:40000 % Maks antal iterationer
    
    v = v + a.*dt; % Hastighedsopdatering
    x = x + v.*dt; % Positionsopdatering
    
    
    % Randbetingelser 
    if x(2) >= boks_sider
        v = [v(1) -v(2)];
    end
    
    if x(1) >= boks_sider
        v = [-v(1) v(2)];
    end
    
    if x(2) <= 0
        v = [v(1) -v(2)];
    end
    
    if x(1) <= 0
        v = [-v(1) v(2)];
    end
    
    
    if x(1) >= boks_sider || x(2) >= boks_sider || x(1) <= 0 || x(2) <= 0
       plot(x(1), x(2), '+g') 
    end
    if mod(i, 50) == 0
        plot(x(1), x(2), 'ob')    
    end
    
    
    if sqrt(v(1)*v(1) + v(2)*v(2)) <= 0.01
        break
    end

    pause(0.0005)

end