{ etude.pas (c) Kenrou Adachi, 1993-2004 }

program etude(input, output);

uses crt, graph;   { for Turbo Pascal 5.5 later }
{ uses crtansi, graph for Turbo Pascal 4.0 }

var
   x, y, xx, yy : array[1..40] of real;

{$I a:\work\turbo\ode2\define.pas}

procedure inputdata(var a, b, eps : real; var m, l : integer);
begin
   writeln('x座標の範囲 [-a,a] と y座標の範囲 [-b,b] を入力して下さい');
   write(' a = ? ');
   readln(a);
   write(' b = ? ');
   readln(b);
   writeln('精度の限界 eps を入力して下さい');
   write('eps= ? ');
   readln(eps);
   writeln('描画回数 n を入力して下さい');
   write(' n = ? ');
   readln(m);
   writeln('計算ステップ数 l を入力して下さい');
   write(' l = ? ');
   readln(l)
end; { inputdata }


procedure axes;
begin
   setcolor(BLUE);
   line(0, 200, 639, 200);
   line(320, 0, 320, 399)
end; { axes }


procedure process(a, b, eps, h, e : real);
var
   j					      : integer;
   x0, y0, r0, x1, y1, r1, wx0, wy0, wx1, wy1 : real;
   u, v, u1, v1, q, q_h			      : real;
begin
   for j := 1 to 40 do
   begin
      x0 := x[j];
      y0 := y[j];
      if (abs(x0) <= 2.0 * a) and (abs(y0) <= 2.0 * b) then
      begin
	 r0 := sqrt(x0 * x0 + y0 * y0);
	 if (r0 >= h) then
	 begin
	    u := fx(x0, y0);
	    v := gy(x0, y0);
	    q := sqrt(u * u + v * v);
	    if (q >= eps) then
	    begin
	       q_h := h / q;
	       x1 := x0 + 0.5 * e * u * q_h;
	       y1 := y0 + 0.5 * e * v * q_h;
	       { r1 := sqrt(x1 * x1 + y1 * y1);
	         used for fx(x, y, r), gy(x, y, r) }
	       u1 := fx(x1, y1);
	       v1 := gy(x1, y1);
	       x[j] := x0 + e * u1 * q_h;
	       y[j] := y0 + e * v1 * q_h;
	       wx0 := 320.0 * x0 / a;
	       wy0 := 200.0 * y0 / b;
	       wx1 := 320.0 * x[j] / a;
	       wy1 := 200.0 * y[j] / b;
	       line(round(wx0 + 320.0), round(200.0 - wy0),
		    round(wx1 + 320.0), round(200.0 - wy1))
	    end
	 end
      end
   end
end; { process }


procedure drawing(a, b, eps, h : real; l :integer);
var
   j, count	  : integer;
   rnd, e1, e2, e : real;
begin
   randomize;
   j := 1;
   e1 = -1.0;
   while (e1 <= 0.9) do
   begin
      e2 := -1.0;
      while (e2 <= 0.9) do
      begin
	 rnd := random;
	 xx[j] := a * (e1 + 0.25 * rnd);
	 rnd := random;
	 yy[j] := b * (e2 + 0.4 * rnd);
	 e2 := e2 + 0.4;
	 j := j + 1
      end;
      e1 := e1 + 0.25
   end;

   e := 1.0;
   repeat
      for j := 1 to 40 do
      begin
	 x[j] := xx[j];
	 y[j] := yy[j]
      end;
      
      if (e > 0) then
	 setcolor(GREEN)
      else
	 setcolor(RED);

      for count := 1 to l do
      begin
       {
	 gotoxy(2, 0);
	 write('      ');
	 gotoxy(2, 0);
	 write(count);
       }
	 process(a, b, eps, h, e)
      end;
      e := -e;
   until (e > 0)
end; { drawing }


{ main }

var
   a, b, eps, h		   : real;
   m, l, n, gdriver, gmode : integer;
begin
   gdriver := DETECT;
   clrscr;

   inputdata(a, b, eps, m, l);

   clrscr;
   initgraph(gdriver, gmode, '');

   h := sqrt((a / 320.0) * (a / 320.0) + (b / 200.0) * (b / 200.0));

   axes;

   for n := 1 to m do
   begin
    {
      gotoxy(0, 0);
      write('      ');
      gotoxy(0, 0);
      write(n);
    }

      drawing(a, b, eps, h, l)
   end;

   writeln('終了しました');
   writeln('リターンキーを押して下さい');
   readln

end. { main }


{ end of etude.pas }