%%function bolgepakkemange

clear all;

pi = 3.141592654;
%
% Antall punkter langs x-aksen:
Nx = 1600;
% Antall tidssteg:
Nt = 500;
% Antall k-verdier:
Nk = 200;
% Konstanter; luft vann grenseflate
g = 9.81;
rho = 1000;
gamma = 0.073;
% Senter-bolgetall
% Hvis k0 = sqrt(rho*g/gamma) = 366.6 m^(-1), 
% er tyngdeleddet og overflateleddet i uttrykket for omega(k) like store.
% Da er bolgelengden ca 17 mm.
% Hvis k0 >> sqrt(rho*g/gamma): Kapillaerbolger, vg = 1.5*vf
% Hvis k0 << sqrt(rho*g/gamma): Tyngdebolger, vg = 0.5*vf
k0 = 0.04*sqrt(rho*g/gamma);
% Senter-bolgelengde
lambda0 = 2*pi/k0;
% Senter-omega
omega0 = sqrt(g*k0 + gamma*k0^3/rho);
% Spredning i k-rommet
sigma = k0/10;
% maksimal x
xmax = 20/sigma;
% Tabeller
k = linspace(0,k0+5*sigma,Nk);
omega = sqrt(g.*k + gamma.*k.^3./rho);
% Fourierkoeffisienter:
y0 = exp(-(k-k0).^2./sigma^2);
ymax = max(y0)*60;
x = linspace(0,xmax,Nx);
y = linspace(0,ymax,Nx);
% Animasjon av overflatebolgepakke
p = plot(x,y,'-','LineWidth',1.5,'EraseMode','xor'); 
axis([0 xmax -ymax ymax])
hold on
xlabel('posisjon (m)');
ylabel('Utsving');
title('Sum av mange harmoniske');
for i=1:Nt
    teller = i
    % Velger tidssteg:
    t = i*10*2*pi/(5.0*sigma*Nt);
    %
    for jx = 1:Nx
        y(jx) = 0;
        for jk = 1:Nk
            y(jx) = y(jx) + y0(jk)*sin(k(jk)*x(jx)-omega(jk)*t);
        end;
    end;
    set(p,'XData',x,'YData',y) 
    % MATLAB help: drawnow flushes the event queue and updates the figure window
    drawnow    
    pause(0.01);  % Forsinker framvisningen i 0.1 s (juster hvis nodvendig)
end
