solar system 7 (update), Pascal

solar system 7 (update)

uses GraphABC;

const
  n = 4;
  tmax = 128000;
  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;
 
begin
   
  // вводим начальные условия для первого тела
  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.00001;
  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]);
  
  // задаём начальные расстояния между телами
  r[1,2] := x[1] - x[2];  r[2,1] := x[2] - x[1];
  r[2,3] := x[2] - x[3];  r[3,2] := x[3] - x[2];
  r[1,3] := x[1] - x[3];  r[3,1] := x[3] - x[1];

a.x[1] := m[2] * (x[1]-x[2]) / r[1,2] + m[3] * (x[1]-x[3]) / r[1,3];
    a.y[1] := m[2] * (y[1]-y[2]) / r[1,2] + m[3] * (y[1]-y[3]) / r[1,3];
 
    a.x[2] := m[1] * (x[2]-x[1]) / r[2,1] + m[3] * (x[2]-x[3]) / r[2,3];
    a.y[2] := m[1] * (y[2]-y[1]) / r[2,1] + m[3] * (y[2]-y[3]) / r[2,3];
    
    a.x[3] := m[1] * (x[3]-x[1]) / r[3,1] + m[2] * (x[3]-x[2]) / r[3,2];
    a.y[3] := m[1] * (y[3]-y[1]) / r[3,1] + m[2] * (y[3]-y[2]) / r[3,2];
 
 
  SetWindowTop (10);
  SetWindowLeft(10);
  SetWindowSize(1000, 1000);
  
  while x[1] < 10000 do
  begin
    x[1] := x[1] + t * v.x[1, 3];
    y[1] := y[1] + t * v.y[1, 3];
    
    v.x[1, 1] := v.x[1, 1] + t * a.x[1];
    v.y[1, 1] := v.y[1, 1] + t * a.y[1];
    
    v.x[1, 2] := v.x[1, 2] + (t+t ) * a.x[1];
    v.y[1, 2] := v.y[1, 2] + (t+t ) * a.y[1];
    
    v.x[1, 3] := 0.5 * (v.x[1, 1] + v.x[1, 2] );
    v.y[1, 3] := 0.5 * (v.y[1, 1] + v.y[1, 2] );
    
    a.x[1] := m[2] * (x[1]-x[2]) / r[1,2] + m[3] * (x[1]-x[3]) / r[1,3];
    a.y[1] := m[2] * (y[1]-y[2]) / r[1,2] + m[3] * (y[1]-y[3]) / r[1,3];
  
    r[1,2] := sqr(x[1] - x[2]) + sqr(y[1] - y[2]);
    r[1,2] := sqrt(r[1,2]);
    r[1,2] := -r[1,2] * r[1,2] * r[1,2];
    
    r[1,3] := sqr(x[1] - x[3]) + sqr(y[1] - y[3]);
    r[1,3] := sqrt(r[1,3]);
    r[1,3] := -r[1,3] * r[1,3] * r[1,3];
 
 
    x[2] := x[2] + t * v.x[2, 3];
    y[2] := y[2] + t * v.y[2, 3];
    
    v.x[2, 1] := v.x[2, 1] + t * a.x[2];
    v.y[2, 1] := v.y[2, 1] + t * a.y[2];
    
    v.x[2, 2] := v.x[2, 2] + (t+t ) * a.x[2];
    v.y[2, 2] := v.y[2, 2] + (t+t ) * a.y[2];
    
    v.x[2, 3] := 0.5 * (v.x[2, 1] + v.x[2, 2] );
    v.y[2, 3] := 0.5 * (v.y[2, 1] + v.y[2, 2] );
    
    a.x[2] := m[1] * (x[2]-x[1]) / r[2,1] + m[3] * (x[2]-x[3]) / r[2,3];
    a.y[2] := m[1] * (y[2]-y[1]) / r[2,1] + m[3] * (y[2]-y[3]) / r[2,3];
    
    r[2,1] := sqr(x[2] - x[1]) + sqr(y[2] - y[1]);
    r[2,1] := sqrt(r[2,1]);
    r[2,1] := -r[2,1] * r[2,1] * r[2,1];
   
    r[2,3] := sqr(x[2] - x[3]) + sqr(y[2] - y[3]);
    r[2,3] := sqrt(r[2,3]);
    r[2,3] := -r[2,3] * r[2,3] * r[2,3];

x[3] := x[3] + t * v.x[3, 3];
    y[3] := y[3] + t * v.y[3, 3];
    
    v.x[3, 1] := v.x[3, 1] + t * a.x[3];
    v.y[3, 1] := v.y[3, 1] + t * a.y[3];
    
    v.x[3, 2] := v.x[3, 2] + (t+t ) * a.x[3];
    v.y[3, 2] := v.y[3, 2] + (t+t ) * a.y[3];
    
    v.x[3, 3] := 0.5 * (v.x[3, 1] + v.x[3, 2] );
    v.y[3, 3] := 0.5 * (v.y[3, 1] + v.y[3, 2] );
    
    a.x[3] := m[1] * (x[3]-x[1]) / r[3,1] + m[2] * (x[3]-x[2]) / r[3,2];
    a.y[3] := m[1] * (y[3]-y[1]) / r[3,1] + m[2] * (y[3]-y[2]) / r[3,2];
  
    r[3,1] := sqr(x[3] - x[1]) + sqr(y[3] - y[1]);
    r[3,1] := sqrt(r[3,1]);
    r[3,1] := -r[3,1] * r[3,1] * r[3,1];
    
    r[3,2] := sqr(x[3] - x[2]) + sqr(y[3] - y[2]);
    r[3,2] := sqrt(r[3,2]);
    r[3,2] := -r[3,2] * r[3,2] * r[3,2];
 
    ClearWindow;
    LockDrawing;
    
    line(100, 500, 800, 500);
    line(500, 100, 500, 800);
    circle(round(500 + 100 * x[1]), round(500 + 100 * y[1]), round(m[1] * 10 ) );
    circle(round(500 + 100 * x[2]), round(500 + 100 * y[2]), round({m[2] *} 4 ) );
    circle(round(500 + 100 * x[3]), round(500 + 100 * y[3]), round({m[3] *} 2 ) );
    line(round(500 + 100 * x[1]), round(500 + 100 * y[1]), round(500 + 100 * x[2]), round(500 + 100 * y[2]));
    line(round(500 + 100 * x[2]), round(500 + 100 * y[2]), round(500 + 100 * x[3]), round(500 + 100 * y[3]));
    line(round(500 + 100 * x[1]), round(500 + 100 * y[1]), round(500 + 100 * x[3]), round(500 + 100 * y[3]));
    
    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 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