program INTDUPLA;
uses printer;
var
  HY,SOMAY,SOMA2Y,SOMA4Y,Y,EXP1,EXP2  :real;
  N,I,J,K1,K2  :integer;
  NOVO :string[1];
  TITULO :string[70];
function SIMPSON(Y:real;N:integer) :real;
var
  A,B,HX,SOMAX,SOMA2X,SOMA4X,X :real;
  I :integer;
function FXY(X,Y :real):real;
 begin
  FXY :=  X * Y
 end;
function F1Y(Y,EXPN:real):real;
 begin
  if Y = 0 then F1Y := 0
           else F1Y := EXP( EXPN * LN(Y))
 end;
begin
  SOMA2X := 0; SOMA4X := 0;
  X := F1Y(Y,EXP1);
  HX := (F1Y(Y,EXP2) - X)/N; SOMAX := FXY(X,Y);
  for I := 1 to N div 2 do
    begin
      X := X + HX; SOMA4X := SOMA4X + FXY(X,Y);
      X := X + HX; SOMA2X := SOMA2X + FXY(X,Y)
    end;
  SIMPSON := ( SOMAX + 4 * SOMA4X + 2 * SOMA2X - FXY(X,Y) ) * HX / 3
end;

begin
 writeln('Indique abaixo a integral dupla a ser calculada');
 readln(TITULO);writeln(Lst,#12,TITULO);writeln(Lst,' ');
 repeat
  SOMA2Y := 0; SOMA4Y := 0;
  writeln('Forneca o numero de divisoes e os inversos dos expoentes');
  write(' N = ');readln(N);
  write(' K1 = ');readln(K1);EXP1 := 1 / K1;
  write(' K2 = ');readln(K2);EXP2 := 1 / K2;
    HY := ( 1 - 0 ) / N;
  Y := 0;
  SOMAY := SIMPSON(Y,N);
  for J := 1 to N div 2 do
    begin
      Y := Y + HY; SOMA4Y := SOMA4Y + SIMPSON(Y,N);
      Y := Y + HY; SOMA2Y := SOMA2Y + SIMPSON(Y,N)
    end;
  SOMAY := ( SOMAY + 4 * SOMA4Y + 2 * SOMA2Y - SIMPSON(Y,N)) * HY / 3;
  writeln(Lst,'Num. de div. = ',N:4,
    '  K1 = ',K1:2,'  K2 = ',K2:2,'  Integral -> ',SOMAY:15:12);
  write('Tecle S se deseja novo calculo ');
  readln(NOVO);
 until NOVO <> 'S'
end. { fim da INTDUPLA }