A particle physics simulation program that models ideal gas behavior. Uses object-oriented programming to simulate particles with velocity, acceleration, and forces. Features graphical visualization using BGI graphics, particle interactions, statistical analysis, and real-time rendering.
program Gas;
uses
Crt, Graph;
const
{Константы - шаг времени, максимальное число частиц и измерений}
dt = 5e-3;
maxCount= 1000;
maxStat = 100;
type
{Тип данных - вектор и указатель на вектор (т.к. функции возвращают только простые типы)}
pTVector = ^TVector;
TVector = record
x,y : Double;
end;
{Объект - частица}
pTParticleObj = ^TParticleObj;
TParticleObj = object
acc,
vel,
cur : TVector;
radius,
sforce : Double;
underF : Boolean;
isFree : Boolean;
number : Integer;
constructor Init(_cur, _vel: TVector; _radius, _sforce: Double; _number : Integer);
procedure IntegrateFirst(dt: Double);
procedure IntegrateSecond(dt: Double);
procedure Draw;virtual;
destructor Done;
private
function DistTo(That: pTParticleObj): Double;
function GetForce(That: pTParticleObj): pTVector;
end;
{Объект - частица}
pTPartWlkmObj = ^TPartWlkmObj;
TPartWlkmObj = object(TParticleObj)
ch : Char;
constructor Init(_cur, _vel: TVector; _radius, _sforce: Double; _number : Integer; _ch: Char);
procedure Draw;virtual;
end;
{Тип данных - запись, информация о контейнере}
TContainer = record
part : array[1..maxCount] of pTParticleObj;
firstcount, count: Integer;
radius: Integer;
sforce: Integer;
minX, minY, maxX, maxY : Double;
minVel, maxVel: Integer;
enableGraph: Boolean;
needDispose: Boolean;
skipFrames: Integer;
end;
{Тип данных - запись, информация о статистике}
TStatistic = record
statData : array[1..maxStat] of Integer;
fromStat, toStat : Integer;
curStat: Integer;
statCount: Integer;
statGaps: Integer;
end;
var
Xmax,
Ymax,
ViewXmax,
ViewYmax : Integer;
TextHi, TextWi : Integer;
GraphDriver, GraphMode : Integer;
CurPage : Integer;
EnableStat, StatChanged : Boolean;
container: TContainer;
statistic: TStatistic;
zero: TVector;
{Конструктор объекта - частицы}
constructor TParticleObj.Init(_cur, _vel: TVector; _radius, _sforce: Double; _number : Integer);
begin
cur := _cur;
vel := _vel;
radius := _radius;
sforce := _sforce;
number := _number;
end;
{метод - расстояние до другой частицы}
function TParticleObj.DistTo(That : pTParticleObj): Double;
var res: Double;
begin
with That^.cur do
res := Sqrt((cur.x - x) * (cur.x - x) + (cur.y - y) * (cur.y - y));
if res = 0 then DistTo := 1
else DistTo := res;
end;
{метод, возвращающий указатель на вектор силы}
function TParticleObj.GetForce(That : pTParticleObj): pTVector;
var resS: Double;
resV: TVector;
presV: pTVector;
dist: Double;
begin
dist := DistTo(That);
resS := (That^.radius + radius) - dist;
if resS > 0 then resS := resS * sforce
else resS := 0;
with That^.cur do
begin
resV.x := (cur.x - x) / dist * resS;
resV.y := (cur.y - y) / dist * resS;
end;
GetForce := @resV;
end;
{Интегрирование уравнения движения: вычисление сила действующий на частицу, скорости}
procedure TParticleObj.IntegrateFirst(dt: Double);
var i: Integer;
force: TVector;
curWall: TVector;
wall: pTParticleObj;
begin
acc.x := 0;
acc.y := 0;
curWall.x := container.minX;
curWall.y := cur.y;
New(wall, Init(curWall, zero, 0, 0, 0));
force := GetForce(wall)^;
Dispose(wall);
acc.x := acc.x + force.x;
acc.y := acc.y + force.y;
curWall.x := container.maxX;
curWall.y := cur.y;
New(wall, Init(curWall, zero, 0, 0, 0));
force := GetForce(wall)^;
Dispose(wall);
acc.x := acc.x + force.x;
acc.y := acc.y + force.y;
curWall.x := cur.x;
curWall.y := container.minY;
New(wall, Init(curWall, zero, 0, 0, 0));
force := GetForce(wall)^;
Dispose(wall);
acc.x := acc.x + force.x;
acc.y := acc.y + force.y;
curWall.x := cur.x;
curWall.y := container.maxY;
New(wall, Init(curWall, zero, 0, 0, 0));
force := GetForce(wall)^;
Dispose(wall);
acc.x := acc.x + force.x;
acc.y := acc.y + force.y;
for i:=1 to container.count do
if i<>number then
begin
force := GetForce(container.part[i])^;
acc.x := acc.x + force.x;
acc.y := acc.y + force.y;
end;
if (acc.x <> 0) or (acc.y <> 0) then
underF := True
else underF := False;
EnableStat := EnableStat and not underF; {если на частицу не действуют силы, то можно измерять скорость}
vel.x := vel.x + acc.x * dt;
vel.y := vel.y + acc.y * dt;
end;
{Интегрирование уравнения движения: вычисление новой координаты}
procedure TParticleObj.IntegrateSecond(dt: Double);
begin
cur.x := cur.x + vel.x * dt;
cur.y := cur.y + vel.y * dt;
with container, cur do
isFree := (x > maxX) or (x < minX) or (y > maxY) or (y < minY);
end;
{рисуем частицу}
procedure TParticleObj.Draw;
begin
if underF then
begin
SetColor(4); {красным, если на нее действует сила}
SetFillStyle(SolidFill, 4);
FillEllipse(Trunc(cur.x), Trunc(cur.y), Trunc(radius), Trunc(radius));
end
else
begin
SetColor(9); {иначе - синим}
SetFillStyle(SolidFill, 9);
FillEllipse(Trunc(cur.x), Trunc(cur.y), Trunc(radius), Trunc(radius));
end;
end;
{не нужен}
destructor TParticleObj.Done;
begin
end;
constructor TPartWlkmObj.Init(_cur, _vel: TVector; _radius, _sforce: Double; _number : Integer; _ch: Char);
begin
TParticleObj.Init(_cur, _vel, _radius, _sforce, _number);
ch := _ch;
end;
procedure TPartWlkmObj.Draw;
begin
SetColor(12);
SetFillStyle(SolidFill, 4);
OutTextXY(Trunc(cur.x - TextWi/2), Trunc(cur.y - TextHi/2), ch);
end;
{из целого в строку}
function IntToStr(I: Real): String;
var
S: string[16];
begin
Str(I:8:0, S);
IntToStr := S;
end;
{меняем страницу}
procedure ChangePage;
begin
if CurPage = 0 then
begin
CurPage := 1;
SetActivePage(1);
SetVisualPage(0);
end
else
begin
CurPage := 0;
SetActivePage(0);
SetVisualPage(1);
end;
end;
{рисование на всем доступном экране}
procedure FullPort;
begin
end;
{рама}
procedure DrawFrame;
var i: integer;
begin
ChangePage;
SetViewPort(0, 0, Xmax, Ymax, ClipOn);
SetFillStyle(SolidFill, 0);
Bar(0, 0, Xmax, Ymax);
SetColor(7);
SetFillStyle(SolidFill, 1);
Bar(0, 0, Xmax, Ymax - (TextHi + 4) - 2);
Rectangle(0, 0, Xmax, Ymax - (TextHi + 4) - 2);
SetFillStyle(SolidFill, 8);
Bar(0, Ymax - (TextHi + 4), Xmax, Ymax);
Rectangle(0, Ymax - (TextHi + 4), Xmax, Ymax);
ChangePage;
SetViewPort(0, 0, Xmax, Ymax, ClipOn);
SetFillStyle(SolidFill, 0);
Bar(0, 0, Xmax, Ymax);
SetColor(7);
SetFillStyle(SolidFill, 1);
Bar(0, 0, Xmax, Ymax - (TextHi + 4) - 2);
Rectangle(0, 0, Xmax, Ymax - (TextHi + 4) - 2);
SetFillStyle(SolidFill, 8);
Bar(0, Ymax - (TextHi + 4), Xmax, Ymax);
Rectangle(0, Ymax - (TextHi + 4), Xmax, Ymax);
{draw on top field only}
SetViewPort(1, 1, Xmax - 1, Ymax - (TextHi + 4) - 3, ClipOn);
end;
{строка состояния}
procedure MessageFrame(Msg:string);
begin
ChangePage;
SetViewPort(1, Ymax - (TextHi + 4) + 1, Xmax - 1, Ymax - 1, ClipOn);
SetFillStyle(SolidFill, 8);
Bar(1, 1, Xmax - 1, TextHi + 2);
SetColor(15);
SetTextStyle(DefaultFont, HorizDir, 1);
SetTextJustify(CenterText, TopText);
OutTextXY(Xmax div 2, 2, Msg);
ChangePage;
SetViewPort(1, Ymax - (TextHi + 4) + 1, Xmax - 1, Ymax - 1, ClipOn);
SetFillStyle(SolidFill, 8);
Bar(1, 1, Xmax - 1, TextHi + 2);
SetColor(15);
SetTextStyle(DefaultFont, HorizDir, 1);
SetTextJustify(CenterText, TopText);
OutTextXY(Xmax div 2, 2, Msg);
{draw on top field only}
SetViewPort(1, 1, Xmax - 1, Ymax - (TextHi + 4) - 3, ClipOn);
end;
{добавляем измерение}
procedure AddMeasure;
var i,k: Integer;
cvel: Double;
gapW: Double;
begin
if (statistic.curStat >= statistic.fromStat) and
(statistic.curStat <= statistic.toStat) then
begin
gapW := 3 * container.maxVel / statistic.statGaps; {промежуток скоростей}
with statistic do
for i:=1 to statGaps do
for k:=1 to container.count do
with container.part[k]^ do
begin
if isFree then continue;
cvel := Sqrt(vel.x * vel.x + vel.y * vel.y);
if (cvel > gapW * (i - 1)) and (cvel < gapW * i) then
Inc(statData[i]);
end;
end;
Inc(statistic.curStat);
MessageFrame(Concat('Current measure #', IntToStr(statistic.curStat)));
end;
{рисуем гистограмму}
procedure DrawStatPlot;
var i, maxStat: Integer;
gapW : Double;
begin
with statistic, container do
begin
gapW := 3 * maxVel / statGaps;
maxStat := 0;
for i:=1 to statGaps do
if maxStat < statData[i] then maxStat := statData[i];
ChangePage;
SetFillStyle(SolidFill, 1);
Bar(0, 0, Xmax, Ymax - (TextHi + 4) - 2);
SetColor(9);
Graph.Rectangle(Trunc(minX), Trunc(minY), Trunc(maxX), Trunc(maxY));
SetColor(11);
for i:=1 to statGaps do
begin
Graph.Rectangle(
Trunc(minX + (maxX - minX) / statGaps * (i - 1)),
Trunc(maxY),
Trunc(minX + (maxX - minX) / statGaps * i ),
Trunc(maxY - (maxY - minY) / maxStat * statData[i]));
OutTextXY(
Trunc(minX + (maxX - minX) / statGaps * (i - 1)),
Trunc(maxY + (i mod 3) * TextHi + 3),
IntToStr(GapW*i));
end;
for i:=1 to Trunc((maxY-minY)/(TextHi+3)) do
OutTextXY(
Trunc(minX - 4*(TextWi + 2)),
Trunc(maxY - (TextHi+3)*i),
IntToStr(Trunc(maxStat/(maxY-minY)*(TextHi+3)*i)));
ChangePage;
end;
end;
{рисуем молекулы}
procedure DrawPoints;
var
i : Integer;
begin
SetFillStyle(SolidFill, 1);
Bar(0, 0, Xmax, Ymax - (TextHi + 4) - 2);
SetColor(9);
with container do
Graph.Rectangle(Trunc(minX), Trunc(minY), Trunc(maxX), Trunc(maxY));
for i:=1 to container.count do
container.part[i]^.Draw;
ChangePage;
end;
{procedure WaitOrExit(Msg: String);forward;}
procedure Remove(That : pTParticleObj);
var i: Integer;
begin
MessageFrame('Particle #' + IntToStr(That^.number) + ' has gone far...');
Dec(container.count);
for i:= That^.number to container.count do
begin
container.part[i] := container.part[i + 1];
container.part[i]^.number := i;
end
container.part[container.firstcount] := That;
{ Dispose(That);}
end;
{цикл метода Эйлера, также проверяет, можно ли производить измерение}
procedure MainCycle(dt: Double);
var i,k : Integer;
begin
EnableStat := True;
for i:=1 to container.count do
begin
container.part[i]^.IntegrateFirst(dt);
end;
for i:=1 to container.count do
begin
container.part[i]^.IntegrateSecond(dt);
end;
for i:=1 to container.count do
begin
if container.part[i]^.isFree then Remove(container.part[i]);
end;
StatChanged := StatChanged or not EnableStat;
if EnableStat and StatChanged then
begin
AddMeasure;
StatChanged := False;
end;
end;
{ввод целого параметра}
procedure InputParam(var param: Integer; Msg: String);
var Ch : Char;
Done: Boolean;
StrVal: String;
i: Integer;
begin
Done := False;
Str(param, StrVal);
repeat
MessageFrame(Concat(Msg, ' : ', StrVal));
Ch:=ReadKey;
case Ch of
#13: Done:=True;
'0'..'9' : StrVal := Concat(StrVal, Ch);{только цифры}
#8: Delete(StrVal, Length(StrVal), 1);{enter - конец ввода}
else
MessageFrame('Must be a numeric value!');
repeat until KeyPressed;
ReadKey;
end;
until Done;
Val(StrVal, param, i);
end;
{ввод логического параметра}
procedure InputParamB(var param: Boolean; Msg: String);
var Ch : Char;
Done: Boolean;
StrVal: String;
i: Integer;
begin
Done := False;
if param then StrVal := 'Y'
else StrVal := 'N';
repeat
MessageFrame(Concat(Msg, ' [Y/N]: ', StrVal));
Ch:=ReadKey;
case Ch of
#13: Done:=True;
'y','Y','n','N' : StrVal := Ch;
#8: Delete(StrVal, Length(StrVal), 1);
else
MessageFrame('Must be "y" or "n"!');
repeat until KeyPressed;
ReadKey;
end;
until Done and (Length(StrVal) = 1);
param := (StrVal = 'y') or (StrVal = 'Y');
end;
{создает контейнер с частицами}
procedure MakeContainer;
var i: Integer;
cur, vel: TVector;
absVel, angle : Double;
begin
if container.needDispose then
for i:=1 to container.firstcount do
Dispose(container.part[i]);
repeat
InputParam(container.firstcount, 'Input quantity of particles (1..1000)');
if (container.count < 1) or (container.count > MaxCount) then
begin
MessageFrame('Wrong input!');
repeat until KeyPressed;
ReadKey;
end;
until not ((container.count < 1) or (container.count > MaxCount));
container.count := container.firstcount;
repeat
InputParam(container.radius, 'Input radius of particles (1..20)');
if (container.radius < 1) or (container.radius > 20) then
begin
MessageFrame('Wrong input!');
repeat until KeyPressed;
ReadKey;
end;
until not ((container.radius < 1) or (container.radius > 20));
repeat
InputParam(container.sforce, 'Input spring force of particles (0..100000)');
if (container.sforce > 100000) then
begin
MessageFrame('Wrong input!');
repeat until KeyPressed;
ReadKey;
end;
until not ((container.sforce > 100000));
repeat
InputParam(container.minVel, 'Input minimal velocity of particles');
InputParam(container.maxVel, 'Input maximal velocity of particles');
if container.minVel > container.minVel then
begin
MessageFrame('Wrong input!');
repeat until KeyPressed;
ReadKey;
end;
until not (container.minVel > container.maxVel);
InputParamB(container.enableGraph, 'Enable graphics');
if (container.enableGraph) then
repeat
InputParam(container.skipFrames, 'Amount of skipped frames per calculation (0..50)');
if (container.skipFrames < 0) or (container.skipFrames > 50) then
begin
MessageFrame('Wrong input!');
repeat until KeyPressed;
ReadKey;
end;
until not ((container.skipFrames < 0) or (container.skipFrames > 50));
for i:=1 to container.count do
begin
cur.x := (container.maxX - container.minX) * Random + container.minX;
cur.y := (container.maxY - container.minY) * Random + container.minY;
absVel := container.minVel + (container.maxVel - container.minVel)*Random;
angle := 2 * Pi * Random;
vel.x := absVel * Cos(angle);
vel.y := absVel * Sin(angle);
New(container.part[i], Init(cur, vel, container.radius, container.sforce, i));
end;
container.needDispose := True;
end;
{создает статистику}
procedure MakeStatistic;
var i : Integer;
begin
repeat
InputParam(statistic.statCount, 'Input amount of measurements (<100)');
if statistic.statCount > maxStat then
begin
MessageFrame('Wrong input!');
repeat until KeyPressed;
ReadKey;
end;
until not (statistic.statCount > maxStat);
repeat
InputParam(statistic.fromStat, 'Input first measure of statistics');
InputParam(statistic.toStat, 'Input last measure of statistics');
if (statistic.fromStat > statistic.toStat) or
(statistic.fromStat > statistic.statCount) or
(statistic.toStat > statistic.statCount) then
begin
MessageFrame('Wrong input!');
repeat until KeyPressed;
ReadKey;
end;
until not ((statistic.fromStat > statistic.toStat) or
(statistic.fromStat > statistic.statCount) or
(statistic.toStat > statistic.statCount));
repeat
InputParam(statistic.statGaps, 'Input number of columns on histogram (<100)');
if statistic.statGaps > 100 then
begin
MessageFrame('Wrong input!');
repeat until KeyPressed;
ReadKey;
end;
until not (statistic.statGaps > 100);
statistic.curStat:=0;
for i:=1 to statistic.statCount do statistic.statData[i] := 0;
end;
{ожидает нажатия клавиши и выходит из программы}
procedure WaitOrExit(Msg: String);
var
Ch : char;
WantExit: Boolean;
begin
MessageFrame(Msg);
repeat until KeyPressed;
Ch := ReadKey;
if Ch = #27 then
begin
InputParamB(WantExit, 'Do you really want to exit?');
if WantExit then
begin
CloseGraph;
RestoreCrtMode;
Writeln('Have a nice day!');
end;
end;
end;
{проверяет ввод пользователя}
procedure CheckForUserInput;
var Ch: Char;
begin
if KeyPressed then
begin
Ch := ReadKey;
case Ch of
'n','N': begin
MakeContainer;
MakeStatistic;
end;
's','S': AddMeasure;
else
WaitOrExit('Press any key to continue, Esc quits');
end;
end;
end;
procedure Demo;
const
WelcomeStr : Array[1..4] of String =
('WELCOME!!!',
' ',
'This program illustrates',
'Maxwell distribution for Ideal Gases'
);
var wpart: pTPartWlkmObj;
cur, vel: TVector;
i,k : Integer;
begin
MessageFrame('LOADING...');
begin
container.firstcount := 49;
container.count := 0;
container.radius := 4;
container.sforce := 10000;
container.minVel := 75;
container.maxVel := 300;
container.minX := ViewXmax/6*1;
container.minY := ViewYmax/6*1;
container.maxX := ViewXmax/6*5;
container.maxY := ViewYmax/6*5;
container.needDispose := True;
container.enableGraph := True;
container.skipFrames := 0;
statistic.fromStat:=0;
statistic.toStat:=100;
statistic.statGaps:=30;
statistic.statCount:=100;
end;
with container do
for i:=1 to 4 do
for k:=1 to Length(WelcomeStr[i]) do
begin
if WelcomeStr[i,k] = ' ' then continue;
cur.x := (maxX + minX - Length(WelcomeStr[i]) * (TextWi + 1)) /2 + k * (TextWi + 1);
cur.y := (maxY + minY - 8 * (TextHi + 2)) / 2 + i * (TextHi + 2);
vel.x := 300 * (0.5 - Random);
vel.y := 300 * (0.5 - Random);
New(wpart, Init(cur, vel, 0, sforce, 0, WelcomeStr[i,k]));
Inc(count);
part[count] := wpart;
end;
with container do
begin
for i:=1 to 100 do MainCycle(dt);
MessageFrame('Displaying Demo!');
for i:=1 to 100 do
begin
if (i mod 4 = 0) then DrawPoints;
MainCycle(-dt);
end;
DrawPoints;
ChangePage;
SetTextStyle(DefaultFont, HorizDir, 1);
SetTextJustify(CenterText, TopText);
SetColor(2);
OutTextXY(
Trunc((maxX + minX) / 2),
Trunc((maxY + minY) / 2 + 3 * (TextHi + 2)),
'Arkhipenko Kirill, group #212.');
SetColor(2);
OutTextXY(
Trunc((maxX + minX) / 2),
Trunc((maxY + minY) / 2 + 4 * (TextHi + 2)),
'Term paper.');
SetColor(10);
OutTextXY(
Trunc((maxX + minX) / 2),
Trunc((maxY + minY) / 2 + 6 * (TextHi + 2)),
'Visit http://www.physfac.nm.ru/ to get more!');
ChangePage;
WaitOrExit('Press Any key to continue or Esc to exit.');
end;
end;
{главная функция}
procedure DoMagic;
var i : Integer;
counter: Integer;
begin
MakeContainer;
MakeStatistic;
MessageFrame('Press "n", "s" or "p", Esc quits, any other key pauses.');
counter := container.skipFrames;
repeat
MainCycle(dt);
if container.enableGraph then
if counter > container.skipFrames then
begin
DrawPoints;
counter := 0;
end
else
counter := counter + 1;
CheckForUserInput;
until statistic.curStat = statistic.statCount;
end;
procedure TestGraphError(GraphErr: integer);
begin
if GraphErr <> grOk then begin
Writeln('Graphics error: ', GraphErrorMsg(GraphErr));
repeat until KeyPressed;
ReadKey;
Halt(1);
end;
end;
{инициализация графики}
procedure Init;
begin
GraphDriver := Detect;
DetectGraph(GraphDriver, GraphMode);
TestGraphError(GraphResult);
case GraphDriver of
CGA : begin
GraphDriver := CGA;
GraphMode := CGAC1;
end;
MCGA : begin
case GraphMode of
MCGAMed, MCGAHi: GraphMode := MCGAC1;
end;
end;
EGA : begin
GraphMode := EGAHi;
end;
EGA64 : begin
GraphMode := EGA64Hi;
end;
PC3270 : begin
GraphDriver := CGA;
GraphMode := CGAC1;
end;
ATT400 : case GraphMode of
ATT400C1,
ATT400C2,
ATT400Med,
ATT400Hi :
begin
GraphMode := ATT400C1;
end;
end;
end;
InitGraph(GraphDriver, GraphMode, '');
TestGraphError(GraphResult);
SetGraphMode(1);
SetTextStyle(DefaultFont, HorizDir, 1);
SetTextJustify(CenterText, TopText);
Xmax := GetMaxX;
Ymax := GetMaxY;
TextHi := TextHeight('M');
TextWi := TextWidth('M');
ViewXmax := Xmax - 2;
ViewYmax := (Ymax - (TextHi + 4) - 1) - 2;
end;
{}
begin
Init;
DrawFrame;
repeat
Demo;
DoMagic;
DrawStatPlot;
WaitOrExit('Press any key to start again or Esc to exit.');
until False;
end.