solar_system_array_n_planets, Pascal

solar_system_array_n_planets

uses GraphABC;

const
  n = 6;
  t = 0.01;

type
  velocity = record
    x: array [1..n, 1..3] of real;
    y: array [1..n, 1..3] of real;
  end;
  
  acceleration = record
    x: array [1..n] of real;
    y: array [1..n] of real;
  end;

var
  x: array [1..n] of real;
  y: array [1..n] of real;
  r: array [1..n,1..n] of real;
  m: array [1..n] of real;

var
  v: velocity;
  a: acceleration;
  flag : integer;
 
var i, j : integer;

begin
  flag := 1; 
  // вводим начальные условия для первого тела
  m[1] := 3;
  x[1] := 0;         y[1] := 0;
  v.x[1, 1] := 0.0;   v.y[1, 1] := 0.0;
  v.x[1, 2] := 0.0;   v.y[1, 2] := 0;
  
  v.x[1, 3] := 0.0;   v.y[1, 3] := 0.5 * (v.y[1, 1] + v.y[1, 2]);
  
  // вводим начальные условия для второго тела
  m[2] := 0.01;
  x[2] := 4;          y[2] := 0;
  v.x[2, 1] := 0.0;   v.y[2, 1] := 0.0;
  v.x[2, 2] := 0;     v.y[2, 2] := 1.8;
  
  v.x[2, 3] := 0.0;   v.y[2, 3] := 0.5 * (v.y[2, 1] + v.y[2, 2]);
  
  // вводим начальные условия для третьего тела
  m[3] := 0.000001;
  x[3] := 4.1;          y[3] := 0;
  v.x[3, 1] := 0.0;   v.y[3, 1] := 0.0;
  v.x[3, 2] := 0;     v.y[3, 2] := 1.0;
  
  v.x[3, 3] := 0.0;   v.y[3, 3] := 0.5 * (v.y[3, 1] + v.y[3, 2]);
   
// вводим начальные условия для четвёртого тела
  m[4] := 0.000001;
  x[4] := 4.3;          y[4] := 0;
  v.x[4, 1] := 0.0;   v.y[4, 1] := 0.0;
  v.x[4, 2] := 0;     v.y[4, 2] := 1.3;
  
  v.x[4, 3] := 0.0;   v.y[4, 3] := 0.5 * (v.y[4, 1] + v.y[4, 2]);
 
 // вводим начальные условия для пятого тела
  m[5] := 0.000001;
  x[5] := 4.5;          y[5] := 0;
  v.x[5, 1] := 0.0;   v.y[5, 1] := 0.0;
  v.x[5, 2] := 0;     v.y[5, 2] := 1.6;
  
  v.x[5, 3] := 0.0;   v.y[5, 3] := 0.5 * (v.y[5, 1] + v.y[5, 2]);

// вводим начальные условия для шестого тела
  m[6] := 0.000001;
  x[6] := 4.7;          y[6] := 0;
  v.x[6, 1] := 0.0;   v.y[6, 1] := 0.0;
  v.x[6, 2] := -1;     v.y[6, 2] := 2;
  
  v.x[6, 3] := 0.0;   v.y[6, 3] := 0.5 * (v.y[6, 1] + v.y[6, 2]);

for i := 1 to n do
  begin
    a.x[i] := 0;
    a.y[i] := 0;
  end;
 
 // задаём начальные расстояния между телами
  
  for i := 1 to n do
    for j := 1 to n do
      if i <> j then
      begin
        r[i,j] := -sqrt (sqr(x[i] - x[j]) + sqr(y[i] - y[j]))**3;
        a.x[i] := a.x[i] + m[j] * (x[i]-x[j]) / r[i,j];
        a.y[i] := a.y[i]  + m[j] * (y[i]-y[j]) / r[i,j];
      end;
  
  SetWindowTop (10);
  SetWindowLeft(10);
  SetWindowSize(1000, 1000);
  
  while x[1] < 10000 do
  begin
    for i := 1 to n do
    begin
      
      x[i] := x[i] + t * v.x[i, 3];
      y[i] := y[i] + t * v.y[i, 3];
    
      v.x[i, 1] := v.x[i, 1] + t * a.x[i];
      v.y[i, 1] := v.y[i, 1] + t * a.y[i];
    
      v.x[i, 2] := v.x[i, 2] + (t+t ) * a.x[i];
      v.y[i, 2] := v.y[i, 2] + (t+t ) * a.y[i];
    
      v.x[i, 3] := 0.5 * (v.x[i, 1] + v.x[i, 2] );
      v.y[i, 3] := 0.5 * (v.y[i, 1] + v.y[i, 2] );
     end;
     
     for i := 1 to n do
     begin
        a.x[i] := 0;
        a.y[i] := 0;  
        for j :=  1 to n do
          if i <> j then
          begin
            r[i,j] := -sqrt (sqr(x[i] - x[j]) + sqr(y[i] - y[j]))**3;
            a.x[i] := a.x[i] + m[j] * (x[i]-x[j]) / r[i,j];
            a.y[i] := a.y[i] + m[j] * (y[i]-y[j]) / r[i,j];
          end;
     end;    
     
    ClearWindow;
    LockDrawing;
    
    line(100, 500, 800, 500);
    line(500, 100, 500, 800);
    for i := 1 to n do
    begin
      if m[i] < 0.01 then 
        circle(round(500 + 100 * x[i]), round(500 + 100 * y[i]), round( 2 ) )
      else if m[i] < 0.1 then
        circle(round(500 + 100 * x[i]), round(500 + 100 * y[i]), round( 4 ) )
      else 
        circle(round(500 + 100 * x[i]), round(500 + 100 * y[i]), round( m[i]*10 ) );
       
    end;
    for i := 1 to n - 1 do
      for j := i + 1 to n do
      begin
        line(round(500 + 100 * x[i]), round(500 + 100 * y[i]), round(500 + 100 * x[j]), round(500 + 100 * y[j]));
      end;
    
    // circle(round(500 + 100 * (m[1]*x[1]+m[2]*x[2]+m[3]*x[3])/(m[1]+m[2]+m[3])), round(500 + 100 * (m[1]*y[1]+m[2]*y[2]+m[3]*y[3])/(m[1]+m[2]+m[3])), round( 3 ) );
    
    Redraw;
    if flag = 1 then begin readln; flag := 0; end;
    if random(101)=5 then 
    begin
     // m[1] := 1.05 * m[1];
     // m[2] := 0.95 * m[2];
    end;
  end;
end.


run \| edit \| history \| help	0