PROCEDURE scrsho;
(* Programs using routine SCRSHO must externally define a
function fx(x:real):real which is to be plotted. *)
LABEL 1,99;
CONST
   iscr=60;
   jscr=21;
   blank=' ';
   zero='-';
   yy='l';
   xx='-';
   ff='x';
VAR
   jz,j,i: integer;
   ysml,ybig,x2,x1,x,dyj,dx: real;
   y: ARRAY [1..iscr] OF real;
   scr: ARRAY [1..iscr,1..jscr] OF char;
BEGIN
1:   writeln('Enter x1 x2 (x1=x2 to stop): '); readln(x1,x2);
   IF (x1 = x2) THEN GOTO 99;
   FOR j := 1 to jscr DO BEGIN
      scr[1,j] := yy;
      scr[iscr,j] := yy
   END;
   FOR i := 2 to iscr-1 DO BEGIN
      scr[i,1] := xx;
      scr[i,jscr] := xx;
      FOR j := 2 to jscr-1 DO BEGIN
         scr[i,j] := blank
      END
   END;
   dx := (x2-x1)/(iscr-1);
   x := x1;
   ybig := 0.0;
   ysml := ybig;
   FOR i := 1 to iscr DO BEGIN
      y[i] := fx(x);
      IF (y[i] < ysml) THEN  ysml := y[i];
      IF (y[i] > ybig) THEN  ybig := y[i];
      x := x+dx
   END;
   IF (ybig = ysml) THEN  ybig := ysml+1.0;
   dyj := (jscr-1)/(ybig-ysml);
   jz := 1-trunc(ysml*dyj);
   FOR i := 1 to iscr DO BEGIN
      scr[i,jz] := zero;
      j := 1+trunc((y[i]-ysml)*dyj);
      scr[i,j] := ff
   END;
   write(' ',ybig:10:3,' ');
   FOR i := 1 to iscr DO write(scr[i,jscr]);
   writeln;
   FOR j := jscr-1 DOWNTO 2 DO BEGIN
      write(' ':12);
      FOR i := 1 to iscr DO write(scr[i,j]);
      writeln
   END;
   write(' ',ysml:10:3,' ');
   FOR i := 1 to iscr DO write(scr[i,1]);
   writeln;
   writeln(' ':8,x1:10:3,' ':44,x2:10:3);
   GOTO 1;
99:   END;
PROGRAM sfroid(input,output);
LABEL 99;
CONST
   ne=3;
   m=41;
   nb=1;
   nsi=ne;
   nyj=ne;
   nyk=m;
   nci=ne;
   ncj=3;      (* ncj=ne-nb+1 *)
   nck=42;      (* nck=m+1 *)
   nsj=7;      (* nsj=2*ne+1 *)
TYPE
   glyarray = ARRAY [1..ne,1..m] OF real;
   glcarray = ARRAY [1..nci,1..ncj,1..nck] OF real;
   glsarray = ARRAY [1..nsi,1..nsj] OF real;
   glscalv = ARRAY [1..ne] OF real;
   glindex = ARRAY [1..ne] OF integer;
VAR
   mm,n: integer;
   i,itmax,k: integer;
   h,c2,anorm: real;
   conv,deriv,fac1,fac2: real;
   q1,slowc: real;
   scalv: glscalv;
   indexv: glindex;
   y: glyarray;
   c: glcarray;
   s: glsarray;
   x: ARRAY [1..m] OF real;
(*$I MODFILE.PAS *)
(*$I PLGNDR.PAS *)
(*$I DIFEQ.PAS *)
(*$I RED.PAS *)
(*$I PINVS.PAS *)
(*$I BKSUB.PAS *)
(*$I SOLVDE.PAS *)
BEGIN
   itmax := 100;
   conv := 5.0e-6;
   slowc := 1.0;
   h := 1.0/(m-1);
   writeln('Enter m n');
   readln(mm,n);
   IF (((n+mm) MOD 2) = 1) THEN BEGIN
      indexv[1] := 1;
      indexv[2] := 2;
      indexv[3] := 3
   END ELSE BEGIN
      indexv[1] := 2;
      indexv[2] := 1;
      indexv[3] := 3
   END;
   anorm := 1.0;
   IF (mm <> 0) THEN BEGIN
      q1 := n;
      FOR i := 1 to mm DO BEGIN
         anorm := -0.5*anorm*(n+i)*(q1/i);
         q1 := q1-1.0
      END
   END;
   FOR k := 1 to (m-1) DO BEGIN
      x[k] := (k-1)*h;
      fac1 := 1.0-sqr(x[k]);
      fac2 := exp((-mm/2.0)*ln(fac1));
      y[1,k] := plgndr(n,mm,x[k])*fac2;
      deriv := -((n-mm+1)*plgndr(n+1,mm,x[k])-
         (n+1)*x[k]*plgndr(n,mm,x[k]))/fac1;
      y[2,k] := mm*x[k]*y[1,k]/fac1+deriv*fac2;
      y[3,k] := n*(n+1)-mm*(mm+1)
   END;
   x[m] := 1.0;
   y[1,m] := anorm;
   y[3,m] := n*(n+1)-mm*(mm+1);
   y[2,m] := (y[3,m]-c2)*y[1,m]/(2.0*(mm+1.0));
   scalv[1] := abs(anorm);
   IF (y[2,m] > abs(anorm)) THEN scalv[2] := y[2,m] ELSE scalv[2] := abs(anorm);
   IF (y[3,m] > 1.0) THEN scalv[3] := y[3,m] ELSE scalv[3] := 1.0;
   WHILE true DO BEGIN
      writeln('Enter c**2 or 999 to end.');
      readln(c2);
      IF (c2 = 999) THEN GOTO 99;
      solvde(itmax,conv,slowc,scalv,indexv,ne,nb,m,y,nyj,nyk,
         c,nci,ncj,nck,s,nsi,nsj);
      writeln;
      writeln('m = ',mm:2,'  n = ',n:2,'  c**2 = ',c2:7:3,
         ' lam = ',y[3,1]+mm*(mm+1):10:6);
   END;
99:   END.
PROCEDURE shell(n: integer; VAR arr: glnarray);
(* Programs using routine SHELL must define the type
TYPE
   glnarray = ARRAY [1..np] OF real;
in the main routine, with np >= n.   *)
LABEL 3;
CONST
   aln2i=1.442695022;
   tiny=1.0e-5;
VAR
   nn,m,lognb2,l,k,j,i: integer;
   t: real;
BEGIN
   lognb2 := trunc(ln(n)*aln2i+tiny);
   m := n;
   FOR nn := 1 to lognb2 DO BEGIN
      m := m DIV 2;
      k := n-m;
      FOR j := 1 to k DO BEGIN
         i := j;
3:         l := i+m;
         IF (arr[l] < arr[i]) THEN BEGIN
            t := arr[i];
            arr[i] := arr[l];
            arr[l] := t;
            i := i-m;
            IF (i >= 1) THEN GOTO 3
         END
      END
   END
END;
PROCEDURE shoot(nvar: integer; VAR v: gln2array; delv: gln2array;
      n2: integer; x1,x2,eps,h1,hmin: real; VAR f,dv: glnvar);
(* Programs using routine SHOOT must define the types
TYPE
   gln2array = ARRAY [1..n2] OF real;
   glnvar = ARRAY [1..nvar] OF real;
   gln2byn2 = ARRAY [1..n2,1..n2] OF real;
   glinvar = ARRAY [1..nvar] OF integer;
   glnpbynp = gln2byn2;
in the main routine. *)
VAR
   nok,nbad,iv,i: integer;
   sav,det: real;
   y: glnvar;
   dfdv: gln2byn2;
   indx: glinvar;
BEGIN
   load(x1,v,y);
   odeint(y,nvar,x1,x2,eps,h1,hmin,nok,nbad);
   score(x2,y,f);
   FOR iv := 1 to n2 DO BEGIN
      sav := v[iv];
      v[iv] := v[iv]+delv[iv];
      load(x1,v,y);
      odeint(y,nvar,x1,x2,eps,h1,hmin,nok,nbad);
      score(x2,y,dv);
      FOR i := 1 to n2 DO BEGIN
         dfdv[i,iv] := (dv[i]-f[i])/delv[iv]
      END;
      v[iv] := sav
   END;
   FOR iv := 1 to n2 DO BEGIN
      dv[iv] := -f[iv]
   END;
   ludcmp(dfdv,n2,nvar,indx,det);
   lubksb(dfdv,n2,nvar,indx,dv);
   FOR iv := 1 to n2 DO BEGIN
      v[iv] := v[iv] + dv[iv]
   END
END;
PROCEDURE shootf(nvar: integer; VAR v1: gln2array; VAR v2: gln1array;
       delv1: gln2array; delv2: gln1array; n1,n2: integer;
       x1,x2,xf,eps,h1,hmin: real; VAR f: glnvar;
       VAR dv1: gln2array; VAR dv2: gln1array);
(* Programs using routine SHOOTF must define the types
TYPE
   gln1array = ARRAY [1..n1] OF real;
   gln2array = ARRAY [1..n2] OF real;
   glnvar = ARRAY [1..nvar] OF real;
   glnvarbynvar = ARRAY [1..nvar,1..nvar];
   glindx = ARRAY [1..nvar] OF integer;
   glnpbynp = glnvarbynvar;
in the main routine. *)
VAR
   nok,nbad,j,iv,i: integer;
   sav,det: real;
   y,f1,f2: glnvar;
   dfdv: glnvarbynvar;
   indx: glindx;
BEGIN
   load1(x1,v1,y);
   odeint(y,nvar,x1,xf,eps,h1,hmin,nok,nbad);
   score(xf,y,f1);
   load2(x2,v2,y);
   odeint(y,nvar,x2,xf,eps,h1,hmin,nok,nbad);
   score(xf,y,f2);
   j := 0;
   FOR iv := 1 to n2 DO BEGIN
      j := j+1;
      sav := v1[iv];
      v1[iv] := v1[iv]+delv1[iv];
      load1(x1,v1,y);
      odeint(y,nvar,x1,xf,eps,h1,hmin,nok,nbad);
      score(xf,y,f);
      FOR i := 1 to nvar DO BEGIN
         dfdv[i,j] := (f[i]-f1[i])/delv1[iv]
      END;
      v1[iv] := sav
   END;
   FOR iv := 1 to n1 DO BEGIN
      j := j+1;
      sav := v2[iv];
      v2[iv] := v2[iv]+delv2[iv];
      load2(x2,v2,y);
      odeint(y,nvar,x2,xf,eps,h1,hmin,nok,nbad);
      score(xf,y,f);
      FOR i := 1 to nvar DO BEGIN
         dfdv[i,j] := (f2[i]-f[i])/delv2[iv]
      END;
      v2[iv] := sav
   END;
   FOR i := 1 to nvar DO BEGIN
      f[i] := f1[i]-f2[i];
      f1[i] := -f[i]
   END;
   ludcmp(dfdv,nvar,nvar,indx,det);
   lubksb(dfdv,nvar,nvar,indx,f1);
   j := 0;
   FOR iv := 1 to n2 DO BEGIN
      j := j+1;
      v1[iv] := v1[iv]+f1[j];
      dv1[iv] := f1[j]
   END;
   FOR iv := 1 to n1 DO BEGIN
      j := j+1;
      v2[iv] := v2[iv]+f1[j];
      dv2[iv] := f1[j]
   END
END;
PROCEDURE simp1(a: glmpbynp; mp,np,mm: integer;
      ll: glnparray; nll,iabf: integer;
      VAR kp: integer; VAR bmax: real);
(* Programs using routine SIMP1 must define the types
TYPE
   glmpbynp = ARRAY [1..mp,1..np] OF real;
   glnparray = ARRAY [1..np] OF integer;
in the main routine. *)
LABEL 99;
VAR
   k: integer;
   test: real;
BEGIN
   kp := ll[1];
   bmax := a[mm+1,kp+1];
   IF (nll < 2) THEN GOTO 99;
   FOR k := 2 to nll DO BEGIN
      IF (iabf = 0) THEN BEGIN
         test := a[mm+1,ll[k]+1]-bmax
      END ELSE BEGIN
         test := abs(a[mm+1,ll[k]+1])-abs(bmax)
      END;
      IF (test > 0.0) THEN BEGIN
         bmax := a[mm+1,ll[k]+1];
         kp := ll[k]
      END
   END;
99:   END;
PROCEDURE simp2(a: glmpbynp; m,n,mp,np: integer;
      l2: glmparray; nl2: integer; VAR ip: integer;
      kp: integer; VAR q1: real);
(* Programs using routine SIMP2 must define the types
TYPE
   glmpbynp = ARRAY [1..mp,1..np] OF real;
   glmparray = ARRAY [1..mp] OF integer;
in the main routine. *)
LABEL 2,6,99;
VAR
   k,ii,i: integer;
   qp,q0,q: real;
BEGIN
   ip := 0;
   IF (nl2 < 1) THEN GOTO 99;
   FOR i := 1 to nl2 DO BEGIN
      IF (a[l2[i]+1,kp+1] < 0.0) THEN GOTO 2
   END;
   GOTO 99;
2:   q1 := -a[l2[i]+1,1]/a[l2[i]+1,kp+1];
   ip := l2[i];
   IF ((i+1) > nl2) THEN GOTO 99;
   FOR i := i+1 to nl2 DO BEGIN
      ii := l2[i];
      IF (a[ii+1,kp+1] < 0.0) THEN BEGIN
         q := -a[ii+1,1]/a[ii+1,kp+1];
         IF (q < q1) THEN BEGIN
            ip := ii;
            q1 := q
         END ELSE IF (q = q1) THEN BEGIN
            FOR k := 1 to n DO BEGIN
               qp := -a[ip+1,k+1]/a[ip+1,kp+1];
               q0 := -a[ii+1,k+1]/a[ii+1,kp+1];
               IF (q0 <> qp) THEN GOTO 6
            END;
6:            IF (q0 < qp) THEN ip := ii
         END
      END
   END;
99:   END;
PROCEDURE simp3(VAR a: glmpbynp; mp,np,i1,k1,ip,kp: integer);
(* Programs using routine SIMP3 must define the type
TYPE
   glmpbynp = ARRAY [1..mp,1..np] OF real;
in the main routine. *)
VAR
   kk,ii: integer;
   piv: real;
BEGIN
   piv := 1.0/a[ip+1,kp+1];
   IF (i1 >= 0) THEN BEGIN
      FOR ii := 1 to (i1+1) DO BEGIN
         IF ((ii-1) <> ip) THEN BEGIN
            a[ii,kp+1] := a[ii,kp+1]*piv;
            FOR kk := 1 to k1+1 DO BEGIN
               IF ((kk-1) <> kp) THEN BEGIN
                  a[ii,kk] := a[ii,kk]
                   -a[ip+1,kk]*a[ii,kp+1]
               END
            END
         END
      END
   END;
   FOR kk := 1 to k1+1 DO BEGIN
      IF ((kk-1) <> kp) THEN a[ip+1,kk] := -a[ip+1,kk]*piv
   END;
   a[ip+1,kp+1] := piv
END;
PROCEDURE simplx(VAR a: glmpbynp; m,n,mp,np,m1,m2,m3: integer;
       VAR icase: integer; VAR izrov: glnarray;
       VAR iposv: glmarray);
(* Programs using routine SIMPLX must define the types
TYPE
   glmpbynp = ARRAY [1..mp,1..np] OF real;
   glnarray = ARRAY [1..n] OF integer;
   glmarray = ARRAY [1..m] OF integer;
   glmparray = ARRAY [1..mp] OF integer;
   glnparray = ARRAY [1..np] OF integer;
in the main routine. *)
LABEL 1,2,10,20,30,99;
CONST eps=1.0e-6;
VAR
   nl2,nl1,m12,kp,kh,k,is,ir,ip,i: integer;
   q1,bmax: real;
   l1: glnparray;
   l2,l3: glmparray;
BEGIN
   IF (m <> (m1+m2+m3)) THEN BEGIN
      writeln('pause in routine SIMPLX');
      writeln('bad input constraint counts'); readln
   END;
   nl1 := n;
   FOR k := 1 to n DO BEGIN
      l1[k] := k;
      izrov[k] := k
   END;
   nl2 := m;
   FOR i := 1 to m DO BEGIN
      IF (a[i+1,1] < 0.0) THEN BEGIN
         writeln('pause in routine SIMPLX');
         writeln('bad input tableau'); readln
      END;
      l2[i] := i;
      iposv[i] := n+i
   END;
   FOR i := 1 to m2 DO BEGIN
      l3[i] := 1
   END;
   ir := 0;
   IF ((m2+m3) = 0) THEN GOTO 30;
   ir := 1;
   FOR k := 1 to n+1 DO BEGIN
      q1 := 0.0;
      FOR i := m1+1 to m DO BEGIN
         q1 := q1+a[i+1,k]
      END;
      a[m+2,k] := -q1
   END;
10:   simp1(a,mp,np,m+1,l1,nl1,0,kp,bmax);
   IF ((bmax <= eps) AND (a[m+2,1] < -eps)) THEN BEGIN
      icase := -1; GOTO 99 END
   ELSE IF ((bmax <= eps) AND (a[m+2,1] <= eps)) THEN BEGIN
      m12 := m1+m2+1;
      IF (m12 <= m) THEN BEGIN
         FOR ip := m12 to m DO BEGIN
            IF (iposv[ip] = (ip+n)) THEN BEGIN
               simp1(a,mp,np,ip,l1,nl1,1,kp,bmax);
               IF (bmax > 0.0) THEN GOTO 1
            END
         END
      END;
      ir := 0;
      m12 := m12-1;
      IF ((m1+1) > m12) THEN GOTO 30;
      FOR i := m1+1 to m12 DO BEGIN
         IF (l3[i-m1] = 1) THEN BEGIN
            FOR k := 1 to n+1 DO BEGIN
               a[i+1,k] := -a[i+1,k]
            END
         END
      END;
      GOTO 30
   END;
   simp2(a,m,n,mp,np,l2,nl2,ip,kp,q1);
   IF (ip = 0) THEN BEGIN
      icase := -1; GOTO 99
   END;
1:   simp3(a,mp,np,m+1,n,ip,kp);
   IF (iposv[ip] >= (n+m1+m2+1)) THEN BEGIN
      FOR k := 1 to nl1 DO BEGIN
         IF (l1[k] = kp) THEN GOTO 2
      END;
2:      nl1 := nl1-1;
      FOR is := k to nl1 DO BEGIN
         l1[is] := l1[is+1]
      END
   END ELSE BEGIN
      IF (iposv[ip] < (n+m1+1)) THEN GOTO 20;
      kh := iposv[ip]-m1-n;
      IF (l3[kh] = 0) THEN GOTO 20;
      l3[kh] := 0
   END;
   a[m+2,kp+1] := a[m+2,kp+1]+1.0;
   FOR i := 1 to m+2 DO BEGIN
      a[i,kp+1] := -a[i,kp+1]
   END;
20:   is := izrov[kp];
   izrov[kp] := iposv[ip];
   iposv[ip] := is;
   IF (ir <> 0) THEN GOTO 10;
30:   simp1(a,mp,np,0,l1,nl1,0,kp,bmax);
   IF (bmax <= 0.0) THEN BEGIN
      icase := 0; GOTO 99
   END;
   simp2(a,m,n,mp,np,l2,nl2,ip,kp,q1);
   IF (ip = 0) THEN BEGIN
      icase := 1; GOTO 99
   END;
   simp3(a,mp,np,m,n,ip,kp);
   GOTO 20;
99:   END;
PROCEDURE sinft(VAR y: glyarray; n: integer);
(* Programs using routine SINFT must define the type
TYPE
   glyarray = ARRAY [1..n] OF real;
where n is the dimension of the input data. *)
VAR
   jj,j,m,n2: integer;
   sum,y1,y2: real;
   theta,wi,wr,wpi,wpr,wtemp: double;
BEGIN
   theta := 3.14159265358979/n;
   wr := 1.0;
   wi := 0.0;
   wpr := -2.0*sqr(sin(0.5*theta));
   wpi := sin(theta);
   y[1] := 0.0;
   m := n DIV 2;
   n2 := n+2;
   FOR j := 2 to (m+1) DO BEGIN
      wtemp := wr;
      wr := wr*wpr-wi*wpi+wr;
      wi := wi*wpr+wtemp*wpi+wi;
      y1 := sngl(wi)*(y[j]+y[n2-j]);
      y2 := 0.5*(y[j]-y[n2-j]);
      y[j] := y1+y2;
      y[n2-j] := y1-y2
   END;
   realft(y,m,+1);
   sum := 0.0;
   y[1] := 0.5*y[1];
   y[2] := 0.0;
   FOR jj := 0 to (m-1) DO BEGIN
      j := 2*jj+1;
      sum := sum+y[j];
      y[j] := y[j+1];
      y[j+1] := sum
   END
END;
PROCEDURE smooft(VAR y: glyarray; n: integer; pts: real);
(* Programs using routine SMOOFT must define the types
TYPE
   glyarray = ARRAY [1..n] OF real;
in the main routine, *)
VAR
   nmin,m,mo2,k,j: integer;
   yn,y1,rn1,fac,cnst: real;
BEGIN
   m := 2;
   nmin := n+round(2.0*pts);
   WHILE (m < nmin) DO m := 2*m;
   cnst := sqr(pts/m);
   y1 := y[1];
   yn := y[n];
   rn1 := 1.0/(n-1.0);
   FOR j := 1 to n DO BEGIN
      y[j] := y[j]-rn1*(y1*(n-j)+yn*(j-1))
   END;
   IF (n+1 <= m) THEN BEGIN
      FOR j := n+1 to m DO BEGIN
         y[j] := 0.0
      END
   END;
   mo2 := m DIV 2;
   realft(y,mo2,1);
   y[1] := y[1]/mo2;
   fac := 1.0;
   FOR j := 1 to (mo2-1) DO BEGIN
      k := 2*j+1;
      IF (fac <> 0.0) THEN BEGIN
         fac := (1.0-cnst*j*j)/mo2;
         IF (fac < 0.0) THEN fac := 0.0;
         y[k] := fac*y[k];
         y[k+1] := fac*y[k+1]
      END ELSE BEGIN
         y[k] := 0.0;
         y[k+1] := 0.0
      END
   END;
   fac := (1.0-0.25*pts*pts)/mo2;
   IF (fac < 0.0) THEN fac := 0.0;
   y[2] := fac*y[2];
   realft(y,mo2,-1);
   FOR j := 1 to n DO BEGIN
      y[j] := rn1*(y1*(n-j)+yn*(j-1))+y[j]
   END
END;

PROCEDURE sncndn(uu,emmc: real; VAR sn,cn,dn: real);
LABEL 1;
CONST
   ca=0.0003;
VAR
   a,b,c,d,emc,u: real;
   i,ii,l: integer;
   bo: boolean;
   em,en: ARRAY [1..13] OF real;
FUNCTION cosh(u: real): real;
   BEGIN cosh := 0.5*(exp(u)+exp(-u)) END;
FUNCTION tanh(u: real): real;
   VAR
      u2,epu,emu: real;
   BEGIN
      epu := exp(u);
      emu := 1.0/epu;
      IF (abs(u)<0.3) THEN BEGIN
         u2 := u*u;
         tanh := 2*u*(1+u2/6*(1+u2/20*(1+u2/42*(1+u2/72))))/(epu+emu) END
      ELSE BEGIN tanh := (epu-emu)/(epu+emu) END
   END;
BEGIN
   emc := emmc;
   u := uu;
   IF (emc <> 0.0) THEN BEGIN
      bo := (emc < 0.0);
      IF (bo) THEN BEGIN
         d := 1.0-emc;
         emc := -emc/d;
         d := sqrt(d);
         u := d*u
      END;
      a := 1.0;
      dn := 1.0;
      FOR i := 1 to 13 DO BEGIN
         l := i;
         em[i] := a;
         emc := sqrt(emc);
         en[i] := emc;
         c := 0.5*(a+emc);
         IF (abs(a-emc) <= ca*a) THEN GOTO 1;
         emc := a*emc;
         a := c
      END;
1:      u := c*u;
      sn := sin(u);
      cn := cos(u);
      IF (sn <> 0.0) THEN BEGIN
         a := cn/sn;
         c := a*c;
         FOR ii := l DOWNTO 1 DO BEGIN
            b := em[ii];
            a := c*a;
            c := dn*c;
            dn := (en[ii]+a)/(b+a);
            a := c/b
         END;
         a := 1.0/sqrt(sqr(c)+1.0);
         IF (sn < 0.0) THEN sn := -a
         ELSE sn := a;
         cn := c*sn
      END;
      IF (bo) THEN BEGIN
         a := dn;
         dn := cn;
         cn := a;
         sn := sn/d
      END;
   END ELSE BEGIN
      cn := 1.0/cosh(u);
      dn := cn;
      sn := tanh(u)
   END
END;
PROCEDURE solvde(itmax: integer; conv,slowc: real; scalv: glscalv;
       indexv: glindex; ne,nb,m: integer; VAR y: glyarray;
       nyj,nyk: integer; VAR c: glcarray; nci,ncj,nck: integer;
       VAR s: glsarray; nsi,nsj: integer);
(* Programs using routine SOLVDE must define the types
TYPE
   glindex = ARRAY [1..nyj] OF integer;
   glscalv = ARRAY [1..nyj] OF real;
   glyarray = ARRAY [1..nyj,1..nyk] OF real;
   glcarray = ARRAY [1..nci,1..ncj,1..nck] OF real;
   glsarray = ARRAY [1..nsi,1..nsj] OF real;
in the main routine. *)
LABEL 99;
CONST
   nmax=10;
VAR
   err,errj,fac,vmax,vz: real;
   ic1,ic2,ic3,ic4,it: integer;
   j,j1,j2,j3,j4,j5,j6,j7,j8,j9: integer;
   jc1,jcf,jv,k,k1,k2,km,kp,nvars: integer;
   ermax: ARRAY [1..nmax] OF real;
   kmax: ARRAY [1..nmax] OF integer;
BEGIN
   k1 := 1;
   k2 := m;
   nvars := ne*m;
   j1 := 1;
   j2 := nb;
   j3 := nb+1;
   j4 := ne;
   j5 := j4+j1;
   j6 := j4+j2;
   j7 := j4+j3;
   j8 := j4+j4;
   j9 := j8+j1;
   ic1 := 1;
   ic2 := ne-nb;
   ic3 := ic2+1;
   ic4 := ne;
   jc1 := 1;
   jcf := ic3;
   FOR it := 1 to itmax DO BEGIN
      k := k1;
      difeq(k,k1,k2,j9,ic3,ic4,indexv,ne,s,nsi,nsj,y,nyj,nyk);
      pinvs(ic3,ic4,j5,j9,jc1,k1,c,nci,ncj,nck,s,nsi,nsj);
      FOR k := k1+1 to k2 DO BEGIN
         kp := k-1;
         difeq(k,k1,k2,j9,ic1,ic4,indexv,ne,s,nsi,nsj,y,nyj,nyk);
         red(ic1,ic4,j1,j2,j3,j4,j9,ic3,jc1,jcf,kp,c,nci,ncj,nck,s,nsi,nsj);
         pinvs(ic1,ic4,j3,j9,jc1,k,c,nci,ncj,nck,s,nsi,nsj)
      END;
      k := k2+1;
      difeq(k,k1,k2,j9,ic1,ic2,indexv,ne,
         s,nsi,nsj,y,nyj,nyk);
      red(ic1,ic2,j5,j6,j7,j8,j9,ic3,jc1,jcf,k2,
         c,nci,ncj,nck,s,nsi,nsj);
      pinvs(ic1,ic2,j7,j9,jcf,k2+1,
         c,nci,ncj,nck,s,nsi,nsj);
      bksub(ne,nb,jcf,k1,k2,c,nci,ncj,nck);
      err := 0.0;
      FOR j := 1 to ne DO BEGIN
         jv := indexv[j];
         ermax[j] := 0.0;
         errj := 0.0;
         kmax[j] := 0;
         vmax := 0.0;
         FOR k := k1 to k2 DO BEGIN
            vz := abs(c[j,1,k]);
            IF (vz > vmax)  THEN BEGIN
                vmax := vz;
                km := k
            END;
            errj := errj+vz
         END;
         err := err+errj/scalv[jv];
         ermax[j] := c[j,1,km]/scalv[jv];
         kmax[j] := km
      END;
      err := err/nvars;
      fac := 1.0;
      IF (err > slowc) THEN fac := slowc/err;
      FOR jv := 1 to ne DO BEGIN
         j := indexv[jv];
         FOR k := k1 to k2 DO BEGIN
            y[j,k] := y[j,k]-fac*c[jv,1,k]
         END
      END;
      writeln;
      writeln('Iter.':8,'Error':9,'FAC':9);
      writeln(it:6,err:12:6,fac:11:6);
      writeln('Var.':8,'Kmax':8,'Max. Error':14);
      FOR j := 1 to ne DO writeln(indexv[j]:6,
         kmax[j]:9,ermax[j]:14:6);
      IF (err < conv) THEN GOTO 99
   END;
   writeln('pause in routine SOLVDE');
   writeln('too many iterations'); readln;
99:   END;
PROCEDURE sor(a,b,c,d,e,f: gljmax; VAR u: gljmax;
         jmax: integer; rjac: double);
(* Programs using routine SOR must define the type
TYPE
   gljmax = ARRAY [1..jmax,1..jmax] OF double;
in the main routine. *)
LABEL 99;
CONST
   maxits=1000;
   eps=1.0e-5;
   zero=0.0;
   half=0.5;
   qtr=0.25;
   one=1.0;
VAR
   n,l,j: integer;
   resid,omega,anormf,anorm: double;
BEGIN
   anormf := zero;
   FOR j := 2 to jmax-1 DO BEGIN
      FOR l := 2 to jmax-1 DO BEGIN
         anormf := anormf+abs(f[j,l])
      END
   END;
   omega := one;
   FOR n := 1 to maxits DO BEGIN
      anorm := zero;
      FOR j := 2 to (jmax-1) DO BEGIN
         FOR l := 2 to (jmax-1) DO BEGIN
            IF (((j+l) MOD 2) = (n MOD 2)) THEN BEGIN
               resid := a[j,l]*u[j+1,l]+b[j,l]*u[j-1,l]
                  +c[j,l]*u[j,l+1]+d[j,l]*u[j,l-1]
                  +e[j,l]*u[j,l]-f[j,l];
               anorm := anorm+abs(resid);
               u[j,l] := u[j,l]-omega*resid/e[j,l]
            END
         END
      END;
      IF (n = 1) THEN BEGIN
         omega := one/(one-half*sqr(rjac))
      END ELSE BEGIN
         omega := one/(one-qtr*sqr(rjac)*omega)
      END;
      IF ((n > 1) AND (anorm < (eps*anormf))) THEN GOTO 99
   END;
   writeln('pause in routine SOR');
   writeln('too many iterations'); readln;
99:   END;
PROCEDURE sort(n: integer; VAR ra: glsarray);
(* Program using routine SORT must define the type
TYPE
   glsarray = ARRAY [1..np] OF real;
in the main routine, with np >= n.   *)
LABEL 99;
VAR
   l,j,ir,i: integer;
   rra: real;
BEGIN
   l := (n DIV 2)+1;
   ir := n;
   WHILE true DO BEGIN
      IF (l > 1) THEN BEGIN
         l := l-1;
         rra := ra[l]
      END ELSE BEGIN
         rra := ra[ir];
         ra[ir] := ra[1];
         ir := ir-1;
         IF (ir = 1) THEN BEGIN
            ra[1] := rra;
            GOTO 99
         END
      END;
      i := l;
      j := l+l;
      WHILE (j <= ir) DO BEGIN
         IF (j < ir) THEN
            IF (ra[j] < ra[j+1]) THEN j := j+1;
         IF (rra < ra[j]) THEN BEGIN
            ra[i] := ra[j];
            i := j;
            j := j+j
         END ELSE
            j := ir+1
      END;
      ra[i] := rra
   END;
99:   END;
PROCEDURE sort2(n: integer;VAR ra,rb: glsarray);
(* Programs using routine SORT2 must define type
TYPE
   glsarray = ARRAY [1..np] OF real;
in the main routine, with np >= n.   *)
LABEL 99;
VAR
   l,j,ir,i: integer;
   rrb,rra: real;
BEGIN
   l := (n DIV 2)+1;
   ir := n;
   WHILE true DO BEGIN
      IF (l > 1) THEN BEGIN
         l := l-1;
         rra := ra[l];
         rrb := rb[l]
      END ELSE BEGIN
         rra := ra[ir];
         rrb := rb[ir];
         ra[ir] := ra[1];
         rb[ir] := rb[1];
         ir := ir-1;
         IF (ir = 1) THEN BEGIN
            ra[1] := rra;
            rb[1] := rrb;
            GOTO 99
         END
      END;
      i := l;
      j := l+l;
      WHILE (j <= ir) DO BEGIN
         IF (j < ir) THEN
            IF (ra[j] < ra[j+1]) THEN j := j+1;
         IF (rra < ra[j]) THEN BEGIN
            ra[i] := ra[j];
            rb[i] := rb[j];
            i := j;
            j := j+j
         END ELSE j := ir+1
      END;
      ra[i] := rra;
      rb[i] := rrb
   END;
99:   END;
PROCEDURE sort3(n: integer; VAR ra,rb,rc,wksp: glsarray;
      VAR iwksp: gliarray);
(* Programs using routine SORT3 must define the types
TYPE
   glsarray = ARRAY [1..np] OF real;
   gliarray = ARRAY [1..np] OF integer;
in the main routine, with np >= n.   *)
VAR
   j: integer;
BEGIN
   indexx(n,ra,iwksp);
   FOR j := 1 to n DO BEGIN
      wksp[j] := ra[j]
   END;
   FOR j := 1 to n DO BEGIN
      ra[j] := wksp[iwksp[j]]
   END;
   FOR j := 1 to n DO BEGIN
      wksp[j] := rb[j]
   END;
   FOR j := 1 to n DO BEGIN
      rb[j] := wksp[iwksp[j]]
   END;
   FOR j := 1 to n DO BEGIN
      wksp[j] := rc[j]
   END;
   FOR j := 1 to n DO BEGIN
      rc[j] := wksp[iwksp[j]]
   END
END;
PROCEDURE sparse(b: glnarray; n: integer; VAR x: glnarray; VAR rsq: real);
(* Programs using routine SPARSE must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine. They must also provide two routines,
PROCEDURE asub(x: glnarray; VAR y: glnarray; n: integer);
and
PROCEDURE atsub(x: glnarray; VAR z: glnarray; n: integer);
which calculate A*x and (A transpose)*x *)
LABEL 1,99;
CONST
   eps=1.0e-6;
VAR
   j,iter,irst: integer;
   rp,gg,gam,eps2,dgg,bsq,anum,aden: real;
   g,h,xi,xj: glnarray;
BEGIN
   eps2 := n*sqr(eps);
   irst := 0;
1:   irst := irst+1;
   asub(x,xi,n);
   rp := 0.0;
   bsq := 0.0;
   FOR j := 1 to n DO BEGIN
      bsq := bsq+sqr(b[j]);
      xi[j] := xi[j]-b[j];
      rp := rp+sqr(xi[j])
   END;
   atsub(xi,g,n);
   FOR j := 1 to n DO BEGIN
      g[j] := -g[j];
      h[j] := g[j]
   END;
   FOR iter := 1 to 10*n DO BEGIN
      asub(h,xi,n);
      anum := 0.0;
      aden := 0.0;
      FOR j := 1 to n DO BEGIN
         anum := anum+g[j]*h[j];
         aden := aden+sqr(xi[j])
      END;
      IF (aden = 0.0) THEN BEGIN
         writeln('pause in routine SPARSE');
         writeln('very singular matrix'); readln
      END;
      anum := anum/aden;
      FOR j := 1 to n DO BEGIN
         xi[j] := x[j];
         x[j] := x[j]+anum*h[j]
      END;
      asub(x,xj,n);
      rsq := 0.0;
      FOR j := 1 to n DO BEGIN
         xj[j] := xj[j]-b[j];
         rsq := rsq+sqr(xj[j])
      END;
      IF ((rsq = rp) OR (rsq <= bsq*eps2)) THEN GOTO 99;
      IF (rsq > rp) THEN BEGIN
         FOR j := 1 to n DO BEGIN
            x[j] := xi[j]
         END;
         IF (irst >= 3) THEN GOTO 99;
         GOTO 1
      END;
      rp := rsq;
      atsub(xj,xi,n);
      gg := 0.0;
      dgg := 0.0;
      FOR j := 1 to n DO BEGIN
         gg := gg+sqr(g[j]);
         dgg := dgg+(xi[j]+g[j])*xi[j]
      END;
      IF (gg = 0.0) THEN GOTO 99;
      gam := dgg/gg;
      FOR j := 1 to n DO BEGIN
         g[j] := -xi[j];
         h[j] := g[j]+gam*h[j]
      END
   END;
   writeln('pause in routine SPARSE');
   writeln('too many iterations'); readln;
99:   END;
PROCEDURE spctrm(VAR p: glmarray; m,k: integer; ovrlap: boolean;
       VAR w1: gl4marray; VAR w2: glmarray);
(* Programs using routine SPCTRM must include 'dfile' as
a main routine parameter, and declare
VAR
   dfile: text;
They must open and close this file appropriately.
Also they must define data types
TYPE
   glmarray = ARRAY [1..m] OF real;
   gl4marray = ARRAY [1..4*m] OF real;
and must declare two workspace matrices w1,w2 with types as shown
in the arguments above. *)
VAR
   mm,m44,m43,m4,kk,joffn,joff,j2,j,jj: integer;
   w,sumw,facp,facm,den: real;
FUNCTION window(j: integer; facm,facp: real): real;
   BEGIN
      window := (1.0-abs(((j-1)-facm)*facp))      (* Parzen *)
(*      window := 1.0            *)   (* Square *)
(*      window := (1.0-sqr(((j-1)-facm)*facp))    *)   (* Welch  *)
   END;
BEGIN
   mm := m+m;
   m4 := mm+mm;
   m44 := m4+4;
   m43 := m4+3;
   den := 0.0;
   facm := m-0.5;
   facp := 1.0/(m+0.5);
   sumw := 0.0;
   FOR j := 1 to mm DO sumw := sumw+sqr(window(j,facm,facp));
   FOR j := 1 to m DO p[j] := 0.0;
   IF (ovrlap) THEN BEGIN
      FOR j := 1 to m DO read(dfile,w2[j])
   END;
   FOR kk := 1 to k DO BEGIN
      FOR joff := -1 to 0 DO BEGIN
         IF (ovrlap) THEN BEGIN
            FOR j := 1 to m DO w1[joff+j+j] := w2[j];
            FOR j := 1 to m DO read(dfile,w2[j]);
            joffn := joff+mm;
            FOR j := 1 to m DO w1[joffn+j+j] := w2[j] END
         ELSE BEGIN
            FOR jj := 0 to ((m4-joff-2) DIV 2) DO BEGIN
               j := joff+2+2*jj;
               read(dfile,w1[j])
            END
         END
      END;
      FOR j := 1 to mm DO BEGIN
         j2 := j+j;
         w := window(j,facm,facp);
         w1[j2] := w1[j2]*w;
         w1[j2-1] := w1[j2-1]*w
      END;
      four1(w1,mm,1);
      p[1] := p[1]+sqr(w1[1])+sqr(w1[2]);
      FOR j := 2 to m DO BEGIN
         j2 := j+j;
         p[j] := p[j]+sqr(w1[j2])+sqr(w1[j2-1])
            +sqr(w1[m44-j2])+sqr(w1[m43-j2])
      END;
      den := den+sumw
   END;
   den := m4*den;
   FOR j := 1 to m DO p[j] := p[j]/den
END;
PROCEDURE spear(data1,data2: narray; n: integer;
      VAR wksp1,wksp2: narray;
      VAR d,zd,probd,rs,probrs: real);
(* Programs using routine SPEAR must define types
TYPE
   narray = ARRAY [1..n] OF real;
   glsarray = narray;
in the calling routine *)
VAR
   j: integer;
   vard,t,sg,sf,fac,en3n,en,df,aved: real;
BEGIN
   FOR j := 1 to n DO BEGIN
      wksp1[j] := data1[j];
      wksp2[j] := data2[j]
   END;
   sort2(n,wksp1,wksp2);
   crank(n,wksp1,sf);
   sort2(n,wksp2,wksp1);
   crank(n,wksp2,sg);
   d := 0.0;
   FOR j := 1 to n DO d := d+sqr(wksp1[j]-wksp2[j]);
   en := n;
   en3n := en*en*en-en;
   aved := en3n/6.0-(sf+sg)/12.0;
   fac := (1.0-sf/en3n)*(1.0-sg/en3n);
   vard := ((en-1.0)*sqr(en)*sqr(en+1.0)/36.0)*fac;
   zd := (d-aved)/sqrt(vard);
   probd := erfcc(abs(zd)/1.4142136);
   rs := (1.0-(6.0/en3n)*(d+0.5*(sf+sg)))/fac;
   t := rs*sqrt((en-2.0)/((1.0+rs)*(1.0-rs)));
   df := en-2.0;
   probrs := betai(0.5*df,0.5,df/(df+sqr(t)))
END;
PROCEDURE splie2(x1a,x2a: glnarray; ya: glmbyn;
       m,n: integer; VAR y2a: glmbyn);
(* Programs using routine SPLIE2 must define the types
TYPE
   glnarray = ARRAY [1..nn] OF real;
   glmbyn = ARRAY [1..m,1..n] OF real;
in the main routine. The dimension nn of glnarray must be set to the larger
of n and m. *)
VAR
   k,j: integer;
   ytmp,y2tmp: glnarray;
BEGIN
   FOR j := 1 to m DO BEGIN
      FOR k := 1 to n DO BEGIN
         ytmp[k] := ya[j,k]
      END;
      spline(x2a,ytmp,n,1.0e30,1.0e30,y2tmp);
      FOR k := 1 to n DO BEGIN
         y2a[j,k] := y2tmp[k]
      END
   END
END;
PROCEDURE splin2(x1a,x2a: glnarray; ya,y2a: glmbyn;
       m,n: integer; x1,x2: real; VAR y: real);
(* Programs using routine SPLIN2 must define the types
TYPE
   glnarray = ARRAY [1..nn] OF real;
   glmbyn = ARRAY [1..m,1..n] OF real;
in the main routine. The dimension nn of glnarray must be set
to the larger of n and m *)
VAR
   k,j: integer;
   ytmp,y2tmp,yytmp: glnarray;
BEGIN
   FOR j := 1 to m DO BEGIN
      FOR k := 1 to n DO BEGIN
         ytmp[k] := ya[j,k];
         y2tmp[k] := y2a[j,k]
      END;
      splint(x2a,ytmp,y2tmp,n,x2,yytmp[j])
   END;
   spline(x1a,yytmp,m,1.0e30,1.0e30,y2tmp);
   splint(x1a,yytmp,y2tmp,m,x1,y)
END;
PROCEDURE spline(x,y: glnarray; n: integer; yp1,ypn: real;
       VAR y2: glnarray);
(* Programs using routine SPLINE must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine. *)
VAR
   i,k: integer;
   p,qn,sig,un: real;
   u: glnarray;
BEGIN
   IF (yp1 > 0.99e30) THEN BEGIN
      y2[1] := 0.0;
      u[1] := 0.0
   END ELSE BEGIN
      y2[1] := -0.5;
      u[1] := (3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1)
   END;
   FOR i := 2 to n-1 DO BEGIN
      sig := (x[i]-x[i-1])/(x[i+1]-x[i-1]);
      p := sig*y2[i-1]+2.0;
      y2[i] := (sig-1.0)/p;
      u[i] := (y[i+1]-y[i])/(x[i+1]-x[i])
         -(y[i]-y[i-1])/(x[i]-x[i-1]);
      u[i] := (6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p
   END;
   IF (ypn > 0.99e30) THEN BEGIN
      qn := 0.0;
      un := 0.0
   END ELSE BEGIN
      qn := 0.5;
      un := (3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]))
   END;
   y2[n] := (un-qn*u[n-1])/(qn*y2[n-1]+1.0);
   FOR k := n-1 DOWNTO 1 DO BEGIN
      y2[k] := y2[k]*y2[k+1]+u[k]
   END
END;
PROCEDURE splint(xa,ya,y2a: glnarray; n: integer;
       x: real; VAR y: real);
(* Programs using routine SPLINT must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine. *)
VAR
   klo,khi,k: integer;
   h,b,a: real;
BEGIN
   klo := 1;
   khi := n;
   WHILE (khi-klo > 1) DO BEGIN
      k := (khi+klo) DIV 2;
      IF (xa[k] > x) THEN khi := k ELSE klo := k
   END;
   h := xa[khi]-xa[klo];
   IF (h = 0.0) THEN BEGIN
      writeln ('pause in routine SPLINT');
      writeln (' ... bad XA input'); readln END;
   a := (xa[khi]-x)/h;
   b := (x-xa[klo])/h;
   y := a*ya[klo]+b*ya[khi]+
      ((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0
END;
PROCEDURE svbksb(u: glmpbynp; w: glnparray; v: glnpbynp;
       m,n,mp,np: integer; b: glmparray; VAR x: glnparray);
(* Programs using SVBKSB must define the types
TYPE
   glnparray = ARRAY [1..np] OF real;
   glmparray = ARRAY [1..mp] OF real;
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glmpbynp = ARRAY [1..mp,1..np] OF real;
in the main routine. *)
VAR
   jj,j,i: integer;
   s: real;
   tmp: glnparray;
BEGIN
   FOR j := 1 to n DO BEGIN
      s := 0.0;
      IF (w[j] <> 0.0) THEN BEGIN
         FOR i := 1 to m DO BEGIN
            s := s+u[i,j]*b[i]
         END;
         s := s/w[j]
      END;
      tmp[j] := s
   END;
   FOR j := 1 to n DO BEGIN
      s := 0.0;
      FOR jj := 1 to n DO BEGIN
         s := s+v[j,jj]*tmp[jj];
      END;
      x[j] := s
   END
END;
PROCEDURE svdcmp(VAR a: glmpbynp; m,n,mp,np: integer;
       VAR w: glnparray; VAR v: glnpbynp);
(* Programs using routine SVDCMP must define the types
TYPE
   glnparray = ARRAY [1..np] OF real;
   glmpbynp = ARRAY [1..mp,1..np] OF real;
   glnpbynp = ARRAY [1..np,1..np] OF real;
in the main routine. *)
LABEL 1,2,3;
CONST
   nmax=100;
VAR
   nm,l,k,j,its,i: integer;
   z,y,x,scale,s,h,g,f,c,anorm: real;
   rv1: ARRAY [1..nmax] OF real;
FUNCTION sign(a,b: real): real;
   BEGIN
      IF (b >= 0.0) THEN sign := abs(a) ELSE sign := -abs(a)
   END;
FUNCTION max(a,b: real): real;
   BEGIN
      IF (a > b) THEN max := a ELSE max := b
   END;
BEGIN
   g := 0.0;
   scale := 0.0;
   anorm := 0.0;
   FOR i := 1 to n DO BEGIN
      l := i+1;
      rv1[i] := scale*g;
      g := 0.0;
      s := 0.0;
      scale := 0.0;
      IF (i <= m) THEN BEGIN
         FOR k := i to m DO BEGIN
            scale := scale+abs(a[k,i])
         END;
         IF (scale <> 0.0) THEN BEGIN
            FOR k := i to m DO BEGIN
               a[k,i] := a[k,i]/scale;
               s := s+a[k,i]*a[k,i]
            END;
            f := a[i,i];
            g := -sign(sqrt(s),f);
            h := f*g-s;
            a[i,i] := f-g;
            IF (i <> n) THEN BEGIN
               FOR j := l to n DO BEGIN
                  s := 0.0;
                  FOR k := i to m DO BEGIN
                     s := s+a[k,i]*a[k,j]
                  END;
                  f := s/h;
                  FOR k := i to m DO BEGIN
                     a[k,j] := a[k,j]+
                        f*a[k,i]
                  END
               END
            END;
            FOR k := i to m DO BEGIN
               a[k,i] := scale*a[k,i]
            END
         END
      END;
      w[i] := scale*g;
      g := 0.0;
      s := 0.0;
      scale := 0.0;
      IF ((i <= m) AND (i <> n)) THEN BEGIN
         FOR k := l to n DO BEGIN
            scale := scale+abs(a[i,k])
         END;
         IF (scale <> 0.0) THEN BEGIN
            FOR k := l to n DO BEGIN
               a[i,k] := a[i,k]/scale;
               s := s+a[i,k]*a[i,k]
            END;
            f := a[i,l];
            g := -sign(sqrt(s),f);
            h := f*g-s;
            a[i,l] := f-g;
            FOR k := l to n DO BEGIN
               rv1[k] := a[i,k]/h
            END;
            IF (i <> m) THEN BEGIN
               FOR j := l to m DO BEGIN
                  s := 0.0;
                  FOR k := l to n DO BEGIN
                     s := s+a[j,k]*a[i,k]
                  END;
                  FOR k := l to n DO BEGIN
                     a[j,k] := a[j,k]
                        +s*rv1[k]
                  END
               END
            END;
            FOR k := l to n DO BEGIN
               a[i,k] := scale*a[i,k]
            END
         END
      END;
      anorm := max(anorm,(abs(w[i])+abs(rv1[i])))
   END;
   FOR i := n DOWNTO 1 DO BEGIN
      IF (i < n) THEN BEGIN
         IF (g <> 0.0) THEN BEGIN
            FOR j := l to n DO BEGIN
               v[j,i] := (a[i,j]/a[i,l])/g
            END;
            FOR j := l to n DO BEGIN
               s := 0.0;
               FOR k := l to n DO BEGIN
                  s := s+a[i,k]*v[k,j]
               END;
               FOR k := l to n DO BEGIN
                  v[k,j] := v[k,j]+s*v[k,i]
               END
            END
         END;
         FOR j := l to n DO BEGIN
            v[i,j] := 0.0;
            v[j,i] := 0.0
         END
      END;
      v[i,i] := 1.0;
      g := rv1[i];
      l := i
   END;
   FOR i := n DOWNTO 1 DO BEGIN
      l := i+1;
      g := w[i];
      IF (i < n) THEN BEGIN
         FOR j := l to n DO BEGIN
            a[i,j] := 0.0
         END
      END;
      IF (g <> 0.0) THEN BEGIN
         g := 1.0/g;
         IF (i <> n) THEN BEGIN
            FOR j := l to n DO BEGIN
               s := 0.0;
               FOR k := l to m DO BEGIN
                  s := s+a[k,i]*a[k,j]
               END;
               f := (s/a[i,i])*g;
               FOR k := i to m DO BEGIN
                  a[k,j] := a[k,j]+f*a[k,i]
               END
            END
         END;
         FOR j := i to m DO BEGIN
            a[j,i] := a[j,i]*g
         END
      END ELSE BEGIN
         FOR j := i to m DO BEGIN
            a[j,i] := 0.0
         END
      END;
      a[i,i] := a[i,i]+1.0
   END;
   FOR k := n DOWNTO 1 DO BEGIN
      FOR its := 1 to 30 DO BEGIN
         FOR l := k DOWNTO 1 DO BEGIN
            nm := l-1;
            IF ((abs(rv1[l])+anorm) = anorm) THEN GOTO 2;
            IF ((abs(w[nm])+anorm) = anorm) THEN GOTO 1
         END;
1:         c := 0.0;
         s := 1.0;
         FOR i := l to k DO BEGIN
            f := s*rv1[i];
            IF ((abs(f)+anorm) <> anorm) THEN BEGIN
               g := w[i];
               h := sqrt(f*f+g*g);
               w[i] := h;
               h := 1.0/h;
               c := (g*h);
               s := -(f*h);
               FOR j := 1 to m DO BEGIN
                  y := a[j,nm];
                  z := a[j,i];
                  a[j,nm] := (y*c)+(z*s);
                  a[j,i] := -(y*s)+(z*c)
               END
            END
         END;
2:         z := w[k];
         IF (l = k) THEN BEGIN
            IF (z < 0.0) THEN BEGIN
               w[k] := -z;
               FOR j := 1 to n DO BEGIN
               v[j,k] := -v[j,k]
            END
         END;
         GOTO 3
         END;
         IF (its = 30) THEN BEGIN
            writeln ('no convergence in 30 SVDCMP iterations'); readln
         END;
         x := w[l];
         nm := k-1;
         y := w[nm];
         g := rv1[nm];
         h := rv1[k];
         f := ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
         g := sqrt(f*f+1.0);
         f := ((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x;
         c := 1.0;
         s := 1.0;
         FOR j := l to nm DO BEGIN
            i := j+1;
            g := rv1[i];
            y := w[i];
            h := s*g;
            g := c*g;
            z := sqrt(f*f+h*h);
            rv1[j] := z;
            c := f/z;
            s := h/z;
            f := (x*c)+(g*s);
            g := -(x*s)+(g*c);
            h := y*s;
            y := y*c;
            FOR nm := 1 to n DO BEGIN
               x := v[nm,j];
               z := v[nm,i];
               v[nm,j] := (x*c)+(z*s);
               v[nm,i] := -(x*s)+(z*c)
            END;
            z := sqrt(f*f+h*h);
            w[j] := z;
            IF (z <> 0.0) THEN BEGIN
               z := 1.0/z;
               c := f*z;
               s := h*z
            END;
            f := (c*g)+(s*y);
            x := -(s*g)+(c*y);
            FOR nm := 1 to m DO BEGIN
               y := a[nm,j];
               z := a[nm,i];
               a[nm,j] := (y*c)+(z*s);
               a[nm,i] := -(y*s)+(z*c)
            END
         END;
         rv1[l] := 0.0;
         rv1[k] := f;
         w[k] := x
      END;
3:   END
END;
PROCEDURE svdfit(x,y,sig: glndata; ndata: integer; VAR a: glmma; mma: integer;
      VAR u: glmpbynp; VAR v: glnpbynp; VAR w: glnparray; mp,np:
      integer; VAR chisq: real);
(* Programs using routine SVDFIT must define the types
TYPE
   glndata = ARRAY [1..ndata] OF real;
   glmma = ARRAY [1..mma] OF real;
   glmpbynp = ARRAY [1..mp,1..np] OF real;
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glnparray = ARRAY [1..np] OF real;
in the main routine. Implementations without conformant arrays require mma=np
and the declaration
   glnparray = glmma;
and also mp=ndata, with the declaration
   glmparray = glndata;
in the main routine for use by the procedure SVBKSB. All user programs must also
include a routine that computes the desired basis functions with the declaration
PROCEDURE func(x: real: VAR p: glmma; mma: integer);
which should return the values of first mma basis functions at x in p.
FPOLY and FLEG (renamed to FUNC) are possible choices. *)
CONST
   tol=1.0e-5;
VAR
   j,i: integer;
   wmax,tmp,thresh,sum: real;
   b: glndata;
   afunc: glmma;
BEGIN
   FOR i := 1 to ndata DO BEGIN
      func(x[i],afunc,mma);
      tmp := 1.0/sig[i];
      FOR j := 1 to mma DO u[i,j] := afunc[j]*tmp;
      b[i] := y[i]*tmp
   END;
   svdcmp(u,ndata,mma,mp,np,w,v);
   wmax := 0.0;
   FOR j := 1 to mma DO IF (w[j] > wmax) THEN wmax := w[j];
   thresh := tol*wmax;
   FOR j := 1 to mma DO IF (w[j] < thresh) THEN w[j] := 0.0;
   svbksb(u,w,v,ndata,mma,mp,np,b,a);
   chisq := 0.0;
   FOR i := 1 to ndata DO BEGIN
      func(x[i],afunc,mma);
      sum := 0.0;
      FOR j := 1 to mma DO sum := sum+a[j]*afunc[j];
      chisq := chisq+sqr((y[i]-sum)/sig[i])
   END
END;
PROCEDURE svdvar(v: glnpbynp; ma,np: integer; w: glnparray;
       VAR cvm: glcvm; ncvm: integer);
(* Programs using routine SVDVAR must define the types
TYPE
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glnparray = ARRAY [1..np] OF real;
   glcvm = ARRAY [1..ncvm,1..ncvm] OF real;
in the calling program. *)
VAR
   k,j,i: integer;
   sum: real;
   wti: glnparray;
BEGIN
   FOR i := 1 to ma DO BEGIN
      wti[i] := 0.0;
      IF (w[i] <> 0.0) THEN  wti[i] := 1.0/(w[i]*w[i])
   END;
   FOR i := 1 to ma DO BEGIN
      FOR j := 1 to i DO BEGIN
         sum := 0.0;
         FOR k := 1 to ma DO BEGIN
            sum := sum+v[i,k]*v[j,k]*wti[k]
         END;
         cvm[i,j] := sum;
         cvm[j,i] := sum
      END
   END
END;
