Дифракция электромагнитного поля миллиметрового диапазона на плоских объектах. Струков И.Ф. - 39 стр.

  for j1:=1 to N do
   for j2:=1 to N do begin
     fe1:=2*pi*j1/N;
     ro1:=ro0*j2/N;
     Icos:=Icos+cos(2*pi/z*ro1*ro*cos(fe1))*ro1;
     Isin:=Isin+sin(2*pi/z*ro1*ro*cos(fe1))*ro1;
    end;
  U:=(sqr(Icos)+sqr(Isin)) *sqr(cos(arctan(ro/z)/2));
end;
{     Ф ункци яU1(ro) - вы ч и слени еД Н ч ерезполи номы }
function U1 (ro:real):real;
  function f(x:real):real;
  begin f:=0.79788456+0.00000156*x+0.01659667*x*x+0.00017105*x*x*x
      -
      0.00249511*x*x*x*x+0.00113653*x*x*x*x*x+0.00020033*x*x*x*x*x*
      x;
  end;
  function t(x:real):real;
  begin t:=-2.35619449+0.12499612*x+0.00005650*x*x-0.00637879*x*x*x

      +0.000743448*x*x*x*x+0.00079824*x*x*x*x*x+0.00029166*x*x*x*x*
      x*x;
  end;
  function g(x:real):real;
  begin g:=-0.56249985*x+0.21093573*x*x-0.03954289*x*x*x
      +0.00443319*x*x*x*x-
      0.00031761*x*x*x*x*x+0.00001109*x*x*x*x*x*x;
  end;
  function J(x:real):real;
  begin
        if x>=3 then J:=1/x*1/sqrt(x)*f(3/x)*cos(x+t(3/x))
        else J:=(1/2+g(sqr(x/3)));
  end;
begin
  if ro>0 then U1:=sqr(J(ro/ro0*1.61)) else U1:=sqr(J(-ro/ro0*1.61));
end;
{     П роцедураи ни ци али заци и граф и ки }
procedure init;
begin GraphDriver:=Detect;
InitGraph (GraphDriver, GraphMode, '');
ErrorCode:=GraphResult;
if ErrorCode<>grOk then begin
writeln ('Error: ', GraphErrorMsg(ErrorCode)); readln; end; end;
{ П роцедурапостроени якоорди натной сетки }
procedure XYplot; begin Line (Xm, 0, Xm, round(Ym*0.9));
Line (0, Round(0.9*Ym), Xm*2, Round(0.9*Ym));