PROGRAM d17r2(input,output);
(* driver for routine ADI *)
LABEL 1;
CONST
   jmax=11;
   pi=3.1415926;
TYPE
   double = real;
   gljmax = ARRAY [1..jmax,1..jmax] OF double;
VAR
   alim,alpha,beta,eps : double;
   i,j,k,mid,twotok : integer;
   a,b,c,d,e,f,g,u : gljmax;

(*$I ADI.PAS *)

BEGIN
   FOR i := 1 to jmax DO BEGIN
      FOR j := 1 to jmax DO BEGIN
         a[i,j] := -1.0;
         b[i,j] := 2.0;
         c[i,j] := -1.0;
         d[i,j] := -1.0;
         e[i,j] := 2.0;
         f[i,j] := -1.0;
         g[i,j] := 0.0;
         u[i,j] := 0.0
      END
   END;
   mid := (jmax DIV 2)+1;
   g[mid,mid] := 2.0;
   alpha := 2.0*(1.0-cos(pi/jmax));
   beta := 2.0*(1.0-cos((jmax-1)*pi/jmax));
   alim := ln(4.0*jmax/pi);
   k := 0;
   twotok := 1;
   REPEAT
      k := k+1;
      twotok := 2*twotok;
   UNTIL twotok >= alim;
   eps := 1.0e-4;
   adi(a,b,c,d,e,f,g,u,jmax,k,alpha,beta,eps);
   writeln('ADI Solution:');
   FOR i := 1 to jmax DO BEGIN
      FOR j := 1 to jmax DO write(u[i,j]:7:2);
      writeln
   END;
   writeln;
   writeln('Test that solution satisfies difference eqns:');
   FOR i := 2 to (jmax-1) DO BEGIN
      FOR j := 2 to (jmax-1) DO BEGIN
         g[i,j] := -4.0*u[i,j]+u[i+1,j]
            +u[i-1,j]+u[i,j-1]+u[i,j+1]
      END;
      write(' ':7);
      FOR j := 2 to (jmax-1) DO write(g[i,j]:7:2);
      writeln
   END
END.
PROGRAM d10r5(input,output);
(* driver for routine AMOEBA *)
CONST
   np=3;
   mp=4;
   ftol=1.0e-6;
TYPE
   glmpnp = ARRAY [1..mp,1..np] OF real;
   glmp = ARRAY [1..mp] OF real;
   glnp = ARRAY [1..np] OF real;
VAR
   i,iter,j,ndim : integer;
   x : glnp;
   y : glmp;
   p : glmpnp;

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

FUNCTION func(x: glnp): real;
(* calling function must define type
TYPE
   glnp = ARRAY [1..np] OF real;
where np is the physical dimension of the argument x. *)
BEGIN
   func := 0.6-bessj0(sqr(x[1]-0.5)+sqr(x[2]-0.6)+sqr(x[3]-0.7))
END;

(*$I AMOEBA.PAS *)

BEGIN
   p[1,1] := 0.0; p[1,2] := 0.0; p[1,3] := 0.0;
   p[2,1] := 1.0; p[2,2] := 0.0; p[2,3] := 0.0;
   p[3,1] := 0.0; p[3,2] := 1.0; p[3,3] := 0.0;
   p[4,1] := 0.0; p[4,2] := 0.0; p[4,3] := 1.0;
   ndim := np;
   FOR i := 1 to mp DO BEGIN
      FOR j := 1 to np DO BEGIN
         x[j] := p[i,j]
      END;
      y[i] := func(x)
   END;
   amoeba(p,y,ndim,ftol,iter);
   writeln;
   writeln('Iterations: ',iter:3);
   writeln('Vertices of final 3-d simplex and');
   writeln('function values at the vertices:');
   writeln;
   writeln('i':3,
      'x[i]':10,'y[i]':12,'z[i]':12,'function':14);
   writeln;
   FOR i := 1 to mp DO BEGIN
      write(i:3);
      FOR j := 1 to np DO write(p[i,j]:12:6);
      writeln(y[i]:12:6)
   END;
   writeln;
   writeln('True minimum is at (0.5,0.6,0.7)')
END.
PROGRAM D10R13(input,output);
CONST
        ncity=10;
TYPE
        cityarray = ARRAY [1..ncity] OF real;
        iarray = ARRAY [1..ncity] OF integer;
VAR
        glinext,glinextp,gljdum,idum,i,ii : integer;
        glma : ARRAY [1..55] OF real;
        x,y : cityarray;
        iorder : iarray;

(*$I MODFILE.PAS*)

(*$I RAN3.PAS*)

(*$I IRBIT1.PAS*)

(*$I ANNEAL.PAS*)
 
BEGIN
        gljdum := 1;
        idum := -1;
        FOR i := 1 TO ncity DO BEGIN
                x[i] := ran3(idum);
                y[i] := ran3(idum);
                iorder[i] := i;
        END;
        anneal(x,y,iorder,ncity);
        writeln('*** System Frozen ***');
        writeln('Final path:');
        writeln(' ':3,'city',' ':6,'x',' ':9,'y');
        FOR i := 1 TO ncity DO BEGIN
                ii := iorder[i];
                writeln(ii:4,x[ii]:10:4,y[ii]:10:4)
        END
END.
PROGRAM d13r4(input,output);
(* driver for routine AVEVAR *)
CONST
   npts=1000;
   eps=0.1;
TYPE
   glnparray = ARRAY [1..npts] OF real;
VAR
   gliset : integer;
   glgset : real;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   i,idum,j : integer;
   ave,shift,vrnce : real;
   data : glnparray;

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

(*$I GASDEV.PAS *)

(*$I AVEVAR.PAS *)

BEGIN
(* generate gaussian distributed data *)
   gliset := 0;
   idum := -5;
   writeln('shift':9,'average':11,'variance':12);
   FOR i := 1 to 11 DO BEGIN
      shift := (i-1)*eps;
      FOR j := 1 to npts DO
         data[j] := shift+i*gasdev(idum);
      avevar(data,npts,ave,vrnce);
      writeln(shift:8:2,ave:11:2,vrnce:12:2)
   END
END.
