unit zplot;

Interface

uses
Graph,zMath; (* NOTE ZMATH IS FURTHER DOWN THIS FILE *)
Type
Dimension = (Xaxis,Yaxis);
InBlock = array[1..8192] of Single ;

const
Margin_L: integer = 85;
Margin_B: integer = 25;

(* Draw semilog Axis *)
Procedure AxisLog(Dir : Dimension;
Mmin,Mmax :Single ;
Plot_L,Plot_R,Plot_T,Plot_B: integer;
AxisLabel : string);

(* Draw a LINEAR axis *)
Procedure AxisLin(Dir : Dimension;
Mmin,Mmax :Single ;
Plot_L,Plot_R,Plot_T,Plot_B: integer;
AxisLabel : string);

(* Plot data on a Linear/Linear graph *)
Procedure PlotLin(VAR Data: inblock; (* data to be plotted *)
N: integer; (* number of data points *)
var Ymin, Ymax: Single ; (* info for drawing axis *)
Plot_L,Plot_R,Plot_T,Plot_B: integer);

(* Plot data on a Log/Linear graph *)
Procedure PlotSemi(Data : InBlock;
dX : Single;
N : Integer;
Plot_L,Plot_R,Plot_T,Plot_B : integer);

Procedure EndGraph;

Implementation

Procedure EndGraph;
Begin
CloseGraph
End;

(* Draw semilog Axis *)
Procedure AxisLog(Dir : Dimension;
Mmin,Mmax :Single ;
Plot_L,Plot_R,Plot_T,Plot_B: integer;
AxisLabel : string);
(* Draw semilog Axis *)

(* This proceedure will draw a semilog axis to fit the points of a
function begining at Mmin and ending at Mmax, according to
plot specifications. *)

(* Plot_L : integer; left edge in viewport
Plot_R : integer; right edge in viewport
Plot_T : integer; top edge in viewport
Plot_B : integer; Bottom edge in viewport
Mmin : Single ; smallest axis point, must be positive
Mmax : Single ; largest axis point, also >=0
Dir : enumerated On which dimension to generate the axis *)

var
LastTick : Single ; (* Right-most tick on the plot *)
Tick : Single ; (* Plot tick marks *)
DeltaTick : Single ; (* space between tick marks *)
DeltaTick15 : Single ; (* Sense new decade *)
Scale : Single ; (* scaler used in mapping to Plot dim *)
LnMmax : Single ; (* Log of last data point *)
TicMap : integer;(* Tick mapped to Plot *)
TickLabel : string[9]; (* label at tick mark*)

Begin

plot_l := plot_l + Margin_L; (* leave room for label *)
plot_b := plot_b - Margin_B; (* leave room for label *)
(* Computing the nearest tick mark for the last point of the plot.
Tick marks are at 1 2 3 4 5 6 7 8 9 10 times powers of 10.
Last tick mark is computed as follows:
shift decimal point between first and second significant digit
take the integer part. *)

(* Get the nearest power of 10 less than Mmax, initial DeltaTick *)
DeltaTick := Exp10(Floor(Log10(Mmax)));
(* Divide into Mmax to get nearest integer < *)
LastTick := Int(Mmax/DeltaTick);
If LastTick = 1.0 then
Begin (* If last tick at 1 then must be a solid line *)
SetLineStyle(SolidLn,0,NormWidth);
LastTick := LastTick*DeltaTick; (* restore power of 10 *)
DeltaTick := 0.1*DeltaTick; (* DeltaTick is for next decade *)
End
else
Begin
SetLineStyle(DottedLn,0,NormWidth);
LastTick := LastTick*DeltaTick; (* restore power of 10 *)
end;
LnMmax := Ln(Mmax);(* part of map to plot, save the calc in loop *)

(* When scaling the axis points to the plot dimensions, the
following is required *)

If Dir = Xaxis then
Scale := (Plot_R - Plot_L)/(Ln(Mmax/Mmin))
else
Scale := (Plot_B - Plot_T)/(Ln(Mmax/Mmin));

(* Generate tic marks from LastTick to Mmin *)
Tick := LastTick;
(* DeltaTick15 is 1.5 times DeltaTick, used to determine if tick
has gone to the next decade *)
DeltaTick15 := 1.5*DeltaTick;
While Tick >= Mmin Do (* Generate tick marks from LastTick to Mmin *)
Begin
If Dir = Xaxis then
Begin
(* Tick point mapped to Plot dimensions *)
TicMap := Plot_R-Round((LnMmax-Ln(Tick))*Scale);

(* Draw line at Tick mark
Assumes graphics has been initialized.
Note the conversion from real to integer above *)

Line(TicMap,Plot_T,TicMap,Plot_B);
Str(Tick:10:2,TickLabel);
OutTextXY(plot_l - Margin_L,ticmap,TickLabel);
End
Else (* Y axis plot *)
Begin
(* Tick point mapped to Plot dimensions *)
TicMap := Plot_T+Round((LnMmax-Ln(Tick))*Scale);

(* Draw line at Tick mark
Assumes graphics has been initialized.
Note the conversion from real to integer above *)

Line(Plot_L,TicMap,Plot_R,TicMap);
Str(Tick:10:2,TickLabel);
OutTextXY(plot_l - Margin_L,ticmap,TickLabel);
End;
Tick := Tick - DeltaTick; (* Decrement tick mark *)
If Tick < DeltaTick15 then (* In next decade ? *)
Begin (* Yes *)
DeltaTick := 0.1*DeltaTick; (* New DeltaTick *)
DeltaTick15 := 1.5*DeltaTick; (* New DeltaTick15 *)
SetLineStyle(SolidLn,0,NormWidth); (* solid line *)
End
else
SetLineStyle(DottedLn,0,NormWidth);(* dotted line *)
End;
(* Now draw at the min and max points with a thick linetype *)
SetLineStyle(SolidLn,0,ThickWidth); (* thick solid line *)
If Dir = Xaxis then
Begin
Line(Plot_L,Plot_T,Plot_L,Plot_B);
Line(Plot_R,Plot_T,Plot_R,Plot_B);
End
else (* Y axis *)
Begin
Line(Plot_L,Plot_T,Plot_R,Plot_T);
Line(Plot_L,Plot_B,Plot_R,Plot_B);
End;

(* Now add axis labels *)
If Dir = Xaxis then
Begin
OutTextXY(Plot_L,Plot_B + 15,AxisLabel);
End
else (* Y axis *)
Begin
OutTextXY(0,0,AxisLabel);
End;
end; (* AxisLog *)

Procedure AxisLin(Dir : Dimension;
Mmin,Mmax : Single ;
Plot_L,Plot_R,Plot_T,Plot_B : integer;
AxisLabel : string);
(* Draw a LINEAR axis *)

(* This proceedure will draw a linear axis to fit the points of a
function begining at Mmin and ending at Mmax, according to
plot specifications. *)

(* Plot_L : integer; left edge in viewport
Plot_R : integer; right edge in viewport
Plot_T : integer; top edge in viewport
Plot_B : integer; Bottom edge in viewport
Mmin : Single ; smallest axis point
Mmax : Single ; largest axis point
Dir : enumerated On which dimension to generate the axis *)

var
FirstTick : Single ; (* Left most tick mark *)
LastTick : Single ; (* Right-most tick mark *)
Tick : Single ; (* Plot tick marks *)
DeltaTick : Single ; (* space between tick marks *)
Scale : Single ; (* scaler used in mapping to Plot dim *)
TicMap : integer;(* Tick mapped to Plot *)
TickLabel : string[8];(* label at tick mark *)
Begin

plot_l := plot_l + Margin_L; (* leave room for label *)
plot_b := plot_b - Margin_B; (* leave room for label *)
(* Computing the nearest tick mark for the last point of the plot.
Tick marks are at 0 1 2 3 4 5 6 7 8 9
Last tick mark is computed as follows:
Find x, the power of 10 of the span of the plot.
Tick marks are at the units of this power *)

(* Get the nearest power of 10 of the span, initial DeltaTick *)
DeltaTick := Exp10(Floor(Log10(Mmax-Mmin)));

SetLineStyle(SolidLn,0,NormWidth); (* Linear plot tic marks use solid lines *)
If ((Mmax-Mmin)/DeltaTick) < 1.5 then (* assure at least 2 ticks *)
DeltaTick := 0.1*DeltaTick; (* DeltaTick is finer resolution *)

LastTick := Int(Mmax/DeltaTick)*DeltaTick;

(* When scaling the axis points to the plot dimensions, the
following is required *)

If Dir = Xaxis then
Scale := (Plot_R - Plot_L)/(Mmax-Mmin)
else
Scale := (Plot_B - Plot_T)/(Mmax-Mmin);

(* Generate tic marks from LastTick to Mmin *)
Tick := LastTick;
While Tick >= Mmin Do
Begin
If Dir = Xaxis then
Begin
(* Tick point mapped to Plot dimensions *)
TicMap := Plot_R-Round((Mmax-Tick)*Scale);

(* Draw line at Tick mark and label it
Assumes graphics has been initialized.
Note the conversion from real to integer above *)

Line(TicMap,Plot_T,TicMap,Plot_B);
Str(Tick:10:2,TickLabel);
OutTextXY(ticmap - Margin_B,plot_b + 5, TickLabel);
End
Else (* Y axis plot *)
Begin
(* Tick point mapped to Plot dimensions *)
TicMap := Plot_T+Round((Mmax-Tick)*Scale);

(* Draw line at Tick mark and label it
Assumes graphics has been initialized.
Note the conversion from real to integer above *)

Line(Plot_L,TicMap,Plot_R,TicMap);
Str(Tick:10:2,TickLabel);
OutTextXY(plot_l - Margin_L,ticmap,TickLabel);
End;
Tick := Tick - DeltaTick; (* Decrement tick mark *)
End; (* while *)

(* Now draw at the min and max points with a thick linetype *)
SetLineStyle(SolidLn,0,ThickWidth); (* thick solid line *)
If Dir = Xaxis then
Begin
Line(Plot_L,Plot_T,Plot_L,Plot_B);
Line(Plot_R,Plot_T,Plot_R,Plot_B);
End
else (* Y axis *)
Begin
Line(Plot_L,Plot_T,Plot_R,Plot_T);
Line(Plot_L,Plot_B,Plot_R,Plot_B);
End;

(* Now add axis labels *)
If Dir = Xaxis then
Begin
OutTextXY(Plot_L, Plot_B + 15, AxisLabel);
End
else (* Y axis *)
Begin
OutTextXY(Plot_L - Margin_L, Plot_T - 15, AxisLabel);
End;
End; (* AxisLin *)

Procedure PlotLin(VAR Data: inblock; (* data to be plotted *)
N: integer; (* number of data points *)
var Ymin, Ymax: Single ; (* return to calling program *)
(* info for drawing axis *)
Plot_L,Plot_R,Plot_T,Plot_B: integer);
(* Plot data on a linear linear graph *)
(* This procedure will plot data contained in Data on a graph with x and
y axis both linear. *)

var X,
Vertsize, HorzIncr : Single ;
X1,Y1,X2,Y2,I: integer;
Bot : integer;
begin
SetLineStyle(SolidLn,0,NormWidth);
(* artificially set max and min Y and refine later *)
Ymax := Data[1]; (* maximum value on Y axis *)
Ymin := Data[1]; (* minimum value on Y axis *)
For i := 2 to N Do
begin
if Ymax < Data[i] then Ymax := Data[i];
if Ymin > Data[i] then Ymin := Data[i];
end; (* endfor *)
If Ymin = Ymax Then (* avoid divide by zero *)
Begin
Ymax := Ymax + 1.0 ; (* data will be centered when plotted *)
Ymin := Ymin - 1.0 ;
End;
Vertsize := (Plot_B - Plot_T - Margin_B)/(Ymax - Ymin);
HorzIncr := (Plot_R - Plot_L - Margin_L)/(N - 1.0);
Bot := Plot_B - Margin_B;
Y1 := Bot - round((data[1]-Ymin)*Vertsize);
X1 := Margin_L + Plot_L;
X := 1.0 * X1;
For i := 2 to N Do
begin
Y2 := Bot - round((data[i] - Ymin) * Vertsize);
X := X + HorzIncr;
X2 := Round(X);
line(X1, Y1, X2, Y2);
X1 := X2;
Y1 := Y2;
end; {endfor}
end; {plotlin}

Procedure PlotSemi(Data: inblock; dX: Single; N: integer;
Plot_L,Plot_R,Plot_T,Plot_B: integer);
(* Plot data on a log linear graph *)
(* This procedure will plot data contained in Data on a graph with x linear
and y axis log. X axis goes from 0 to dX * N/2. *)

var maxY,minY,maxX,minX: Single;
Vertsize,Horzsize: Single;
X1,Y1,X2,Y2,i: integer;

begin
SetLineStyle(SolidLn,0,NormWidth);
(* artificially set max and min Y and refine later *)
maxY := Data[1]; (* maximum value on Y axis *)
minY := Data[1]; (* minimum value on Y axis *)
maxX := dX * (N/2);
minX := 0;
For i := 2 to N div 2 do
begin
if maxY < Data[i] then maxY := Data[i];
if minY > Data[i] then minY := Data[i];
end; (* endfor *)
If minY = maxY Then (* avoid divide by zero *)
Begin
maxY := maxY * 1.1 ; (* data will be centered when plotted *)
minY := minY * 0.9 ;
End;
Vertsize := (Plot_B - Plot_T - Margin_B)/Ln(maxY / minY);
Horzsize := (Plot_R - Plot_L - Margin_L)/(maxX - minX);
Y1 := (Plot_B - Plot_T - Margin_B) - round((Ln(data[1])-Ln(minY))*Vertsize);
X1 := Margin_L;
For i := 2 to N div 2 do
begin
Y2 := (Plot_B - Plot_T - Margin_B) - round((Ln(data[i]) - Ln(minY)) * Vertsize);
X2 := Margin_L + round((i-1) * dX * Horzsize);
line(X1, Y1, X2, Y2);
X1 := X2;
Y1 := Y2;
end; {endfor}
end; {plotsemi}

Var
Graphmode : Integer;
Graphdriver : integer;

Begin (* Initialize Graphics *)
Graphdriver := Detect;
InitGraph(Graphdriver, Graphmode, '');
End.

{ What follows is unit zMath, which is just a shortened version of Mathunit }

(* written by Andrew L. Hamm, IBM Corporation, Manassas Virginia 22110 *)

UNIT zMath;
INTERFACE
FUNCTION Floor(X : Single) : Single; (* The next lowest integer *)
FUNCTION Ceiling(X : Single) : Single; (* The next highest integer *)
FUNCTION X_intercept(P1x,P1y,P2x,P2y:Single):Single; (* x intrcpt from P1&P2 *)
FUNCTION Y_intercept(P1x,P1y,P2x,P2y:Single):Single; (* y intrcpt from P1&P2 *)
FUNCTION Log10(X : Single) : Single; (* Base 10 log of x *)
FUNCTION Exp10(X : Single) : Single; (* 10 raised to the x power *)
FUNCTION PwrXY(X, Y : Single) : Single; (* raise x to the power y *)
FUNCTION QDBV(Volts : Single) : Single; (* convert from volts to dB *)
FUNCTION QDBW(Watts : Single) : Single; (* convert from watts to dB *)
FUNCTION QWATTS(DB : Single) : Single; (* convert from dB to watts *)
FUNCTION QVOLTS(DB : Single) : Single; (* convert from dB to volts *)
FUNCTION SINH(X : Single) : Single; (* hyperbolic sine *)
FUNCTION COSH(X : Single) : Single; (* hyperbolic cosine *)
FUNCTION TANH(X : Single) : Single; (* hyperbolic tangent *)
FUNCTION ISINH(X : Single) : Single; (* arc hyperbolic sine *)
FUNCTION ICOSH(X : Single) : Single; (* arc hyperbolic cosine *)
FUNCTION ITANH(X : Single) : Single; (* arc hyperbolic tangent *)
FUNCTION ARCSIN(X : Single) : Single; (* arc sine using TP arctan function *)
FUNCTION ARCCOS(X : Single) : Single; (* inverse cosine using TP arctan *)
FUNCTION ATAN2(X, Y : Single) : Single; (* arctan function with quadrant check
X is Single axis value or denominator, Y is imaginary axis value or numerator *)
FUNCTION TAN(X : Single) : Single; (* tangent of X *)
FUNCTION GAUSS (Mean, StdDev : Single) : Single; (* gaussian random number *)
FUNCTION RESIDUE(Radix, Number : Single) : Single; (* remainder of number/radix
*)
FUNCTION MINIMUM(A, B : Single) : Single; (* the minimum of a and b *)
FUNCTION MAXIMUM(A, B : Single) : Single; (* the maximum of a and b *)

IMPLEMENTATION

VAR
Ln10 : Single;

FUNCTION Floor (X : Single) : Single; (* The next lowest integer *)

(* note that INT will not work when X is a negative number *)
BEGIN
IF X >= 0.0 THEN
Floor := Int(X)
ELSE (* Use int for positive x *)
IF Frac(X) = 0.0 THEN
Floor := X
ELSE (* no shift on exact integer *)
Floor := - Int(1 - X) (* Round away from zero *)
END; (* Floor *)

FUNCTION Ceiling (X : Single) : Single; (* The next highest integer *)
BEGIN
IF X <= 0.0 THEN
Ceiling := Int(X)
ELSE (* Use int for negative x *)
IF Frac(X) = 0.0 THEN
Ceiling := X
ELSE (* no shift on exact integer *)
Ceiling := 1 - Int(- X) (* Shift x to negative *)
END; (* Ceiling *)

FUNCTION X_intercept(P1x,P1y,P2x,P2y:Single):Single; (* x intrcpt from P1&P2 *)
var Minv : Single; (* inverse slope *)
BEGIN
Minv := (P2x - P1x)/(P2y - P1y);
X_intercept := P2x - Minv * P2y
END;

FUNCTION Y_intercept(P1x,P1y,P2x,P2y:Single):Single; (* Y intrcpt from P1&P2 *)
var M : Single; (* slope *)
BEGIN
M := (P2y - P1y)/(P2x - P1x);
Y_intercept := P2y - M * P2x
END;

FUNCTION Log10 (X : Single) : Single; (* Base 10 log of x *)
BEGIN (* Ln(10) supplied for speed *)
Log10 := Ln(X) / Ln10 (* Easily derived *)
END; (* Log10 *)

FUNCTION Exp10 (X : Single) : Single; (* 10 raised to the x power *)
BEGIN (* Ln(10) supplied for speed *)
Exp10 := Exp(X * Ln10) (* easily derived *)
END; (* Exp10 *)

FUNCTION PwrXY (X, Y : Single) : Single; (* raise x to the power y *)
BEGIN
IF (Y = 0.0) THEN
PwrXY := 1.0
ELSE
IF (X <= 0.0) AND (Frac(Y) = 0.0) THEN
IF (Frac(Y / 2)) = 0.0 THEN
PwrXY := Exp(Y * Ln(Abs(X)))
ELSE
PwrXY := - Exp(Y * Ln(Abs(X)))
ELSE
PwrXY := Exp(Y * Ln(X));
END; (* PwrXY *)

FUNCTION QDBV (Volts : Single) : Single; (* convert from volts to dB *)
BEGIN
QDBV := 20.0 * Log10(Volts)
END; (* QDBV *)

FUNCTION QDBW (Watts : Single) : Single; (* convert from watts to dB *)
BEGIN
QDBW := 10.0 * Log10(Watts)
END; (* QDBW *)

FUNCTION QWATTS (DB : Single) : Single; (* convert from dB to watts *)
BEGIN
QWATTS := Exp10(DB / 10.0);
END; (* QWATTS *)

FUNCTION QVOLTS (DB : Single) : Single; (* convert from dB to volts *)
BEGIN
QVOLTS := Exp10(DB / 20.0);
END; (* QVOLTS *)

FUNCTION SINH (X : Single) : Single; (* hyperbolic sine *)
BEGIN
SINH := 0.5 * (Exp(X) - Exp(- X))
END; (* SINH *)

FUNCTION COSH (X : Single) : Single; (* hyperbolic cosine *)
BEGIN
COSH := 0.5 * (Exp(X) + Exp(- X))
END; (* COSH *)

FUNCTION TANH (X : Single) : Single; (* hyperbolic tangent *)
BEGIN
X := Exp(2.0 * X);
TANH := (X - 1.0) / (X + 1.0)
END; (* TANH *)

FUNCTION ISINH (X : Single) : Single; (* arc hyperbolic sine *)
BEGIN
ISINH := Ln(Sqrt(1.0 + X * X) + X)
END; (* ISINH *)

FUNCTION ICOSH (X : Single) : Single; (* arc hyperbolic cosine *)
BEGIN
ICOSH := Ln(X + Sqrt(X * X - 1.0))
END; (* ICOSH *)

FUNCTION ITANH (X : Single) : Single; (* arc hyperbolic tangent *)
BEGIN
ITANH := Ln((1.0 + X) / (1.0 - X))
END; (* ITANH *)

FUNCTION ARCSIN (X : Single) : Single; (* arc sine using TP arctan function *)
BEGIN (* answer returned in radians *)
IF X = 1.0 THEN
ARCSIN := Pi / 2.0
ELSE
IF X = - 1.0 THEN
ARCSIN := Pi / - 2.0
ELSE
ARCSIN := Arctan(X / Sqrt(1.0 - Sqr(X)))
END; (* ARCSIN *)

FUNCTION ARCCOS (X : Single) : Single; (* inverse cosine using TP arctan *)
BEGIN (* answer returned in radians *)
IF X = 0.0 THEN
ARCCOS := Pi / 2.0
ELSE
IF X < 0.0 THEN
ARCCOS := Pi - Arctan(Sqrt(1.0 - Sqr(X)) / Abs(X))
ELSE
ARCCOS := Arctan(Sqrt(1.0 - Sqr(X)) / Abs(X))
END; (* ARCCOS *)

FUNCTION ATAN2(X, Y : Single) : Single; (* arctan function with quadrant check
X is Single axis value or denominator, Y is imaginary axis value or numerator *)
BEGIN (* answer returned in radians *)
IF Y <> 0.0 THEN
IF X <> 0.0 THEN (* point not on an axis *)
IF X > 0.0 THEN (* 1st or 4th quadrant use std arctan *)
ATAN2 := Arctan(Y / X)
ELSE
IF Y > 0.0 THEN (* 2nd quadrant *)
ATAN2 := Pi + Arctan(Y / X)
ELSE
ATAN2 := Arctan(Y / X) - Pi (* 3rd quadrant *)
ELSE (* point on the Y axis *)
IF Y >= 0.0 THEN
ATAN2 := Pi / 2.0 (* positive Y axis *)
ELSE
ATAN2 := - Pi / 2.0 (* negative Y axis *)
ELSE (* point on the X axis *)
IF X >= 0.0 THEN
ATAN2 := 0.0 (* positive X axis *)
ELSE
ATAN2 := - Pi (* negative X axis *)
END; (* ATAN2 *)

FUNCTION TAN (X : Single) : Single; (* tangent of X *)
BEGIN
TAN := Sin(X) / Cos(X)
END; (* TAN *)

FUNCTION GAUSS (Mean, StdDev : Single) : Single; (* gaussian random number *)
VAR
I : BYTE; (* index for loop *)
T : Single; (* temporary variable *)

BEGIN (* Based on the central limit theorem *)
T := - 6.0; (* maximum deviation is 6 sigma, remove
the mean first *)
FOR I := 1 TO 12 DO
T := T + Random; (* 12 uniform over 0 to 1 *)
GAUSS := Mean + T * StdDev (* adjust mean and standard deviation *)
END; (* GAUSS *)

FUNCTION RESIDUE (Radix, Number : Single) : Single; (* remainder of
radix/number *)

(* uses APL residue definition *)
BEGIN
RESIDUE := Number - Radix * Floor(Number / Radix)
END; (* RESIDUE *)

FUNCTION MINIMUM (A, B : Single) : Single; (* the minimum of a and b *)
BEGIN
IF A < B THEN
MINIMUM := A
ELSE
MINIMUM := B
END; (* MINIMUM *)

FUNCTION MAXIMUM (A, B : Single) : Single; (* the maximum of a and b *)
BEGIN
IF A < B THEN
MAXIMUM := B
ELSE
MAXIMUM := A
END; (* MAXIMUM *)

BEGIN
Randomize;
Ln10 := Ln(10.0)
END.