//Solar system simulation

clc
clear; 
xdel(winsid()); //winsid()-return the list of graphics windows
//xdel() — delete a graphics window

t0=0
t_fin=10*365.25*24*60*60; // years
step=10*60*60;
t=t0:step:t_fin

G=6.67408e-11;
Nbodies = 10;

function dy=f(t,y)
  dy=zeros(y);
  dy(1:Nbodies*3)=y(Nbodies*3+1:$);
  for i=1:Nbodies
    ii=Nbodies*3+1+(i-1)*3;
    for j=1:Nbodies
      jj=Nbodies*3+1+(j-1)*3;
      if jj<>ii
        dR=y(1+(j-1)*3:j*3)-y(1+(i-1)*3:i*3);
        dy(ii:ii+2)=dy(ii:ii+2)+G*ObjMass(j).*dR./norm(dR)^3; 
      end
    end
  end 
endfunction

// Solar System obects mass, radius, initial position and velocity
// ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/de436/
// http://ilrs.gsfc.nasa.gov/docs/2014/196C.pdf

ObjName=["Sun","Mercury","Venus","Earth","Mars","Jupiter",...
         "Saturn","Uranus","Neptune","Pluto"];
         
// Object colors
ObjColor=[7, 2, 3, 4, 5, 6, 1, 8, 9, 10];

// Planetary masses are given as "mass parameters" 
// A mass parameter of 1 = 6.72006*10^33 kg

ObjMass=6.72006e+33*[...//units=[kg]
         0.295912208285591100D-03,... //mass of Sun
         0.491248045036476000D-10,... // mass of Mercury
         0.724345233264412000D-09,... // mass of Venus
         0.899701139019987100D-09,... // mass of Earth/Moon system
         0.954954869555077000D-10,... // mass of Mars
         0.282534584083387000D-06,... // mass of Jupiter
         0.845970607324503000D-07,... // mass of Saturn
         0.129202482578296000D-07,... // mass of Uranus
         0.152435734788511000D-07,... // mass of Neptune
         0.217844105197418000D-11]    // mass of Pluto

// Object Size
ObjRadius=log([900000 4850 12000 12000 6790 142980 120536 51118 49500 2368]./1.00001);

// Object position and velocity
// the epoch is 1969-Jun-28 00:00:00.0000 UTC = Julian day
// 2440400.500000000 = Unix second -16156800 = Unix day -187
// Positions are in the ICRF J2000 reference frame in astronomical units
//(1 AU = 149597870.700 km)

ObjPos=149597870700*[...//units=[m]
0.450250878464055400D-02;0.767076427091007100D-03;0.266057917766977600D-03;...//position of Sun at epoch
0.361762716560282000D+00;-0.907819721567660000D-01;-0.857149725627511600D-01;...//position of Mercury at epoch
0.612751940835072200D+00;-0.348365369033622200D+00;-0.195278286675943800D+00;...//position of Venus at epoch
0.120517414101384700D+00;-0.925838474769148700D+00;-0.401540226453152200D+00;...//position of Earth-Moon barycenter at epoch
-0.110186077148798200D+00;-0.132759945030298300D+01;-0.605889140484291400D+00;...//position of Mars at epoch
-0.537970676855393700D+01;-0.830481326563397800D+00;-0.224828874426565500D+00;...//position of Jupiter at epoch
0.789439068290953100D+01;0.459647805517127300D+01;0.155869584283190000D+01;...//position of Saturn at epoch
-0.182654022538723600D+02;-0.116195541867587000D+01;-0.250106057721338000D+00;...//position of Uranus at epoch
-0.160550357802333700D+02;-0.239421915598547100D+02;-0.940015796880239500D+01;...//position of Neptune at epoch
-0.304833137671838400D+02;-0.872405556841050100D+00;0.891157617249955100D+01];   //position of Pluto at epoch

//Velocity units (au/day) (1 AU = 149597870700/(24*60*60) m/s)
ObjVel=149597870700/24/60/60*[...//units=[m/s]
-0.351749536075523200D-06;0.517762640983340500D-05;0.222910217891202900D-05;...//vel. of Sun at epoch
0.336749397200575900D-02;0.248945205576834300D-01;0.129463004097040900D-01;...//vel. of Mercury at epoch
0.109520684235282300D-01;0.156176842678676800D-01;0.633110570297786400D-02;...//vel. of Venus at epoch
0.168112683097837900D-01;0.174830923073434400D-02;0.758202897383129100D-03;...//vel. of Earth-Moon barycenter at epoch
0.144816530570475700D-01;0.242463076836468600D-03;-0.281520727924338800D-03;...//vel. of Mars at epoch
0.109201259423733800D-02;-0.651811661280738400D-02;-0.282078276229867900D-02;...//vel. of Jupiter at epoch
-0.321755651650091600D-02;0.433581034174662600D-02;0.192864631686015500D-02;...//vel. of Saturn at epoch
0.221190391015614700D-03;-0.376247500810884400D-02;-0.165101502742995000D-02;...//vel. of Uranus at epoch
0.264276984798005500D-02;-0.149831255054097800D-02;-0.679041960802913300D-03;...//vel. of Neptune at epoch
0.322207373497780800D-03;-0.314357639364532900D-02;-0.107794975959731300D-02];  //vel. of Pluto at epoch


y0 = [ObjPos;ObjVel];
result = ode(y0, t0, t, f);
[m,n] = size(result);
minP=min(result(2:Nbodies*3,:))
maxP=max(result(2:Nbodies*3,:))
result = [t(1:n);result];

//configurar ventana
figure(1)
clf(1)
fig=gcf();
fig.figure_position=[100,100]
fig.axes_size=[800,600]
ax=gca();
ax.auto_scale="off"
ax.data_bounds=[minP,minP,minP;maxP,maxP,maxP]

//Body settings

for i=1:Nbodies
  param3d1(result((i-1)*3+2,1),result((i-1)*3+3,1),result((i-1)*3+4),1)
  SurfHandles(i)=gce()
  SurfHandles(i).mark_mode="on"
  SurfHandles(i).mark_size_unit="point"
  SurfHandles(i).mark_size=ObjRadius(i)
  SurfHandles(i).mark_foreground=ObjColor(i)
end 

for i=1:Nbodies
  param3d1(result((i-1)*3+2,:)',result((i-1)*3+3,:)',result((i-1)*3+4,:)');
end 

//Animation Loop
for i=1:length(result(1,:))
//  sleep(1)
  for j=1:Nbodies
    SurfHandles(j).data(1)=result((j-1)*3+2,i);
    SurfHandles(j).data(2)=result((j-1)*3+3,i);
    SurfHandles(j).data(3)=result((j-1)*3+4,i);
  end
end