Category : Pascal Source Code
Archive   : 3DF.ZIP
Filename : 3DF.PAS

{ 3dF --- Simple viewer for real functions of form F(x,y). }
{ Parse an expression to form a syntax tree. Walk the tree }
{ to evaluate the expression on a grid over the x-y plane. }
{ Draw orthographic projection with hidden lines removed, }
{ view angle, etc. adjustable by user with cursor controls. }
{ Disclaimer: This code is a mess. It was written in a }
{ very short time to prove that expensive hardware was not }
{ required for a specific purpose. Let the user beware. }

{ Compile with TP 5.0. }

Program p3dF;
Uses Crt, Parse;

{ Assembly language assist. According to tprof, most
opportunity for getting faster is in DrawLine,
especially Plot (on my 4.77Mhz 8088 w/o NCP). }
procedure SetVideo(mode: Byte); external;
procedure ClearScreen; external;
procedure SetPallette(pal: Byte); external;
procedure SetBackGround(color: Byte); external;
procedure MovePen(x, y: Word); external;
procedure DrawLine(x, y: Word); external;
procedure InitCloud; external;
procedure UpdateCloud; external;
procedure ClearImage; external;
procedure ShowImage; external;
{$ifdef profile}
procedure initx; external;
procedure inity; external;
procedure initXinc; external;
procedure initYinc; external;
procedure plot; external;
{$endif}
{$L 3df.obj}

const
nGrids = 15; { make me 10 for faster displays of fewer grids }
DevP = 192; { vertical device size in pixels }
StepsInQuad = 18;
IntScale = 32;
ScaleBits = 5;
ScaleBitsM1 = 4; { scale bits minus 1 (accounts for 2/1 aspect of pixels) }
pi = 3.1415926535897932385; { Save fun call time. }
var
eData: array[-nGrids..nGrids, -nGrids..nGrids] of Real;
eMin, eMax: Real;
yData: array[0..3, -nGrids..nGrids, -nGrids..nGrids] of Integer;
scrImg: array[0..199, 0..79] of Byte;
baseXTbl,
baseYTbl,
deltaXiTbl,
deltaYiTbl,
deltaXjTbl, { how much to inc x for step in j }
deltaYjTbl: array[0..StepsInQuad] of Integer; { how much to inc y for step in j }

procedure FillGrid(size, exag: Real; phi: Integer);
var
e, de, i, j: Integer;
sPhi, cPhi,
imageHt, imageWd,
unitsToDev, eScale,
theta, dTheta, sTheta, cTheta, sizeP, scaleFac: Real;
begin
sPhi := Sin(phi*pi/180);
cPhi := Cos(phi*pi/180);
{ 1.5 is really sqrt(2)*fudge to _ensure_ images lie in screen. }
imageHt := sPhi*size*1.5+cPhi*(eMax-eMin)*exag;
imageWd := size*1.5;
if imageHt > imageWd/1.6 then
unitsToDev := DevP/imageHt
else
unitsToDev := DevP*1.6/imageWd;
sizeP := size*unitsToDev;
eScale := exag*unitsToDev*cPhi*IntScale;
for i := -nGrids to nGrids do
for j := -nGrids to nGrids do begin
e := Round((eData[i,j]-eMin)*eScale);
yData[0, i, j] := e;
yData[1, j,-i] := e;
yData[2, -i,-j] := e;
yData[3, -j, i] := e;
end;
theta := 0;
dTheta := (pi/2)/StepsInQuad;
scaleFac := sizeP/(nGrids+nGrids)*IntScale;
for i := 0 to stepsInQuad-1 do begin
sTheta := scaleFac*Sin(theta);
cTheta := scaleFac*Cos(theta);
deltaXiTbl[i] := Round(sTheta);
deltaYiTbl[i] := Round(cTheta*sPhi);
deltaXjTbl[i] := Round(cTheta);
deltaYjTbl[i] := Round(sTheta*sPhi);
baseXtbl[i] := Round((sizeP*0.707*Cos(theta+pi/4)+160)*IntScale);
baseYtbl[i] := Round((DevP-sizeP*0.707*(1-Sin(theta+pi/4))*sPhi)*IntScale);
theta := theta+dTheta;
end;
end;

procedure DrawIt(cx, cy, size: Real);
const
initStep = 4;
initDstep = 1;
initQuad = 1;
initPhi = 20;
initExag = 1;
var
exag: Real;
i, j, phi,
step, quad, dStep, x00, y00, x0, y0, x, y, tx, ty,
px, py0, py, lastPx, lastPy0, dxi, dyi, dxj, dyj: Integer;
lastCutX, lastCutY: array[-nGrids..nGrids] of Integer;
ch: Char;

procedure WritePhi;
begin
GoToXY(19, 25);
Write(phi:2);
end;

procedure WriteExag;
begin
GoToXY(36, 25);
Write(exag:3:1);
end;

begin
DirectVideo := False;
dStep := initDstep;
step := initStep;
quad := initQuad;
phi := initPhi;
exag := initExag;
FillGrid(size, exag, phi);
SetVideo(6); { 640 x 200 }
ClearScreen;
GoToXY(1, 25);
Write(
'eXit '#27#26'rot '#24#25'elev(__) PgUp/Dn exag(___) ',
'x ',cx:0:2, ' y ', cy:0:2, ' sz ', size:0:2,
' f ', eMin:0:2, '/', eMax:0:2);
WritePhi;
WriteExag;
repeat
ClearImage;
InitCloud;
x0 := baseXtbl[step];
x00 := x0;
y0 := baseYtbl[step];
y00 := y0;
dxi := deltaXiTbl[step];
dyi := deltaYiTbl[step];
dxj := deltaXjTbl[step];
dyj := deltaYjTbl[step];
{ draw first cut }
y := y0;
x := x0;
px := x shr ScaleBitsM1;
py0 := y shr ScaleBits;
py := (y-yData[quad, -nGrids,-nGrids]) shr ScaleBits;
MovePen(px, py0);
DrawLine(px, py);
lastPx := px;
lastPy0 := py0;
lastCutX[-nGrids] := px;
lastCutY[-nGrids] := py;
for j := -nGrids+1 to nGrids do begin
x := x - dxj;
y := y - dyj;
px := x shr ScaleBitsM1;
py0 := y shr ScaleBits;
py := (y-yData[quad, -nGrids, j]) shr ScaleBits;
DrawLine(px, py);
DrawLine(px, py0);
MovePen(lastPx, lastPy0);
DrawLine(px, py0);
if j = 0 then
if Odd(quad) then begin
tx := x-dxi*3;
ty := y+dyi*3;
DrawLine(tx shr ScaleBitsM1, ty shr ScaleBits);
if quad = 1 then begin
DrawLine((tx-dxi+dxj) shr ScaleBitsM1,
(ty+dyi+dyj) shr ScaleBits);
MovePen(tx shr ScaleBitsM1, ty shr ScaleBits);
DrawLine((tx-dxi-dxj) shr ScaleBitsM1,
(ty+dyi-dyj) shr ScaleBits);
end
else begin
DrawLine((tx+dxi+dxj) shr ScaleBitsM1,
(ty-dyi+dyj) shr ScaleBits);
MovePen(tx shr ScaleBitsM1, ty shr ScaleBits);
DrawLine((tx+dxi-dxj) shr ScaleBitsM1,
(ty-dyi-dyj) shr ScaleBits);
end;
end;
MovePen(px, py);
lastPx := px;
lastPy0 := py0;
lastCutX[j] := px;
lastCutY[j] := py;
end;
UpdateCloud;
lastPx := x00 shr ScaleBitsM1;
lastPy0 := y00 shr ScaleBits;
for i := -nGrids+1 to nGrids do begin
x0 := x0+dxi;
y0 := y0-dyi;
y := y0;
x := x0;
px := x shr ScaleBitsM1;
py0 := y shr ScaleBits;
py := (y-yData[quad, i,-nGrids]) shr ScaleBits;
MovePen(lastPx, lastPy0);
DrawLine(px, py0);
lastPx := px;
lastPy0 := py0;
if (i = 0) and (quad and 1 = 0) then begin
tx := x+dxj*3;
ty := y+dyj*3;
MovePen(tx shr ScaleBitsM1, ty shr ScaleBits);
if quad = 0 then begin
DrawLine((tx+dxj-dxi) shr ScaleBitsM1,
(ty+dyj+dyi) shr ScaleBits);
MovePen(tx shr ScaleBitsM1, ty shr ScaleBits);
DrawLine((tx+dxj+dxi) shr ScaleBitsM1,
(ty+dyj-dyi) shr ScaleBits);
end
else begin
DrawLine((tx-dxj-dxi) shr ScaleBitsM1,
(ty-dyj+dyi) shr ScaleBits);
MovePen(tx shr ScaleBitsM1, ty shr ScaleBits);
DrawLine((tx-dxj+dxi) shr ScaleBitsM1,
(ty-dyj-dyi) shr ScaleBits);
end;
MovePen(tx shr ScaleBitsM1, ty shr ScaleBits);
DrawLine(px, py0);
end;
DrawLine(px, py);
DrawLine(lastCutX[-nGrids], lastCutY[-nGrids]);
UpdateCloud;
MovePen(px, py);
lastCutX[-nGrids] := px;
lastCutY[-nGrids] := py;
for j := -nGrids+1 to nGrids do begin
x := x - dxj;
y := y - dyj;
px := x shr ScaleBitsM1;
py := (y-yData[quad, i, j]) shr ScaleBits;
DrawLine(px, py);
DrawLine(lastCutX[j], lastCutY[j]);
MovePen(px, py);
lastCutX[j] := px;
lastCutY[j] := py;
UpdateCloud;
end;
end;
if KeyPressed then begin
ch := ReadKey;
if ch = #0 then begin
ch := ReadKey;
case ch of
#77: if dStep > -3 then { right arrow }
dStep := dStep-1;
#75: if dStep < 3 then { left arrow }
dStep := dStep+1;
#73: if exag < 9.5 then begin { pg up }
if exag < 1 then
exag := exag+0.1
else
exag := exag+0.5;
FillGrid(size, exag, phi);
WriteExag;
end;
#81: if exag > 0.1 then begin { pg dn }
if exag > 1 then
exag := exag-0.5
else
exag := exag-0.1;
FillGrid(size, exag, phi);
WriteExag;
end;
#72: if phi < 60 then begin { up arrow }
phi := phi + 5;
FillGrid(size, exag, phi);
WritePhi;
end;
#80: if phi > 5 then begin { dn arrow }
phi := phi - 5;
FillGrid(size, exag, phi);
WritePhi;
end;
#71: dStep := initDstep; { home }
#79: begin { end }
phi := initPhi;
exag := initExag;
FillGrid(size, exag, phi);
WriteExag;
WritePhi;
end;
end;
end;
end
else ch := ' ';
while KeyPressed do ch := ReadKey; { purge keyboard buffer }
step := step+dStep;
if step >= stepsInQuad then begin
step := step mod stepsInQuad;
quad := quad+1;
if quad > 3 then
quad := 0;
end;
if step < 0 then begin
step := stepsInQuad+step mod stepsInQuad;
quad := quad-1;
if quad < 0 then
quad := 3;
end;
ShowImage;
until ch in [#27, 'X', 'x'];
SetVideo(CO80);
end;

procedure ComputeFun(var expr: String; cx, cy, size: Real);
var
tr: Tree;
i, j, n: Integer;
x, y, halfSize, e: Real;
begin
tr := ParseExpr(expr);
halfSize := size/2;
eMax := -1e30; eMin := 1e30;
n := (nGrids + nGrids + 1);
n := n*n;
Write('Computing F (');
for j := -nGrids to nGrids do begin
Write(n:4, ' left).'#8#8#8#8#8#8#8#8#8#8#8);
x := cX + j/nGrids*halfSize;
for i := -nGrids to nGrids do begin
y := cY + i/nGrids*halfSize;
e := Eval(tr, x, y);
if e < eMin then eMin := e;
if e > eMax then eMax := e;
eData[i,j] := e;
end;
Dec(n, nGrids+nGrids+1);
end;
WriteLn(n:4);
end;

procedure Usage(msg: String);
const
nLines= 4;
verbiage: array[1..nLines] of string = (
'Usage: 3df [[options] expression]'#13#10+
' options:'#13#10+
' /x -- x of domain center (default 0)'#13#10+
' /y -- y of domain center (default 0)'#13#10+
' /s -- x and y size of domain (default 1)'#13#10+
' /? -- This help.'#13#10,
' expression: Optionally "quoted" function of x and y.'#13#10+
' Operators by precedence:'#13#10+
' a^b - a to b power, b >= 0 or int < 0 right associative'#13#10+
' -a - negate a'#13#10+
' a*b - a times b a/b - a divided by b left associative'#13#10,
' a+b - a plus b a-b - a minus b left associative'#13#10+
' Functions:'#13#10+
' Sin(a), Cos(a), Ln(a), Exp(a), Sqrt(a), Atan(a), Abs(a),'#13#10+
' Min(a,b), Max(a,b). Case is ignored.'#13#10+
'Example:'#13#10+
' 3df /x1 /y1.2 /s2 min(2,sqrt(abs(1 - x^2-y^2)))'#13#10#10,
'This is freeware. I ask only for credit and bug reports.'#13#10+
' MAJ Gene Ressler ([email protected])'#13#10+
' 124 Pine Tree Rd, Ithaca, NY 14580'#13#10);
var
i: Integer;
begin
WriteLn(msg+'.');
for i := 1 to nLines do Write(verbiage[i]);
Halt;
end;

function GetRealParam(name: String; var i: Integer): Real;
var
x: Real;
code: Integer;
param: String;
begin
param := Copy(ParamStr(i), 3, 255);
if Length(param) > 0 then
Val(param, x, code)
else begin
inc(i); if i > ParamCount then Usage('Missing '+name);
Val(ParamStr(i), x, code);
end;
if code <> 0 then Usage('Bad '+name);
GetRealParam := x;
end;

procedure GetParams(var expr: String; var cx, cy, size: Real);
var
i, n, code: Integer;
param: String[2];
begin
n := ParamCount;

{ Demo case. }
if n = 0 then begin
expr := '0.25*cos(28*x^2 + 20*y^2) / (9*x^2 + 6*y^2 + 1)';
Writeln('/? for help. Running demo F: ', expr);
cx := 0.3;
cy := 0.2;
size := 1;
Exit;
end;

expr := '';
cx := 0;
cy := 0;
size := 1;

i := 1;
repeat
param := Copy(ParamStr(i), 1, 2);
if (param[1] = '/') then
if Length(param) > 1 then
case param[2] of
'?': Usage('Help..');
's': begin
size := GetRealParam('size', i);
if size <= 0 then Usage('Non-positive size');
end;
'x': cx := GetRealParam('x', i);
'y': cy := GetRealParam('y', i);
else Usage('Unknown switch '+param);
end
else Usage('Missing switch character')
else begin
repeat
if Length(expr) > 0 then expr := expr+' ';
expr := expr + ParamStr(i);
inc(i);
until i > n;
if expr[1] = '"' then begin
expr := Copy(expr, 2, 255);
if expr[Length(expr)] = '"' then
dec(expr[0])
else
Usage('Missing "');
end;
Exit; { Normal exit. }
end;
inc(i);
until i > n;
Usage('Missing expression');
end;

procedure Main;
var
cx, cy, size: Real;
expr: String;
begin
WriteLn('3dF v0.1 MAJ Gene Ressler April 91');
GetParams(expr, cx, cy, size);
ComputeFun(expr, cx, cy, size);
DrawIt(cx, cy, size);
end;

begin
Main;
end.