program vdemo;
(*  compile this demo program from the command-line as
    DCC32 -u"LIB4D" VDEMO.DPR
    or open as a project in the Delphi IDE
    and add LIB4D to the search path in
    Project/Options
    For multi-core support, use "LIB7M" instead of "LIB4D"

    Copyright 1996-2008 OptiCode - Dr. Martin Sander Software Dev.
*)

{$N+} {$D+}
uses SysUtils, WinTypes, WinProcs, Messages, Classes, Graphics, Controls,
     Forms, Dialogs, ExtCtrls, StdCtrls,
     VecLib, VFstd, VFmath, VCFstd, VCFmath, VIstd, Vgraph;

type
      TForm1 = class(TForm)
        PaintBox1: TPaintBox;
        BNextView: TButton;
        procedure ShowView(Sender: TObject);
        procedure NextView(Sender: TObject);
    end;

var
   Form1:                              TForm1;
   X1, X2, Y1, Y2, Y3, Spc, Freq, Win: fVector ;
   CX1, CX2:                           cfVector;
   I1:                                 iVector;
   deltat, fOscill, fNyquist:          Single;
   vsize, spcsize:                     UInt;
   tStart, tStop:                      TDateTime;
   vview:                              Integer;
const
   ProgName         = 'OptiVec Demo';
   {$R *.DFM}

function StopTime:Double;
var  tBegin, tEnd:         double;
     Hour, Min, Sec, MSec: word;
begin
   tStop := Now;
   DecodeTime( tStart, Hour, Min, Sec, MSec );
   tBegin := (Hour*60 + Min) * 60 + Sec + MSec * 0.001;
   DecodeTime( tStop, Hour, Min, Sec, MSec );
   tEnd := (Hour*60 + Min) * 60 + Sec + MSec * 0.001;
   StopTime := tEnd - tBegin;
end;


procedure TForm1.NextView(Sender: TObject);
begin
       inc(vview);
       PaintBox1.Invalidate;
end;

procedure TForm1.ShowView(Sender: TObject);
var   StrDum:       array[0..80] of char;
      ExTime:       array[0..10] of char;
      i, j:         Word;
      xt:           Double;
      X1Dum, X2Dum: fVector;
      DC:           HDC;
const A:Single =  1.21;
      B:Single = -0.13;
      C:Single =  0.85;  (*any arbitrary values! *)
begin
   V_initPlot( PaintBox1 );
   DC := PaintBox1.Canvas.Handle;
   case vview  of
       0: begin TextOut( DC, 10, 10,
                'This is a series of graphs illustrating OptiVec functions.', 58 );
                TextOut( DC, 10, 40,
                'Press [Alt] [F4] to end the demonstration.', 41 );
          end;
       1: begin TextOut( DC, 0, 10, 'Sine wave f = 100 Hz', 20 );
                deltat   := 0.001;  (* sampling interval 1 millisecond *)
                fNyquist := 0.5 / deltat;   (* Nyquist frequency = 500 Hz *)
                fOscill  := 100;  (* say we are sampling a 100 Hz oscillation *)
                VF_ramp( X1, vsize, 0, deltat );  (* Time axis *)
                VF_ramp( Freq, spcsize+1, 0, fNyquist / spcsize );
                VFx_sin( Y1, X1, vsize, 2*PI*fOscill, 0, 1 ); (* sine wave( omega*t) *)
                VF_xyAutoPlot( X1, Y1, 101, PS_SOLID, LIGHTRED );
                       (* show only 10 waves *)
          end;
       2: begin TextOut( DC, 0, 10, 'Portion of an asymmetric square wave, f = 100 Hz', 48 );
                VF_cmp_gtC( Y2, Y1, vsize, 0.6 );
                    (*  transform the sine into an asymmetric square wave *)
                V_drawAxes( 0, 100, -0.3, 1.2 );
                VF_yDataPlot( Y2, 101, PS_SOLID, BLUE );
                       (* show again only 10 waves *)
           end;
       3: begin TextOut( DC, 0, 8, 'Frequency spectrum of the square wave from the last page', 56 );
                TextOut( DC, 0, 23, 'Note the higher harmonics of the 100 Hz ground freq. and the asymmetry-induced f=0 term', 87 );
                VF_Hanning( Win, 2*spcsize );   (* or VF_Welch, VF_Parzen *)
                VF_Pelement( Spc, spcsize )^ :=
                VF_spectrum( Spc, spcsize, Y2, vsize, Win );
                    (*  analyse the frequency spectrum of the square wave *)
                VF_xyAutoPlot( Freq, Spc, spcsize+1, PS_SOLID or SY_CROSS, GREEN );
           end;
       4: begin TextOut( DC, 0,  8, 'Trigonometric functions:', 24 );
                TextOut( DC, 0,  23, 'Red: sine (*10, +20),     green: tangent,     blue:  cosecant (-20)', 67 );
                VI_ramp( I1, vsize, 0, 1 );
                VF_ramp( X1, vsize, 0, 4.0 /(vsize-1) );
                VF_sinrpi2( Y1, I1, vsize, (vsize-1)div 2);
                    (* get a sine wave, magnify it by a factor of 10
                       and displace it vertically by +20 units. *)
                VFx_equV( Y1, Y1, vsize, 10.0, 20.0 );
                V_setFPErrorHandling (fperrIgnore);
                 (* avoid error messages for the SING errors
                    in the following function calls  *)
                VF_tanrpi2( Y2, I1, vsize, (vsize-1)div 2);
                VF_limit( Y2, Y2, vsize, -40, 40 );
                   (*  limit the range of tangent values  *)
                VF_cosecrpi2( Y3, I1, vsize, (vsize-1)div 2);
                VF_limit( Y3, Y3, vsize, -40, 40 );
                    (*  the same for cosecant values  *)
                VF_subC( Y3, Y3, vsize, 20.0 );
                     (*  displace the cosecant by -20 units. *)
                V_setFPErrorHandling (fperrDefaultHandling);
                VF_xy2AutoPlot( X1, Y2, vsize, PS_SOLID, LIGHTGREEN,
                                X1, Y3, vsize, PS_SOLID, LIGHTBLUE );
                VF_xyDataPlot( X1, Y1, vsize, PS_SOLID, LIGHTRED );
           end;
       5: begin TextOut( DC, 0,  8, 'Playing with functions in the complex plane:', 44 );
                TextOut( DC, 0,  23, 'Red:   complex sine,  green: complex cosine',  43 );
                CX1^.Re := PI;   CX1^.Im := 0;
                CX2^.Re := 0.1; CX2^.Im := 0.002;
                VCF_ramp( CX1, vsize, CX1^, CX2^ );
                VCF_cos( CX2, CX1, vsize );
                VCF_sin( CX1, CX1, vsize );
                 (*  Try also with other complex functions! *)
                VCF_2AutoPlot( CX1, vsize, PS_SOLID, LIGHTRED,
                               CX2, vsize, PS_SOLID, LIGHTGREEN );
           end;
       6: begin TextOut( DC, 0, 10, 'To end the demo, let us make a crude speed comparison', 53 );
                TextOut( DC, 0, 30, 'between compiled code and OptiVec code.', 39 );
                VF_ramp( X1, vsize, -10, 0.005 );
                TextOut( DC, 0, 55, 'E.g., evaluate the exponential function for a whole vector', 58 );
                TextOut( DC, 20, 80, 'Y[i] = c * exp[ a * X[i] + b],       i=0,..,size-1', 50 );
                TextOut( DC, 0,  105, '(1000 runs of 16384 elements each).   Please wait...', 50 );
                tStart := Now;
                for i:=1 to 1000 do
                begin
                    X1Dum := X1;  X2Dum := X2;
                    for j:=1 to vsize do
                    begin
                      X2Dum^ := C * exp( A * X1Dum^ + B );
                      inc( VRec( X1Dum ).Off, sizeof( Single ));
                      inc( VRec( X2Dum ).Off, sizeof( Single ));
                        (* at least, use fast pointers! *)
                    end;
                end;
                xt := StopTime;
                Str( xt:5:2, ExTime );
                StrCopy( StrDum, 'Execution Time for classic loop version:  ');
                StrCat( StrDum, Extime ); StrCat( StrDum, ' sec' );
                TextOut( DC, 0, 130, StrDum, StrLen( StrDum ) );
                tStart := Now;
                for i:=1 to 1000 do
                       VFx_exp( X2, X1, vsize, A, B, C );
                xt := StopTime;
                Str( xt:5:2, ExTime );
                StrCopy( StrDum, 'Execution Time for OptiVec version:     ');
                StrCat( StrDum, Extime ); StrCat( StrDum, ' sec' );
                TextOut( DC, 0, 155, StrDum, StrLen( StrDum ) );
                TextOut( DC, 0, 190, 'Generally, the OptiVec version is about 2 times faster than', 59 );
                TextOut( DC, 0, 210, 'compiled loops; sometimes, up to a factor of 8 can be gained!', 61 );
                vview := -1;
          end;
   end;
end;

{$R *.RES}

begin
    vsize  := 16384;
    spcsize := 1024;
    vview  :=    0;
    V_initMT( 2 );
       (* multi-core initialization for 2 processors; no effect for single-thread *)
    I1  := VI_vector( vsize );
    X1  := VF_vector( vsize );
    X2  := VF_vector( vsize );
    Y1  := VF_vector( vsize );
    Y2  := VF_vector( vsize );
    Y3  := VF_vector( vsize );
    Spc := VF_vector0( spcsize+1 );
    Freq:= VF_vector( spcsize+1 );
    Win := VF_vector( 2*spcsize );
    CX1 := VCF_vector( vsize );
    CX2 := VCF_vector( vsize );
    Application.Initialize;
    Application.CreateForm(TForm1, Form1);
    Application.Run;
    V_freeAll;
end.

