{ etude.pas (c) Kenrou Adachi, 2007/3/1 }
{ to compile: 
  gpc -o etude etude.p g.c \
  -I/usr/X11R6/include -I/usr/local/include \
  -L/usr/X11R6/lib -L/usr/local/lib \
  -lm -lX11 -leggx }



program etude(input, output);

const
   RED	   = 2;
   GREEN   = 3;
   BLUE	   = 4;

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


{ eggx library (in g.c) }
   
procedure ginit;                          external name 'ginit';
procedure gfinish;                        external name 'gfinish';
procedure gsetcolor(pc:  integer);        external name 'gsetcolor';
procedure gline(x0, y0, x1, y1 : real);   external name 'gline';


{ example fuctions }

function fx(x, y :  real) : real;
begin
   fx := cos(x) - x * sin(y)
end; { fx }


function gy(x, y :  real) : real;
begin
   gy := x * sin(y) + y * sin(x)
end; { gy }


{ etpro.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
   gsetcolor(BLUE);
   gline(0.0, 200.0, 639.0, 200.0);
   gline(320.0, 0.0, 320.0, 399.0)
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;
	       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;
	       gline(wx0 + 320.0, wy0 + 200.0, wx1 + 320.0, wy1 + 200.0)
	    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
	 gsetcolor(GREEN)
      else
	 gsetcolor(RED);
      for count := 1 to l do
	 process(a, b, eps, h, e);
      e := -e
   until (e > 0)
end; { drawing }


{ main }

var
   a, b, eps, h	: real;
   m, l, n	: integer;
   dummy	:  char;
begin
   inputdata(a, b, eps, m, l);
   ginit;

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

   axes;

   for n := 1 to m do
      drawing(a, b, eps, h, l);

   writeln('終了しました');
   writeln('任意のキーに続いてリターンキーを押して下さい');
   readln(dummy);
   gfinish;
end.