UNIT math;

{ ******************************************
  ****   Scientific Subroutine Library  ****
  ****         for Turbo Pascal         ****
  ******************************************

  The following programs were written by Allen Miller and appear in the
  book entitled "Pascal Programs For Scientists And Engineers" which is
  published by Sybex, 1981.
  They were originally typed and submitted to MTPUG in Oct. 1982
    Juergen Loewner
    Hoher Heckenweg 3
    D-4400 Muenster
  They have had minor corrections and adaptations for Turbo Pascal by
    Jeff Weiss
    1572 Peacock Ave.
    Sunnyvale, CA 94087. }

INTERFACE

PROCEDURE gaussj(VAR b: ary2s; y: arys;	VAR  coef: arys;
                 ncol: Integer;	VAR error: Boolean);
PROCEDURE linfit(x, y: ary; VAR y_calc: ary;
	         VAR resid: ary; VAR coef: arys; VAR sig: arys;
		 nrow: Integer;	VAR ncol: Integer);
PROCEDURE meanstd(x: ary; length: Integer; VAR mean: Real;
                  VAR std_dev: Real);
PROCEDURE newton1(VAR x: Real);
PROCEDURE newton2(VAR x: Real);
PROCEDURE plot(x, y, ycalc: ary; m : Integer);
FUNCTION randg(mean,sigma: Real): Real;
FUNCTION random(dummy: Integer): Real;
PROCEDURE simps(lower,upper,tol	: Real; VAR sum: Real);
PROCEDURE bubble_sort1(VAR a: ary; n: Integer);
PROCEDURE bubble_sort2(VAR a: ary; n: Integer);
PROCEDURE quick_sort1(VAR x: ary; n: Integer);
PROCEDURE quick_sort2(VAR x: ary; n: Integer);
PROCEDURE shell_sort(VAR a: ary; n: Integer);
PROCEDURE square(x: ary2; y: ary; VAR a: ary2s; VAR g: arys;
                 nrow,ncol: Integer);
PROCEDURE trapez(lower,upper,tol: Real;	VAR sum: Real);
FUNCTION gamma(x: Real): Real;
FUNCTION bessj(x,n: Real): Real;
PROCEDURE romb(lower,upper,tol: Real; VAR ans: Real);

FUNCTION julian(dd,mm,yyyy: Integer): Integer;

{ ************************************************************************ }

IMPLEMENTATION

PROCEDURE gaussj
	   (VAR	b: ary2s;	{ square matrix of coefficients }
		y: arys;	{ constant vector }
	VAR  coef: arys;	{ solution vector }
	     ncol: Integer;	{ order of matrix }
	VAR error: Boolean);	{ true if matrix singular }

{ Gauss Jordan matrix inversion and solution }

{ B(n,n) coefficient matrix becomes inverse }
{ Y(n) original constant vector }
{ W(n,m) constant vector(s) become solution vector }
{ DETERM is the determinant }
{ ERROR=1 if singular }
{ INDEX(n,3) }
{ NV is number of constant vectors }

LABEL 99;
VAR	w	: ARRAY[1..maxc,1..maxc] OF Real;
	index	: ARRAY[1..maxc,1..3] OF Integer;
	i,j,k,l,nv,irow,icol,n,l1 : Integer;
	determ,pivot,hold,sum,t,ab,big : Real;

PROCEDURE swap(VAR a,b: Real);
VAR hold : Real;

BEGIN	{ swap }
  hold := a;
  a := b;
  b := hold
END;		{ procedure swap }


PROCEDURE gausj2;
label	98;
VAR	i,j,k,l,l1	: Integer;


PROCEDURE gausj3;

VAR	l	: Integer;

BEGIN	{ procedure gausj3 }
	{ interchange rows to put pivot on diagonal }
  IF irow<>icol THEN
    BEGIN
      determ := -determ;
      FOR l:=1 to n DO
	swap(b[irow,l],b[icol,l]);
      IF nv>0 THEN
        FOR l:=1 to nv DO
	  swap(w[irow,l],w[icol,l])
  END	{ if iroe<>icol }
END;	{ gausj3 }

BEGIN	{ procedure gausj2 }
	{ actual start of gaussj }
  error := false;
  nv := 1;	{ single constant vector }
  n := ncol;
  FOR i:=1 to n DO
    BEGIN
      w[i,1] := y[i];	{ copy constant vector }
      index[i,3] := 0
    END;
  determ := 1.0;
  FOR i:=1 to n DO
    BEGIN
    { search for largest element }
    big := 0.0;
    FOR j:=1 to n DO
      BEGIN
	IF index[j,3]<>1 THEN
	  BEGIN
	    FOR k:=1 to n DO
	      BEGIN
		IF index[k,3]>1 THEN
		  BEGIN
		    writeln('ERROR: matrix is singular');
		    error := true;
		    goto 98		{ abort }
                  END;
               IF index[k,3]<1 THEN
		  IF abs(b[j,k])>big THEN
		    BEGIN
		      irow := j;
		      icol := k;
		      big := abs(b[j,k])
		    END
	     END	{ k-loop }
	END
  END;		{ j-loop }

  index[icol,3] := index[icol,3]+1;
  index[i,1] := irow;
  index[i,2] := icol;

  gausj3;	{ further subdivision of gaussj }
		{ divide pivot row by pivot column }
  pivot := b[icol,icol];
  determ := determ*pivot;
  b[icol,icol] := 1.0;

  FOR l:=1 to n DO
    b[icol,l] := b[icol,l]/pivot;
  IF nv>0 THEN
    FOR l:=1 to nv DO
      w[icol,l] := w[icol,l]/pivot;

	{ reduce nonpivot rows }

  FOR l1:=1 to n DO
    BEGIN
      IF l1<>icol THEN
	BEGIN
	  t := b[l1,icol];
	  b[l1,icol] := 0.0;
	  FOR l:=1 to n DO
	    b[l1,l] := b[l1,l]-b[icol,l]*t;
	  IF nv>0 THEN
	    FOR l:=1 to nv DO
	      w[l1,l] := w[l1,l]-w[icol,l]*t;
      END	{ if l1<>icol }
     END
    END;	{ i-loop }
98:
END;		{ gausj2 }

BEGIN		{ gaus-jordan main program }
  gausj2;	{ first half of gaussj }
  IF error THEN goto 99;
     { interchange columns }
  FOR i:=1 to n DO
    BEGIN
      l := n-i+1;
      IF index[l,1]<>index[l,2] THEN
	BEGIN
	  irow := index[l,1];
	  icol := index[l,2];
	  FOR k:=1 to n DO
	    swap(b[k,irow],b[k,icol])
      END	{ if index }
  END;		{ i-loop }
FOR k:=1 to n DO
  IF index[k,3]<>1 THEN
    BEGIN
      writeln(chr(7),'ERROR: matrix is singular');
      error := true;
      goto 99		{ abort }
    END;
  FOR i:=1 to n DO
    coef[i] := w[i,1];
99:
END;	{ procedure gaussj }

{ ************************************************************************ }

PROCEDURE linfit(X,	{ independent variable }
		 y	: ary;	{ dependent variable }
	VAR y_calc	: ary;	{ calculated dep. variable }
	VAR resid	: ary;	{ array of residuals }
	VAR coef	: arys;	{ coefficients }
	VAR sig		: arys;	{ error on coefficients }
		nrow	: Integer;	{ length of ary }
	VAR ncol	: Integer);	{ number of terms }

{ least-squares fit to nrow sets of x and y pairs of points }
{ Seperate  procedure needed:
	SQUARE -> form square coefficient matrix
	GAUSSJ -> Gauss-Jordan elimination }

VAR	xmatr		: ary2;		{ data matrix }
	a		: ary2s;	{ coefficient matrix }
	g		: arys;		{ constant vector }
	error		: Boolean;
	i,j,nm		: Integer;
	xi,yi,yc,srs,see,
	sum_y,sum_y2	: Real;

BEGIN		{ procedure linfit }
  ncol := 3;	{ number of terms }
  FOR i:=1 to nrow DO
    BEGIN	{ setup x matrix }
	xi := x[i];
	xmatr[i,1] := 1.0;	{ first column }
	xmatr[i,2] := xi;		{ second column }
	xmatr[i,3] := 1.0/sqr(xi)	{ third column }
    END;
  square(xmatr,y,a,g,nrow,ncol);
  gaussj(a,g,coef,ncol,error);
  sum_y := 0.0;
  sum_y2 := 0.0;
  srs := 0.0;
  FOR i:=1 to nrow DO
    BEGIN
      yi := y[i];
      yc := 0.0;
      FOR j:=1 to ncol DO
	yc := yc+coef[j]*xmatr[i,j];
      y_calc[i] := yc;
      resid[i] := yc-yi;
      srs := srs+sqr(resid[i]);
      sum_y := sum_y+yi;
      sum_y2 := sum_y2+yi*yi
    END;
  correl_coef := sqrt(1.0-srs/(sum_y2-sqr(sum_y)/nrow));
  IF nrow=ncol THEN nm := 1
  ELSE nm := nrow-ncol;
  see := sqrt(srs/nm);
  FOR i:=1 to ncol DO		{ errors on solution }
    sig[i] := see*sqrt(a[i,i])
END;		{ LINFIT }

{ *********************************************************************** }


FUNCTION julian(dd,mm,yyyy: Integer): Integer;
VAR jm,jy: Integer;
BEGIN
  case mm of
    1,2	: BEGIN
		jy := yyyy-1;
		jm := mm+13
	  END
    ELSE
	  BEGIN
	  	jy := yyyy;
		jm := mm
	  END;
  julian := trunc(365.25*jy)+trunc(30.6001*jm)+dd+1720982
END;	{ julian }

{ ************************************************************************ }

PROCEDURE meanstd
                (x              : ary;  {array of values}
                 length         : Integer;
                 VAR mean       : Real;
                 VAR std_dev    : Real);
 
VAR             i               : Integer;
                sum_x,sum_sq    : Real;
BEGIN {main}
 sum_x := 0;
 sum_sq := 0;
 FOR i:=1 to length DO
        BEGIN
         sum_x := sum_x+x[i];
         sum_sq := sum_sq+x[i]*x[i]
        END;
 mean := sum_x/length;
 std_dev := sqrt((sum_sq-sqr(sum_x)/length)/(length-1))
END     { procedure meanstd };

{ ************************************************************************ }

PROCEDURE newton1(VAR x: Real);
CONST	tol		= 1.0E-6;
	max		= 20;
VAR	fx,dfx,dx,x1	: Real;
	i		: Integer;

BEGIN	{ newton }
  error := false;
  i := 0;
  REPEAT
    i := i+1;
    x1 := x;
    func(x,fx,dfx);
    IF dfx=0.0 THEN
      BEGIN
	error := true;
	x := 1.0;
	writeln(chr(7),'ERROR: slope zero')
      END
    ELSE
      BEGIN
	dx := fx/dfx;
	x := x1-dx;
	writeln('x=',x,'  fx=',fx,'  dfx=',dfx)
      END
  UNTIL
    error OR
      (i>max) OR
	(abs(dx)<=abs(tol*x));
  IF i>max THEN
    BEGIN
      writeln(chr(7),'ERROR: no convergence in ',max,' loops');
      error := true
  END
END;		{ newton }

{ ************************************************************************ }

PROCEDURE newton2(VAR x: Real);
CONST	tol		= 1.0E-6;
VAR	fx,dfx,dx,x1	: Real;

BEGIN	{ newton }
  error := false;
  REPEAT
    x1 := x;
    func(x,fx,dfx);
    IF dfx=0.0 THEN
      BEGIN
	error := true;
	x := 1.0;
	writeln(chr(7),'ERROR: slope zero')
      END
    ELSE
      BEGIN
	dx := fx/dfx;
	x := x1-dx;
	writeln('x=',x,'  fx=',fx,'  dfx=',dfx)
      END
  UNTIL
    error OR
      (abs(dx)<=abs(tol*x))
END;		{ newton }

{ ************************************************************************ }

PROCEDURE plot(		{ with arrays }
		x,		{ as independant variable }
		y,		{ as dependant variable }
		ycalc		{ as fitted curve }
		: ary;
	{ and }	m : Integer	{ number of points });

{ plot y and ycalc as a function of x for m points }
{ if m is negative, only x and y are plotted }

CONST	blank	= ' ';
	linel	= 51;

VAR
	ylabel		: ARRAY[1..6] of Real;
	out		: ARRAY[1..linel] of char;
	lines,i,j,jp,l,n: Integer;
	iskip,yonly	: Boolean;

	xlow,xhigh,xnext,xlabel,xscale,signxs,
	ymin,ymax,change,yscale,ys10		: Real;

function pscale(p: Real): Integer;
BEGIN
  pscale := trunc((p-ymin)/yscale+1)
END;	{ pscale}

procedure outlin(xname: Real);
{ output a line }

VAR	i,max	: Integer;

BEGIN
  write(xname:8:2,blank);	{ line label }
  max := linel+1;
  REPEAT		{ skip blanks on end of line }
    max := max-1
  UNTIL (out[max]<>blank) OR (max=1);
  FOR i:=1 to max DO
    write(out[i]);
  writeln;
  FOR i:=1 to max DO
    out[i] := blank	{ blank next line }
END;	{ outlin}

procedure setup(index: Integer);
{ setup the plus and asterisk for printing }

CONST	star = '*';
	plus = '+';

VAR	i    : Integer;

BEGIN
  i := pscale(y[index]);
  out[i] := plus;
  IF not yonly THEN
    BEGIN		{ add ycalc too }
      i := pscale(ycalc[index]);
      out[i] := star
    END
END;		{ setup }


BEGIN		{ body of plot }
  IF m>0 THEN		{ plot y and ycalc vs x }
    BEGIN
      n := m;
      yonly := false
    END
  ELSE		{ plot only y vs x }
    BEGIN
      n := -m;
      yonly := true
    END;
  { space out alternate lines }
  lines := 2*(n-1)+1;
  writeln;
  xlow := x[1];
  xhigh := x[n];
  ymax := y[1];
  ymin := ymax;
  xscale := (xhigh-xlow)/(lines-1);
  signxs := 1.0;
  IF xscale<0.0 THEN signxs := -1.0;
  FOR i:=1 to n DO
    BEGIN
      IF y[i]<ymin THEN ymin := y[i];
      IF y[i]>ymax THEN ymax := y[i];
      IF not yonly THEN
	BEGIN
	  IF ycalc[i]<ymin THEN ymin := ycalc[i];
	  IF ycalc[i]>ymax THEN ymax := ycalc[i]
	END	{ if yonly }
  END;
  yscale := (ymax-ymin)/(linel-1);
  ys10 := yscale*10;
  ylabel[1] := ymin;	{ y axis }
  FOR i:=1 to 4 DO
    ylabel[i+1] := ylabel[i]+ys10;
  ylabel[6] := ymax;
  FOR i:=1 to linel DO
    out[i] := blank;	{ blank line }
  setup(1);
  l := 1;
  xlabel := xlow;
  iskip := false;

  FOR i:=2 to lines DO		{ set up a line }
    BEGIN
      xnext := xlow+xscale*(i-1);
      IF iskip THEN writeln(' -')
      ELSE
	BEGIN
	  l := l+1;
	  WHILE
	    (x[l]-(xnext-0.5*xscale))*signxs<=0.0 DO
	      BEGIN
		setup(l);	{ setup print line }
		l := l+1
	   END;		{ while }
	outlin(xlabel);	{ print a line }
	FOR j:=1 to linel DO
	  out[j] := blank		{ blank line }
      END;		{ if skip }
    IF (x[l]-(xnext+0.5*xscale))*signxs>0.0 THEN iskip := true
      ELSE
	BEGIN
	  iskip := false;
	  xlabel := xnext;
	  setup(l)		{ setup print line }
	END
    END;		{ for-loop }
  outlin(xhigh);	{ last line }
  write('    ');
  FOR i:=1 to 6 DO
    write('     ^    ');
  writeln;
  write('   ');
  FOR i:=1 to 6 DO
    write(ylabel[i]:9:1,blank);
  writeln;
  writeln
END;			{ PLOT }

{ ************************************************************************ }

function randg(mean,sigma: Real): Real;

{ produce random numbers with a gaussian distribution }
{ MEAN and SIGMA are supplied by calling program }
{ function RANDOM is required !!! }

VAR
 i	: Integer;
 sum	: Real;

BEGIN
  sum := 0.0;
  FOR i:=1 to 12 DO
	sum := sum+random(0);
  randg := (sum-6)*sigma+mean
END;	{ function randg }

{ ************************************************************************ }

FUNCTION random(dummy: Integer): Real;	{ --> 29}
{ random number 0-1 }
		{ DEFINE SEED=4.0 AS GLOBAL !!!!!!!! }
{ adapted from HP-35 applications programs }

CONST
 pi	= 3.14159;

VAR
 x	: Real;
 i	: Integer;

BEGIN	{ RANDOM }
 x := seed+pi;
 x := exp(5.0*ln(x));
 seed := x-trunc(x);
 random := seed
END;		{ RANDOM }

{ ************************************************************************ }

PROCEDURE simps(
		lower,upper,tol	: Real;
		VAR sum		: Real);
 
{ numerical integration by Simpson's rule }
{ function is fx, limits are lower and upper }
{ with number of regions equal to pieces }
{ partition is delta_x, answer is sum }
 
VAR	i		: Integer;
	x,delta_x,even_sum,
	odd_sum,end_sum,
	end_cor,sum1	: Real;
	pieces		: Integer;
 
BEGIN
  pieces := 2;
  delta_x := (upper-lower)/pieces;
  odd_sum := fx(lower+delta_x);
  even_sum := 0.0;
  end_sum := fx(lower)+fx(upper);
  end_cor := dfx(lower)-dfx(upper);
  sum := (end_sum+4.0*odd_sum)*delta_x/3.0;
  REPEAT
    pieces := pieces*2;
    sum1 := sum;
    delta_x := (upper-lower)/pieces;
    even_sum := even_sum+odd_sum;
    odd_sum := 0.0;
    FOR i:=1 to pieces div 2 DO
      BEGIN
	x := lower+delta_x*(2.0*i-1.0);
	odd_sum := odd_sum+fx(x)
      END;
    sum := (7.0*end_sum+14.0*even_sum+16.00*odd_sum
				+end_cor*delta_x)*delta_x/15.0;
  UNTIL (sum<>sum1) AND (abs(sum-sum1)<=abs(tol*sum))
END;	{ simps }
 
{ ************************************************************************ }

PROCEDURE bubble_sort1(VAR a: ary; n: Integer);
VAR	i,j	: Integer;
	hold	: Real;

BEGIN		{ procedure sort }
  FOR i:=1 to n-1 DO
    FOR j := i+1 to n DO
      BEGIN
	IF a[i]>a[j] THEN
	  BEGIN
	    hold := a[i];
	    a[i] := a[j];
	    a[j] := hold
	  END
    END		{ for }
END;		{ procedure sort }

{ ************************************************************************ }

PROCEDURE bubble_sort2(VAR a: ary; n: Integer);
{ adapted from 'Introduction to PASCAL',
	R.Zaks, Sybex, 1980 }

VAR	no_change	: Boolean;
	j		: Integer;

procedure swap(p,q: Real);
VAR	hold	: Real;
BEGIN
  hold := p;
  p := q;
  q := hold
END;		{ swap }

BEGIN		{ procedure sort }
  REPEAT
    no_change := true;
    FOR j:=1 to n-1 DO
      BEGIN
	IF a[j]>a[j+1] THEN
	  BEGIN
	    swap(a[j],a[j+1]);
	    no_change := false
	  END
	END	{ for }
  UNTIL no_change
END;	{ procedure sort }

{ ************************************************************************ }

PROCEDURE quick_sort(VAR x: ary; n: Integer);
{ a NONRECURSIVE quicksort routine }
{ Adapted from 'Software-Tools',
	B.Kernighan, Addison Wesley, 1976 }

VAR	left,right	: ARRAY[1..20] of Integer;
	i,j,sp,mid	: Integer;
	pivot		: Real;

procedure swap(VAR p,q: Real);
VAR	hold	: Real;

BEGIN
  hold := p;
  p := q;
  q := hold
END;		{ swap }


BEGIN
  left[1] := 1;
  right[1] := n;
  sp := 1;
  WHILE sp>0 DO
    BEGIN
      IF left[sp]>=right[sp] THEN sp := sp-1
      ELSE
	BEGIN
	  i := left[sp];
	  j := right[sp];
	  pivot := x[j];
	  mid := (i+j)div 2;
	  IF (j-i)>5 THEN
	    IF ((x[mid]<pivot)AND(x[mid]>x[i]))
	      OR
		((x[mid]>pivot)AND(x[mid]<x[i]))
		  THEN swap(x[mid],x[j])
	    ELSE
	      IF((x[i]<x[mid])AND(x[i]>pivot))
		OR ((x[i]>x[mid])AND(x[i]<pivot))
		  THEN swap(x[i],x[j]);
    pivot := x[j];
    WHILE i<j DO
      BEGIN
	WHILE x[i]<pivot DO
	  i := i+1;
	j := j-1;
	WHILE (i<j)AND(pivot<x[j]) DO
	  j := j-1;
	IF i<j THEN swap(x[i],x[j])
    END;	{ while }
    j := right[sp];	{ pivot to i }
    swap(x[i],x[j]);
    IF i-left[sp]>=right[sp]-i THEN
      BEGIN		{ put shorter part first }
	left[sp+1] := left[sp];
	right[sp+1] := i-1;
	left[sp] := i+1
      END
    ELSE
      BEGIN
	left[sp+1] := i+1;
	right[sp+1] := right[sp];
	right[sp] := i-1
      END;
    sp := sp+1		{ push stack }
  END		{ if }
 END		{ while }
END;	{ QUICK SORT }

{ ************************************************************************ }

PROCEDURE quick_sort2(VAR x: ary; n: Integer);
{ a RECURSIVE sorting routine }
{ Adapted from 'The design of Well-Structured and Correct Programs',
	S. Alagic, Springer-Verlag, 1978 }


procedure qsort(VAR x: ary; m,n: Integer);
VAR	i,j	: Integer;


procedure partit(VAR a: ary; VAR i,j: Integer; left,right: Integer);
VAR	pivot	: Real;

procedure swap(VAR p,q: Real);
VAR	hold	: Real;
BEGIN
  hold := p;
  p := q;
  q := hold
END;		{ swap }

BEGIN
  pivot := a[(left+right)div 2];
  i := left;
  j := right;
  WHILE i<=j DO
    BEGIN
      WHILE a[i]<pivot DO
	i := i+1;
      WHILE pivot<a[j] DO
	j := j-1;
      IF i<=j THEN
	BEGIN
	  swap(a[i],a[j]);
	  i := i+1;
	  j := j-1
	END
      END	{ while }
  END;	{ partit }

BEGIN		{ q-sort }
  IF m<n THEN
    BEGIN
      partit(x,i,j,m,n);	{ divide in two }
      qsort(x,m,j);		{ sort left part }
      qsort(x,i,n)		{ sort right part }
    END
END;		{ QSORT }

BEGIN	{ sort }
  qsort(x,1,n)
END;		{ SORT }

{ ************************************************************************ }

PROCEDURE shell_sort(VAR a: ary; n: Integer);
{ Shell-Metzner sort }
{ Adapted from 'Programming in pascal',
	P. Grogono, Addison-Wesley, 1980 }

VAR	done	: Boolean;
	jump,i,j: Integer;

procedure swap(VAR p,q: Real);
VAR	hold	: Real;

BEGIN
  hold := p;
  p := q;
  q := hold
END;		{ swap }

BEGIN
  jump := n;
  WHILE jump>1 DO
    BEGIN
      jump := jump div 2;
      REPEAT
	done := true;
	FOR j:=1 to n DO
	  BEGIN
	    i := j+jump;
	    IF a[j]>a[i] THEN
	      BEGIN
		swap(a[j],a[i]);
		done := false
	    END		{ if }
	END		{ for }
      UNTIL done
  END			{ while }
END;		{ SORT }

{ ************************************************************************ }

PROCEDURE square(x: ary2;
		 y: ary;
	     VAR a: ary2s;
	     VAR g: arys;
	 nrow,ncol: Integer);

{ matrix multiplication routine }
{ a= transpose x times x }
{ g= y times x }

VAR
 i,k,l	: Integer;

BEGIN		{ square }
  FOR k:=1 to ncol DO
    BEGIN
	FOR l:=1 to k DO
	  BEGIN
	    a[k,l] := 0.0;
	    FOR i:=1 to nrow DO
	    BEGIN
		  a[k,l] := a[k,l]+x[i,l]*x[i,k];
		  IF k<>l THEN a[l,k] := a[k,l]
	    END
	 END;		{ l-loop }
	 g[k] := 0.0;
	 FOR i:=1 to nrow DO
	 g[k] := g[k]+y[i]*x[i,k]
  END	{ k-loop }
END;	{ square }

{ ************************************************************************ }

PROCEDURE trapez(
		lower,upper,tol: Real;
		VAR sum		: Real);

{ numerical integration by the trapezoid method }
{ function is f (as parameter), limits are LOWER and UPPER }
{ with number of regions equal to PIECES }
{ fixed partition is DELTA_X, answer is SUM }

VAR	pieces,i			: Integer;
	x,delta_x,end_sum,mid_sum,
	end_cor,sum1			: Real;

BEGIN
  pieces := 1;
  delta_x := (upper-lower)/pieces;
  end_sum := fx(lower)+fx(upper);
  end_cor := (dfx(upper)-dfx(lower))/12.0;
  sum := end_sum*delta_x/2.0;
  writeln('    1',sum);
  mid_sum := 0.0;
  REPEAT
    pieces := pieces*2;
    sum1 := sum;
    delta_x := (upper-lower)/pieces;
    FOR i:=1 to pieces div 2 DO
    BEGIN
      x := lower+delta_x*(2.0*i-1.0);
      mid_sum := mid_sum+fx(x)
    END;
  sum := (end_sum+2.0*mid_sum)*delta_x*0.5-sqr(delta_x)*end_cor;
  writeln(pieces:5,sum)
  UNTIL abs(sum-sum1)<=abs(tol*sum)
END;		{ TRAPEZ }

{ ************************************************************************ }

FUNCTION gamma(x: Real): Real;
CONST	pi	= 3.1415926;

VAR	i,j	: Integer;
	y,gam	: Real;

BEGIN		{ gamma function }
  IF x>=0.0 THEN
    BEGIN
      y := x+2.0;
      gam := sqrt(2*pi/y)*exp(y*ln(y)+(1-1/(30*y*y))/(12*y)-y);
      gamma := gam/(x*(x+1))
    END
  ELSE		{ x<0 }
    BEGIN
      j := 0;
      y := x;
      REPEAT
	j := j+1;
	y := y+1.0
      UNTIL y>0.0;
      gam := gamma(y);		{ recursive call }
      FOR i:=0 to j-1 DO
	gam := gam/(x+1);
      gamma := gam
    END	{ x<0 }
END;		{ gamma function }

{ ************************************************************************ }

FUNCTION bessj(x,n: Real): Real;
{ cylindrical Bessel function of the first kind }
{ the gamma function is required }

CONST	tol	= 1.0E-4;
	pi	= 3.1415926;

VAR	i		: Integer;
	term,new_term,
	sum,x2		: Real;

BEGIN	{ bessj }
  x2 := x*x;
  IF (x=0.0) AND (N=1.0) THEN bessj := 0.0
  ELSE IF x>15 THEN { asymptotic expansion }
    bessj := sqrt(2/(pi*x))*cos(x-pi/4-n*pi/2)
  ELSE
    BEGIN
      IF n=0.0 THEN sum := 1.0
      ELSE sum := exp(n*ln(x/2))/gamma(n+1.0);
      new_term := sum;
      i := 0;
  REPEAT
    i := i+1;
    term := new_term;
    new_term := -term*x2*0.25/(i*(n+1));
    sum := sum+new_term
  UNTIL abs(new_term)<=abs(sum*tol);
  bessj := sum
 END	{ if}
END;	{ bessj }

{ ************************************************************************ }

PROCEDURE romb(lower, upper, tol: Real; VAR ans: Real);
{ numerical integration by the Romberg method }
{ function fx, name cannot be passed by Turbo Pascal}
VAR
	nx			: ARRAY[1..16] of Integer;
	t			: ARRAY[1..136] of Real;
	done,error		: Boolean;
	pieces,nt,i,ii,n,nn,
	l,ntra,k,m,j		: Integer ;
	delta_x,c,sum,fotom,x	: Real ;
BEGIN
  done := false;
  error := false;
  pieces := 1;
  nx[1] := 1;
  delta_x := (upper-lower)/pieces;
  c := (fx(lower)+fx(upper))*0.5;
  t[1] := delta_x*c;
  n := 1;
  nn := 2;
  sum := c;
  REPEAT
    n := n+1;
    fotom := 4.0;
    nx[n] := nn;
    pieces := pieces*2;
    l := pieces-1;
    delta_x := (upper-lower)/pieces;
    { compute trapezoidal sum for 2^(n-1)+1 points }
    FOR ii:=1 to (l+1) div 2 DO
      BEGIN
	i := ii*2-1;
	x := lower+i*delta_x;
	sum := sum+fx(x)
      END;
    t[nn] := delta_x*sum;
    write(pieces:5,t[nn]);
    ntra := nx[n-1];
    k := n-1;
    { compute n-th row of T array }
    FOR m:=1 to k DO
      BEGIN
	j := nn+m;
	nt := nx[n-1]+m-1;
	t[j] := (fotom*t[j-1]-t[nt])/(fotom-1.0);
	fotom := fotom*4.0
      END;
    writeln(j:4,t[j]);
    IF n>4 THEN
      BEGIN
	IF t[nn+1]<>0.0 THEN
	  IF (abs(t[ntra+1]-t[nn+1])<=abs(t[nn+1]*tol))
	    OR (abs(t[nn-1]-t[j])<=abs(t[j]*tol)) THEN
	      done := true
	  ELSE IF n>15 THEN
	    BEGIN
	      done := true;
	      error := true
	    END
	END;	{ if n>4 }
    nn := j+1
  UNTIL done;
  ans := t[j]
END;		{ ROMBERG }

{ ************************************************************************ }

{ ************************************************************************ }

{ ************************************************************************ }

{ ************************************************************************ }

{ ************************************************************************ }

BEGIN
END.