- /
- 
        sun
        on 23 Oct 2022
        
        
 
    - 19
- 105
- 4
- 1
- 280
%credits to Sebastian Kraemer Mini-hack2021
figure(Color='k')
hold;
% function shortcuts
m=@(a,c)max(min(a,c),-c);
s=@vecnorm;
n=@(p)p./s(p);
% random uniformly distributed values on the unit sphere
P=randn(3,1e3);
% iterate gravitational force exchange
A=hot(580);
for k=1:540
    Q=P;
    % distance of all points to each other
    D=reshape(P,3,1,[])-P;
    % calculate gravitational exchange force
    U=sum(m(D./s(D).^2,10),3)/1e4-1e-3;
    % add bounded, cumulative displacement and project result to unit sphere
    P=n(P+n(U).*m(s(U),.01));   
    % plot imitating more points than actually calculated
    h=scatter(P(1,:),P(2,:),m(s(P-Q).^-2/2e4,2e5/(2+k)),A(40+k,:),'f'); % 
    alpha(h,.03)
end
axis off equal


 

 
             
             
             
               
               
              