---------------------------------------------------------------------------
--  Champ de vecteurs / Vector Field Plot
---------------------------------------------------------------------------

with Graph;                             use Graph;
with Ada.Numerics;                      use Ada.Numerics;
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;

procedure Champ_vt is

  function F (x, y : Float) return vector2 is
    r2 : constant Float := x * x + y * y;
    r : constant Float := Sqrt (r2);
    a : constant := 0.05;
  begin
    return (a * ((r - 1.0) * (r - 1.0) / r) * x - y,
            a * ((r - 1.0) * (r - 1.0) / r) * y + x);
    --  return (x, 0);
  end F;

  procedure Dessine_Champ (xa, ya, xb, yb : Float; points_x, points_y : Positive) is
    xab : constant Float := xb - xa;
    yab : constant Float := yb - ya;
    ipx : constant Float := 1.0 / Float (points_x + 1);
    ipy : constant Float := 1.0 / Float (points_y + 1);
    carr : constant Float := Float'Max (xab * ipx, yab * ipy);
    x, y, cal : Float;

  begin
    Set_math_plane (xa, ya, xb, yb);
    Draw_axes (units => True, grid_points => True, scale => True);
    cal := 0.0;
    for calibrate in reverse False .. True loop
      for i in 1 .. points_x loop
        for j in 1 .. points_y loop
          x := xa + Float (i) * ipx * xab;
          y := ya + Float (j) * ipy * yab;
          if calibrate then
            cal := Float'Max (cal, Norm2 (F (x, y)));
          else
            if cal > 0.0 then DrawVector (x, y, F (x, y) * (carr / cal)); end if;
          end if;
        end loop;
      end loop;
    end loop;
  end Dessine_Champ;

  procedure Dessin (d : device_type) is
  begin
    current_device := d;
    InitGraph (d, "Champ_Vt.ps");
    Circle (0.0, 0.0, 1.0);
    Dessine_Champ (-3.0, -3.0, 3.0, 3.0, 30, 20);
    CloseGraph;
  end Dessin;

begin
  Dessin (PostScript);
end Champ_vt;