PROGRAM d2r7(input,output);
(* driver for routine TOEPLZ *)
CONST
   n=5;
   twon=10;    (* twon=2*n *)
TYPE
   glnarray = ARRAY [1..n] OF real;
   gltwon = ARRAY [1..twon] OF real;
VAR
   i,j : integer;
   sum : real;
   r : gltwon;
   x,y : glnarray;

(*$I MODFILE.PAS *)
(*$I TOEPLZ.PAS *)

BEGIN
   FOR i := 1 to n DO BEGIN
      y[i] := 0.1*i
   END;
   FOR i := 1 to 2*n-1 DO BEGIN
      r[i] := 0.1*i
   END;
   toeplz(r,x,y,n);
   writeln('Solution vector:');
   FOR i := 1 to n DO BEGIN
      writeln('x[':7,i:1,'] :=',x[i]:13)
   END;
   writeln;
   writeln('Test of solution:');
   writeln('mtrx*soln':13,'original':12);
   FOR i := 1 to n DO BEGIN
      sum := 0.0;
      FOR j := 1 to n DO BEGIN
         sum := sum+r[n+i-j]*x[j]
      END;
      writeln(sum:12:4,y[i]:12:4)
   END
END.
PROGRAM d13r6(input,output);
(* driver for routine TPTEST *)
(* compare two correlated distributions vs. two *)
(* uncorrelated distributions *)
CONST
   npts=500;
   eps=0.01;
   nshft=10;
   anoise=0.3;
TYPE
   glnparray = ARRAY [1..npts] OF real;
VAR
   gliset : integer;
   glgset : real;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   ave1,ave2,ave3,gauss : real;
   offset,prob1,prob2,shift,t1,t2: real;
   var1,var2,var3 : real;
   i,idum,j : integer;
   data1,data2,data3 : glnparray;

(*$I MODFILE.PAS *)
(*$I RAN3.PAS *)

(*$I GASDEV.PAS *)

(*$I GAMMLN.PAS *)

(*$I BETACF.PAS *)

(*$I BETAI.PAS *)

(*$I AVEVAR.PAS *)

(*$I TPTEST.PAS *)

BEGIN
   gliset := 0;
   idum := -5;
   writeln('Correlated:':29,'Uncorrelated:':30);
   writeln('Shift':7,'t':11,'Probability':17,'t':11,'Probability':17);
   offset := (nshft DIV 2)*eps;
   FOR j := 1 to npts DO BEGIN
      gauss := gasdev(idum);
      data1[j] := gauss;
      data2[j] := gauss+anoise*gasdev(idum);
      data3[j] := gasdev(idum)+anoise*gasdev(idum)
   END;
   avevar(data1,npts,ave1,var1);
   avevar(data2,npts,ave2,var2);
   avevar(data3,npts,ave3,var3);
   FOR j := 1 to npts DO BEGIN
      data1[j] := data1[j]-ave1+offset;
      data2[j] := data2[j]-ave2;
      data3[j] := data3[j]-ave3
   END;
   FOR i := 1 to nshft DO BEGIN
      shift := i*eps;
      FOR j := 1 to npts DO BEGIN
         data2[j] := data2[j]+eps;
         data3[j] := data3[j]+eps
      END;
      tptest(data1,data2,npts,t1,prob1);
      tptest(data1,data3,npts,t2,prob2);
      writeln(shift:6:2,t1:14:4,prob1:12:4,t2:16:4,prob2:12:4)
   END
END.
PROGRAM d11r4(input,output);
(* driver for routine TQLI *)
CONST
   np=10;
   tiny=1.0e-6;
TYPE
   glnp = ARRAY [1..np] OF real;
   glnpnp = ARRAY [1..np,1..np] OF real;
VAR
   i,j,k : integer;
   a,c : glnpnp;
   d,e,f : glnp;

(*$I MODFILE.PAS *)
(*$I TRED2.PAS *)

(*$I TQLI.PAS *)

BEGIN
   c[1,1] := 5.0; c[1,2] := 4.0; c[1,3] := 3.0;
   c[1,4] := 2.0; c[1,5] := 1.0; c[1,6] := 0.0;
   c[1,7] := -1.0; c[1,8] := -2.0; c[1,9] := -3.0;
   c[1,10] := -4.0;
   c[2,1] := 4.0; c[2,2] := 5.0; c[2,3] := 4.0;
   c[2,4] := 3.0; c[2,5] := 2.0; c[2,6] := 1.0;
   c[2,7] := 0.0; c[2,8] := -1.0; c[2,9] := -2.0;
   c[2,10] := -3.0;
   c[3,1] := 3.0; c[3,2] := 4.0; c[3,3] := 5.0;
   c[3,4] := 4.0; c[3,5] := 3.0; c[3,6] := 2.0;
   c[3,7] := 1.0; c[3,8] := 0.0; c[3,9] := -1.0;
   c[3,10] := -2.0;
   c[4,1] := 2.0; c[4,2] := 3.0; c[4,3] := 4.0;
   c[4,4] := 5.0; c[4,5] := 4.0; c[4,6] := 3.0;
   c[4,7] := 2.0; c[4,8] := 1.0; c[4,9] := 0.0;
   c[4,10] := -1.0;
   c[5,1] := 1.0; c[5,2] := 2.0; c[5,3] := 3.0;
   c[5,4] := 4.0; c[5,5] := 5.0; c[5,6] := 4.0;
   c[5,7] := 3.0; c[5,8] := 2.0; c[5,9] := 1.0;
   c[5,10] := 0.0;
   c[6,1] := 0.0; c[6,2] := 1.0; c[6,3] := 2.0;
   c[6,4] := 3.0; c[6,5] := 4.0; c[6,6] := 5.0;
   c[6,7] := 4.0; c[6,8] := 3.0; c[6,9] := 2.0;
   c[6,10] := 1.0;
   c[7,1] := -1.0; c[7,2] := 0.0; c[7,3] := 1.0;
   c[7,4] := 2.0; c[7,5] := 3.0; c[7,6] := 4.0;
   c[7,7] := 5.0; c[7,8] := 4.0; c[7,9] := 3.0;
   c[7,10] := 2.0;
   c[8,1] := -2.0; c[8,2] := -1.0; c[8,3] := 0.0;
   c[8,4] := 1.0; c[8,5] := 2.0; c[8,6] := 3.0;
   c[8,7] := 4.0; c[8,8] := 5.0; c[8,9] := 4.0;
   c[8,10] := 3.0;
   c[9,1] := -3.0; c[9,2] := -2.0; c[9,3] := -1.0;
   c[9,4] := 0.0; c[9,5] := 1.0; c[9,6] := 2.0;
   c[9,7] := 3.0; c[9,8] := 4.0; c[9,9] := 5.0;
   c[9,10] := 4.0;
   c[10,1] := -4.0; c[10,2] := -3.0; c[10,3] := -2.0;
   c[10,4] := -1.0; c[10,5] := 0.0; c[10,6] := 1.0;
   c[10,7] := 2.0; c[10,8] := 3.0; c[10,9] := 4.0;
   c[10,10] := 5.0;
   FOR i := 1 to np DO BEGIN
      FOR j := 1 to np DO BEGIN
         a[i,j] := c[i,j]
      END
   END;
   tred2(a,np,d,e);
   tqli(d,e,np,a);
   writeln('eigenvectors for a real symmetric matrix');
   FOR i := 1 to np DO BEGIN
      FOR j := 1 to np DO BEGIN
         f[j] := 0.0;
         FOR k := 1 to np DO BEGIN
            f[j] := f[j]+c[j,k]*a[k,i]
         END
      END;
      writeln('eigenvalue',i:3,'  := ',d[i]:10:6);
      writeln('vector':11,'mtrx*vect.':14,'ratio':9);
      FOR j := 1 to np DO BEGIN
         IF (abs(a[j,i]) < tiny) THEN BEGIN
            writeln(a[j,i]:12:6,f[j]:12:6,'div. by 0':12)
         END ELSE BEGIN
            writeln(a[j,i]:12:6,f[j]:12:6,f[j]/a[j,i]:12:6)
         END
      END;
   writeln('press enter to continue...');
   readln
   END
END.
PROGRAM d4r1(input,output);
(* driver for routine trapzd *)
CONST
   nmax=14;
   pio2=1.5707963;
VAR
   glit : integer;
   i : integer;
   a,b,s : real;

FUNCTION func(x : real) : real;
(* Test function *)
BEGIN
   func := sqr(x)*(sqr(x)-2.0)*sin(x)
END;

FUNCTION fint(x : real) : real;
(* Integral of test function *)
BEGIN
   fint := 4.0*x*(sqr(x)-7.0)*sin(x)-
         (sqr(sqr(x))-14.0*sqr(x)+28.0)*cos(x);
END;

(*$I MODFILE.PAS *)
(*$I TRAPZD.PAS *)

BEGIN
   a := 0.0;
   b := pio2;
   writeln ('integral of func with 2^(n-1) points');
   writeln ('actual value of integral is',fint(b)-fint(a):12:6);
   writeln ('n':6,'approx. integral':24);
   FOR i := 1 to nmax DO BEGIN
      trapzd(a,b,s,i);
      writeln (i:6,s:20:6)
   END
 END.
PROGRAM d11r3(input,output);
(* driver for routine TRED2 *)
CONST
   np=10;
TYPE
   glnp = ARRAY [1..np] OF real;
   glnpnp = ARRAY [1..np,1..np] OF real;
VAR
   i,j,k,l,m : integer;
   a,c,f : glnpnp;
   d,e : glnp;

(*$I MODFILE.PAS *)
(*$I TRED2.PAS *)

BEGIN
   c[1,1] := 5.0; c[1,2] := 4.0; c[1,3] := 3.0;
   c[1,4] := 2.0; c[1,5] := 1.0; c[1,6] := 0.0;
   c[1,7] := -1.0; c[1,8] := -2.0; c[1,9] := -3.0;
   c[1,10] := -4.0;
   c[2,1] := 4.0; c[2,2] := 5.0; c[2,3] := 4.0;
   c[2,4] := 3.0; c[2,5] := 2.0; c[2,6] := 1.0;
   c[2,7] := 0.0; c[2,8] := -1.0; c[2,9] := -2.0;
   c[2,10] := -3.0;
   c[3,1] := 3.0; c[3,2] := 4.0; c[3,3] := 5.0;
   c[3,4] := 4.0; c[3,5] := 3.0; c[3,6] := 2.0;
   c[3,7] := 1.0; c[3,8] := 0.0; c[3,9] := -1.0;
   c[3,10] := -2.0;
   c[4,1] := 2.0; c[4,2] := 3.0; c[4,3] := 4.0;
   c[4,4] := 5.0; c[4,5] := 4.0; c[4,6] := 3.0;
   c[4,7] := 2.0; c[4,8] := 1.0; c[4,9] := 0.0;
   c[4,10] := -1.0;
   c[5,1] := 1.0; c[5,2] := 2.0; c[5,3] := 3.0;
   c[5,4] := 4.0; c[5,5] := 5.0; c[5,6] := 4.0;
   c[5,7] := 3.0; c[5,8] := 2.0; c[5,9] := 1.0;
   c[5,10] := 0.0;
   c[6,1] := 0.0; c[6,2] := 1.0; c[6,3] := 2.0;
   c[6,4] := 3.0; c[6,5] := 4.0; c[6,6] := 5.0;
   c[6,7] := 4.0; c[6,8] := 3.0; c[6,9] := 2.0;
   c[6,10] := 1.0;
   c[7,1] := -1.0; c[7,2] := 0.0; c[7,3] := 1.0;
   c[7,4] := 2.0; c[7,5] := 3.0; c[7,6] := 4.0;
   c[7,7] := 5.0; c[7,8] := 4.0; c[7,9] := 3.0;
   c[7,10] := 2.0;
   c[8,1] := -2.0; c[8,2] := -1.0; c[8,3] := 0.0;
   c[8,4] := 1.0; c[8,5] := 2.0; c[8,6] := 3.0;
   c[8,7] := 4.0; c[8,8] := 5.0; c[8,9] := 4.0;
   c[8,10] := 3.0;
   c[9,1] := -3.0; c[9,2] := -2.0; c[9,3] := -1.0;
   c[9,4] := 0.0; c[9,5] := 1.0; c[9,6] := 2.0;
   c[9,7] := 3.0; c[9,8] := 4.0; c[9,9] := 5.0;
   c[9,10] := 4.0;
   c[10,1] := -4.0; c[10,2] := -3.0; c[10,3] := -2.0;
   c[10,4] := -1.0; c[10,5] := 0.0; c[10,6] := 1.0;
   c[10,7] := 2.0; c[10,8] := 3.0; c[10,9] := 4.0;
   c[10,10] := 5.0;
   FOR i := 1 to np DO BEGIN
      FOR j := 1 to np DO BEGIN
         a[i,j] := c[i,j]
      END
   END;
   tred2(a,np,d,e);
   writeln('diagonal elements');
   FOR i := 1 to np DO BEGIN
      write(d[i]:12:6);
      IF ((i MOD 5) = 0) THEN writeln
   END;
   writeln('off-diagonal elements');
   FOR i := 2 to np DO BEGIN
      write(e[i]:12:6);
      IF ((i MOD 5) = 0) THEN writeln
   END;
(* check transformation matrix *)
   FOR j := 1 to np DO BEGIN
      FOR k := 1 to np DO BEGIN
      f[j,k] := 0.0;
         FOR l := 1 to np DO BEGIN
            FOR m := 1 to np DO BEGIN
               f[j,k] := f[j,k]
               +a[l,j]*c[l,m]*a[m,k]
            END
         END
      END
   END;
(* how does it look? *)
   writeln('tridiagonal matrix');
   FOR i := 1 to np DO BEGIN
      FOR j := 1 to np DO write(f[i,j]:7:2);
      writeln
   END
END.
PROGRAM d2r4(input,output,dfile);
(* driver for routine TRIDAG *)
LABEL 10,99;
CONST
   np=20;
TYPE
   glnarray = ARRAY [1..np] OF real;
VAR
   k,n: integer;
   diag,superd,subd,rhs,u: glnarray;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I TRIDAG.PAS *)

BEGIN
   glopen(dfile,'matrx2.dat');
10:   readln(dfile);
   readln(dfile);
   readln(dfile,n);
   readln(dfile);
   FOR k := 1 to n-1 DO read(dfile,diag[k]);
   readln(dfile,diag[n]);
   readln(dfile);
   FOR k := 1 to n-2 DO read(dfile,superd[k]);
   readln(dfile,superd[n-1]);
   readln(dfile);
   FOR k := 2 to n-1 DO read(dfile,subd[k]);
   readln(dfile,subd[n]);
   readln(dfile);
   FOR k := 1 to n-1 DO read(dfile,rhs[k]);
   readln(dfile,rhs[n]);
(* carry out solution *)
   tridag(subd,diag,superd,rhs,u,n);
   writeln ('the solution vector is:');
   FOR k := 1 to n-1 DO write(u[k]:12:6);
   writeln(u[n]:12:6);
(* test solution *)
   writeln ('(matrix)*(sol''n vector) should be:');
   FOR k := 1 to n-1 DO write(rhs[k]:12:6);
   writeln(rhs[n]:12:6);
   writeln ('actual result is:');
   FOR k := 1 to n DO BEGIN
      IF (k = 1) THEN BEGIN
         rhs[k] := diag[1]*u[1] + superd[1]*u[2]
      END ELSE IF (k = n) THEN BEGIN
         rhs[k] := subd[n]*u[n-1] + diag[n]*u[n]
      END ELSE BEGIN
         rhs[k] := subd[k]*u[k-1] + diag[k]*u[k]
         + superd[k]*u[k+1]
      END
   END;
   FOR k := 1 to n-1 DO write(rhs[k]:12:6);
   writeln(rhs[n]:12:6);
   writeln ('***********************************');
   IF eof(dfile) THEN GOTO 99;
   writeln ('press return for next problem:');
   readln;
   GOTO 10;
99:   close(dfile)
END.
PROGRAM d13r3(input,output);
(* driver for routine TTEST *)
(* generate gaussian distributed data *)
CONST
   npts=1024;   (* npts must be > mpts for TTEST.PAS *)
   mpts=512;
   eps=0.03;
   nshft=10;
TYPE
   glnarray = ARRAY [1..npts] OF real;
   glmarray = ARRAY [1..mpts] OF real;
   glnparray = glnarray;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   data1 : glnarray;
   data2 : glmarray;
   i,idum,j : integer;
   prob,t : real;

(*$I MODFILE.PAS *)
(*$I AVEVAR.PAS *)

(*$I GAMMLN.PAS *)

(*$I BETACF.PAS *)

(*$I BETAI.PAS *)

(*$I RAN3.PAS *)

(*$I GASDEV.PAS *)

(*$I TTEST.PAS *)

BEGIN
   gliset := 0;
   idum := -11;
   FOR i := 1 to npts DO data1[i] := gasdev(idum);
   FOR i := 1 to mpts DO data2[i] := (nshft DIV 2)*eps+gasdev(idum);
   writeln('shift':6,'t':8,'probability':16);
   FOR i := 1 to (nshft+1) DO BEGIN
      ttest(data1,npts,data2,mpts,t,prob);
      writeln((i-1)*eps:6:2,t:10:2,prob:10:2);
      FOR j := 1 to npts DO data1[j] := data1[j]+eps
   END
END.
PROGRAM d13r5(input,output);
(* driver for routine TUTEST *)
CONST
   npts=3000;   (* npts must be > mpts for TUTEST *)
   mpts=600;
   eps=0.03;
   var1=1.0;
   var2=4.0;
   nshft=10;
TYPE
   glnarray = ARRAY [1..npts] OF real;
   glmarray = ARRAY [1..mpts] OF real;
   glnparray = glnarray;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   fctr1,fctr2,prob,t : real;
   i,idum,j : integer;
   data1 : glnarray;
   data2 : glmarray;

(*$I MODFILE.PAS *)
(*$I AVEVAR.PAS *)

(*$I GAMMLN.PAS *)

(*$I BETACF.PAS *)

(*$I BETAI.PAS *)

(*$I RAN3.PAS *)

(*$I GASDEV.PAS *)

(*$I TUTEST.PAS *)

BEGIN
(* generate two gaussian distributions of different variance *)
   gliset := 0;
   idum := -1773;
   fctr1 := sqrt(var1);
   FOR i := 1 to npts DO data1[i] := fctr1*gasdev(idum);
   fctr2 := sqrt(var2);
   FOR i := 1 to mpts DO data2[i] := (nshft DIV 2)*eps+fctr2*gasdev(idum);
   writeln;
   writeln('Distribution #1 : variance = ',var1:6:2);
   writeln('Distribution #2 : variance = ',var2:6:2);
   writeln;
   writeln('shift':7,'t':8,'probability':16);
   FOR i := 1 to (nshft+1) DO BEGIN
      tutest(data1,npts,data2,mpts,t,prob);
      writeln((i-1)*eps:6:2,t:10:2,prob:11:2);
      FOR j := 1 to npts DO data1[j] := data1[j]+eps
   END
END.
PROGRAM d12r2(input,output);
(* driver for routine TWOFFT *)
CONST
   n=32;
   n2=64;   (* n2=2*n *)
   per=8;
   pi=3.1415926;
TYPE
   glnarray = ARRAY [1..n] OF real;
   gl2narray = ARRAY [1..n2] OF real;
   gldarray = gl2narray;
VAR
   i,isign : integer;
   data1,data2 : glnarray;
   fft1,fft2 : gl2narray;

PROCEDURE prntft(data : gldarray; nn : integer);
VAR
   ii,mm,n : integer;
BEGIN
   writeln('n':4,'real(n)':13,'imag.(n)':13,'real(N-n)':12,'imag.(N-n)':13);
   writeln (0:4,data[1]:14:6,data[2]:12:6,data[1]:12:6,data[2]:12:6);
   mm := nn DIV 2;
   FOR ii := 1 to mm DO BEGIN
      n := 2*ii+1;
      writeln (((n-1) DIV 2):4,data[n]:14:6,data[n+1]:12:6,
         data[2*nn+2-n]:12:6,data[2*nn+3-n]:12:6)
   END;
   writeln (' press return to continue ...');
   readln
END;

(*$I MODFILE.PAS *)
(*$I FOUR1.PAS *)

(*$I TWOFFT.PAS *)

BEGIN
   FOR i := 1 to n DO BEGIN
      data1[i] := round(cos(i*2.0*pi/per));
      data2[i] := round(sin(i*2.0*pi/per));
   END;
   twofft(data1,data2,fft1,fft2,n);
   writeln ('fourier transform of first function:');
   prntft(fft1,n);
   writeln ('fourier transform of second function:');
   prntft(fft2,n);
(* invert transform *)
   isign := -1;
   four1(fft1,n,isign);
   writeln ('inverted transform  =  first function:');
   prntft(fft1,n);
   four1(fft2,n,isign);
   writeln ('inverted transform  =  second function:');
   prntft(fft2,n)
END.
