Построение линий уровня(pascal)

[cherkasenok]

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

uses graph,crt;

  const
   n=220;
   m=90;
  // iterk=10000;
   //eps=0.0001;
   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.

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

[Anatole]

Линий уровня именно вашей функции имеют вид прямых линий наклоненных под 45градусов. Это можно доказать математически. поэтому их рисование можно очень легко запрограмировать. А это не интересно. Посему я написал небольшую програмку для этой цели, намного изменив саму функцию и диапазон изменемия аргументов.

program LineUrov;
Uses Graph;
Type
TFun = Function (_x,_y : extended):extended;

Function Fun(x,y : extended):extended;  far;
begin
  Fun:=sin(x+y)+cos(x*y);
end;

Procedure InitGraphMode;
var Gd, Gm, errCode : Integer;
begin
  Gd := Detect;
  InitGraph(Gd, Gm, '..\bp\bgi');
  errCode :=  GraphResult;
  if errCode  <> grOk then
   begin
    case errCode of
     1:writeln('ЋиЁЎЄ*. ѓа*дЁЄ* *Ґ Ё*ЁжЁ*«Ё§Ёа®ў***.');
     2:writeln('ЋиЁЎЄ*. ЌҐ ®Ў**а㦥*® Ја*дЁзҐбЄ®Ґ гбва®©бвў®.');
     3:writeln('ЋиЁЎЄ*. ЌҐ **©¤Ґ* ¤а*©ўҐа гбва®©бвў*.');
     4:writeln('ЋиЁЎЄ*. ЌҐўҐа*л© ¤а*©ўҐа.');
     5:writeln('ЋиЁЎЄ*. ЌҐ еў*в*Ґв Ї*¬пвЁ ¤«п §*Јаг§ЄЁ ¤а*©ўҐа*.');
     6:writeln('ЋиЁЎЄ*. ЌҐ еў*в*Ґв Ї*¬пвЁ ¤«п §*Ї®«*Ґ*Ёп бЄ**ЁагойЁ¬Ё Є®¤*¬Ё.');
     7:writeln('ЋиЁЎЄ*. ЌҐ еў*в*Ґв Ї*¬пвЁ ¤«п §*Ї®«*Ґ*Ёп.');
     8:writeln('ЋиЁЎЄ*. ЌҐ **©¤Ґ* иаЁдв.');
     9:writeln('ЋиЁЎЄ*. ЌҐ еў*в*Ґв Ї*¬пвЁ ¤«п §*Јаг§ЄЁ иаЁдв*.');
     10:writeln('ЋиЁЎЄ*. ‚лЎа***л© ¤а*©ўҐа *Ґ ¬®¦Ґв а*Ў®в*вм ў §*¤***®¬ ०Ё¬Ґ.');
     11:writeln('ѓа*дЁзҐбЄ*п ®иЁЎЄ*.');
     12:writeln('ЋиЁЎЄ* Ја*дЁзҐбЄ®Ј® ўў®¤*/ўлў®¤*.');
     13:writeln('ЋиЁЎЄ*. ЌҐўҐа*л© д*©« иаЁдв*.');
     14:writeln('ЋиЁЎЄ*.ЌҐўҐа*л© *®¬Ґа иаЁдв*.');
     else writeln('ЌҐЁ*ЁжЁ*«Ё§Ёа®ў****п ®иЁЎЄ*.');
    end;
     Halt(1);
   end;
end;

Procedure BuildLineUrov(F:TFun; xmin,ymin,xmax,ymax:extended);
Const dz=0.05;
Var
 i,j,k,xMaxPixel,yMaxPixel:word;
 x,y,z,hx,hy,hz,zmin,zmax : extended;
 Pal:PaletteType;
begin
  InitGraphMode;
  GetPalette(Pal);
  xMaxPixel := GetMaxX;    yMaxPixel:= GetMaxY;
  hx := (xmax-xmin)/(xMaxPixel+1); hy := (ymax-ymin)/(yMaxPixel+1);
  x := xmin; y := ymin;
  zmin := f(x,y); zmax := zmin;
  for i:= 0 to xMaxPixel do
    begin
    y:=ymin;
    for j := 0 to yMaxPixel do
      begin
      z:=f(x,y);
      if z<zmin then zmin:=z;
      if z>zmax then zmax:=z;
      y:=y+hy;
      end;
    x:=x+hx;
    end;
  hz := (zmax-zmin)/pal.size;
  zmax := hz*dz;
  x := xmin;
for i:= 0 to xMaxPixel do
    begin
    y:=ymin;
    for j := 0 to yMaxPixel do
      begin
      z:=f(x,y);
      for k:=0 to Pred(Pal.Size) do
        if (z<zmin+k*hz+zmax) and (z>zmin+k*hz-zmax) then
                  putPixel(i,j,pal.Colors[k]);
      y:=y+hy;
      end;
    x:=x+hx;
    end;
Readln;
CloseGraph;
end;

Begin
BuildLineUrov(Fun, 0{xmin},0{ymin},2.2{xmax},10.9{ymax});
End.

В ней вам необходимо указать ваши пути к файлам драйверов BGI, изменить функцию для построения и нарисовать сетку. Последнее вы умеете делать. Только не делайте сетку слышком яркой, это не главное в вашей програме