PROGRAM d1r4(input,output,dfile);
(* driver for routine CALDAT *)
TYPE
   month = PACKED ARRAY [1..10] OF char;
VAR
   i,id,idd,im,imm,iy,iyy,j,n : integer;
   name : ARRAY [1..12] OF month;
   dfile : text;

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

(*$I CALDAT.PAS *)

BEGIN
(* check whether caldat properly undoes the operation of julday *)
   name[1] := 'january   ';
   name[2] := 'february  ';
   name[3] := 'march     ';
   name[4] := 'april     ';
   name[5] := 'may       ';
   name[6] := 'june      ';
   name[7] := 'july      ';
   name[8] := 'august    ';
   name[9] := 'september ';
   name[10] := 'october   ';
   name[11] := 'november  ';
   name[12] := 'december  ';
   glopen(dfile,'dates1.dat');
   readln(dfile);
   readln(dfile,n);
   writeln;
   writeln('original date:','reconstructed date:':43);
   writeln('month':5,'day':8,'year':6,'julian day':15,
      'month':9,'day':8,'year':6);
   FOR i := 1 to n DO BEGIN
      readln(dfile,im,id,iy);
      j := julday(im,id,iy);
      caldat(j,imm,idd,iyy);
      writeln(name[im],id:3,iy:6,j:13,name[imm]:16,idd:3,iyy:6)
   END
END.
PROGRAM d6r29(input,output);
(* driver for routine CEL *)
CONST
   pio2=1.5707963;
VAR
   glit : integer;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gla,glb,glp,glakc : real;
   ago,astop,s : real;
   i,idum : integer;

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

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

(*$I TRAPZD.PAS *)

(*$I QSIMP.PAS *)

(*$I CEL.PAS *)

BEGIN
   writeln('complete elliptic integral');
   writeln('kc':7,'p':10,'a':10,'b':10,'cel':11,'integral':12);
   idum := -55;
   ago := 0.0;
   astop := pio2;
   FOR i := 1 to 20 DO BEGIN
      glakc := 0.1+ran3(idum);
      gla := 10.0*ran3(idum);
      glb := 10.0*ran3(idum);
      glp := 0.1+ran3(idum);
      qsimp(ago,astop,s);
      writeln(glakc:10:6,glp:10:6,gla:10:6,glb:10:6,
         cel(glakc,glp,gla,glb):10:6,s:10:6)
   END
END.
PROGRAM d5r7(input,output);
(* driver for routine CHDER *)
LABEL 10,99;
CONST
   nval=40;
   pio2=1.5707963;
TYPE
   glcarray=ARRAY [1..nval] OF real;
VAR
   a,b,x : real;
   i,mval : integer;
   c,cder : glcarray;

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

FUNCTION fder(x: real): real;
BEGIN
   fder := 4.0*x*(sqr(x)-1.0)*sin(x)+sqr(x)*(sqr(x)-2.0)*cos(x)
END;

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

(*$I CHEBFT.PAS *)

(*$I CHDER.PAS *)

BEGIN
   a := -pio2;
   b := pio2;
   chebft(a,b,c,nval);
(* test derivative *)
10:   writeln;
   writeln('How many terms in Chebyshev evaluation?');
   write('Enter n between 6 and ',nval:2,
         '. (n := 0 to end).  ');
   readln (mval);
   IF ((mval <= 0) OR (mval > nval)) THEN GOTO 99;
   chder(a,b,c,cder,mval);
   writeln;
   writeln('x':9,'actual':14,'Cheby. deriv.':16);
   FOR i := -8 to 8 DO BEGIN
      x := i*pio2/10.0;
      writeln(x:12:6,fder(x):12:6,chebev(a,b,cder,mval,x):12:6)
   END;
   GOTO 10;
99:
END.
PROGRAM d5r5(input,output);
(* driver for routine CHEBEV *)
LABEL 10,99;
CONST
   nval=40;
   pio2=1.5707963;
TYPE
   glcarray=ARRAY [1..nval] OF real;
VAR
   a,b,x : real;
   i,mval : integer;
   c : glcarray;

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

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

(*$I CHEBEV.PAS *)

BEGIN
   a := -pio2;
   b := pio2;
   chebft(a,b,c,nval);
(* test chebyshev evaluation routine *)
10:   writeln;
   writeln('How many terms in Chebyshev evaluation?');
   write('Enter n between 6 and ',nval:2,
         '. (n := 0 to end).  ');
   readln(mval);
   IF ((mval <= 0) OR (mval > nval)) THEN GOTO 99;
   writeln;
   writeln('x':9,'actual':14,'chebyshev fit':16);
   FOR i := -8 to 8 DO BEGIN
      x := i*pio2/10.0;
      writeln(x:12:6,func(x):12:6,chebev(a,b,c,mval,x):12:6)
   END;
   GOTO 10;
99:
END.
PROGRAM d5r4(input,output);
(* driver for routine CHEBFT *)
LABEL 10,99;
CONST
   nval=40;
   pio2=1.5707963;
   eqs=1.0e-6;
TYPE
   glcarray=ARRAY [1..nval] OF real;
VAR
   a,b,dum,f : real;
   t0,t1,term,x,y : real;
   i,j,mval : integer;
   c : glcarray;

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

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

BEGIN
   a := -pio2;
   b := pio2;
   chebft(a,b,c,nval);
(* test result *)
10:   writeln;
   writeln('How many terms in Chebyshev evaluation?');
   write('Enter n between 6 and ',nval:2,
         '. (n := 0 to end).  ');
   readln(mval);
   IF ((mval <= 0) OR (mval > nval)) THEN GOTO 99;
   writeln;
   writeln('x':9,'actual':14,'chebyshev fit':16);
   FOR i := -8 to 8 DO BEGIN
      x := i*pio2/10.0;
      y := (x-0.5*(b+a))/(0.5*(b-a));
(* evaluate chebyshev polynomial without using routine chebev *);
      t0 := 1.0;
      t1 := y;
      f := c[2]*t1+c[1]*0.5;
      FOR j := 3 to mval DO BEGIN
         dum := t1;
         t1 := 2.0*y*t1-t0;
         t0 := dum;
         term := c[j]*t1;
         f := f+term
      END;
      writeln(x:12:6,func(x):12:6,f:12:6)
   END;
   GOTO 10;
99:
END.
PROGRAM d5r8(input,output);
(* driver for routine CHEBPC *)
LABEL 10,99;
CONST
   nval=40;
   pio2=1.5707963;
TYPE
   glcarray= ARRAY [1..nval] OF real;
VAR
   a,b,poly,x,y : real;
   i,j,mval : integer;
   c,d : glcarray;

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

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

(*$I CHEBPC.PAS *)

BEGIN
   a := -pio2;
   b := pio2;
   chebft(a,b,c,nval);
10:   writeln;
   writeln('How many terms in Chebyshev evaluation?');
   write('Enter n between 6 and ',nval:2,
         '. (n := 0 to end).  ');
   readln(mval);
   IF ((mval <= 0) OR (mval > nval)) THEN GOTO 99;
   chebpc(c,d,mval);
(* test polynomial *)
   writeln;
   writeln('x':9,'actual':14,'polynomial':14);
   FOR i := -8 to 8 DO BEGIN
      x := i*pio2/10.0;
      y := (x-(0.5*(b+a)))/(0.5*(b-a));
      poly := d[mval];
      FOR j := mval-1 DOWNTO 1 DO BEGIN
         poly := poly*y+d[j]
      END;
      writeln(x:12:6,func(x):12:6,poly:12:6)
   END;
   GOTO 10;
99:
END.
PROGRAM d5r6(input,output);
(* driver for routine CHINT *)
LABEL 10,99;
CONST
   nval=40;
   pio2=1.5707963;
TYPE
   glcarray=ARRAY [1..nval] OF real;
VAR
   a,b,x : real;
   i,mval : integer;
   c,cint : glcarray;

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

FUNCTION fint(x: real): real;
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 CHEBFT.PAS *)

(*$I CHEBEV.PAS *)

(*$I CHINT.PAS *)

BEGIN
   a := -pio2;
   b := pio2;
   chebft(a,b,c,nval);
(* test integral *)
10:   writeln;
   writeln('How many terms in Chebyshev evaluation?');
   write('Enter n between 6 and ',nval:2,
         '. (n := 0 to end).  ');
   readln(mval);
   IF ((mval <= 0) OR (mval > nval)) THEN GOTO 99;
   chint(a,b,c,cint,mval);
   writeln;
   writeln('x':9,'actual':14,'Cheby. integ.':16);
   FOR i := -8 to 8 DO BEGIN
      x := i*pio2/10.0;
      writeln(x:12:6,fint(x)-fint(-pio2):12:6,
            chebev(a,b,cint,mval,x):12:6)
   END;
   GOTO 10;
99:
END.
PROGRAM d13r8(input,output);
(* driver for routine CHSONE *)
CONST
   nbins=10;
   npts=2000;
TYPE
   barray = ARRAY [0..nbins] OF real;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   chsq,df,prob,x : real;
   i,ibin,idum,j : integer;
   bins,ebins : barray;

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

(*$I GCF.PAS *)

(*$I GSER.PAS *)

(*$I GAMMQ.PAS *)

(*$I RAN3.PAS *)

(*$I EXPDEV.PAS *)

(*$I CHSONE.PAS *)

BEGIN
   idum := -15;
   FOR j := 1 to nbins DO BEGIN
      bins[j] := 0.0
   END;
   FOR i := 1 to npts DO BEGIN
      x := expdev(idum);
      ibin := trunc(x*nbins/3.0)+1;
      IF (ibin <= nbins) THEN
         bins[ibin] := bins[ibin]+1.0
   END;
   FOR i := 1 to nbins DO BEGIN
      ebins[i] := 3.0*npts/nbins*exp(-3.0*(i-0.5)/nbins)
   END;
   chsone(bins,ebins,nbins,-1,df,chsq,prob);
   writeln('expected':15,'observed':15);
   FOR i := 1 to nbins DO BEGIN
      writeln(ebins[i]:14:2,bins[i]:15:2)
   END;
   writeln;
   writeln('chi-squared:':19,chsq:10:4);
   writeln('probability:':19,prob:10:4)
END.
PROGRAM d13r9(input,output);
(* driver for routine CHSTWO *)
CONST
   nbins=10;
   npts=2000;
TYPE
   barray = ARRAY [1..nbins] OF real;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   chsq,df,prob,x : real;
   i,ibin,idum,j : integer;
   bins1,bins2 : barray;

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

(*$I GSER.PAS *)

(*$I GCF.PAS *)

(*$I GAMMQ.PAS *)

(*$I RAN3.PAS *)

(*$I EXPDEV.PAS *)

(*$I CHSTWO.PAS *)

BEGIN
   idum := -18;
   FOR j := 1 to nbins DO BEGIN
      bins1[j] := 0.0;
      bins2[j] := 0.0
   END;
   FOR i := 1 to npts DO BEGIN
      x := expdev(idum);
      ibin := trunc(x*nbins/3.0)+1;
      IF (ibin <= nbins) THEN bins1[ibin] := bins1[ibin]+1.0;
      x := expdev(idum);
      ibin := trunc(x*nbins/3.0)+1;
      IF (ibin <= nbins) THEN bins2[ibin] := bins2[ibin]+1.0
   END;
   chstwo(bins1,bins2,nbins,-1,df,chsq,prob);
   writeln;
   writeln('dataset 1':15,'dataset 2':15);
   FOR i := 1 to nbins DO BEGIN
      writeln(bins1[i]:13:2,bins2[i]:15:2)
   END;
   writeln;
   writeln('chi-squared:':18,chsq:12:4);
   writeln('probability:':18,prob:12:4)
END.
PROGRAM d13r13(input,output,dfile);
(* driver for routine CNTAB1 *)
(* contingency table in file TABLE1.DAT *)
CONST
   ndat=9;
   nmon=12;
TYPE
   narray = ARRAY [1..ndat,1..nmon] OF integer;
   chr15 = string[15];
   chr5 = string[5];
VAR
   ccc,chisq,cramrv,df,prob : real;
   i,j : integer;
   nmbr : narray;
   fate : ARRAY [1..ndat] OF chr15;
   mon : ARRAY [1..nmon] OF chr5;
   txt : string[64];
   dfile : text;

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

(*$I GCF.PAS *)

(*$I GSER.PAS *)

(*$I GAMMQ.PAS *)

(*$I CNTAB1.PAS *)

BEGIN
   glopen(dfile,'table1.dat');
   readln(dfile);
   readln(dfile,txt);
   read(dfile,fate[1]);
   FOR i := 1 to 12 DO read(dfile,mon[i]);
   readln(dfile);
   readln(dfile);
   FOR i := 1 to ndat DO BEGIN
      read(dfile,fate[i]);
      FOR j := 1 to 12 DO read(dfile,nmbr[i,j]);
      readln(dfile)
   END;
   close(dfile);
   writeln;
   writeln(txt);
   writeln;
   write(' ':15);
   FOR i := 1 to 12 DO write(mon[i]:5);
   writeln;
   FOR i := 1 to ndat DO BEGIN
      write(fate[i]);
      FOR j := 1 to 12 DO write(nmbr[i,j]:5);
      writeln
   END;
   cntab1(nmbr,ndat,nmon,chisq,df,prob,cramrv,ccc);
   writeln;
   writeln('chi-squared       ',chisq:20:2);
   writeln('degrees of freedom',df:20:2);
   writeln('probability       ',prob:20:4);
   writeln('cramer-v          ',cramrv:20:4);
   writeln('contingency coeff.',ccc:20:4)
END.
PROGRAM d13r14(input,output,dfile);
(* driver for routine CNTAB2 *)
(* contingency table in file TABLE1.DAT *)
CONST
   ni=9;
   nmon=12;
TYPE
   narray = ARRAY [1..ni,1..nmon] OF integer;
   chr5 = string[5];
   chr15 = string[15];
VAR
   h,hx,hxgy,hy,hygx : real;
   uxgy,uxy,uygx : real;
   i,j : integer;
   nmbr : narray;
   fate : ARRAY [1..ni] OF chr15;
   mon : ARRAY [1..nmon] OF chr5;
   txt : string[64];
   dfile : text;

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

BEGIN
   glopen(dfile,'table1.dat');
   readln(dfile);
   readln(dfile,txt);
   read(dfile,fate[1]);
   FOR i := 1 to 12 DO read(dfile,mon[i]);
   readln(dfile);
   readln(dfile);
   FOR i := 1 to ni DO BEGIN
      read(dfile,fate[i]);
      FOR j := 1 to 12 DO read(dfile,nmbr[i,j]);
      readln(dfile)
   END;
   close(dfile);
   writeln;
   writeln(txt);
   writeln;
   write(' ':15);
   FOR i := 1 to 12 DO write(mon[i]:5);
   writeln;
   FOR i := 1 to ni DO BEGIN
      write(fate[i]);
      FOR j := 1 to 12 DO write(nmbr[i,j]:5);
      writeln
   END;
   cntab2(nmbr,ni,nmon,h,hx,hy,hygx,hxgy,uygx,uxgy,uxy);
   writeln;
   writeln('entropy of table           ',h:10:4);
   writeln('entropy of x-distribution  ',hx:10:4);
   writeln('entropy of y-distribution  ',hy:10:4);
   writeln('entropy of y given x       ',hygx:10:4);
   writeln('entropy of x given y       ',hxgy:10:4);
   writeln('dependency of y on x       ',uygx:10:4);
   writeln('dependency of x on y       ',uxgy:10:4);
   writeln('symmetrical dependency     ',uxy:10:4)
END.
PROGRAM d12r6(input,output);
(* driver for routine CONVLV *)
CONST
   n=16;      (* data array size *)
   m=9;      (* response function dimension - must be odd *)
   n2p2=34;      (* n2=2*n+2 *)
   pi=3.14159265;
TYPE
   glnarray = ARRAY [1..n] OF real;
   gln2array = ARRAY [1..n2p2] OF real;
   gl2narray = gln2array;
   gldarray = gln2array;
VAR
   i,isign,j : integer;
   cmp : real;
   ans : gln2array;
   data,respns,resp : glnarray;

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

(*$I TWOFFT.PAS *)

(*$I REALFT.PAS *)

(*$I CONVLV.PAS *)

BEGIN
   FOR i := 1 to n DO BEGIN
      data[i] := 0.0;
      IF ((i >= ((n DIV 2)-(n DIV 8))) AND
         (i <= ((n DIV 2)+(n DIV 8)))) THEN
         data[i] := 1.0
   END;
   FOR i := 1 to m DO BEGIN
      respns[i] := 0.0;
      IF ((i > 2) AND (i < 7)) THEN respns[i] := 1.0;
      resp[i] := respns[i]
   END;
   isign := 1;
   convlv(data,n,resp,m,isign,ans);
(* compare with a direct convolution *)
   writeln ('i':3,'CONVLV':14,'Expected':13);
   FOR i := 1 to n DO BEGIN
      cmp := 0.0;
      FOR j := 1 to (m DIV 2) DO BEGIN
         cmp := cmp+data[((i-j-1+n) MOD n)+1]*respns[j+1];
         cmp := cmp+data[((i+j-1) MOD n)+1]*respns[m-j+1]
      END;
      cmp := cmp+data[i]*respns[1];
      writeln (i:3,ans[i]:15:6,cmp:12:6)
   END
END.
PROGRAM d12r7(input,output);
(* driver for routine CORREL *)
CONST
   n=64;
   n2=128;      (* n2=2*nn *)
   pi=3.1415927;
TYPE
   glnarray = ARRAY [1..n] OF real;
   gl2narray = ARRAY [1..n2] OF real;
   gldarray = gl2narray;
VAR
   i,j : integer;
   cmp : real;
   data1,data2 : glnarray;
   ans : gl2narray;

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

(*$I TWOFFT.PAS *)

(*$I REALFT.PAS *)

(*$I CORREL.PAS *)

BEGIN
   FOR i := 1 to n DO BEGIN
      data1[i] := 0.0;
      IF ((i > ((n DIV 2)-(n DIV 8))) AND 
         (i < ((n DIV 2)+(n DIV 8)))) THEN
         data1[i] := 1.0;
      data2[i] := data1[i]
   END;
   correl(data1,data2,n,ans);
(* calculate directly *)
   writeln ('n':3,'CORREL':14,'direct calc.':18);
   FOR i := 0 to 16 DO BEGIN
      cmp := 0.0;
      FOR j := 1 to n DO BEGIN
         cmp := cmp+data1[((i+j-1) MOD n)+1]*data2[j]
      END;
      writeln (i:3,ans[i+1]:15:6,cmp:15:6)
   END
END.
PROGRAM d12r5(input,output);
(* driver for routine COSFT *)
LABEL 1,99;
CONST
   eps=1.0e-3;
   np=16;
   np2=18;   (* np2=np+2 *)
   width=30.0;
   pi=3.1415926;
TYPE
   gldarray=ARRAY [1..np2] OF real;
   glyarray=gldarray;
VAR
   big,per,scal,small : real;
   i,jj,j,nlim : integer;
   data,size : gldarray;

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

(*$I REALFT.PAS *)

(*$I COSFT.PAS *)

BEGIN
1:   writeln('period of cosine in channels (2-',np:2,')');
   readln(per);
   IF (per <= 0.0) THEN GOTO 99;
   FOR i := 1 to np DO BEGIN
      data[i] := cos(2.0*pi*(i-1)/per)
   END;
   cosft(data,np,+1);
   big := -1.0e10;
   small := 1.0e10;
   FOR i := 1 to np DO BEGIN
      IF  (data[i] < small) THEN small := data[i];
      IF  (data[i] > big) THEN big := data[i]
   END;
   scal := width/(big-small);
   FOR i := 1 to np DO BEGIN
      nlim := round(scal*(data[i]-small)+eps);
      write(i:4,' ');
      FOR j := 1 to nlim+1 DO write('*');
      writeln
   END;
   writeln('press RETURN to continue ...');
   readln;
   cosft(data,np,-1);
   big := -1.0e10;
   small := 1.0e10;
   FOR i := 1 to np DO BEGIN
      IF (data[i] < small) THEN small := data[i];
      IF (data[i] > big) THEN big := data[i]
   END;
   scal := width/(big-small);
   FOR i := 1 to np DO BEGIN
      nlim := round(scal*(data[i]-small)+eps);
      write(i:4,' ');
      FOR j := 1 to nlim+1 DO write('*');
      writeln
   END;
   GOTO 1;
99:
END.
PROGRAM d14r3(input,output);
(* driver for routine COVSRT *)
CONST
   ma=10;
   mfit=5;
TYPE
   glcovar = ARRAY [1..ma,1..ma] OF real;
   gllista = ARRAY [1..mfit] OF integer;
VAR
   i,j : integer;
   covar : glcovar;
   lista : gllista;

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

BEGIN
   FOR i := 1 to ma DO BEGIN
      FOR j := 1 to ma DO BEGIN
         covar[i,j] := 0.0;
         IF ((i <= 5) AND (j <= 5)) THEN BEGIN
            covar[i,j] := i+j-1
         END
      END
   END;
   writeln;
   writeln('original matrix');
   FOR i := 1 to ma DO BEGIN
      FOR j := 1 to ma DO write(covar[i,j]:4:1);
      writeln
   END;
   writeln(' press RETURN to continue...');
   readln;
(* test 1 - spread by 2 *)
   writeln;
   writeln('test #1 - spread by two');
   FOR i := 1 to mfit DO BEGIN
      lista[i] := 2*i
   END;
   covsrt(covar,ma,ma,lista,mfit);
   FOR i := 1 to ma DO BEGIN
      FOR j := 1 to ma DO write(covar[i,j]:4:1);
      writeln
   END;
   writeln(' press RETURN to continue...');
   readln;
(* test 2 - reverse *)
   writeln;
   writeln('test #2 - reverse');
   FOR i := 1 to ma DO BEGIN
      FOR j := 1 to ma DO BEGIN
         covar[i,j] := 0.0;
         IF  ((i <= 5) AND (j <= 5)) THEN BEGIN
            covar[i,j] := i+j-1
         END
      END
   END;
   FOR i := 1 to mfit DO BEGIN
      lista[i] := mfit+1-i
   END;
   covsrt(covar,ma,ma,lista,mfit);
   FOR i := 1 to ma DO BEGIN
      FOR j := 1 to ma DO write(covar[i,j]:4:1);
      writeln
   END;
   writeln(' press RETURN to continue...');
   readln;
(* test 3 - spread and reverse *)
   writeln;
   writeln('test #3 - spread and reverse');
   FOR i := 1 to ma DO BEGIN
      FOR j := 1 to ma DO BEGIN
         covar[i,j] := 0.0;
         IF ((i <= 5) AND (j <= 5)) THEN BEGIN
            covar[i,j] := i+j-1
         END
      END
   END;
   FOR i := 1 to mfit DO BEGIN
      lista[i] := ma+2-2*i
   END;
   covsrt(covar,ma,ma,lista,mfit);
   FOR i := 1 to ma DO BEGIN
      FOR j := 1 to ma DO write(covar[i,j]:4:1);
      writeln
   END
END.
PROGRAM d13r17(input,output,dfile);
(* driver for routine CRANK *)
CONST
   ndat=20;
   nmon=12;
TYPE
   narray = ARRAY [1..ndat] OF real;
   glsarray = narray;
   cityname = string[15];
   monthname = string[4];
VAR
   i,j : integer;
   data,order,s : narray;
   rays : ARRAY [1..ndat,1..nmon] OF real;
   city : ARRAY [1..ndat] OF cityname;
   mon : ARRAY [1..12] OF monthname;
   txt : string[64];
   dfile : text;

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

(*$I CRANK.PAS *)

BEGIN
   glopen(dfile,'table2.dat');
   readln(dfile);
   readln(dfile,txt);
   read(dfile,city[1]);
   FOR i := 1 to 12 DO read(dfile,mon[i]);
   readln(dfile);
   readln(dfile);
   FOR i := 1 to ndat DO BEGIN
      read(dfile,city[i]);
      FOR j := 1 to 12 DO read(dfile,rays[i,j]);
      readln(dfile)
   END;
   close(dfile);
   writeln(txt);
   write(' ':15);
   FOR i := 1 to 12 DO write(mon[i]:4);
   writeln;
   FOR i := 1 to ndat DO BEGIN
      write(city[i]);
      FOR j := 1 to 12 DO write(round(rays[i,j]):4);
      writeln
   END;
   writeln(' press return to continue ...');
   readln;
(* replace solar flux in each column by rank order *)
   FOR j := 1 to 12 DO BEGIN
      FOR i := 1 to ndat DO BEGIN
         data[i] := rays[i,j];
         order[i] := i
      END;
      sort2(ndat,data,order);
      crank(ndat,data,s[i]);
      FOR i := 1 to ndat DO BEGIN
         rays[round(order[i]),j] := data[i]
      END
   END;
   write(' ':15);
   FOR i := 1 to 12 DO write(mon[i]:4);
   writeln;
   FOR i := 1 to ndat DO BEGIN
      write(city[i]);
      FOR j := 1 to 12 DO BEGIN
         write(round(rays[i,j]):4)
      END;
      writeln
   END
END.
