Форум программистов [Powered by Invision Power Board]

Код

Program KRAI;
Uses CRT;
Type
   mix=Array[1..50,1..50] of Real;
   intmx=Array[1..50] of Integer;
Var k,i,j:Byte;
    A,B,R,h,x,y,m,n,d,d1,ep,e:Real;
    u,u2,u3:mix;
    indx:intmx;

Function Uh(x,y:Real):Real;
Begin
   Uh:=sqr(x)*sqr(y);
End;

Function way(x2,y2,x1,y1:Real):Real;
Begin
   way:=sqrt(sqr(x2-x1)+sqr(y2-y1));
End;

Function d_x(A,B,R,x,y:Real):Real;
Var dx:Real;
Begin
   If (R-B*y*y)/A>=0 Then
     Begin
        dx:=sqrt((R-B*y*y)/A);
        d_x:=way(x,y,dx,y);
     End
   Else
     d_x:=h+1;
End;

Function d_y(A,B,R,x,y:Real):Real;
Var dy:Real;
Begin
   If (R-B*x*x)/A>=0 Then
     Begin
        dy:=sqrt((R-A*x*x)/B);
        d_y:=way(x,y,x,dy);
     End
   Else
     d_y:=h+1;
End;

Procedure min_d(Var m,n:Real; A,B,x,y:Real);
Var min:Real;
Begin
   If d_x(A,B,R,x,y)<h Then
      min:=d_x(A,B,R,x,y);
   If d_y(A,B,R,x,y)<h Then
     If d_y(A,B,R,x,y)<min Then
       Begin
          n:=sqrt((R-A*sqr(x))/B);
          m:=x;
       End
     Else
       Begin
          n:=y;
          m:=sqrt((R-B*y*y)/A);
       End;
End;



Procedure Calc;
Begin
   For i:=k DownTo 1 Do
     Begin
        j:=1;
        Repeat
          Begin
             x:=(j-1)*h;
             y:=(i-1)*h;
             If ((way(x,y,0,0)>=d) and (way(x,y,0,0)<d1)) Then
               If (d_x(A,B,R,x,y)<h) or (d_y(A,B,R,x,y)<h) Then
                 Begin
                    If i>k-sqrt(R)/2 Then
                      Begin
                         min_d(m,n,A,B,x,y);
                         u[i][j]:=Uh(m,n);
                         u[j][i]:=u[i][j];
                      End;
                 End
               Else
                 Begin
                   indx[i]:=j-1;
                   Continue
                 End
             Else
               Begin
                  If ((x=0) and (y<k+1) and (y>0)) Then
                    u[i][j]:=1/4*(u[i+1][j]+u[i-1][j]+2*u[i][j+1]);
                  If ((y=0) and (x<k+1) and (x>0)) Then
                    Begin
                       u[i][j]:=1/4*(2*u[i+1][j]+u[i][j+1]+u[i][j-1]);
                       u[j][i]:=u[i][j];
                    End;
                  If ((x=0) and (y=0)) Then
                    u[i][j]:=1/4*u[i+1][j];
                  If ((x>0) and (y>0)) Then
                    u[i][j]:=1/4*(u[i+1][j]+u[i-1][j]+u[i][j+1]+u[i][j-1]);
               End;
             j:=j+1;
          End;
        Until(indx[i]>0)
     End;
End;




Procedure init;
Begin
   For i:=1 To 50 Do
     For j:=1 To 50 Do
       u[i][j]:=0;
   For i:=1 To 50 Do
       indx[i]:=0;
End;
Begin ClrScr;
   WriteLn('Ax^2+By^2=R^2');
   WriteLn('vvedite A, B, R');
   Read(A); Read(B); ReadLn(R);
   WriteLn('Vvedite shag setki H'); ReadLn(h);
   WriteLn('Vvedite tochnost E'); ReadLn(e);
   init;
   k:=round(sqrt(R)/h)+1;
   i:=0; j:=k;
   d:=sqrt(R)-(sqrt(2*sqr(h))/7);
   d1:=sqrt(R)+2*h-h/30;
   Calc;
   Repeat
     Begin
        ep:=0;
        For i:=k DownTo 1 Do
           Begin
             For j:=1 to indx[i] Do
               Begin
                  x:=(j-1)*h;
                  y:=(i-1)*h;
                  If ((way(x,y,0,0)>=d) and (way(x,y,0,0)<d1)) Then
                    Begin
                       If (d_x(A,B,R,x,y)<=h) or (d_y(A,B,R,x,y)<=h) Then
                         Begin
                            If i>k-sqrt(R)/2 Then
                              Begin
                                 min_d(m,n,A,B,x,y);
                                 u2[i][j]:=Uh(m,n);
                                 u[j][i]:=u2[i][j];
                              End;
                         End;
                    End
                  Else
                    Begin
                       If ((x=0) and (y<k+1) and (y>0)) Then
                         u2[i][j]:=1/4*(u[i+1][j]+u[i-1][j]+2*u[i][j+1]);
                       If ((y=0) and (x<k+1) and (x>0)) Then
                         Begin
                            u2[i][j]:=1/4*(2*u[i+1][j]+u[i][j+1]+u[i][j-1]);
                            u[j][i]:=u2[i][j];
                         End;
                       If ((x=0) and (y=0)) Then
                         u2[i][j]:=1/4*4*u[i+1][j];
                       If ((x>0) and (y>0)) Then
                         u2[i][j]:=1/4*(u[i+1][j]+u[i-1][j]+u[i][j+1]+u[i][j-1]);
                    End;
                  If u2[i][j]-u[i][j]>ep Then
                    ep:=u2[i][j]-u[i][j];
                  u[i][j]:=u2[i][j];
                End;
           End;
     End;
   Until (ep<e);

   For i:=k DownTo 1 Do
     Begin
     WriteLn;
       For j:=1 to indx[i] Do
             Write(' ',u[i][j]:2:2);
     End;
ReadLn;
End.