Графика.Изолинии(Паскаль)., Построение изолиний для функции z=f(x,y) Опции

[cherkasenok]

Добрый день!
Есть задача: нужно построить линии уровня для функции двух переменных f(x,y)=z=Const, заданной в узлах прямоугольной сетки {xi,yi}.
вот что я написала, но она не хочет работать.Рисует только сетку.

Код

uses graph,crt;

  const
   n=220;
   m=90;
   niz=10;

var

x0e,xke,y0e,yke        :integer;
  x0,xk,y0,yk            :real;
  z,zn,hz,zk             :real;
  hxe,hye                :real;
  alpha,beta,gamma,delta :real;
  hx,hy                  :real;
  xr,yr :array[1..4] of integer;
  p:array [1..n,1..m] of real;
  i,j,ii                 :integer;
  grdriver,grmode        :integer;
  max,min                :real;
  imax,jmax,imin,jmin    :integer;
   kk   :integer;
   x,y      :real;
   lx,ly    :real;


BEGIN
clrscr;

  grDriver := Detect;
  initGraph(grDriver, grMode,'');

  //setlinestyle(solidln,1,thickwidth);
  setcolor(yellow);

x0e:=1; y0e:=1;
  xke:=560;yke:=560;
  hxe:=(xke-x0e)/n;hye:=(yke-y0e)/m;
  x0:=0.0;y0:=0.0;xk:=2.2;yk:=0.9;
  hx:=(xk-x0)/n; hy:=(yk-y0)/m;
  alpha:=(x0e-xke)/(x0-xk);
  beta:=(x0*xke-xk*x0e)/(x0-xk);
  gamma:=(y0e-yke)/(yk-y0);
  delta:=(yk*yke-y0*y0e)/(yk-y0);

  x0:=1;y0:=1;
  for i:=1 to n do
   begin
    x:=x0+(i-1)*hx;
     for j:=1 to m do
       begin
        y:=y0+(j-1)*hy;
        p[i,j]:=sin(x+y);
       end;
   end;

{построение сетки}
  for i := 1 to n do
   begin
    x:=i * hxe;
    MoveTo(round(x), round(y0e));
    LineTo(round(x), round(yke));
   end;
  for j := 1 to m do
   begin
    y := j * hye;
    MoveTo(round(x0e), round(y));
    LineTo(round(xke), round(y));
   end;

    max:=p[1,1];imax:=1;jmax:=1;
   min:=p[1,1];imin:=1;jmin:=1;
   for i:=1 to n do
    begin
     for j:=1 to m do
      begin
       if max<p[i,j] then
        begin
         max:=p[i,j];
         imax:=i;
         jmax:=j;
        end;
       if min>p[i,j] then
        begin
         min:=p[i,j];
         imin:=i;
         jmin:=j;
        end;

      end;
    end;
   zn:=min;zk:=max;hz:=(max-min)/niz;

z:=zn;
  while z<=zk-0.000001*hz do
   begin

for i:=1 to n-1 do
     begin
      lx:=x0+(i)*hx;
      for j:=1 to m-1 do

         begin

         ly:=y0+(j)*hy;

kk:=0;
        if ((p[i,j] <= z) and (z < p[i+1,j])) or
           ((p[i,j] >= z) and (z > p[i+1,j])) then
         begin
          x:=lx+(hx*(z-p[i,j]))/(p[i+1,j]-p[i,j]);
          kk:=kk+1;
          xr[kk]:=round(alpha*x+beta);
          yr[kk]:=round(gamma*ly+delta);
         end;

        if ((p[i+1,j] <= z) and (z < p[i+1,j+1])) or
           ((p[i+1,j] >= z) and (z > p[i+1,j+1])) then
         begin
          y:=ly+(hy*(z-p[i+1,j]))/(p[i+1,j+1]-p[i+1,j]);
          kk:=kk+1;
          xr[kk]:=round(alpha*(lx+hx)+beta);
          yr[kk]:=round(gamma*y+delta);

          if kk=2 then
           begin
            MoveTo(xr[2],yr[2]);
            LineTo(xr[1],yr[1]);
            kk:=0;
           end;
         end;

        if ((p[i,j+1] <= z) and (z < p[i+1,j+1])) or
           ((p[i,j+1] >= z) and (z > p[i+1,j+1])) then
         begin
          x:=lx+(hx*(z-p[i,j+1]))/(p[i+1,j+1]-p[i,j+1]);
          kk:=kk+1;
          xr[kk]:=round(alpha*x+beta);
          yr[kk]:=round(gamma*(ly+hy)+delta);

          if kk=2 then
           begin
            MoveTo(xr[2],yr[2]);
            LineTo(xr[1],yr[1]);
            kk:=0;
           end;
         end;

if ((p[i,j] <= z) and (z < p[i,j+1])) or
           ((p[i,j] >= z) and (z > p[i,j+1])) then
         begin
          y:=ly+(hy*(z-p[i,j]))/(p[i,j+1]-p[i,j]);
          kk:=kk+1;
          xr[kk]:=round(alpha*lx+beta);
          yr[kk]:=round(gamma*y+delta);

          if kk=2 then
           begin
            MoveTo(xr[2],yr[2]);
            LineTo(xr[1],yr[1]);
            kk:=0;
           end;
         end;
       end;
      end;
z:=z+hz;

   end;
readln; readln;
closegraph;
END.

Помогите,Пожалуйста!
Буду очень благодарна!