{$A+,B-,D+,E+,F+,G-,I+,L+,N+,O-,P-,Q+,R+,S+,T+,V+,X+,Y+}
{$M 16384,16384,655360}

uses rkf, graph;

type real = double;

var G: rkf.T;
    GD, GM: integer;

{ Math hacks: }

function cub(x: real): real;
  begin cub := x * sqr(x) end;

{ Problem definition: }

const
  n = 2;
  NProblems = 6;
  W4 = 20.0;
  S5 = 2.0;
  E6 = 0.05;

var
  Problem: integer;

procedure TrueSolution(t: REAL; n: integer; s: PVector);
  var ct, st, et, A, rho: real;
  begin
    case Problem of
      1:  begin { quartic }
            s^[0] := sqr(sqr(t));
            s^[1] := 4*cub(t);
          end;
      2:  begin { cubic}
            s^[0] := cub(t);
            s^[1] := 3*sqr(t);
          end;
      3:  begin { parabola }
            s^[0] := sqr(t);
            s^[1] := 2*t;
          end;
      4:  begin { sinusoid }
            ct := cos(W4*t);
            st := sin(W4*t);
            s^[0] := ct;
            s^[1] := - W4*st;
          end;
      5:  begin { gaussian }
            et := exp(-sqr(S5*t));
            s^[0] := et;
            s^[1] := -2.0*sqr(S5)*t*et;
          end;
      6:  begin { hyperbola }
            et := sqrt(sqr(t) + sqr(E6));
            s^[0] := et;
            s^[1] := t/et;
          end;
    end
  end;

PROCEDURE Derivative(t: REAL; n: integer; s, v: PVector);
  begin
    case Problem of
      1:  begin
            v^[0] := s^[1];
            v^[1] := 12.0*sqr(t);
          end;
      2:  begin
            v^[0] := s^[1];
            v^[1] := 6.0*t;
          end;
      3:  begin
            v^[0] := s^[1];
            v^[1] := 2.0;
          end;
      4:  begin
            v^[0] := s^[1];
            v^[1] := -sqr(W4)*s^[0];
          end;
      5:  begin
            v^[0] := s^[1];
            v^[1] := s^[0]*(sqr(s^[1]/s^[0]) - 2.0*sqr(S5));
          end;
      6:  begin
            v^[0] := s^[1];
            v^[1] := sqr(E6)/cub(s^[0]);
          end;
    end
  end;

function Boring(t0, t1: real; n: integer; s0, s1: PVector): real;
  begin
    Boring := t1
  end;

{ Integrator test: }

var
  PlotColor: integer;
  PlotTMin, PlotTMax, PlotSMin, PlotSMax: real;
  PlotDotSize: integer;

procedure PlotState(t: real; n: integer; s: PVector);
  var dx, dy: integer;
  begin
    for dx := -PlotDotSize to +PlotDotSize do
      for dy := -PlotDotSize to +PlotDotSize do
        putpixel(
          1 + dx + Round(638.0*(t-PlotTMin)/(PlotTMax-PlotTMin)),
          1 + dy + Round(478.0*(PlotSMax-s^[0])/(PlotSMax-PlotSMin)),
          PlotColor
        );
  end;

procedure TestIntegrator(
    t0, t1: real;
    dtmin, dtmax: real;
    tol: real;
    dotSize: integer;
    color: integer
  );

  var t, tstop: real;
      s0, s: PVector;
      i: integer;
  begin
    PlotColor := color;
    PlotTMin := t0;
    PlotTMax := t1;
    PlotSMin := -2.0;
    PlotSMax := +2.0;
    PlotDotSize := dotSize;

    if tol = 0.0 then begin
      { Plot true solution: }
      s := NEWVECTOR(n);
      for i := 0 to 639 do begin
        t := t0 + (i/639.0)*(t1-t0);
        TrueSolution(t, n, s);
        PlotState(t, n, s)
      end;
      tstop := t1;
      GETRIDOFVECTOR(n, s)
    end else begin
      { Plot numerical solution: }
      s0 := NEWVECTOR(n);
      TrueSolution(t0, n, s0);
      s := NEWVECTOR(n);
      CopyState(n, s0, s);
      tstop := G.integrate(
        t0, t1, dtmin, dtmax, tol, s,
        Derivative, Boring, PlotState
      );
      GETRIDOFVECTOR(n, s0);
      GETRIDOFVECTOR(n, s)
    end;

    { Wait for user OK: }
    readln;
    if tstop <> t1 then writeln('*** tstop = ', tstop);
  end;

begin
  GD:=detect;
  InitGraph(GD, GM, '');

  G.init(n);
  for Problem := 1 to NProblems do begin
    ClearDevice;
    TestIntegrator(-1.0, 1.0,  0.0001, 1.0,  0.0,   0, green);
    TestIntegrator(-1.0, 1.0,  0.0001, 1.0,  0.1,   1, yellow);
    TestIntegrator(-1.0, 1.0,  0.0001, 1.0,  0.01,  1, red);
    TestIntegrator(-1.0, 1.0,  0.0001, 1.0,  0.001, 1, cyan);
  end;

  CloseGraph;
end.