% testprogramm Molekulardynamik
clf;hold off; % ggf. Plotfenster freigeben

L=6;  % Kantenlaenge der Box

% Anfangskonfiguration:
N=5*5;  % Teilchenzahl
x=zeros(N,1); y=zeros(N,1); % reservieren
vx=zeros(N,1); vy=zeros(N,1); % reservieren

i=0;
for yi=1:5
    for xi=1:5
        i=i+1;
        x(i)=xi; y(i)=yi;
    end 
end

% Bild Startkonfiguration:
plot(x,y,'*r');axis([0 L 0 L]);hold on;
pause;

% zufaellige Anfangsimpulse  
vmax=1; % Maximalwert
vx=vmax*(2*rand(N,1)-1);vy=vmax*(2*rand(N,1)-1); % 0 < rand < 1 zufaellig
vx=vx-mean(vx);vy=vy-mean(vy); % Gesamtimpuls Null

tf=20;              % Zielzeit der Integration
tau=0.005;          % Schrittweite
nstep=fix(tf/tau);  % Anzahl Schritte

[X,Y,VX,VY]=leapfrog(x,y,vx,vy,L,tau,nstep);

% plot aller Trajektorien mit duennen Punkten
for i=1:N
    plot(X(i,:),Y(i,:),'MarkerSize',1,'Marker','.','LineStyle','none');
end
pause
%print -depsc molecular.eps

% zeitliche Folge von Konfigurationen
strob=fix((tf/50)/tau); % Haeufigkeit Plot
for n=1:strob:nstep+1
    clf
    plot(X(:,n),Y(:,n),'.');axis([0 L 0 L])
    pause(0.2)
end
%print -depsc mol_final.eps
% Parameter und letzte Konfiguration speichern
x=X(:,nstep+1);y=Y(:,nstep+1); vx=VX(:,nstep+1);vy=VY(:,nstep+1);
save gas.mat N L x y vx vy