function vannbolger
%Bolgepakke med smal fordeling i k-rommet
clear all;
% Tyngdens akselerasjon
g = 9.81;
% Overflatespenning
gamma = 0.073;
% Massetetthet
rho = 1000;
% Parameteren a angir spredningen i bolgetall, slik at
% kmin = k0*(1-a) og kmax = k0*(1+a)
a = 0.10;
% Senter-bolgelengde (m)
lambda0 = 1.0;
% Senter-bolgetall (1/m)
k0 = 2*pi/lambda0;
% Senter-vinkelfrekvens
w0 = sqrt(g*k0+gamma*k0*k0*k0/rho);
% Fasehastighet vf og gruppehastighet vg
vf = w0/k0
vg = (g+3*gamma*k0*k0/rho)/(2*w0)
% Antall fourierkomponenter (oddetall)
nk = 25;
for i = 1:nk
    k(i) = k0*(1+a*(2*i-nk-1)/(nk-1));
    w(i) = sqrt(g*k(i)+gamma*k(i)*k(i)*k(i)/rho);
    % Gaussisk form ("klokkeform") i k-rommet
    b(i) = exp(-4*(2*i-nk-1)^2/(nk-1)^2);
end;

% Bestem skjermstorrelsen
scrsz = get(0,'ScreenSize');
% Naa er scrsz(1)=venstre posisjon til skjermbildet = 1
%        scrsz(2)=nedre posisjon til skjermbildet = 1
%        scrsz(3)=skjermens bredde = 1280 (for eksempel)
%        scrsz(4)=skjermens hoyde = 1024 (for eksempel)
% Lag forste figur overst til venstre
figure('Position',[10 scrsz(4)/1.9 scrsz(3)/2.6 scrsz(4)/2.6]);
% Plott fordelingen i k-rommet
plot(k,b,'-*');
xlabel('k (1/m)');
ylabel('andel');
title('Fordeling i k-rommet');
pause(1.00);

kmin = k0*(1-a);
kmax = k0*(1+a);
wmin = sqrt(g*kmin+gamma*kmin*kmin*kmin/rho);
wmax = sqrt(g*kmax+gamma*kmax*kmax*kmax/rho);
N = 400;
xmax = round(8*pi/(kmax-kmin));
ymax = nk*0.5;
tmax = round(xmax/vf)
x = linspace(-xmax,xmax,N);

% Lager bolgepakken ved t=0, dvs y(x,0)
t = 0.0;
for i = 1:N
  y(i) = 0.0;
  for j = 1:nk
    y(i)=y(i)+b(j)*sin(k(j)*x(i)-w(j)*t);
  end;
end;
figure('Position',[scrsz(3)/2 scrsz(4)/1.9 scrsz(3)/2.6 scrsz(4)/2.6]);
plot(x,y,'-b');
hold on;
xlabel('Posisjon (m)');
ylabel('Amplitude (vilkårlig enhet)');
title('Bolgepakke ved t=0');
pause(1.00);

figure('Position',[10 30 scrsz(3)/1.05 scrsz(4)/2.6]);
%Animasjon av bølgen
p = plot(x,y,'-','EraseMode','xor'); 
axis([0 xmax -ymax ymax])
xlabel('Posisjon (m)');
ylabel('Amplitude (vilkårlig enhet)');
title('Animasjon av bolgepakken');
hold on;
for i=1:tmax*5
    t=i*0.2;
    for j=1:N
        y(j)=0.0;
        for m = 1:nk
            y(j)=y(j)+b(m)*sin(k(m)*x(j)-w(m)*t);
        end; 
    end;
    set(p,'XData',x,'YData',y); 
    drawnow;
    plot(t*vf,ymax*0.9,'.r');   % Plotter rød stav for å vise fasehastigheten 
    plot(t*vg,-ymax*0.9,'.g');  % Plotter grønn stav for gruppehastighet
    pause(0.25);                % Forsink animasjonen (juster denne for din pc)
end;
pause(1.00);

figure('Position',[scrsz(3)/2 30 scrsz(3)/2.6 scrsz(4)/2.6]);
% Lager bolgepakken ved t=10*tmax
t = 10.0*tmax;
x = linspace(3.75*xmax,5.75*xmax,N);
for i = 1:N
  y(i) = 0.0;
  for j = 1:nk
    y(i)=y(i)+b(j)*sin(k(j)*x(i)-w(j)*t);
  end;
end;
plot(x,y,'-b');
hold on;
xlabel('Posisjon (m)');
ylabel('Amplitude (vilkårlig enhet)');
title('Bolgepakke ved t=10tmax');
