← Back to branch summary

~ubuntu-branches/ubuntu/utopic/mricron/utopic

« back to all changes in this revision

Viewing changes to fpmath/uwinplot.pas

  • Committer: Bazaar Package Importer
  • Author(s): Michael Hanke
  • Date: 2010-07-29 22:07:43 UTC
  • Revision ID: james.westby@ubuntu.com-20100729220743-q621ts2zj806gu0n
Tags: upstream-0.20100725.1~dfsg.1
Import upstream version 0.20100725.1~dfsg.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
fpmath/uwinplot.pas
 
1
{ ******************************************************************
 
2
  Plotting routines for Delphi
 
3
  ****************************************************************** }
 
4
 
 
5
unit uwinplot;
 
6
 
 
7
interface
 
8
 
 
9
uses
 
10
  Classes, Graphics,
 
11
  utypes, umath, uround, uinterv, ustrings;
 
12
 
 
13
function InitGraphics(Canvas : TCanvas; Width, Height : Integer) : Boolean;
 
14
 
 
15
{ ------------------------------------------------------------------
 
16
  Enters graphic mode.
 
17
  ------------------------------------------------------------------
 
18
  The parameters Width and Height refer to the object on which the
 
19
  graphic is plotted.
 
20
 
 
21
  Examples:
 
22
 
 
23
  To draw on a TImage object:
 
24
  InitGraph(Image1.Canvas, Image1.Width, Image1.Height)
 
25
 
 
26
  To print the graphic:
 
27
  InitGraph(Printer.Canvas, Printer.PageWidth, Printer.PageHeight)
 
28
  ------------------------------------------------------------------ }
 
29
 
 
30
procedure SetWindow(Canvas         : TCanvas;
 
31
                    X1, X2, Y1, Y2 : Integer;
 
32
                    GraphBorder    : Boolean);
 
33
{ ------------------------------------------------------------------
 
34
  Sets the graphic window
 
35
 
 
36
  X1, X2, Y1, Y2 : Window coordinates in % of maximum
 
37
  GraphBorder    : Flag for drawing the window border
 
38
  ------------------------------------------------------------------ }
 
39
 
 
40
procedure AutoScale(X                     : PVector;
 
41
                    Lb, Ub                : Integer;
 
42
                    Scale                 : TScale;
 
43
                    var XMin, XMax, XStep : Float);
 
44
{ ------------------------------------------------------------------
 
45
  Finds an appropriate scale for plotting the data in X[Lb..Ub]
 
46
  ------------------------------------------------------------------ }
 
47
 
 
48
procedure SetOxScale(Scale                : TScale;
 
49
                     OxMin, OxMax, OxStep : Float);
 
50
{ ------------------------------------------------------------------
 
51
  Sets the scale on the Ox axis
 
52
  ------------------------------------------------------------------ }
 
53
 
 
54
procedure SetOyScale(Scale                : TScale;
 
55
                     OyMin, OyMax, OyStep : Float);
 
56
{ ------------------------------------------------------------------
 
57
  Sets the scale on the Oy axis
 
58
  ------------------------------------------------------------------ }
 
59
 
 
60
procedure SetGraphTitle(Title : String);
 
61
{ ------------------------------------------------------------------
 
62
  Sets the title for the graph
 
63
  ------------------------------------------------------------------ }
 
64
 
 
65
procedure SetOxTitle(Title : String);
 
66
{ ------------------------------------------------------------------
 
67
  Sets the title for the Ox axis
 
68
  ------------------------------------------------------------------ }
 
69
 
 
70
procedure SetOyTitle(Title : String);
 
71
{ ------------------------------------------------------------------
 
72
  Sets the title for the Oy axis
 
73
  ------------------------------------------------------------------ }
 
74
 
 
75
procedure PlotOxAxis(Canvas : TCanvas);
 
76
{ ------------------------------------------------------------------
 
77
  Plots the horizontal axis
 
78
  ------------------------------------------------------------------ }
 
79
 
 
80
procedure PlotOyAxis(Canvas : TCanvas);
 
81
{ ------------------------------------------------------------------
 
82
  Plots the vertical axis
 
83
  ------------------------------------------------------------------ }
 
84
 
 
85
procedure PlotGrid(Canvas : TCanvas; Grid : TGrid);
 
86
{ ------------------------------------------------------------------
 
87
  Plots a grid on the graph
 
88
  ------------------------------------------------------------------ }
 
89
 
 
90
procedure WriteGraphTitle(Canvas : TCanvas);
 
91
{ ------------------------------------------------------------------
 
92
  Writes the title of the graph
 
93
  ------------------------------------------------------------------ }
 
94
 
 
95
procedure SetMaxCurv(NCurv : Byte);
 
96
{ ------------------------------------------------------------------
 
97
  Sets the maximum number of curves and re-initializes their
 
98
  parameters
 
99
  ------------------------------------------------------------------ }
 
100
 
 
101
procedure SetPointParam(CurvIndex, Symbol, Size : Integer;
 
102
                        Color                   : TColor);
 
103
{ ------------------------------------------------------------------
 
104
  Sets the point parameters for curve # CurvIndex
 
105
  ------------------------------------------------------------------ }
 
106
 
 
107
procedure SetLineParam(CurvIndex : Integer;
 
108
                       Style     : TPenStyle;
 
109
                       Width     : Integer;
 
110
                       Color     : TColor);
 
111
{ ------------------------------------------------------------------
 
112
  Sets the line parameters for curve # CurvIndex
 
113
  ------------------------------------------------------------------ }
 
114
 
 
115
procedure SetCurvLegend(CurvIndex : Integer; Legend : String);
 
116
{ ------------------------------------------------------------------
 
117
  Sets the legend for curve # CurvIndex
 
118
  ------------------------------------------------------------------ }
 
119
 
 
120
procedure SetCurvStep(CurvIndex, Step : Integer);
 
121
{ ------------------------------------------------------------------
 
122
  Sets the step for curve # CurvIndex
 
123
  ------------------------------------------------------------------ }
 
124
 
 
125
procedure PlotPoint(Canvas    : TCanvas;
 
126
                    X, Y      : Float;
 
127
                    CurvIndex : Integer);
 
128
{ ------------------------------------------------------------------
 
129
  Plots a point on the screen
 
130
  ------------------------------------------------------------------
 
131
  Input parameters : X, Y      = point coordinates
 
132
                     CurvIndex = index of curve parameters
 
133
                                  (Symbol, Size, Color)
 
134
  ------------------------------------------------------------------ }
 
135
 
 
136
procedure PlotCurve(Canvas    : TCanvas;
 
137
                    X, Y      : PVector;
 
138
                    Lb, Ub,
 
139
                    CurvIndex : Integer);
 
140
{ ------------------------------------------------------------------
 
141
  Plots a curve
 
142
  ------------------------------------------------------------------
 
143
  Input parameters : X, Y      = point coordinates
 
144
                     Lb, Ub    = indices of first and last points
 
145
                     CurvIndex = index of curve parameters
 
146
  ------------------------------------------------------------------ }
 
147
 
 
148
procedure PlotCurveWithErrorBars(Canvas    : TCanvas;
 
149
                                 X, Y, S   : PVector;
 
150
                                 Ns, Lb, Ub,
 
151
                                 CurvIndex : Integer);
 
152
{ ------------------------------------------------------------------
 
153
  Plots a curve with error bars
 
154
  ------------------------------------------------------------------
 
155
  Input parameters : X, Y      = point coordinates
 
156
                     S         = errors
 
157
                     Ns        = number of SD to be plotted
 
158
                     Lb, Ub    = indices of first and last points
 
159
                     CurvIndex = index of curve parameters
 
160
  ------------------------------------------------------------------ }
 
161
 
 
162
procedure PlotFunc(Canvas     : TCanvas;
 
163
                   Func       : TFunc;
 
164
                   Xmin, Xmax : Float;
 
165
                   Npt,
 
166
                   CurvIndex  : Integer);
 
167
{ ------------------------------------------------------------------
 
168
  Plots a function
 
169
  ------------------------------------------------------------------
 
170
  Input parameters:
 
171
    Func       = function to be plotted
 
172
    Xmin, Xmax = abscissae of 1st and last point to plot
 
173
    Npt        = number of points
 
174
    CurvIndex  = index of curve parameters (Width, Style, Color)
 
175
  ------------------------------------------------------------------
 
176
  The function must be programmed as :
 
177
  function Func(X : Float) : Float;
 
178
  ------------------------------------------------------------------ }
 
179
 
 
180
procedure WriteLegend(Canvas    : TCanvas;
 
181
                      NCurv     : Integer;
 
182
                      ShowPoints,
 
183
                      ShowLines : Boolean);
 
184
{ ------------------------------------------------------------------
 
185
  Writes the legends for the plotted curves
 
186
  ------------------------------------------------------------------
 
187
  NCurv      : number of curves (1 to MaxCurv)
 
188
  ShowPoints : for displaying points
 
189
  ShowLines  : for displaying lines
 
190
  ------------------------------------------------------------------ }
 
191
 
 
192
procedure ConRec(Canvas     : TCanvas;
 
193
                 Nx, Ny, Nc : Integer;
 
194
                 X, Y, Z    : PVector;
 
195
                 F          : PMatrix);
 
196
{ ------------------------------------------------------------------
 
197
  Contour plot
 
198
  Adapted from Paul Bourke, Byte, June 1987
 
199
  http://astronomy.swin.edu.au/~pbourke/projection/conrec/
 
200
  ------------------------------------------------------------------
 
201
  Input parameters:
 
202
  Nx, Ny             = number of steps on Ox and Oy
 
203
  Nc                 = number of contour levels
 
204
  X[0..Nx], Y[0..Ny] = point coordinates in pixels
 
205
  Z[0..(Nc - 1)]     = contour levels in increasing order
 
206
  F[0..Nx, 0..Ny]    = function values, such that F[I,J] is the
 
207
                       function value at (X[I], Y[I])
 
208
  ------------------------------------------------------------------ }
 
209
 
 
210
function Xpixel(X : Float) : Integer;
 
211
{ ------------------------------------------------------------------
 
212
  Converts user abscissa X to screen coordinate
 
213
  ------------------------------------------------------------------ }
 
214
 
 
215
function Ypixel(Y : Float) : Integer;
 
216
{ ------------------------------------------------------------------
 
217
  Converts user ordinate Y to screen coordinate
 
218
  ------------------------------------------------------------------ }
 
219
 
 
220
function Xuser(X : Integer) : Float;
 
221
{ ------------------------------------------------------------------
 
222
  Converts screen coordinate X to user abscissa
 
223
  ------------------------------------------------------------------ }
 
224
 
 
225
function Yuser(Y : Integer) : Float;
 
226
{ ------------------------------------------------------------------
 
227
  Converts screen coordinate Y to user ordinate
 
228
  ------------------------------------------------------------------ }
 
229
 
 
230
procedure LeaveGraphics;
 
231
{ ------------------------------------------------------------------
 
232
  Quits graphic mode
 
233
  ------------------------------------------------------------------ }
 
234
 
 
235
 
 
236
implementation
 
237
 
 
238
const
 
239
  MaxSymbol    = 9;          { Max. number of symbols for plotting curves }
 
240
  MaxCurvColor = 9;          { Max. number of colors for curves }
 
241
  Eps          = 1.0E-10;    { Lower limit for an axis label }
 
242
  MaxColor     = $02FFFFFF;  { Max. color value for Delphi }
 
243
 
 
244
const
 
245
  CurvColor : array[1..MaxCurvColor] of TColor =
 
246
    (clRed,
 
247
     clGreen,
 
248
     clBlue,
 
249
     clFuchsia,
 
250
     clAqua,
 
251
     clLime,
 
252
     clNavy,
 
253
     clOlive,
 
254
     clPurple);
 
255
 
 
256
type
 
257
  TAxis = record        { Coordinate axis }
 
258
    Scale : TScale;
 
259
    Min   : Float;
 
260
    Max   : Float;
 
261
    Step  : Float;
 
262
  end;
 
263
 
 
264
  TPointParam = record  { Point parameters                            }
 
265
    Symbol : Integer;   { Symbol: 0: point (.)                        }
 
266
    Size   : Integer;   {         1: solid circle    2: open circle   }
 
267
    Color  : TColor;    {         3: solid square    4: open square   }
 
268
  end;                  {         5: solid triangle  6: open triangle }
 
269
                        {         7: plus (+)        8: multiply (x)  }
 
270
                        {         9: star (* )                        }
 
271
 
 
272
  TLineParam = record   { Line parameters }
 
273
    Style : TPenStyle;
 
274
    Width : Integer;
 
275
    Color : TColor;
 
276
  end;
 
277
 
 
278
  TCurvParam = record          { Curve parameters }
 
279
    PointParam : TPointParam;
 
280
    LineParam  : TLineParam;
 
281
    Legend     : Str30;        { Legend of curve }
 
282
    Step       : Integer;      { Plot 1 point every Step points }
 
283
  end;
 
284
 
 
285
  TCurvParamVector = array[1..255] of TCurvParam;
 
286
  PCurvParamVector = ^TCurvParamVector;
 
287
 
 
288
var
 
289
  Xwin1, Xwin2, Ywin1, Ywin2 : Integer;
 
290
  XminPixel, XmaxPixel       : Integer;
 
291
  YminPixel, YmaxPixel       : Integer;
 
292
  FactX, FactY               : Float;
 
293
  XAxis, YAxis               : TAxis;
 
294
  GraphTitle, XTitle, YTitle : String;
 
295
  MaxCurv                    : Integer;
 
296
  CurvParam                  : PCurvParamVector;
 
297
  GraphWidth, GraphHeight    : Integer;
 
298
  SymbolSizeUnit             : Integer;
 
299
  PenWidth                   : Integer;
 
300
  PenStyle                   : TPenStyle;
 
301
  PenColor, BrushColor       : TColor;
 
302
  BrushStyle                 : TBrushStyle;
 
303
 
 
304
procedure DimCurvParamVector(var CurvParam : PCurvParamVector; Ub : Byte);
 
305
var
 
306
  I : Integer;
 
307
begin
 
308
  { Allocate vector }
 
309
  GetMem(CurvParam, Ub * SizeOf(TCurvParam));
 
310
  if CurvParam = nil then Exit;
 
311
 
 
312
  MaxCurv := Ub;
 
313
 
 
314
  { Initialize curve parameters }
 
315
  for I := 1 to Ub do
 
316
    with CurvParam^[I] do
 
317
      begin
 
318
        PointParam.Symbol := (I - 1) mod MaxSymbol + 1;
 
319
        PointParam.Size := 2;
 
320
        PointParam.Color := CurvColor[(I - 1) mod MaxCurvColor + 1];
 
321
        Legend := 'Curve ' + LTrim(IntStr(I));
 
322
        LineParam.Width := 1;
 
323
        LineParam.Style := psSolid;
 
324
        LineParam.Color := PointParam.Color;
 
325
        Step := 1;
 
326
      end;
 
327
end;
 
328
 
 
329
procedure DelCurvParamVector(var CurvParam : PCurvParamVector; Ub : Byte);
 
330
begin
 
331
  if CurvParam <> nil then
 
332
    begin
 
333
      FreeMem(CurvParam, Ub * SizeOf(TCurvParam));
 
334
      CurvParam := nil;
 
335
      MaxCurv := 0;
 
336
    end;
 
337
end;
 
338
 
 
339
function InitGraphics(Canvas : TCanvas; Width, Height : Integer) : Boolean;
 
340
begin
 
341
  GraphWidth := Width;
 
342
  GraphHeight := Height;
 
343
  SymbolSizeUnit := GraphWidth div 250;
 
344
 
 
345
  XmaxPixel := Width;
 
346
  YmaxPixel := Height;
 
347
 
 
348
  XminPixel := 0;
 
349
  YminPixel := 0;
 
350
 
 
351
  XTitle     := 'X';
 
352
  YTitle     := 'Y';
 
353
  GraphTitle := '';
 
354
 
 
355
  MaxCurv := MaxSymbol;
 
356
  DimCurvParamVector(CurvParam, MaxCurv);
 
357
 
 
358
  InitGraphics := True;
 
359
end;
 
360
 
 
361
procedure SetWindow(Canvas         : TCanvas;
 
362
                    X1, X2, Y1, Y2 : Integer;
 
363
                    GraphBorder    : Boolean);
 
364
var
 
365
  R : Float;
 
366
begin
 
367
  if (X1 >= 0) and (X2 <= 100) and (X1 < X2) then
 
368
    begin
 
369
      Xwin1 := X1;
 
370
      Xwin2 := X2;
 
371
      R := 0.01 * GraphWidth;
 
372
      XminPixel := Round(X1 * R);
 
373
      XmaxPixel := Round(X2 * R);
 
374
    end;
 
375
 
 
376
  if (Y1 >= 0) and (Y2 <= 100) and (Y1 < Y2) then
 
377
    begin
 
378
      Ywin1 := Y1;
 
379
      Ywin2 := Y2;
 
380
      R := 0.01 * GraphHeight;
 
381
      YminPixel := Round(Y1 * R);
 
382
      YmaxPixel := Round(Y2 * R);
 
383
    end;
 
384
 
 
385
  XAxis.Scale := LinScale;
 
386
  XAxis.Min   := 0.0;
 
387
  XAxis.Max   := 1.0;
 
388
  XAxis.Step  := 0.2;
 
389
 
 
390
  YAxis.Scale := LinScale;
 
391
  YAxis.Min   := 0.0;
 
392
  YAxis.Max   := 1.0;
 
393
  YAxis.Step  := 0.2;
 
394
 
 
395
  FactX := (XmaxPixel - XminPixel) / (XAxis.Max - XAxis.Min);
 
396
  FactY := (YmaxPixel - YminPixel) / (YAxis.Max - YAxis.Min);
 
397
 
 
398
  if GraphBorder then
 
399
    Canvas.Rectangle(XminPixel, YminPixel, Succ(XmaxPixel), Succ(YmaxPixel));
 
400
end;
 
401
 
 
402
procedure AutoScale(X : PVector; Lb, Ub : Integer; Scale : TScale;
 
403
                    var XMin, XMax, XStep : Float);
 
404
var
 
405
  I      : Integer;
 
406
  X1, X2 : Float;
 
407
begin
 
408
  { Minimum and maximum of X }
 
409
 
 
410
  X1 := X^[Lb];
 
411
  X2 := X1;
 
412
  for I := Lb to Ub do
 
413
    if X^[I] < X1 then
 
414
      X1 := X^[I]
 
415
    else if X^[I] > X2 then
 
416
      X2 := X^[I];
 
417
 
 
418
  { Linear scale }
 
419
 
 
420
  if Scale = LinScale then
 
421
    begin
 
422
      Interval(X1, X2, 2, 6, XMin, XMax, XStep);
 
423
      Exit;
 
424
    end;
 
425
 
 
426
  { Logarithmic scale }
 
427
 
 
428
  XMin := 1.0E-3;
 
429
  XMax := 1.0E+3;
 
430
  XStep := 10.0;
 
431
 
 
432
  if X1 <= 0.0 then Exit;
 
433
 
 
434
  XMin := Int(Log10(X1)); if X1 < 1.0 then XMin := XMin - 1.0;
 
435
  XMax := Int(Log10(X2)); if X2 > 1.0 then XMax := XMax + 1.0;
 
436
  XMin := Exp10(XMin);
 
437
  XMax := Exp10(XMax);
 
438
end;
 
439
 
 
440
procedure SetOxScale(Scale : TScale; OxMin, OxMax, OxStep : Float);
 
441
begin
 
442
  XAxis.Scale := Scale;
 
443
  case Scale of
 
444
    LinScale :
 
445
      begin
 
446
        if OxMin < OxMax then
 
447
          begin
 
448
            XAxis.Min := OxMin;
 
449
            XAxis.Max := OxMax;
 
450
          end;
 
451
        if OxStep > 0.0 then XAxis.Step := OxStep;
 
452
      end;
 
453
    LogScale :
 
454
      begin
 
455
        if (OxMin > 0.0) and (OxMin < OxMax) then
 
456
          begin
 
457
            XAxis.Min := Floor(Log10(OxMin));
 
458
            XAxis.Max := Ceil(Log10(OxMax));
 
459
          end;
 
460
        XAxis.Step := 1.0;
 
461
      end;
 
462
  end;
 
463
  FactX := (XmaxPixel - XminPixel) / (XAxis.Max - XAxis.Min);
 
464
end;
 
465
 
 
466
procedure SetOyScale(Scale : TScale; OyMin, OyMax, OyStep : Float);
 
467
begin
 
468
  YAxis.Scale := Scale;
 
469
  case Scale of
 
470
    LinScale :
 
471
      begin
 
472
        if OyMin < OyMax then
 
473
          begin
 
474
            YAxis.Min := OyMin;
 
475
            YAxis.Max := OyMax;
 
476
          end;
 
477
        if OyStep > 0.0 then YAxis.Step := OyStep;
 
478
      end;
 
479
    LogScale :
 
480
      begin
 
481
        if (OyMin > 0.0) and (OyMin < OyMax) then
 
482
          begin
 
483
            YAxis.Min := Floor(Log10(OyMin));
 
484
            YAxis.Max := Ceil(Log10(OyMax));
 
485
          end;
 
486
        YAxis.Step := 1.0;
 
487
      end;
 
488
  end;
 
489
  FactY := (YmaxPixel - YminPixel) / (YAxis.Max - YAxis.Min);
 
490
end;
 
491
 
 
492
function Xpixel(X : Float) : Integer;
 
493
var
 
494
  P : Float;
 
495
begin
 
496
  P := FactX * (X - XAxis.Min);
 
497
  if Abs(P) > 30000 then
 
498
    Xpixel := 30000
 
499
  else
 
500
    Xpixel := Round(P) + XminPixel;
 
501
end;
 
502
 
 
503
function Ypixel(Y : Float) : Integer;
 
504
var
 
505
  P : Float;
 
506
begin
 
507
  P := FactY * (YAxis.Max - Y);
 
508
  if Abs(P) > 30000 then
 
509
    Ypixel := 30000
 
510
  else
 
511
    Ypixel := Round(P) + YminPixel;
 
512
end;
 
513
 
 
514
function Xuser(X : Integer) : Float;
 
515
begin
 
516
  Xuser := XAxis.Min + (X - XminPixel) / FactX;
 
517
end;
 
518
 
 
519
function Yuser(Y : Integer) : Float;
 
520
begin
 
521
  Yuser := YAxis.Max - (Y - YminPixel) / FactY;
 
522
end;
 
523
 
 
524
procedure SetGraphTitle(Title : String);
 
525
begin
 
526
  GraphTitle := Title;
 
527
end;
 
528
 
 
529
procedure SetOxTitle(Title : String);
 
530
begin
 
531
  XTitle := Title;
 
532
end;
 
533
 
 
534
procedure SetOyTitle(Title : String);
 
535
begin
 
536
  YTitle := Title;
 
537
end;
 
538
 
 
539
procedure PlotLine(Canvas : TCanvas; X1, Y1, X2, Y2 : Integer);
 
540
begin
 
541
  Canvas.MoveTo(X1, Y1);
 
542
  Canvas.LineTo(X2, Y2);
 
543
end;
 
544
 
 
545
procedure PlotOxAxis(Canvas : TCanvas);
 
546
var
 
547
  W, X, Z          : Float;
 
548
  Wp, Xp, Yp1, Yp2 : Integer;
 
549
  N, I, J          : Integer;
 
550
  TickLength       : Integer;
 
551
  MinorTickLength  : Integer;
 
552
  XLabel           : String;
 
553
begin
 
554
  TickLength := Canvas.TextHeight('M') div 2;
 
555
  MinorTickLength := Round(0.67 * TickLength);
 
556
 
 
557
  PlotLine(Canvas, XminPixel, YmaxPixel, XmaxPixel, YmaxPixel);
 
558
 
 
559
  N := Round((XAxis.Max - XAxis.Min) / XAxis.Step);  { Nb of intervals }
 
560
  X := XAxis.Min;                                    { Tick mark position }
 
561
 
 
562
  Yp1 := YmaxPixel + TickLength;       { End of tick mark }
 
563
  Yp2 := YmaxPixel + MinorTickLength;  { End of minor tick mark (log scale) }
 
564
 
 
565
  for I := 0 to N do  { Label axis }
 
566
    begin
 
567
      if (XAxis.Scale = LinScale) and (Abs(X) < Eps) then X := 0.0;
 
568
 
 
569
      Xp := Xpixel(X);
 
570
 
 
571
      PlotLine(Canvas, Xp, YmaxPixel, Xp, Yp1);
 
572
 
 
573
      if XAxis.Scale = LinScale then Z := X else Z := Exp10(X);
 
574
 
 
575
      XLabel := Trim(FloatStr(Z));
 
576
 
 
577
      Canvas.TextOut(Xp - Canvas.TextWidth(XLabel) div 2, Yp1, XLabel);
 
578
 
 
579
      { Plot minor divisions on logarithmic scale }
 
580
 
 
581
      if (XAxis.Scale = LogScale) and (I < N) then
 
582
        for J := 2 to 9 do
 
583
          begin
 
584
            W := X + Log10(J);
 
585
            Wp := Xpixel(W);
 
586
            PlotLine(Canvas, Wp, YmaxPixel, Wp, Yp2);
 
587
          end;
 
588
 
 
589
      X := X + XAxis.Step;
 
590
    end;
 
591
 
 
592
  { Write axis title }
 
593
 
 
594
  if XTitle <> '' then
 
595
    Canvas.TextOut(XminPixel + (XmaxPixel - XminPixel -
 
596
                   Canvas.TextWidth(XTitle)) div 2,
 
597
                   YmaxPixel + 4 * TickLength, XTitle);
 
598
end;
 
599
 
 
600
procedure PlotOyAxis(Canvas : TCanvas);
 
601
var
 
602
  W, Y, Z          : Float;
 
603
  Wp, Xp1, Xp2, Yp : Integer;
 
604
  N, I, J          : Integer;
 
605
  TickLength       : Integer;
 
606
  MinorTickLength  : Integer;
 
607
  Yoffset          : Integer;
 
608
  YLabel           : String;
 
609
begin
 
610
  TickLength := Canvas.TextWidth('M') div 2;
 
611
  MinorTickLength := Round(0.67 * TickLength);
 
612
  Yoffset := Canvas.TextHeight('M') div 2;
 
613
 
 
614
  PlotLine(Canvas, XminPixel, YminPixel, XminPixel, YmaxPixel);
 
615
 
 
616
  N := Round((YAxis.Max - YAxis.Min) / YAxis.Step);
 
617
  Y := YAxis.Min;
 
618
 
 
619
  Xp1 := XminPixel - TickLength;
 
620
  Xp2 := XminPixel - MinorTickLength;
 
621
 
 
622
  for I := 0 to N do
 
623
    begin
 
624
      if (YAxis.Scale = LinScale) and (Abs(Y) < Eps) then Y := 0.0;
 
625
 
 
626
      Yp := Ypixel(Y);
 
627
 
 
628
      PlotLine(Canvas, XminPixel, Yp, Xp1, Yp);
 
629
 
 
630
      if YAxis.Scale = LinScale then Z := Y else Z := Exp10(Y);
 
631
 
 
632
      YLabel := Trim(FloatStr(Z));
 
633
 
 
634
      Canvas.TextOut(Xp1 - Canvas.TextWidth(YLabel), Yp - Yoffset, YLabel);
 
635
 
 
636
      if (YAxis.Scale = LogScale) and (I < N) then
 
637
        for J := 2 to 9 do
 
638
          begin
 
639
            W := Y + Log10(J);
 
640
            Wp := Ypixel(W);
 
641
            PlotLine(Canvas, XminPixel, Wp, Xp2, Wp);
 
642
          end;
 
643
 
 
644
      Y := Y + YAxis.Step;
 
645
    end;
 
646
 
 
647
  if YTitle <> '' then
 
648
    Canvas.TextOut(XminPixel, YminPixel - 3 * Yoffset, YTitle);
 
649
end;
 
650
 
 
651
procedure PlotGrid(Canvas : TCanvas; Grid : TGrid);
 
652
var
 
653
  X, Y         : Float;
 
654
  I, N, Xp, Yp : Integer;
 
655
 
 
656
var
 
657
  PenStyle : TpenStyle;
 
658
 
 
659
begin
 
660
  PenStyle := Canvas.Pen.Style;
 
661
  Canvas.Pen.Style := psDot;
 
662
 
 
663
  if Grid in [HorizGrid, BothGrid] then  { Horizontal lines }
 
664
    begin
 
665
      N := Round((YAxis.Max - YAxis.Min) / YAxis.Step);  { Nb of intervals }
 
666
      for I := 1 to Pred(N) do
 
667
        begin
 
668
          Y := YAxis.Min + I * YAxis.Step;  { Origin of line }
 
669
          Yp := Ypixel(Y);
 
670
          PlotLine(Canvas, XminPixel, Yp, XmaxPixel, Yp);
 
671
        end;
 
672
    end;
 
673
 
 
674
  if Grid in [VertiGrid, BothGrid] then  { Vertical lines }
 
675
    begin
 
676
      N := Round((XAxis.Max - XAxis.Min) / XAxis.Step);
 
677
      for I := 1 to Pred(N) do
 
678
        begin
 
679
          X := XAxis.Min + I * XAxis.Step;
 
680
          Xp := Xpixel(X);
 
681
          PlotLine(Canvas, Xp, YminPixel, Xp, YmaxPixel);
 
682
        end;
 
683
    end;
 
684
 
 
685
  Canvas.Pen.Style := PenStyle;
 
686
end;
 
687
 
 
688
procedure WriteGraphTitle(Canvas : TCanvas);
 
689
begin
 
690
  if GraphTitle <> '' then
 
691
    with Canvas do
 
692
      TextOut((XminPixel + XmaxPixel - TextWidth(GraphTitle)) div 2,
 
693
               YminPixel - 2 * TextHeight(GraphTitle), GraphTitle);
 
694
end;
 
695
 
 
696
procedure SetMaxCurv(NCurv : Byte);
 
697
begin
 
698
  if NCurv < 1 then Exit;
 
699
  DelCurvParamVector(CurvParam, MaxCurv);
 
700
  MaxCurv := NCurv;
 
701
  DimCurvParamVector(CurvParam, MaxCurv);
 
702
end;
 
703
 
 
704
procedure SetPointParam(CurvIndex, Symbol, Size : Integer;
 
705
                        Color                   : TColor);
 
706
begin
 
707
  if (CurvIndex < 1) or (CurvIndex > MaxCurv) then Exit;
 
708
 
 
709
  if (Symbol >= 0) and (Symbol <= MaxSymbol) then
 
710
    CurvParam^[CurvIndex].PointParam.Symbol := Symbol;
 
711
 
 
712
  if Size > 0 then
 
713
    CurvParam^[CurvIndex].PointParam.Size := Size;
 
714
 
 
715
  if (Color >= 0) and (Color <= MaxColor) then
 
716
    CurvParam^[CurvIndex].PointParam.Color := Color;
 
717
end;
 
718
 
 
719
procedure SetLineParam(CurvIndex : Integer;
 
720
                       Style     : TPenStyle;
 
721
                       Width     : Integer;
 
722
                       Color     : TColor);
 
723
 
 
724
begin
 
725
  if (CurvIndex < 1) or (CurvIndex > MaxCurv) then Exit;
 
726
 
 
727
  CurvParam^[CurvIndex].LineParam.Style := Style;
 
728
 
 
729
  if Width > 0 then
 
730
    CurvParam^[CurvIndex].LineParam.Width := Width;
 
731
 
 
732
  if (Color >= 0) and (Color <= MaxColor) then
 
733
    CurvParam^[CurvIndex].LineParam.Color := Color;
 
734
end;
 
735
 
 
736
procedure SetCurvLegend(CurvIndex : Integer; Legend : String);
 
737
begin
 
738
  if (CurvIndex >= 1) and (CurvIndex <= MaxCurv) then
 
739
    CurvParam^[CurvIndex].Legend := Legend;
 
740
end;
 
741
 
 
742
procedure SetCurvStep(CurvIndex, Step : Integer);
 
743
begin
 
744
  if (CurvIndex >= 1) and (CurvIndex <= MaxCurv) and (Step > 0) then
 
745
    CurvParam^[CurvIndex].Step := Step;
 
746
end;
 
747
 
 
748
function XOutOfBounds(X : Integer) : Boolean;
 
749
{ Checks if an absissa is outside the graphic bounds }
 
750
begin
 
751
  XOutOfBounds := (X < XminPixel) or (X > XmaxPixel);
 
752
end;
 
753
 
 
754
function YOutOfBounds(Y : Integer) : Boolean;
 
755
{ Checks if an ordinate is outside the graphic bounds }
 
756
begin
 
757
  YOutOfBounds := (Y < YminPixel) or (Y > YmaxPixel);
 
758
end;
 
759
 
 
760
function CheckPoint(X, Y       : Float;
 
761
                    var Xp, Yp : Integer) : Boolean;
 
762
{ Computes the pixel coordinates of a point and
 
763
  checks if it is enclosed within the graph limits }
 
764
begin
 
765
  Xp := Xpixel(X);
 
766
  Yp := Ypixel(Y);
 
767
  CheckPoint := not(XOutOfBounds(Xp) or YOutOfBounds(Yp));
 
768
end;
 
769
 
 
770
procedure PlotSymbol(Canvas    : TCanvas;
 
771
                     Xp, Yp,
 
772
                     CurvIndex : Integer);
 
773
var
 
774
  Xp1, Xp2, Yp1, Yp2, Size : Integer;
 
775
begin
 
776
  Size := CurvParam^[CurvIndex].PointParam.Size * SymbolSizeUnit;
 
777
 
 
778
  Xp1 := Xp - Size;
 
779
  Yp1 := Yp - Size;
 
780
  Xp2 := Xp + Size + 1;
 
781
  Yp2 := Yp + Size + 1;
 
782
 
 
783
  with Canvas do
 
784
    case CurvParam^[CurvIndex].PointParam.Symbol of
 
785
      0    : Pixels[Xp, Yp] := Brush.Color;
 
786
      1, 2 : Ellipse(Xp1, Yp1, Xp2, Yp2);          { Circle }
 
787
      3, 4 : Rectangle(Xp1, Yp1, Xp2, Yp2);        { Square }
 
788
      5, 6 : Polygon([Point(Xp1, Yp2 - 1),
 
789
                      Point(Xp2, Yp2 - 1),
 
790
                      Point(Xp, Yp1 - 1)]);        { Triangle }
 
791
      7 : begin                                    { + }
 
792
            PlotLine(Canvas, Xp, Yp1, Xp, Yp2);
 
793
            PlotLine(Canvas, Xp1, Yp, Xp2, Yp);
 
794
          end;
 
795
      8 : begin                                    { x }
 
796
            PlotLine(Canvas, Xp1, Yp1, Xp2, Yp2);
 
797
            PlotLine(Canvas, Xp1, Yp2 - 1, Xp2, Yp1 - 1);
 
798
          end;
 
799
      9 : begin                                    { * }
 
800
            PlotLine(Canvas, Xp, Yp1, Xp, Yp2);
 
801
            PlotLine(Canvas, Xp1, Yp, Xp2, Yp);
 
802
            PlotLine(Canvas, Xp1, Yp1, Xp2, Yp2);
 
803
            PlotLine(Canvas, Xp1, Yp2 - 1, Xp2, Yp1 - 1);
 
804
          end;
 
805
      end;
 
806
end;
 
807
 
 
808
procedure SetGraphSettings(Canvas : TCanvas; CurvIndex : Integer);
 
809
{ Saves the current graphic properties of the Canvas
 
810
  and sets them to the values of curve # CurvIndex }
 
811
begin
 
812
  PenColor   := Canvas.Pen.Color;
 
813
  PenStyle   := Canvas.Pen.Style;
 
814
  PenWidth   := Canvas.Pen.Width;
 
815
  BrushColor := Canvas.Brush.Color;
 
816
  BrushStyle := Canvas.Brush.Style;
 
817
 
 
818
  Canvas.Pen.Color   := CurvParam^[CurvIndex].LineParam.Color;
 
819
  Canvas.Pen.Style   := CurvParam^[CurvIndex].LineParam.Style;
 
820
  Canvas.Pen.Width   := CurvParam^[CurvIndex].LineParam.Width;
 
821
  Canvas.Brush.Color := CurvParam^[CurvIndex].PointParam.Color;
 
822
 
 
823
  if CurvParam^[CurvIndex].PointParam.Symbol in [0, 1, 3, 5] then
 
824
    Canvas.Brush.Style := bsSolid
 
825
  else
 
826
    Canvas.Brush.Style := bsClear;
 
827
end;
 
828
 
 
829
procedure RestoreGraphSettings(Canvas : TCanvas);
 
830
begin
 
831
  Canvas.Pen.Color   := PenColor;
 
832
  Canvas.Pen.Style   := PenStyle;
 
833
  Canvas.Pen.Width   := PenWidth;
 
834
  Canvas.Brush.Color := BrushColor;
 
835
  Canvas.Brush.Style := BrushStyle;
 
836
end;
 
837
 
 
838
procedure PlotPoint(Canvas    : TCanvas;
 
839
                    X, Y      : Float;
 
840
                    CurvIndex : Integer);
 
841
var
 
842
  Xp, Yp : Integer;
 
843
begin
 
844
  if XAxis.Scale = LogScale then X := Log10(X);
 
845
  if YAxis.Scale = LogScale then Y := Log10(Y);
 
846
 
 
847
  if not CheckPoint(X, Y, Xp, Yp) then Exit;
 
848
 
 
849
  SetGraphSettings(Canvas, CurvIndex);
 
850
  PlotSymbol(Canvas, Xp, Yp, CurvIndex);
 
851
  RestoreGraphSettings(Canvas);
 
852
end;
 
853
 
 
854
procedure PlotErrorBar(Canvas       : TCanvas;
 
855
                       Y, S         : Float;
 
856
                       Ns,
 
857
                       Xp, Yp, Size : Integer);
 
858
{ Plots an error bar with the current canvas settings }
 
859
var
 
860
  Delta, Y1 : Float;
 
861
  Yp1       : Integer;
 
862
  PenStyle  : TPenStyle;
 
863
begin
 
864
  Size := Size * SymbolSizeUnit;
 
865
  PenStyle := Canvas.Pen.Style;
 
866
  Canvas.Pen.Style := psSolid;
 
867
 
 
868
  Delta := Ns * S;
 
869
  Y1 := Y - Delta;
 
870
  if YAxis.Scale = LogScale then Y1 := Log10(Y1);
 
871
  Yp1 := Ypixel(Y1);
 
872
 
 
873
  if Yp1 <= YmaxPixel then
 
874
    begin
 
875
      PlotLine(Canvas, Xp - Size, Yp1, Xp + Size + 1, Yp1);
 
876
      PlotLine(Canvas, Xp, Yp, Xp, Yp1);
 
877
    end
 
878
  else
 
879
    PlotLine(Canvas, Xp, Yp, Xp, YmaxPixel);
 
880
 
 
881
  Y1 := Y + Delta;
 
882
  if YAxis.Scale = LogScale then Y1 := Log10(Y1);
 
883
  Yp1 := Ypixel(Y1);
 
884
 
 
885
  if Yp1 >= YminPixel then
 
886
    begin
 
887
      PlotLine(Canvas, Xp - Size, Yp1, Xp + Size + 1, Yp1);
 
888
      PlotLine(Canvas, Xp, Yp, Xp, Yp1);
 
889
    end
 
890
  else
 
891
    PlotLine(Canvas, Xp, Yp, Xp, YminPixel);
 
892
 
 
893
  Canvas.Pen.Style := PenStyle;
 
894
end;
 
895
 
 
896
procedure GenPlotCurve(Canvas    : TCanvas;
 
897
                       X, Y, S   : PVector;
 
898
                       Ns,
 
899
                       Lb, Ub,
 
900
                       CurvIndex : Integer;
 
901
                       ErrorBars : Boolean);
 
902
{ General curve plotting routine }
 
903
var
 
904
  X1, Y1, X2, Y2     : Float;
 
905
  Xp1, Yp1, Xp2, Yp2 : Integer;
 
906
  I                  : Integer;
 
907
  Flag1, Flag2       : Boolean;
 
908
begin
 
909
  SetGraphSettings(Canvas, CurvIndex);
 
910
 
 
911
  { Plot first point }
 
912
 
 
913
  X1 := X^[Lb]; if XAxis.Scale = LogScale then X1 := Log10(X1);
 
914
  Y1 := Y^[Lb]; if YAxis.Scale = LogScale then Y1 := Log10(Y1);
 
915
 
 
916
  Flag1 := CheckPoint(X1, Y1, Xp1, Yp1);
 
917
 
 
918
  if Flag1 then
 
919
    begin
 
920
      PlotSymbol(Canvas, Xp1, Yp1, CurvIndex);
 
921
      if ErrorBars and (S^[Lb] > 0.0) then
 
922
        PlotErrorBar(Canvas, Y^[Lb], S^[Lb], Ns, Xp1, Yp1, CurvIndex);
 
923
    end;
 
924
 
 
925
  { Plot other points and connect them by lines if necessary }
 
926
 
 
927
  I := Lb + CurvParam^[CurvIndex].Step;
 
928
 
 
929
  while I <= Ub do
 
930
    begin
 
931
      X2 := X^[I]; if XAxis.Scale = LogScale then X2 := Log10(X2);
 
932
      Y2 := Y^[I]; if YAxis.Scale = LogScale then Y2 := Log10(Y2);
 
933
 
 
934
      Flag2 := CheckPoint(X2, Y2, Xp2, Yp2);
 
935
 
 
936
      if Flag2 then
 
937
        begin
 
938
          PlotSymbol(Canvas, Xp2, Yp2, CurvIndex);
 
939
          if ErrorBars and (S^[I] > 0.0) then
 
940
            PlotErrorBar(Canvas, Y^[I], S^[I], Ns, Xp2, Yp2, CurvIndex);
 
941
          if (CurvParam^[CurvIndex].LineParam.Style <> psClear) and Flag1 then
 
942
            PlotLine(Canvas, Xp1, Yp1, Xp2, Yp2);
 
943
        end;
 
944
 
 
945
      Xp1 := Xp2;
 
946
      Yp1 := Yp2;
 
947
      Flag1 := Flag2;
 
948
      Inc(I, CurvParam^[CurvIndex].Step);
 
949
    end;
 
950
 
 
951
  RestoreGraphSettings(Canvas);
 
952
end;
 
953
 
 
954
procedure PlotCurve(Canvas    : TCanvas;
 
955
                    X, Y      : PVector;
 
956
                    Lb, Ub    : Integer;
 
957
                    CurvIndex : Integer);
 
958
begin
 
959
  GenPlotCurve(Canvas, X, Y, nil, 0, Lb, Ub, CurvIndex, False);
 
960
end;
 
961
 
 
962
procedure PlotCurveWithErrorBars(Canvas    : TCanvas;
 
963
                                 X, Y, S   : PVector;
 
964
                                 Ns        : Integer;
 
965
                                 Lb, Ub    : Integer;
 
966
                                 CurvIndex : Integer);
 
967
begin
 
968
  GenPlotCurve(Canvas, X, Y, S, Ns, Lb, Ub, CurvIndex, True);
 
969
end;
 
970
 
 
971
procedure PlotFunc(Canvas     : TCanvas;
 
972
                   Func       : TFunc;
 
973
                   Xmin, Xmax : Float;
 
974
                   Npt        : Integer;
 
975
                   CurvIndex  : Integer);
 
976
var
 
977
  X1, Y1, X2, Y2, H  : Float;
 
978
  Xp1, Yp1, Xp2, Yp2 : Integer;
 
979
  Flag1, Flag2       : Boolean;
 
980
  I                  : Integer;
 
981
begin
 
982
  if (Npt < 2) or (CurvParam^[CurvIndex].LineParam.Style = psClear) then Exit;
 
983
 
 
984
  if Xmin >= Xmax then
 
985
    begin
 
986
      Xmin := XAxis.Min;
 
987
      Xmax := XAxis.Max;
 
988
    end;
 
989
 
 
990
  H := (Xmax - Xmin) / Npt;
 
991
 
 
992
  SetGraphSettings(Canvas, CurvIndex);
 
993
 
 
994
  { Check first point }
 
995
  X1 := Xmin;
 
996
  if XAxis.Scale = LinScale then
 
997
    Y1 := Func(X1)
 
998
  else
 
999
    Y1 := Func(Exp10(X1));
 
1000
 
 
1001
  if YAxis.Scale = LogScale then Y1 := Log10(Y1);
 
1002
  Flag1 := CheckPoint(X1, Y1, Xp1, Yp1);
 
1003
 
 
1004
  { Check other points and plot lines if possible }
 
1005
  for I := 1 to Npt do
 
1006
    begin
 
1007
      X2 := X1 + H;
 
1008
      if XAxis.Scale = LinScale then
 
1009
        Y2 := Func(X2)
 
1010
      else
 
1011
        Y2 := Func(Exp10(X2));
 
1012
 
 
1013
      if YAxis.Scale = LogScale then Y2 := Log10(Y2);
 
1014
 
 
1015
      Flag2 := CheckPoint(X2, Y2, Xp2, Yp2);
 
1016
 
 
1017
      if Flag1 and Flag2 then
 
1018
        PlotLine(Canvas, Xp1, Yp1, Xp2, Yp2);
 
1019
 
 
1020
      X1 := X2;
 
1021
      Xp1 := Xp2;
 
1022
      Yp1 := Yp2;
 
1023
      Flag1 := Flag2;
 
1024
    end;
 
1025
 
 
1026
   RestoreGraphSettings(Canvas);
 
1027
end;
 
1028
 
 
1029
procedure WriteLegend(Canvas     : TCanvas;
 
1030
                      NCurv      : Integer;
 
1031
                      ShowPoints,
 
1032
                      ShowLines  : Boolean);
 
1033
 
 
1034
var
 
1035
  CharHeight, I, L, Lmax : Integer;
 
1036
  N, Nmax, Xp, Xl, Y     : Integer;
 
1037
begin
 
1038
  N := 0;     { Nb of legends to be plotted  }
 
1039
  Lmax := 0;  { Length of the longest legend }
 
1040
 
 
1041
  for I := 1 to NCurv do
 
1042
    if CurvParam^[I].Legend <> '' then
 
1043
      begin
 
1044
        Inc(N);
 
1045
        L := Canvas.TextWidth(CurvParam^[I].Legend);
 
1046
        if L > Lmax then Lmax := L;
 
1047
      end;
 
1048
 
 
1049
  if (N = 0) or (Lmax = 0) then Exit;
 
1050
 
 
1051
  { Character height }
 
1052
  CharHeight := Canvas.TextHeight('M');
 
1053
 
 
1054
  { Max. number of legends which may be plotted }
 
1055
  Nmax := Round((YmaxPixel - YminPixel) / CharHeight) - 1;
 
1056
  if N > Nmax then N := Nmax;
 
1057
 
 
1058
  { Draw rectangle around the legends }
 
1059
  Canvas.Rectangle(XmaxPixel + Round(0.02 * GraphWidth), YminPixel,
 
1060
                   XmaxPixel + Round(0.12 * GraphWidth) + Lmax,
 
1061
                   YminPixel + (N + 1) * CharHeight);
 
1062
 
 
1063
  L := Round(0.02 * GraphWidth);  { Half-length of line }
 
1064
  Xp := XmaxPixel + 3 * L;        { Position of symbol  }
 
1065
  Xl := XmaxPixel + 5 * L;        { Position of legend  }
 
1066
 
 
1067
  if NCurv <= Nmax then N := NCurv else N := Nmax;
 
1068
 
 
1069
  for I := 1 to N do
 
1070
    with Canvas do
 
1071
      begin
 
1072
        SetGraphSettings(Canvas, I);
 
1073
 
 
1074
        { Plot point and line }
 
1075
        Y := YminPixel + I * CharHeight;
 
1076
        if ShowPoints then
 
1077
          PlotSymbol(Canvas, Xp, Y, I);
 
1078
        if ShowLines then
 
1079
          PlotLine(Canvas, Xp - L, Y, Xp + L, Y);
 
1080
 
 
1081
        { Write legend }
 
1082
        Brush.Style := bsClear;
 
1083
        Canvas.TextOut(Xl, Y - CharHeight div 2, CurvParam^[I].Legend);
 
1084
      end;
 
1085
 
 
1086
  RestoreGraphSettings(Canvas);
 
1087
end;
 
1088
 
 
1089
procedure ConRec(Canvas     : TCanvas;
 
1090
                 Nx, Ny, Nc : Integer;
 
1091
                 X, Y, Z    : PVector;
 
1092
                 F          : PMatrix);
 
1093
 
 
1094
const
 
1095
  { Mapping from vertex numbers to X offsets }
 
1096
  Im : array[0..3] of Integer = (0, 1, 1, 0);
 
1097
 
 
1098
  { Mapping from vertex numbers to Y offsets }
 
1099
  Jm : array[0..3] of Integer = (0, 0, 1, 1);
 
1100
 
 
1101
  { Case switch table }
 
1102
  CasTab : array[0..2, 0..2, 0..2] of Integer =
 
1103
  (((0,0,8), (0,2,5), (7,6,9)),
 
1104
   ((0,3,4), (1,3,1), (4,3,0)),
 
1105
   ((9,6,7), (5,2,0), (8,0,0)));
 
1106
 
 
1107
var
 
1108
  I, J, K, M, M1, M2, M3 : Integer;
 
1109
  X1, X2, Y1, Y2         : Float;
 
1110
  Fmin, Fmax             : Float;
 
1111
  Xp, Yp                 : PIntVector;
 
1112
  PrmErr                 : Boolean;
 
1113
 
 
1114
var
 
1115
  H   : array[0..4] of Float;    { Relative heights of the box above contour }
 
1116
  Ish : array[0..4] of Integer;  { Sign of H() }
 
1117
  Xh  : array[0..4] of Integer;  { X coordinates of box }
 
1118
  Yh  : array[0..4] of Integer;  { Y coordinates of box }
 
1119
 
 
1120
label
 
1121
  Case0, NoneInTri, NoneInBox;
 
1122
 
 
1123
begin
 
1124
  { Check the input parameters for validity }
 
1125
 
 
1126
  PrmErr := False;
 
1127
  SetErrCode(MatOk);
 
1128
 
 
1129
  if (Nx <= 0) or (Ny <= 0) or (Nc <= 0) then PrmErr := True;
 
1130
 
 
1131
  for K := 1 to Nc - 1 do
 
1132
    if Z^[K] <= Z^[K - 1] then PrmErr := True;
 
1133
 
 
1134
  if PrmErr then
 
1135
    begin
 
1136
      SetErrCode(MatErrDim);
 
1137
      Exit;
 
1138
    end;
 
1139
 
 
1140
  { Convert user coordinates to pixels }
 
1141
 
 
1142
  DimIntVector(Xp, Nx);
 
1143
  DimIntVector(Yp, Ny);
 
1144
 
 
1145
  for I := 0 to Nx do
 
1146
    Xp^[I] := Xpixel(X^[I]);
 
1147
 
 
1148
  for J := 0 to Ny do
 
1149
    Yp^[J] := Ypixel(Y^[J]);
 
1150
 
 
1151
  { Scan the array, top down, left to right }
 
1152
 
 
1153
  for J := Ny - 1 downto 0 do
 
1154
  begin
 
1155
    for I := 0 to Nx - 1 do
 
1156
    begin
 
1157
      { Find the lowest vertex }
 
1158
      if F^[I]^[J] < F^[I]^[J + 1] then
 
1159
        Fmin := F^[I]^[J]
 
1160
      else
 
1161
        Fmin := F^[I]^[J + 1];
 
1162
 
 
1163
      if F^[I + 1]^[J] < Fmin then
 
1164
        Fmin := F^[I + 1]^[J];
 
1165
 
 
1166
      if F^[I + 1]^[J + 1] < Fmin then
 
1167
        Fmin := F^[I + 1]^[J + 1];
 
1168
 
 
1169
      { Find the highest vertex }
 
1170
      if F^[I]^[J] > F^[I]^[J + 1] then
 
1171
        Fmax := F^[I]^[J]
 
1172
      else
 
1173
        Fmax := F^[I]^[J + 1];
 
1174
 
 
1175
      if F^[I + 1]^[J] > Fmax then
 
1176
        Fmax := F^[I + 1]^[J];
 
1177
 
 
1178
      if F^[I + 1]^[J + 1] > Fmax then
 
1179
        Fmax := F^[I + 1]^[J + 1];
 
1180
 
 
1181
      if (Fmax < Z^[0]) or (Fmin > Z^[Nc - 1]) then
 
1182
        goto NoneInBox;
 
1183
 
 
1184
      { Draw each contour within this box }
 
1185
      for K := 0 to Nc - 1 do
 
1186
      begin
 
1187
        if (Z^[K] < Fmin) or (Z^[K] > Fmax) then
 
1188
          goto NoneInTri;
 
1189
 
 
1190
        for M := 4 downto 0 do
 
1191
        begin
 
1192
          if M > 0 then
 
1193
          begin
 
1194
            H[M] := F^[I + Im[M - 1]]^[J + Jm[M - 1]] - Z^[K];
 
1195
            Xh[M] := Xp^[I + Im[M - 1]];
 
1196
            Yh[M] := Yp^[J + Jm[M - 1]];
 
1197
          end;
 
1198
 
 
1199
          if M = 0 then
 
1200
          begin
 
1201
            H[0] := (H[1] + H[2] + H[3] + H[4]) / 4;
 
1202
            Xh[0] := (Xp^[I] + Xp^[I + 1]) div 2;
 
1203
            Yh[0] := (Yp^[J] + Yp^[J + 1]) div 2;
 
1204
          end;
 
1205
 
 
1206
          if H[M] > 0 then Ish[M] := 2;
 
1207
          if H[M] < 0 then Ish[M] := 0;
 
1208
          if H[M] = 0 then Ish[M] := 1;
 
1209
        end; { next M }
 
1210
 
 
1211
        { Scan each triangle in the box }
 
1212
        X1 := 0.0;
 
1213
        X2 := 0.0;
 
1214
        Y1 := 0.0;
 
1215
        Y2 := 0.0;
 
1216
        for M := 1 to 4 do
 
1217
        begin
 
1218
          M1 := M; M2 := 0; M3 := M + 1;
 
1219
          if M3 = 5 then M3 := 1;
 
1220
 
 
1221
          case CasTab[Ish[M1], Ish[M2], Ish[M3]] of
 
1222
            0 :
 
1223
              goto Case0;
 
1224
 
 
1225
            { Line between vertices M1 and M2 }
 
1226
            1 : begin
 
1227
              X1 := Xh[M1];
 
1228
              Y1 := Yh[M1];
 
1229
              X2 := Xh[M2];
 
1230
              Y2 := Yh[M2];
 
1231
            end;
 
1232
 
 
1233
            { Line between vertices M2 and M3 }
 
1234
            2 : begin
 
1235
              X1 := Xh[M2];
 
1236
              Y1 := Yh[M2];
 
1237
              X2 := Xh[M3];
 
1238
              Y2 := Yh[M3];
 
1239
            end;
 
1240
 
 
1241
            { Line between vertices M3 and M1 }
 
1242
            3 : begin
 
1243
              X1 := Xh[M3];
 
1244
              Y1 := Yh[M3];
 
1245
              X2 := Xh[M1];
 
1246
              Y2 := Yh[M1];
 
1247
            end;
 
1248
 
 
1249
            { Line between vertex M1 and side M2-M3 }
 
1250
            4 : begin
 
1251
              X1 := Xh[M1];
 
1252
              Y1 := Yh[M1];
 
1253
              X2 := (H[M3] * Xh[M2] - H[M2] * Xh[M3]) / (H[M3] - H[M2]);
 
1254
              Y2 := (H[M3] * Yh[M2] - H[M2] * Yh[M3]) / (H[M3] - H[M2]);
 
1255
            end;
 
1256
 
 
1257
            { Line between vertex M2 and side M3-M1 }
 
1258
            5 : begin
 
1259
              X1 := Xh[M2];
 
1260
              Y1 := Yh[M2];
 
1261
              X2 := (H[M1] * Xh[M3] - H[M3] * Xh[M1]) / (H[M1] - H[M3]);
 
1262
              Y2 := (H[M1] * Yh[M3] - H[M3] * Yh[M1]) / (H[M1] - H[M3]);
 
1263
            end;
 
1264
 
 
1265
            { Line between vertex M3 and side M1-M2 }
 
1266
            6 : begin
 
1267
              X1 := Xh[M3];
 
1268
              Y1 := Yh[M3];
 
1269
              X2 := (H[M2] * Xh[M1] - H[M1] * Xh[M2]) / (H[M2] - H[M1]);
 
1270
              Y2 := (H[M2] * Yh[M1] - H[M1] * Yh[M2]) / (H[M2] - H[M1]);
 
1271
            end;
 
1272
 
 
1273
            { Line between sides M1-M2 and M2-M3 }
 
1274
            7 : begin
 
1275
              X1 := (H[M2] * Xh[M1] - H[M1] * Xh[M2]) / (H[M2] - H[M1]);
 
1276
              Y1 := (H[M2] * Yh[M1] - H[M1] * Yh[M2]) / (H[M2] - H[M1]);
 
1277
              X2 := (H[M3] * Xh[M2] - H[M2] * Xh[M3]) / (H[M3] - H[M2]);
 
1278
              Y2 := (H[M3] * Yh[M2] - H[M2] * Yh[M3]) / (H[M3] - H[M2]);
 
1279
            end;
 
1280
 
 
1281
            { Line between sides M2-M3 and M3-M1 }
 
1282
            8 : begin
 
1283
              X1 := (H[M3] * Xh[M2] - H[M2] * Xh[M3]) / (H[M3] - H[M2]);
 
1284
              Y1 := (H[M3] * Yh[M2] - H[M2] * Yh[M3]) / (H[M3] - H[M2]);
 
1285
              X2 := (H[M1] * Xh[M3] - H[M3] * Xh[M1]) / (H[M1] - H[M3]);
 
1286
              Y2 := (H[M1] * Yh[M3] - H[M3] * Yh[M1]) / (H[M1] - H[M3]);
 
1287
            end;
 
1288
 
 
1289
            { Line between sides M3-M1 and M1-M2 }
 
1290
            9 : begin
 
1291
              X1 := (H[M1] * Xh[M3] - H[M3] * Xh[M1]) / (H[M1] - H[M3]);
 
1292
              Y1 := (H[M1] * Yh[M3] - H[M3] * Yh[M1]) / (H[M1] - H[M3]);
 
1293
              X2 := (H[M2] * Xh[M1] - H[M1] * Xh[M2]) / (H[M2] - H[M1]);
 
1294
              Y2 := (H[M2] * Yh[M1] - H[M1] * Yh[M2]) / (H[M2] - H[M1]);
 
1295
            end;
 
1296
          end;  { case }
 
1297
 
 
1298
          Canvas.Pen.Color := CurvParam^[K mod MaxCurv + 1].LineParam.Color;
 
1299
          PlotLine(Canvas, Trunc(X1), Trunc(Y1), Trunc(X2), Trunc(Y2));
 
1300
Case0:
 
1301
        end;  { next M }
 
1302
NoneInTri:
 
1303
      end;  { next K }
 
1304
NoneInBox:
 
1305
    end;  { next I }
 
1306
  end;  { next J }
 
1307
end;
 
1308
 
 
1309
procedure LeaveGraphics;
 
1310
begin
 
1311
  DelCurvParamVector(CurvParam, MaxCurv);
 
1312
end;
 
1313
 
 
1314
end.