Интегрирование Ньютон-Котес

[Gonzo]

Ребят, помогите, пожалуйста, с процедурой формирования коэффициентов Ньютона-Котеса, насколько я помню, они не зависят от функции, а только от кол-ва узлов.
Кол-во узлов n=7.
Пожалуйста, помогите разобраться. Буду признателен!

Uses crt;
var main:real;

Function Func(x:real):real;
begin
 Func:=exp(-1*ln(1-sqr(exp(x))))*exp(x);
end;

FUNCTION NEWCOT(C,D:REAL;N:INTEGER):REAL;
TYPE MS=ARRAY[0..7] OF INTEGER;
     MS2=ARRAY [3..7] OF MS;
     FF=FILE OF MS2;

VAR H1,H2:MS2;
    FC1,FC2:TEXT;
    X1,SUM,H,Y1:REAL;
    I,J:Byte;
    name1,name2:string;

 BEGIN
  writeln('Vvedite imia, soderzhaschego chislitili koefficientov formul Niutona-Kotesa: ');
  readln(name1);
  clrscr;
  writeln('Vvedite imia, soderzhaschego znamenateli koefficientov formul Niutona-Kotesa: ');
  readln(name2);
  clrscr;
  ASSIGN (FC1,name1);
  ASSIGN (FC2,name2);
  RESET(FC1);
  RESET(FC2);
  FOR I := 3 TO 7 DO
  FOR J := 0 TO 7 DO READ(FC1,H1[I,J]);
  FOR I := 3 TO 7 DO
  FOR J := 0 TO 7 DO READ(FC2,H2[I,J]);
  CLOSE(FC1);
  CLOSE(FC2);
  H:=(D-C)/N;
  X1:=C;
  SUM:=0.0;
  FOR I:=0 TO 7 DO
   BEGIN
    SUM:=SUM+FUNC(X1)*H1[N,I]/H2[N,I];
    X1:=X1+H;
   END;
  NEWCOT:=(D-C)*SUM;
END;

begin
clrscr;
main:=NEWCOT(0,1,7);
writeln('Znachenie integrala: ',main:6:2);
readln;
end.

[Gonzo]

Вот, понемногу разбираюсь:

Uses crt;
var main:real;
    HK:array[1..7] of real;
    im:byte;
{ ---------------------------------------------------------- }
                { NAHOZHDENIE KOEFFICIENTOV }
procedure koef(w:array of real;n:integer;var e:array of real);
 var t:integer;
 begin
  for t:=1 to n do
  e[t]:=w[t]/(n-t+2);
 end;

procedure mnogochlen(n,i:integer;var c:array of real);
 var k,j:integer;
     d:array[1..100] of real;
 begin
  d[1]:=1;
  for j:=1 to n do
   begin
    d[j+1]:=d[j]*j*(-1);
    if j>1 then for k:=j downto 2 do
    d[k]:=d[k]+d[k-1]*j*(-1);
   end;
  c[1]:=d[1];
  for j:=1 to n+1 do c[j]:=i*c[j-1]+d[j];
  koef(c,n,c);
end;

function integral(w:array of real;n:integer):double;
var
t,p:integer;
c:real;
s:double;
 begin
  s:=0;
  p:=n;
  for t:=0 to p+1 do
  s:=s+w[t]*exp((p-t+2)*ln(p));
  integral:=s;
 end;

function facktorial(n:integer):real;
 var t:integer;
     s:real;
 begin
  s:=1;
  if n=0 then s:=1
  else for t:=1 to n do s:=s*t;
  facktorial:=s;
 end;

procedure hkoef(n:integer;var h:array of real);
 var p,j,d,c,i:integer;
     kq:array[0..20] of real;
     s:array[0..20] of real;
 begin
  p:=n;
  if (p mod 2)=1 then d:=round((p-1)*0.5)
  else d:=round(0.5*p);
  for i:=0 to n do
   begin
    mnogochlen(p,i,kq);
    s[i]:=integral(kq,p);
   end;
  for i:=0 to d do
   begin
    if ((p-i) mod 2) = 0 then c:=1
    else
    c:=(-1);
    h[i]:=(c*s[i])/(facktorial(i)*facktorial(p-i)*p);
    h[p-i]:=h[i];
   end;
 end;
{ ---------------------------------------------------------- }

Function Func(x:real):double;
var a,b:real;
begin
 a:=(1-sqr(exp(x)));
 if a<=0 then a:=0.00000000001;
 b:=exp((-1)*ln(a));
 Func:=b*exp(x);
end;

FUNCTION NEWCOT(C,D:REAL;N:INTEGER):real;
VAR X1,SUM,H,Y1:REAL;
    I:Byte;
 BEGIN
  H:=(D-C)/N;
  X1:=C;
  SUM:=0.0;
  FOR I:=0 TO 7 DO
   BEGIN
    SUM:=SUM+FUNC(X1)*HK[I];
    X1:=X1+H;
   END;
  NEWCOT:=(D-C)*SUM;
 END;

begin
 clrscr;
 hkoef(7,HK);
 for im:=1 to 7 do writeln('i= ',im,' HK= ',HK[im]:5:3);
 readln;
 main:=NEWCOT(0,1,7);
 writeln('Znachenie integrala: ',main:6:2);
 readln;
end.

Такой вопрос: как обойти взятие логарифма от отрицаетльно числа или нуля? Посмотрите, правильно ли считаются коэффициенты?