test3.m

tm=10;
dt=0.01;
ta=0:dt:tm;
Lta=length(ta);

d0=0.01; % distance of stability
ss=0.03; % square size
m=1e-3;
k=1; % spring coefficient
fcr=2*sqrt(k/m); % critical friction coeffcient
f=2*fcr; % friction coeficent


r=zeros(2,0);
g=[0; -10]; % gravitation

R=0.1;
R2=R^2;
ph=-0.1; % plate level
% make circle:
for yt=-R:d0:R
    for xt=-R:d0:R
        if (xt^2+yt^2<=R2)
            r=[r [xt;yt]];
        end
    end
end
hr=plot(r(1,:),r(2,:),'k.');
hold on;
hpl=plot([-0.15 0.15],ph*[1 1],'r--');
axis equal;

v=zeros(size(r)); % velocity
N=size(r,2);
Na=1:N;

bi=false(1,N);
d=NaN(1,N);
Nx=NaN(N,N);
Ny=NaN(N,N);
for tc=2:Lta
    t=ta(tc);
    D=NaN(N,N); % distance matrix
    % NaN when not need to calulate force
    for n1=1:N-1
        rt=r(:,n1);
        ind2=n1+1:N;
        bi=(abs(rt(1)-r(1,:))<=ss)&(abs(rt(2)-r(2,:))<=ss)&(Na>n1); % rest points in square
        d(bi)=sqrt((rt(1)-r(1,bi)).^2+(rt(2)-r(2,bi)).^2);
        D(n1,bi)=d(bi);
        D(bi,n1)=D(n1,bi)';
        d((d==0)&bi)=NaN; % points in the same place make no force between each other
        Nx(n1,bi)=(r(1,bi)-rt(1))./d(bi);
        Nx(bi,n1)=-Nx(n1,bi)';
        Ny(n1,bi)=(r(2,bi)-rt(2))./d(bi);
        Ny(bi,n1)=-Ny(n1,bi)';
    end
    
    %v(1,i1)=v(1,i1)+dt*sum((Nx(i1,:)*(D(i1,:)-d0))/m);
    %Nx(i1,i2)*(v(1,i1)-v(1,i2))
    %Ny(i1,i2)*(v(2,i1)-v(2,i2))
    % pr=Nx(i1,i2)*(v(1,i1)-v(1,i2))+Ny(i1,i2)*(v(2,i1)-v(2,i2))
    % Nx(i1,i2)*pr
    % Ny(i1,i2)*pr
    % v(1,i1)=v(1,i1)+(dt*f/m)*sum(Nx(i1,:).*(Nx(i1,:).*(v(1,i1)-v(1,:))+Ny(i1,i2)*(v(2,i1)-v(2,:))));
    
    dvx=bsxfun(@minus,v(1,:),v(1,:)');
    
    %A=(D-d0)/m; % acelerations from edges
    %v=v+dt*;
end