A variant of the gas simulation with PMOR (possibly parameter modification or optimization) features.
program Gas;
uses
Crt, Graph;
const
Memory = 100;
Windows = 4;
dt = 5e-3;
maxCount= 1000;
maxStat = 100;
type
ResolutionPreference = (Lower, Higher);
ColorList = array [1..Windows] of integer;
pTVector = ^TVector;
TVector = record
x,y : Double;
end;
pTParticleObj = ^TParticleObj;
TParticleObj = object
acc,
vel,
cur : TVector;
radius,
sforce : Double;
underF : Boolean;
number : Integer;
constructor Init(_cur, _vel: TVector; _radius, _sforce: Double; _number : Integer);
procedure IntegrateFirst(dt: Double);
procedure IntegrateSecond(dt: Double);
procedure Draw;
destructor Done;
private
function DistTo(That: pTParticleObj): Double;
function GetForce(That: pTParticleObj): pTVector;
end;
TContainer = record
part : array[1..maxCount] of pTParticleObj;
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;
Ch: char;
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 := 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;
end;
procedure TParticleObj.Draw;
begin
if underF then
begin
SetColor(4);
SetFillStyle(SolidFill, 4);
FillEllipse(Trunc(cur.x), Trunc(cur.y), Trunc(radius/2), Trunc(radius/2));
end
else
begin
SetColor(9);
SetFillStyle(SolidFill, 9);
FillEllipse(Trunc(cur.x), Trunc(cur.y), Trunc(radius/2), Trunc(radius/2));
end;
end;
destructor TParticleObj.Done;
begin
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
SetViewPort(0, 0, Xmax, Ymax, ClipOn);
end;
procedure Frame;
begin
SetViewPort(0, 0, Xmax, Ymax-(TextHeight('M')+4)-1,ClipOn);
SetColor(MaxColors);
Rectangle(0, 0, Xmax-1, (Ymax-(TextHeight('M')+4)-1)-1);
SetViewPort(1, 1, Xmax-2, (Ymax-(TextHeight('M')+4)-1)-2,ClipOn);
end;
procedure MessageFrame(Msg:string);
begin
ChangePage;
FullPort;
SetColor(MaxColors);
SetTextStyle(DefaultFont, HorizDir, 1);
SetTextJustify(CenterText, TopText);
SetLineStyle(SolidLn, 0, NormWidth);
SetFillStyle(EmptyFill, 0);
Bar(0, Ymax-(TextHeight('M')+4), Xmax, Ymax);
Rectangle(0, Ymax-(TextHeight('M')+4), Xmax, Ymax);
OutTextXY(Xmax div 2, Ymax-(TextHeight('M')+2), Msg);
Frame;
ChangePage;
FullPort;
SetColor(MaxColors);
SetTextStyle(DefaultFont, HorizDir, 1);
SetTextJustify(CenterText, TopText);
SetLineStyle(SolidLn, 0, NormWidth);
SetFillStyle(EmptyFill, 0);
Bar(0, Ymax-(TextHeight('M')+4), Xmax, Ymax);
Rectangle(0, Ymax-(TextHeight('M')+4), Xmax, Ymax);
OutTextXY(Xmax div 2, Ymax-(TextHeight('M')+2), Msg);
Frame;
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
begin
with container.part[i]^.vel do
cvel := Sqrt(x * x + y * 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;
begin
with statistic, container do
begin
maxStat := 0;
for i:=1 to statGaps do
if maxStat < statData[i] then maxStat := statData[i];
ChangePage;
Graph.ClearViewPort;
SetColor(10);
Graph.Rectangle(Trunc(minX), Trunc(minY), Trunc(maxX), Trunc(maxY));
SetColor(4);
for i:=1 to statGaps do
Graph.Rectangle(
Trunc(minX + (maxX - minX)/statGaps*(i-1)),
Trunc(maxY),
Trunc(minX + (maxX - minX)/statGaps*i),
Trunc(maxY - (maxY - minY)/maxStat*statData[i]));
ChangePage;
end;
end;
procedure WaitToGo;
var
Ch : char;
begin
MessageFrame('Press any key to continue, Esc quits.');
repeat until KeyPressed;
Ch := ReadKey;
if Ch = #27 then begin
CloseGraph;
Halt(1);
end
else
MessageFrame('Press "n" to change parameters, Esc quits.');
end;
procedure DrawPoints;
var
i : Integer;
begin
ChangePage;
Graph.ClearViewPort;
SetColor(7);
with container do
Graph.Rectangle(Trunc(minX), Trunc(minY), Trunc(maxX), Trunc(maxY));
for i:=1 to container.count do
container.part[i]^.Draw;
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;
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);
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;
begin
if container.needDispose then
for i:=1 to container.count do
Dispose(container.part[i]);
repeat
InputParam(container.count, 'Input numbers 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));
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');
InputParam(container.maxVel, 'Input maximal velocity');
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, 'Number of skipped frame 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;
vel.x := 2 * (container.maxVel) * (0.5 - Random);
vel.y := 2 * (container.maxVel) * (0.5 - Random);
New(container.part[i], Init(cur, vel, container.radius, container.sforce, i));
end;
container.needDispose := True;
MessageFrame('Press "n" to change parameters, Esc quits.');
end;
procedure MakeStatistic;
begin
repeat
InputParam(statistic.statCount, 'Input number 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 then
begin
MessageFrame('Wrong input!');
repeat until KeyPressed;
ReadKey;
end;
until not (statistic.fromStat > statistic.toStat);
repeat
InputParam(statistic.statGaps, 'Input number columns on histogramm (<100)');
if statistic.statGaps > 100 then
begin
MessageFrame('Wrong input!');
repeat until KeyPressed;
ReadKey;
end;
until not (statistic.statGaps > 100);
statistic.curStat:=0;
end;
procedure CheckForUserInput;
begin
if KeyPressed then begin
Ch := ReadKey;
case Ch of
'n','N': begin
MakeContainer;
MakeStatistic;
end;
's','S': AddMeasure;
else
WaitToGo;
end;
end;
end;
procedure DoMagic;
var i : Integer;
counter: Integer;
begin
container.count := 25;
container.radius := 7;
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 := False;
container.enableGraph := True;
container.skipFrames := 0;
MakeContainer;
statistic.fromStat:=0;
statistic.toStat:=100;
statistic.statGaps:=30;
statistic.statCount:=100;
MakeStatistic;
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);
begin
DrawStatPlot;
end;
end;
procedure TestGraphError(GraphErr: integer);
begin
if GraphErr <> grOk then begin
Writeln('Graphics error: ', GraphErrorMsg(GraphErr));
repeat until keypressed;
Ch := 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;
ViewXmax := Xmax - 2;
ViewYmax := (Ymax - (TextHeight('M') + 4) - 1) - 2;
end;
var WantAgain: Boolean;
begin
Init;
Frame;
repeat
DoMagic;
InputParamB(WantAgain, 'Do you want to start again?');
until not WantAgain;
CloseGraph;
RestoreCrtMode;
Writeln('Have a nice day!');
end.