function [fx fy] = LJkraft(x,y)
% file LJkraft.m
%
% Kraft aus Lennard-Jones Potential fuer N Teilchen (2dim.)
N=length(x);
fx=zeros(size(x)); fy=zeros(size(fx));     % reservieren
% loop ueber Paare von verschiedenen Teilchen (jedes Paar 1 mal, i<j !)
%
for i=1:N-1
    for j=i+1:N
    r2=(x(i)-x(j))^2+(y(i)-y(j))^2; % Abstand^2
    ri2=1/r2;
    fac=24*ri2^4*(2*ri2^3-1);
    fx(i)=fx(i)+(x(i)-x(j))*fac;% Beitrag dieser Paarkraft AUF i
    fy(i)=fy(i)+(y(i)-y(j))*fac;%
    fx(j)=fx(j)-(x(i)-x(j))*fac;% Beitrag dieser Paarkraft AUF j
    fy(j)=fy(j)-(y(i)-y(j))*fac;% Vorz: actio=reactio
    end % loop j
end     % loop i; alle Kraefte beruecksichtigt
