%function orbits(posp,masap,velp,masas,vels,poss,dist,paso,it,pau)
%
%Simula de forma aproximada las órbitas de los cuatro planetas interiores 
%alrededor del Sol.
%
%posp: un vector de longitud igual al número de planetas en el que se
%introduce la posición de los planetas en número de veces con respecto a 
%la posición original.
%
%masap: un vector de longitud igual al número de planetas en el que se
%introduce la masa de los planetas en número de veces con respecto a 
%la masa original.
%
%velp: un vector de longitud igual al número de planetas en el que se
%introduce el número de veces que cambia la velocidad inicial de cada planeta 
%en proporción a la velocidad inicial real.
%
%masas: es un vector de longitud dos. El primer valor es el número de veces
%que cambia la masa del sol en proporción al valor real. El segundo premite
%introducir una segunda estrella, el valor indica el número de masas
%solares que tiene la estrella, si vale cero no la introduce.
%
%vels: es un vector de longitud dos correspondiente a cada estrella.
%Es el número de veces que cambia la velocidad de la estrella en proporción a 
%la velocidad original de La Tierra.
%
%poss: define la posición de la segunda estrella, en unidades astronónicas. Con
%respecto a la posición del Sol.
%
%dist: se introduce el número de unidades astronómicas que mide el panel de
%visualización.
%
%paso: modifica el paso de integración. Esta en segundos.
%
%it: Cada cuantas iteraciones pinta.
%
%pau: permite introducir una pausa entre figuras en segundos.
%
%Ejemplos:
%
%Simulación próxima a la situación original:
%orbits([1,1,1,1],[1,1,1,1],[1,1,1,1],[1,0],[0,0],[2,2],3,1e3,5e2,0)
%
%Simulación de una estrella atravezando el Sistema Solar interno:
%orbits([1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1],[0,3],[5,2],5,1e3,5e2,0)
%
%Guido Santos - 31/08/2014

function orbits(posp,masap,velp,masas,vels,poss,dist,paso,it,pau)

close all
    
nombres  = {' Sol',' Mercurio',' Venus',' Tierra',' Marte',' Estrella'};
colors   = {'*y','.m','.y','.b','.r','*y'};

iters    = 1e6;
planetas = length(nombres);
dir      = zeros(planetas,2);
d        = zeros(planetas,1);
f        = zeros(planetas,1);
a        = zeros(planetas,2);
entrada  = zeros(planetas,1);
crash    = zeros(planetas,1);
tiempo   = 0;%s
G        = 6.67384e-11;%m^3*Kg^-1*s^-2

%Masa del objeto
m(1,1)   = 1.989e30*masas(1);%Kg
m(2,1)   = 3.285e23*masap(1);%Kg
m(3,1)   = 4.867e24*masap(2);%Kg
m(4,1)   = 5.972e24*masap(3);%Kg
m(5,1)   = 6.390323*masap(4);%Kg
m(6,1)   = 1.989e30*masas(2);%Kg

%Radio del cuerpo
r(1,1)   = (((m(1,1)/(1.989e30/((4/3)*pi*(6.958e8)^3)))*3)/(4*pi))^(1/3);%m
r(2,1)   = 2.440e6;%m
r(3,1)   = 6.052e6;%m
r(4,1)   = 6.371e6;%m
r(5,1)   = 3.390e6;%m
r(6,1)   = (((m(6,1)/(1.989e30/((4/3)*pi*(6.958e8)^3)))*3)/(4*pi))^(1/3);%m

%Distancia orbitsl media
s(1,:)   = [0,0];%m
s(2,:)   = [posp(1)*5.78943763e10,0];%m
s(3,:)   = [posp(2)*1.0820893e11,0];%m
s(4,:)   = [posp(3)*1.49597870691e11,0];%m
s(5,:)   = [posp(4)*2.27936640e11,0];%m
s(6,:)   = [1.49597870691e11*poss(1),1.49597870691e11*poss(2)];%m

%Velocidad orbitsl media
v(1,:)   = [0,29780*vels(1)];%m*s^-1
v(2,:)   = [0,47872.5*velp(1)];%m*s^-1
v(3,:)   = [0,35021.4*velp(2)];%m*s^-1
v(4,:)   = [0,29780*velp(3)];%m*s^-1
v(5,:)   = [0,24077*velp(4)];%m*s^-1
v(6,:)   = [-29780*vels(2),0];%m*s^-1

x0       = max(abs(dist*s(4,:)));%m

for k = 1:iters,
    
    tiempo = tiempo + paso;
    s0     = s;
    
    for p = 1:planetas,

        for q = 1:planetas,
            
            if (p ~= q && crash(p) == 0),
                
                if (masas(2) == 0 && (p == 6 || q == 6)),
                
                else
                
                    if (sqrt(sum((s(p,:) - s(q,:)).^2)) <= (r(q) + r(p)) && crash(q) == 0),

                        v(q,1)    = (m(p)*v(p,1) + m(q)*v(q,1))/m(q);
                        v(q,2)    = (m(p)*v(p,2) + m(q)*v(q,2))/m(q);
                        m(q)      = m(q) + m(p);
                        m(p)      = 0;
                        crash(p)  = 1;
                        colors{q} = '*r';

                    end
                
                end
                
                dir(p,:)                           = (s(p,:) - s(q,:))./abs(s(p,:) - s(q,:));
                dir(p,abs(s(p,:) - s(q,:)) == 0)   = 0;
                d(p)                               = sqrt(sum((s(p,:) - s(q,:)).^2));%m  
                f(p)                               = (m(q)*G)/(d(p)^2);%m*s^-2
                a(p,1)                             = abs(cos(atan2(s(p,2) - s(q,2),s(p,1) - s(q,1)))*f(p));%m*s^-2
                a(p,2)                             = abs(sin(atan2(s(p,2) - s(q,2),s(p,1) - s(q,1)))*f(p));%m*s^-2
                v(p,1)                             = v(p,1) - dir(p,1)*a(p,1)*paso;
                v(p,2)                             = v(p,2) - dir(p,2)*a(p,2)*paso;        
            
            end
  
        end
        
        s(p,1)                             = s(p,1) + v(p,1)*paso;
        s(p,2)                             = s(p,2) + v(p,2)*paso; 
        
    end
    
    for l = 2:planetas - 1,

        if ((s(l,2)*s0(l,2)) ~= abs(s(l,2)*s0(l,2)) && crash(l) == 0),

            entrada(l) = entrada(l) + 1;

            if entrada(l) == 2,
                
                dias   = tiempo/(60*60*24);
                texto  = strcat('Año en',nombres{l},': %g días\n');
                fprintf(texto,dias)
                
            end

        end
        
    end
    
    if k/it == fix(k/it)

        figure(1)
        plot(0,x0,'k')
        hold on

        for p = 2:planetas - 1,
            
            if crash(p) == 0,
                
                plot(s(p,1),s(p,2),colors{p},'MarkerSize',15)
            
            end

        end
        
        if (masas(1) ~= 0 && crash(1) == 0),
            
            plot(s(1,1),s(1,2),colors{1},'MarkerSize',2*r(1,1)*1e-8)
        
        end
        
        if (masas(2) ~= 0 && crash(planetas) == 0),
            
            plot(s(planetas,1),s(planetas,2),colors{planetas},'MarkerSize',2*r(planetas,1)*1e-8)
        
        end
        
        xlim([-x0 x0])
        ylim([-x0 x0])
        
        hold off
        whitebg(1,'black')
        dias   = fix(tiempo/(60*60*24));
        dias   = num2str(dias);
        titulo = strcat('Días: ',dias);
        title(titulo)    
        drawnow
        pause(pau)

    end
    
end

end