PROCEDURE toeplz(r: gltwon; VAR x: glnarray; y: glnarray; n: integer);
(* Programs using routine TOEPLZ must define the types
   glnarray = ARRAY [1..n] OF real;
   gltwon = ARRAY [1..2*n] OF real;
in the main routine. *)
LABEL 98,99;
VAR
   m2,m1,m,k,j: integer;
   sxn,shn,sgn,sgd,sd,qt2,qt1,qq,pt2,pt1,pp: real;
   g,h: glnarray;
BEGIN
   IF (r[n] = 0.0) THEN  GOTO 99;
   x[1] := y[1]/r[n];
   IF (n = 1) THEN GOTO 99;
   g[1] := r[n-1]/r[n];
   h[1] := r[n+1]/r[n];
   FOR m := 1 to n DO BEGIN
      m1 := m+1;
      sxn := -y[m1];
      sd := -r[n];
      FOR j := 1 to m DO BEGIN
         sxn := sxn+r[n+m1-j]*x[j];
         sd := sd+r[n+m1-j]*g[m-j+1]
      END;
      IF (sd = 0.0) THEN GOTO 98;
      x[m1] := sxn/sd;
      FOR j := 1 to m DO BEGIN
         x[j] := x[j]-x[m1]*g[m-j+1]
      END;
      IF (m1 = n) THEN GOTO 99;
      sgn := -r[n-m1];
      shn := -r[n+m1];
      sgd := -r[n];
      FOR j := 1 to m DO BEGIN
         sgn := sgn+r[n+j-m1]*g[j];
         shn := shn+r[n+m1-j]*h[j];
         sgd := sgd+r[n+j-m1]*h[m-j+1]
      END;
      IF ((sd = 0.0) OR (sgd = 0.0)) THEN GOTO 98;
      g[m1] := sgn/sgd;
      h[m1] := shn/sd;
      k := m;
      m2 := (m+1) DIV 2;
      pp := g[m1];
      qq := h[m1];
      FOR j := 1 to m2 DO BEGIN
         pt1 := g[j];
         pt2 := g[k];
         qt1 := h[j];
         qt2 := h[k];
         g[j] := pt1-pp*qt2;
         g[k] := pt2-pp*qt1;
         h[j] := qt1-qq*pt2;
         h[k] := qt2-qq*pt1;
         k := k-1
      END
   END;
   writeln('pause in TOEPLZ - should not arrive here!'); readln;
   GOTO 99;
98:   writeln('pause in TOEPLZ - Levinson method fails');
   writeln('matrix has a singular principal minor'); readln;
99:   END;
PROCEDURE tptest(data1,data2: glnparray; n: integer; VAR t,prob: real);
(* Programs using routine TPTEST must define type
TYPE
   glnparray := ARRAY [1..np];
in the main routine, with np >= n. *)
VAR
   j: integer;
   var2,var1,sd,df,cov,ave2,ave1: real;
BEGIN
   avevar(data1,n,ave1,var1);
   avevar(data2,n,ave2,var2);
   cov := 0.0;
   FOR j := 1 to n DO BEGIN
      cov := cov+(data1[j]-ave1)*(data2[j]-ave2);
   END;
   df := n-1;
   cov := cov/df;
   sd := sqrt((var1+var2-2.0*cov)/n);
   t := (ave1-ave2)/sd;
   prob := betai(0.5*df,0.5,df/(df+sqr(t)))
END;
PROCEDURE tqli(VAR d,e: glnp; n: integer; VAR z: glnpnp);
(* Programs using routine TQLI must define the types
TYPE
   glnp = ARRAY [1..np] OF real;
   glnpnp = ARRAY [1..np,1..np] OF real;
where np is the physical dimension of the matrix to be analyzed. *)
LABEL 1,2;
VAR
   m,l,iter,i,k: integer;
   s,r,p,g,f,dd,c,b: real;
FUNCTION sign(a,b: real): real;
   BEGIN
      IF (b < 0) THEN sign := -abs(a) ELSE sign := abs(a)
   END;
BEGIN
   IF  (n > 1)  THEN BEGIN
      FOR i := 2 to n DO BEGIN
         e[i-1] := e[i]
      END;
      e[n] := 0.0;
      FOR l := 1 to n DO BEGIN
         iter := 0;
1:         FOR m := l to n-1 DO BEGIN
            dd := abs(d[m])+abs(d[m+1]);
            IF  (abs(e[m])+dd = dd) THEN  GOTO 2
         END;
         m := n;
2:         IF (m <> l) THEN BEGIN
            IF (iter = 30) THEN BEGIN
               writeln('pause in routine TQLI');
               writeln('too many iterations'); readln
            END;
            iter := iter+1;
            g := (d[l+1]-d[l])/(2.0*e[l]);
            r := sqrt(sqr(g)+1.0);
            g := d[m]-d[l]+e[l]/(g+sign(r,g));
            s := 1.0;
            c := 1.0;
            p := 0.0;
            FOR i := m-1 DOWNTO l DO BEGIN
               f := s*e[i];
               b := c*e[i];
               IF (abs(f) >= abs(g)) THEN BEGIN
                  c := g/f;
                  r := sqrt(sqr(c)+1.0);
                  e[i+1] := f*r;
                  s := 1.0/r;
                  c := c*s
               END ELSE BEGIN
                  s := f/g;
                  r := sqrt(sqr(s)+1.0);
                  e[i+1] := g*r;
                  c := 1.0/r;
                  s := s*c
               END;
               g := d[i+1]-p;
               r := (d[i]-g)*s+2.0*c*b;
               p := s*r;
               d[i+1] := g+p;
               g := c*r-b;
            (* Next loop can be omitted if eigenvectors not wanted *)
               FOR k := 1 to n DO BEGIN
                  f := z[k,i+1];
                  z[k,i+1] := s*z[k,i]+c*f;
                  z[k,i] := c*z[k,i]-s*f
               END
            END;
            d[l] := d[l]-p;
            e[l] := g;
            e[m] := 0.0;
            GOTO 1
         END
      END
   END
END;
PROCEDURE trapzd(a,b: real; VAR s: real; n: integer);
(* Programs calling TRAPZD must provide a function
func(x:real):real which is to be integrated. They must
also define the variable
VAR
   glit: integer;
in the main routine. *)
VAR
   j: integer;
   x,tnm,sum,del: real;
BEGIN
   IF (n = 1) THEN BEGIN
      s := 0.5*(b-a)*(func(a)+func(b));
      glit := 1
   END
   ELSE BEGIN
      tnm := glit;
      del := (b-a)/tnm;
      x := a+0.5*del;
      sum := 0.0;
      FOR j := 1 to glit DO BEGIN
         sum := sum+func(x);
         x := x+del
      END;
      s := 0.5*(s+(b-a)*sum/tnm);
      glit := 2*glit
   END
END;
PROCEDURE tred2(VAR a: glnpnp; n: integer; VAR d,e: glnp);
(* Programs using routine TRED2 must define the types
TYPE
   glnp = ARRAY [1..np] OF real;
   glnpnp = ARRAY [1..np,1..np] OF real;
where 'np by np' is the physical dimension of the matrix to be analyzed. *)
VAR
   l,k,j,i: integer;
   scale,hh,h,g,f: real;
FUNCTION sign(a,b: real): real;
   BEGIN
      IF (b < 0) THEN sign := -abs(a) ELSE sign := abs(a)
   END;
BEGIN
   IF (n > 1) THEN BEGIN
      FOR i := n DOWNTO 2 DO BEGIN
         l := i-1;
         h := 0.0;
         scale := 0.0;
         IF (l > 1) THEN BEGIN
            FOR k := 1 to l DO BEGIN
               scale := scale+abs(a[i,k])
            END;
            IF (scale = 0.0) THEN BEGIN
               e[i] := a[i,l]
            END ELSE BEGIN
               FOR k := 1 to l DO BEGIN
                  a[i,k] := a[i,k]/scale;
                  h := h+sqr(a[i,k])
               END;
               f := a[i,l];
               g := -sign(sqrt(h),f);
               e[i] := scale*g;
               h := h-f*g;
               a[i,l] := f-g;
               f := 0.0;
               FOR j := 1 to l DO BEGIN
            (* Next statement can be omitted if eigenvectors not wanted *)
                  a[j,i] := a[i,j]/h;
                  g := 0.0;
                  FOR k := 1 to j DO BEGIN
                     g := g+a[j,k]*a[i,k]
                  END;
                  IF (l > j) THEN FOR k := j+1 to l DO g := g+a[k,j]*a[i,k];
                  e[j] := g/h;
                  f := f+e[j]*a[i,j]
               END;
               hh := f/(h+h);
               FOR j := 1 to l DO BEGIN
                  f := a[i,j];
                  g := e[j]-hh*f;
                  e[j] := g;
                  FOR k := 1 to j DO a[j,k] := a[j,k]-f*e[k]-g*a[i,k]
               END
            END
         END ELSE BEGIN
            e[i] := a[i,l]
         END;
         d[i] := h
      END
   END;
   (* Next statement can be omitted if eigenvectors not wanted *)
   d[1] := 0.0;
   e[1] := 0.0;
   FOR i := 1 to n DO BEGIN
   (* Contents of this loop can be omitted if eigenvectors not wanted,
      except for statement d[i] := a[i,i]; *)
      l := i-1;
      IF (d[i] <> 0.0) THEN BEGIN
         FOR j := 1 to l DO BEGIN
            g := 0.0;
            FOR k := 1 to l DO BEGIN
               g := g+a[i,k]*a[k,j]
            END;
            FOR k := 1 to l DO BEGIN
               a[k,j] := a[k,j]-g*a[k,i]
            END
         END
      END;
      d[i] := a[i,i];
      a[i,i] := 1.0;
      IF (l >= 1) THEN BEGIN
         FOR j := 1 to l DO BEGIN
            a[i,j] := 0.0;
            a[j,i] := 0.0
         END
      END
   END
END;
PROCEDURE tridag(a,b,c,r: glnarray; VAR u: glnarray; n: integer);
(* Programs using routine TRIDAG should define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine. *)
VAR
   j: integer;
   bet: real;
   gam: glnarray;
BEGIN
   IF (b[1] = 0.0) THEN BEGIN writeln('pause 1 in TRIDAG'); readln END;
   bet := b[1];
   u[1] := r[1]/bet;
   FOR j := 2 to n DO BEGIN
      gam[j] := c[j-1]/bet;
      bet := b[j]-a[j]*gam[j];
      IF (bet = 0.0) THEN BEGIN writeln('pause 2 in TRIDAG'); readln END;
      u[j] := (r[j]-a[j]*u[j-1])/bet
   END;
   FOR j := n-1 DOWNTO 1 DO BEGIN
      u[j] := u[j]-gam[j+1]*u[j+1]
   END
END;

PROCEDURE ttest(data1: glnarray; n1: integer;
      data2: glmarray; n2: integer; VAR t,prob: real);
(* Programs using routine TTEST must define the types
TYPE
   glnarray = ARRAY [1..n1] OF real;
   glmarray = ARRAY [1..n2] OF real;
in the main routine, with n2 less than or equal to n1. *)
VAR
   i: integer;
   var2,var1,svar,df,ave2,ave1: real;
BEGIN
   avevar(data1,n1,ave1,var1);
   IF (n2 > n1) THEN BEGIN
      writeln('pause in routine TTEST');
      writeln('The first array should be the larger'); readln
   END;
   FOR i := 1 to n2 DO data1[i] := data2[i];
   avevar(data1,n2,ave2,var2);
   df := n1+n2-2;
   svar := ((n1-1)*var1+(n2-1)*var2)/df;
   t := (ave1-ave2)/sqrt(svar*(1.0/n1+1.0/n2));
   prob := betai(0.5*df,0.5,df/(df+sqr(t)))
END;

PROCEDURE tutest(data1: glnarray; n1: integer; data2: glmarray;
        n2: integer; VAR t,prob: real);
(* Programs using routine TUTEST must define the types
TYPE
   glnarray = ARRAY [1..n1] OF real;
   glmarray = ARRAY [1..n2] OF real;
in the main routine, with n2 less than or equal to n1.   *)
VAR
   var2,var1,df,ave2,ave1: real;
   i: integer;
BEGIN
   IF (n2 > n1) THEN BEGIN
      writeln('pause in routine TUTEST');
      writeln('first array must be the larger'); readln
   END;
   avevar(data1,n1,ave1,var1);
   FOR i := 1 to n2 DO data1[i] := data2[i];
   avevar(data1,n2,ave2,var2);
   t := (ave1-ave2)/sqrt(var1/n1+var2/n2);
   df := sqr(var1/n1+var2/n2)/
      (sqr(var1/n1)/(n1-1)+sqr(var2/n2)/(n2-1));
   prob := betai(0.5*df,0.5,df/(df+sqr(t)))
END;
PROCEDURE twofft(data1,data2: glnarray;
       VAR fft1,fft2: gl2narray; n: integer);
(* Programs using routine TWOFFT must define types
TYPE
   glnarray = ARRAY [1..n] OF real;
   gl2narray = ARRAY [1..2*n] OF real;
where n is the dimension of the real-valued data arrays. *)
VAR
   nn3,nn2,nn,jj,j: integer;
   rep,rem,aip,aim: real;
BEGIN
   nn := n+n;
   nn2 := nn+2;
   nn3 := nn+3;
   FOR j := 1 to n DO BEGIN
      jj := j+j;
      fft1[jj-1] := data1[j];
      fft1[jj] := data2[j]
   END;
   four1(fft1,n,1);
   fft2[1] := fft1[2];
   fft1[2] := 0.0;
   fft2[2] := 0.0;
   FOR jj := 1 to (n DIV 2) DO BEGIN
      j := 2*jj+1;
      rep := 0.5*(fft1[j]+fft1[nn2-j]);
      rem := 0.5*(fft1[j]-fft1[nn2-j]);
      aip := 0.5*(fft1[j+1]+fft1[nn3-j]);
      aim := 0.5*(fft1[j+1]-fft1[nn3-j]);
      fft1[j] := rep;
      fft1[j+1] := aim;
      fft1[nn2-j] := rep;
      fft1[nn3-j] := -aim;
      fft2[j] := aip;
      fft2[j+1] := -rem;
      fft2[nn2-j] := aip;
      fft2[nn3-j] := rem
   END
END;
