Язык программирования паскаль

Предмет вычислительная физика. Задача: Смешивание газов, имевших разные температуры и давления. Определить температуру и давление в различных частях объема, в который помещена система одинаковых частиц. Весь объем поделен на две части, в которых помещено разное количество частиц, температура тоже различная. В некоторый момент времени перегородка убирается, газы перемешиваются, во всем объеме устанавливается общая температура и давление.

Также к задаче приложен файл, который как-то должен помочь с решением этой задачи.

Program LennardJones;

{$I-}

uses math;

type

tarr = array[1..100] of real; //массивы для числа частиц n<=100

var

ff:text; //файл ввода начальных данных

x,y,vx,vy,ax,ay: tarr; //координаты, скорости, ускорения

i,j,k,n,tk: longint;

dt, xnew,ynew, r,f, pe, potential,Lx,Ly, rf1,rf2,ke,p: real;

dx,dy, vxc, vyc,vmax: double;

procedure separation(var dx,dy: double); //правило ближайшего образа

begin

if abs(dx)>0.5*Lx then dx:=dx-sign(dx)*Lx;

if abs(dy)>0.5*Ly then dy:=dy-sign(dy)*Ly;

end;

procedure force(r: real; var f:real); //вычисляется сила

var

ri,ri3,ri6,g: real;

begin

// if abs(r)<0.8 then r:=0.8;

if abs(r)>2.5 then

begin

f:=0;

potential:=0;

end

else

begin

ri:=1/r;

ri3:=ri*ri*ri;

ri6:=ri3*ri3;

g:=24*ri*ri6*(2*ri6-1);

f:=g/r;

potential:=4*ri6*(ri6-1); //определяется потенциал

end;

end;

procedure Periodic(var xtemp,ytemp: real); //Периодические граничные условия

begin

if xtemp>Lx then xtemp:=xtemp-Lx;

if xtemp<0 then xtemp:=xtemp+Lx;

if ytemp>Ly then ytemp:=ytemp-Ly;

if ytemp<0 then ytemp:=ytemp+Ly;

end;

Procedure Acceler(var x,y,ax,ay:tarr); // вычисляется ускорение

var i,j: integer;

begin

for i:=1 to n do

begin

ax[i]:=0;

ay[i]:=0;

end;

for i:=1 to n-1 do

for j:=(i+1) to n do

begin

dx:=x[i]-x[j];

dy:=y[i]-y[j];

separation(dx,dy);

r:=sqrt(dx*dx+dy*dy);

force(r,f);

ax[i]:=ax[i]+f*dx;

ay[i]:=ay[i]+f*dy;

ax[j]:=ax[j]-f*dx;

ay[j]:=ay[j]-f*dy;

pe:=pe+potential; //потенциальная энергия

end;

end;

procedure Verlet(var x,y,vx,vy,ax,ay: tarr;dt: real);

var i,j: integer;

xnew,ynew: real;

begin

for i:=1 to n do

begin

xnew:= x[i]+vx[i]*dt+0.5*ax[i]*sqr(dt);

ynew:= y[i]+vy[i]*dt+0.5*ay[i]*sqr(dt);

Periodic(xnew,ynew);

x[i]:=xnew;

y[i]:=ynew;

end;

for i:=1 to n do

begin

vx[i]:=vx[i]+0.5*ax[i]*dt;

vy[i]:=vy[i]+0.5*ay[i]*dt;

end;

Acceler(x,y,ax,ay); //определение нового ускорения

for i:=1 to n do

begin

vx[i]:=vx[i]+0.5*ax[i]*dt;

vy[i]:=vy[i]+0.5*ay[i]*dt;

ke:=ke+(sqr(vx[i])+sqr(vy[i]))/2;

end;

end;

//основной блок программы

begin

assign(ff,'init.txt');

reset(ff);

randomize;

readln(ff, Lx,Ly);

readln(ff,dt);

Readln(ff,k);

Readln(ff,tk);

readln(ff,n);

readln(ff,vmax);

close(ff);

for i:=1 to n do

begin

x[i]:=(i-0.5)*Lx/N;

y[i]:=(i-0.5)*Ly/N;

// x[i]:= Lx*random;

// y[i]:=Ly*random;

rf1:=random;

rf2:=0;

while rf2=0 do rf2:=random;

vx[i]:=vmax*(cos(2*pi*rf1)*sqrt(-2*LN(rf2)));

vy[i]:=vmax*(sin(2*pi*rf1)*sqrt(-2*LN(rf2)));

ax[i]:=0;

ay[i]:=0;

end;

vxc:=0;

vyc:=0;

for i:=1 to n do

begin

vxc:=vxc+vx[i]; // полная скорость (импульс) центра масс

vyc:=vyc+vy[i];

end;

vxc:=vxc/n;

vyc:=vyc/n;

for i:= 1 to n do

begin

vx[i]:=vx[i]-vxc;

vy[i]:=vy[i]-vyc;

end;

for i:=1 to k do

begin

ke:=0;

pe:=0;

for j:=1 to tk do

Verlet(x,y,vx,vy,ax,ay,dt);

ke:=ke/tk;

ke:=ke/n;

pe:=pe/tk;

pe:=pe/n;

p:=pe+ke;

writeln(ke,' ',pe,' ',p);

end;

end.

// поставить счётчик