PROGRAM d8r10(input,output);
(* driver for routine ECLASS *)
CONST
   m=11;
   n=15;
TYPE
   glnarray = ARRAY [1..n] OF integer;
   glmarray = ARRAY [1..m] OF integer;
VAR
   i,j,k,lclas,nclass : integer;
   lista,listb : glmarray;
   nf,nflag,nsav : glnarray;

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

BEGIN
   lista[1] := 1; lista[2] := 1; lista[3] := 5; lista[4] := 2;
   lista[5] := 6; lista[6] := 2; lista[7] := 7; lista[8] := 11;
   lista[9] := 3; lista[10] := 4; lista[11] := 12;
   listb[1] := 5; listb[2] := 9; listb[3] := 13; listb[4] := 6;
   listb[5] := 10; listb[6] := 14; listb[7] := 3; listb[8] := 7;
   listb[9] := 15; listb[10] := 8; listb[11] := 4;
   eclass(nf,n,lista,listb,m);
   FOR i := 1 to n DO BEGIN
      nflag[i] := 1
   END;
   writeln;
   writeln('Numbers from 1-15 divided according to');
   writeln('their value modulo 4:');
   writeln;
   lclas := 0;
   FOR i := 1 to n DO BEGIN
      nclass := nf[i];
      IF (nflag[nclass] <> 0)  THEN BEGIN
         nflag[nclass] := 0;
         lclas := lclas+1;
         k := 0;
         FOR j := i to n DO BEGIN
            IF  (nf[j] = nf[i])  THEN BEGIN
               k := k+1;
               nsav[k] := j
            END
         END;
         write('Class',lclas:2,':   ');
         FOR j := 1 to k DO write(nsav[j]:3);
         writeln
      END
   END
END.
PROGRAM d8r11(input,output);
(* driver for routine ECLAZZ *)
CONST
   n=15;
TYPE
   glnarray = ARRAY [1..n] OF integer;
VAR
   i,j,k,lclas,nclass : integer;
   nf,nflag,nsav : glnarray;

FUNCTION equiv(i,j: integer): boolean;
BEGIN
   equiv := FALSE;
   IF ((i MOD 4) = (j MOD 4)) THEN  equiv := TRUE
END;

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

BEGIN
   eclazz(nf,n);
   FOR i := 1 to n DO BEGIN
      nflag[i] := 1
   END;
   writeln;
   writeln('Numbers from 1-15 divided according to');
   writeln('their value modulo 4:');
   lclas := 0;
   FOR i := 1 to n DO BEGIN
      nclass := nf[i];
      IF (nflag[nclass] <> 0)  THEN BEGIN
         nflag[nclass] := 0;
         lclas := lclas+1;
         k := 0;
         FOR j := i to n DO BEGIN
            IF (nf[j] = nf[i]) THEN BEGIN
               k := k+1;
               nsav[k] := j
            END
         END;
         write('Class',lclas:2,':   ');
         FOR j := 1 to k DO write(nsav[j]:3);
         writeln
      END
   END
END.
PROGRAM d11r2(input,output);
(* driver for routine EIGSRT *)
CONST
   np=10;
TYPE
   glnpnp = ARRAY [1..np,1..np] OF real;
   glnp = ARRAY [1..np] OF real;
VAR
   i,j,nrot : integer;
   c,v : glnpnp;
   d : glnp;

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

(*$I EIGSRT.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;
   jacobi(c,np,d,v,nrot);
   writeln('unsorted eigenvectors:');
   FOR i := 1 to np DO BEGIN
      writeln('eigenvalue',i:3,'  := ',d[i]:12:6);
      writeln('eigenvector:');
      FOR j := 1 to np DO BEGIN
         write(v[j,i]:12:6);
         IF ((j MOD 5) = 0) THEN writeln
      END;
      writeln
   END;
   writeln;
   writeln('****** sorting ******');
   writeln;
   eigsrt(d,v,np);
   writeln('sorted eigenvectors:');
   FOR i := 1 to np DO BEGIN
      writeln('eigenvalue',i:3,'  := ',d[i]:12:6);
      writeln('eigenvector:');
      FOR j := 1 to np DO BEGIN
         write(v[j,i]:12:6);
         IF ((j MOD 5) = 0) THEN writeln
      END;
      writeln
   END
END.
PROGRAM d6r28(input,output);
(* driver for routine EL2 *)
VAR
   glit : integer;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gla,glb,glakc : real;
   ago,astop,s,x : real;
   i,idum : integer;

FUNCTION func(phi: real): real;
(* Programs using routine FUNC must declare the variables
VAR
   glakc,gla,glb : real;
in the main routine. *)
VAR
   tn,tsq : real;
BEGIN
   tn := sin(phi)/cos(phi);
   tsq := tn*tn;
   func := (gla+glb*tsq)/sqrt((1.0+tsq)*(1.0+glakc*glakc*tsq))
END;

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

(*$I TRAPZD.PAS *)

(*$I QSIMP.PAS *)

(*$I EL2.PAS *)

BEGIN
   writeln('general elliptic integral of second kind');
   writeln('x':7,'kc':10,'a':10,'b':10,'el2':11,'integral':12);
   idum := -55;
   ago := 0.0;
   FOR i := 1 to 20 DO BEGIN
      glakc := 5.0*ran3(idum);
      gla := 10.0*ran3(idum);
      glb := 10.0*ran3(idum);
      x := 10.0*ran3(idum);
      astop := arctan(x);
      qsimp(ago,astop,s);
      writeln(x:10:6,glakc:10:6,gla:10:6,glb:10:6,
         el2(x,glakc,gla,glb):10:6,s:10:6)
   END
END.
PROGRAM d11r6(input,output);
(* driver for routine ELMHES *)
CONST
   np=5;
TYPE
   glnp = ARRAY [1..np] OF real;
   glnpnp = ARRAY [1..np,1..np] OF real;
VAR
   i,j : integer;
   a : glnpnp;
   c,r : glnp;

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

(*$I ELMHES.PAS *)

BEGIN
   a[1,1] := 1.0; a[1,2] := 2.0; a[1,3] := 300.0;
   a[1,4] := 4.0; a[1,5] := 5.0;
   a[2,1] := 2.0; a[2,2] := 3.0; a[2,3] := 400.0;
   a[2,4] := 5.0; a[2,5] := 6.0;
   a[3,1] := 3.0; a[3,2] := 4.0; a[3,3] := 5.0;
   a[3,4] := 6.0; a[3,5] := 7.0;
   a[4,1] := 4.0; a[4,2] := 5.0; a[4,3] := 600.0;
   a[4,4] := 7.0; a[4,5] := 8.0;
   a[5,1] := 5.0; a[5,2] := 6.0; a[5,3] := 700.0;
   a[5,4] := 8.0; a[5,5] := 9.0;
   writeln('***** original matrix *****');
   FOR i := 1 to np DO BEGIN
      FOR j := 1 to np DO write(a[i,j]:12:2);
      writeln
   END;
   writeln('***** balance matrix *****');
   balanc(a,np);
   FOR i := 1 to np DO BEGIN
      FOR j := 1 to np DO write(a[i,j]:12:2);
      writeln
   END;
   writeln('***** reduce to hessenberg form *****');
   elmhes(a,np);
   FOR j := 1 to np-2 DO BEGIN
      FOR i := j+2 to np DO BEGIN
         a[i,j] := 0.0
      END
   END;
   FOR i := 1 to np DO BEGIN
      FOR j := 1 to np DO write('  ',a[i,j]:11);
      writeln
   END
END.
PROGRAM d6r10(input,output,dfile);
(* driver for routine ERF *)
VAR 
   i,nval : integer;
   val,x : real;
   txt : string[14];
   dfile : text;

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

(*$I GSER.PAS *)

(*$I GCF.PAS *)

(*$I GAMMP.PAS *)

(*$I ERF.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Error Function');
   readln(dfile,nval);
   writeln(txt);
   writeln('x':5,'actual':12,'erf(x)':12);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,val);
      writeln(x:6:2,val:12:7,erf(x):12:7)
   END;
   close(dfile)
END.
PROGRAM d6r11(input,output,dfile);
(* driver for routine ERFC *)
VAR
   i,nval : integer;
   val,x : real;
   txt : string[14];
   dfile : text;

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

(*$I GCF.PAS *)

(*$I GSER.PAS *)

(*$I GAMMQ.PAS *)

(*$I GAMMP.PAS *)

(*$I ERFC.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Error Function');
   readln(dfile,nval);
   writeln ('complementary error function');
   writeln ('x':5,'actual':12,'erfc(x)':12);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,val);
      val := 1.0-val;
      writeln(x:6:2,val:12:7,erfc(x):12:7)
   END;
   close(dfile);
END.
PROGRAM d6r12(input,output,dfile);
(* driver for routine erfcc *)
VAR
   i,nval : integer;
   x,val : real;
   txt : string[14];
   dfile : text;

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

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Error Function');
   readln(dfile,nval);
   writeln ('complementary error function');
   writeln ('x':5,'actual':12,'erfcc(x)':13);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,val);
      val := 1.0-val;
      writeln (x:6:2,val:12:7,erfcc(x):12:7)
   END;
   close(dfile)
END.
PROGRAM d5r1(input,output);
(* driver for routine EULSUM *)
LABEL 10,99;
CONST
   nval=40;
VAR
   i,j,mval,glnterm : integer;
   sum,term,x,xpower : real;
   glwksp : ARRAY [1..nval] OF real;

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

BEGIN
(* evaluate ln(1+x) := x-x^2/2+x^3/3-x^4/4 ... FOR -1<x<1 *)
10:   writeln;
   writeln ('How many terms in polynomial?');
   writeln ('Enter n between 1 and ',nval:2,
            '. (n := 0 to END)');
   readln (mval);
   writeln;
   IF ((mval <= 0) OR (mval > nval)) THEN GOTO 99;
   writeln ('x':9,'actual':14,'polynomial':14);
   FOR i := -8 to 8 DO BEGIN
      x := i/10.0;
      sum := 0.0;
      xpower := -1;
      FOR j := 1 to mval DO BEGIN
         xpower := -x*xpower;
         term := xpower/j;
         eulsum(sum,term,j)
      END;
      writeln (x:12:6,ln(1.0+x):12:6,sum:12:6)
   END;
   GOTO 10;
99:
END.
PROGRAM d12r10(input,output,dfile);
(* driver for routine EVLMEM *)
CONST
   n=1000;
   m=10;
   nfdt=16;
TYPE
   glnarray = ARRAY [1..n] OF real;
   glmarray = ARRAY [1..m] OF real;
VAR
   i : integer;
   fdt,pm : real;
   cof,wkm : glmarray;
   data,wk1,wk2 : glnarray;
   dfile : text;

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

(*$I EVLMEM.PAS *)

BEGIN
   glopen(dfile,'spctrl.dat');
   FOR i := 1 to n DO BEGIN
      read(dfile,data[i])
   END;
   close(dfile);
   memcof(data,n,m,pm,cof,wk1,wk2,wkm);
   writeln('Power spectum estimate of data in SPCTRL.DAT');
   writeln(' ':4,'f*delta',' ':7,'power');
   FOR i := 0 to nfdt DO BEGIN
      fdt := 0.5*i/nfdt;
      writeln (fdt:12:6,evlmem(fdt,cof,m,pm):12:6)
   END
END.
PROGRAM d7r5 (input,output);
(* driver for routine EXPDEV *)
CONST
   npts=1000;
   ee=2.718281828;
VAR
   i,idum,j : integer;
   expect,total,y : real;
   trig,x : ARRAY [0..20] OF real;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;

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

(*$I EXPDEV.PAS *)

BEGIN
   FOR i := 0 to 20 DO BEGIN
      trig[i] := i/20.0;
      x[i] := 0.0
   END;
   idum := -1;
   FOR i := 1 to npts DO BEGIN
      y := expdev(idum);
      FOR j := 1 to 20 DO BEGIN
         IF ((y < trig[j]) AND (y > trig[j-1]))  THEN BEGIN
            x[j] := x[j]+1.0
         END
      END
   END;
   total := 0.0;
   FOR i := 1 to 20 DO BEGIN
      total := total+x[i]
   END;
   writeln;
   writeln ('exponential distribution with',npts:7,' points');
   writeln ('   interval','     observed','    expected');
   writeln;
   FOR i := 1 to 20 DO BEGIN
      x[i] := x[i]/total;
      expect := exp(-(trig[i-1]+trig[i])/2.0);
      expect := expect*0.05*ee/(ee-1);
      writeln (trig[i-1]:6:2,trig[i]:6:2,x[i]:12:6,expect:12:6)
   END
END.
