solar sysmtem 6 (with random mass m1), Pascal

solar sysmtem 6 (with random mass m1)

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] of real;
  m: array [1..n] of real;

var
  v: velocity;
  a: acceleration;
 
begin
   
  // вводим начальные условия для первого тела
  m[1] := 3;
  x[1] := -4;
  y[1] := 0;
  v.x[1, 1] := 0.0;
  v.y[1, 1] := 0.0;
  v.x[1, 2] := 0;
  v.y[1, 2] := -1;
  
  v.x[1, 3] := 0.0;
  v.y[1, 3] := 0.5 * (v.y[1, 1] + v.y[1, 2]);
  
  // вводим начальные условия для второго тела
  m[2] := 3;
  x[2] := 2;
  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;
  
  v.x[2, 3] := 0.0;
  v.y[2, 3] := 0.5 * (v.y[2, 1] + v.y[2, 2]);
  
  r[1] := x[1] - x[2];
  r[2] := x[2] - x[1];
  
  
  SetWindowSize(800, 800);
  
  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];
    a.y[1] := m[2] * (y[1]-y[2]) / r[1];
    r[1] := sqr(x[1] - x[2]) + sqr(y[1] - y[2]);
    r[1] := sqrt(r[1]);
    r[1] := -r[1] * r[1] * r[1];
    
    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];
    a.y[2] := m[1] * (y[2]-y[1]) / r[2];
    r[2] := sqr(x[1] - x[2]) + sqr(y[1] - y[2]);
    r[2] := sqrt(r[2]);
    r[2] := -r[2] * r[2] * r[2];
   
    ClearWindow;
    LockDrawing;
    
    line(0, 400, 800, 400);
    line(400, 0, 400, 800);
    circle(round(400 + 100 * x[1]), round(400 + 100 * y[1]), round(m[1] * 10 ) );
    circle(round(400 + 100 * x[2]), round(400 + 100 * y[2]), round(m[2] * 10 ) );
    line(round(400 + 100 * x[1]), round(400 + 100 * y[1]), round(400 + 100 * x[2]), round(400 + 100 * y[2]));
    circle(round(400 + 100 * (m[1]*x[1]+m[2]*x[2])/(m[1]+m[2])), round(400 + 100 * (m[1]*y[1]+m[2]*y[2])/(m[1]+m[2])), round( 3 ) );
    
    Redraw;
    if random(101)=5 then m[1] := 1.05 * m[1];
  end;
end.


run \| edit \| history \| help	0