%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% GRAVITATION %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
clc;
%n=randi([10,15]);
n=3;%n är antalet kroppar, 10-15st
XYM=rand(3,n)*(1000-10);                %En matris med x- och y-koordninater
                                        %och massa för varje kropp.
                                        %(varje kolumn är en kropp)

t=0;                                    %tid, vilket används för
                                        %beräkningar, är i början 0

v0X=zeros(1,n);                         %Begynnelsehastigheter=0
v0Y=zeros(1,n);

G=6.673e-11;                            %Gravitationskonstanten



while n>1


    %Kropparna plottas
    sz=XYM(3,:);
    scatter(XYM(1,:),XYM(2,:),sz,XYM(3,:),'filled')
    axis([-250 1250 -250 1250])
    whitebg('k')
    title(['Gravitationssimulering, antal kroppar: ',num2str(n),'st']);
    xlabel('x','FontSize',14);
    ylabel('y','FontSize',14);
    colormap(jet(size(XYM,1)));
    grid on

    pause(.02)


    %%Avstånden beräknas mellan kropparna, även i x och y



    %dist är absolutavstånden mellan kropparna
    dist = hypot(XYM(1, :) - XYM(1, :)', XYM(2, :) - XYM(2, :)');


    %distX är avstånden i x
    m=n;
    distX=zeros(m,n);
    for i2=1:m
        for j2=1:n

            distX(j2, i2) = XYM(1, i2) - XYM(1,j2);

        end
    end



    %distY är avstånden i y
    distY=zeros(m,n);
    for i3=1:m
        for j3=1:n

            distY(j3, i3) = XYM(2, i3) - XYM(2,j3);

        end
    end


    %Krafterna beräknas med formeln F=G*M1*M2/dist^2. Krafterna i x och y
    %ska även beräknas, vilket är varför distX och distY behövs
    F=zeros(m,n);
    for i4=1:m
        for j4=1:n

            F(j4,i4)=(G*XYM(3,i4)*XYM(3,j4))/(dist(i4,j4).^2);

        end
    end

    F(~isfinite(F))=0;
    F=-F;


    %k är en variabel vars syfte är att göra beräkningar enklare

    %Vinkeln mellan linjen från en kropp till en annan och x-axeln
    %beräknas, för alla kroppar
    k=zeros(m,n);
    vinkel=zeros(m,n);
    for i5=1:m
        for j5=1:m

            if distX(i5,j5)<0 && distY(i5,j5)<0

                k(i5,j5)=distY(i5,j5)./distX(i5,j5);
                vinkel(i5,j5)=atand(k(i5,j5));
                vinkel(i5,j5)=vinkel(i5,j5)+180;

            elseif distX(i5,j5)<0 && distY(i5,j5)>0

                k(i5,j5)=distY(i5,j5)./distX(i5,j5);
                vinkel(i5,j5)=atand(k(i5,j5));
                vinkel(i5,j5)=abs(vinkel(i5,j5))+90;

            elseif distX(i5,j5)==0 && distY(i5,j5)==0

                vinkel(i5,j5)=0;

            else

                k(i5,j5)=distY(i5,j5)./distX(i5,j5);
                vinkel(i5,j5)=atan(k(i5,j5));

            end


        end
    end


    %Krafterna i varje axel beräknas
    FX=F.*cosd(vinkel);
    FY=F.*sind(vinkel);
    FX=sum(FX);
    FY=sum(FY);



    %Acceleration i x och y beräknas, för varje kropp
    aX=FX./XYM(3,:);
    aY=FY./XYM(3,:);

    %Hastigheter i x och y beräknas från accelerationen
    v0X=v0X+aX*t;
    v0Y=v0Y+aY*t;

    %Avstånden som varje kropp färdas efter t tid beräknas från hastigheten
    %och accelerationen
    sX=v0X*t+(aX.*t^2)/2;
    sY=v0Y*t+(aY.*t^2)/2;

    %Avstånden som varje kropp färdas läggs till på deras nuvarande
    %position
    XYM(1,:)=XYM(1,:)+sX;
    XYM(2,:)=XYM(2,:)+sY;


    %Vi begränsar x- och y-koordinaterna till mellan -1000 och 2000
    for i6=1:m

        if XYM(1,i6)>1250

            v0X(1,i6)=0;
            XYM(1,i6)=1250;

        elseif XYM(1,i6)<-250

            v0X(1,i6)=0;
            XYM(1,i6)=-250;

        elseif XYM(2,i6)>1250

            v0Y(1,i6)=0;
            XYM(2,i6)=1250;

        elseif XYM(2,i6)<-250

            v0X(1,i6)=0;
            XYM(2,i6)=-250;

        else

        end

    end



    %%Kropparna kommer adderas om de är like med eller mindre än 10
    %%längdenheter ifrån varandra

    %d är en matris som ger ettor om två kroppar är nära varandra
    d=dist<=10;
    d=d-eye(m,n);

    %c summerar d
    c=sum(sum(d));

    %om c>0 så är två kroppar nära varandra och måste adderas
    if c>0

        %r används för att hitta vilka kroppar som är nära varandra i
        %XYM-matrisen genom att hitta var ettorna finns i matrisen d
        [r,c]=find(d);
        %a är den första kroppen i kollisionen
        a=r(1,1);
        %b är den andra kroppen i kollisionen
        b=r(2,1);

        %De nya hastigheterna efter kollisionen beräknas och ersätter
        %hastigheten för kropp a
        v0X(:,a)=(XYM(3,a).*v0X(:,a)+XYM(3,b).*v0X(:,b))/(XYM(3,a)+XYM(3,b));
        v0Y(:,a)=(XYM(3,a).*v0Y(:,a)+XYM(3,b).*v0Y(:,b))/(XYM(3,a)+XYM(3,b));

        %Hastigheterna för kroppen b tas bort
        v0X(:,b)=[];
        v0Y(:,b)=[];

        %Massorna av a och b adderas och ersätter massan för a
        XYM(3,a)=XYM(3,a)+XYM(3,b);
        %Kropp b tas bort
        XYM(:,b)=[];

        %n, vilket är antalet kroppar, reduceras
        n=n-1;
        %m används för loopar och måste vara lika med n
        m=n;

    end

    %Det går 60 sekunder för varje while-loop
    t=t+60;

end