/********************************************************/
/*  Program M V1.1 -- Term Project for 18.336           */
/*                                                      */
/*  Simulate low Peierls barrier dislocation motion     */
/*  impeded by a random array of point-like obstacles.  */
/*  A successful theoretical analysis neglecting mass   */
/*  and dissipation has been given by Friedel (1956),   */
/*  and later was greatly extended by Schwarz and       */
/*  Labush (J. Appl. Phys. 49, 5174, 1978). A review    */
/*  was given by Ardell, Metal.Trans. A16, 2131(1985).  */
/*                                                      */
/*                            Ju Li, Nov.28, 1998 MIT   */
/********************************************************/

/*
  The eqn of motion for nearly straight dislocation is 
  
  m(d^2y/dt^2)=e(d^2y/dx^2)-B(dy/dt)+sb-F(x,y|c,a,K) 
  
  where m: dislocation mass / length = e/vs^2; (E=mc^2)
        e: dislocation line tension = ub^2/2;
	u: material shear modulus;
	B: phonon drag = kT/vs/b^2 (for T>>Debye temperature);
	s: external applied stress;
	c: obstacle volume concentration
	a: average obstacle size
	K: maximum impeding force the obstacle offers.

  If we use reduced units in which
        length in b; time in b/vs; energy in e*b;
	mass/length in e/vs^2; stress in e/b^2=u/2;
	B in e/b/vs; c in 1/b^2; K in e;
  then
        b = 1; vs = 1; e = 1;
        m = 1; u = 2; B = kT/e/b;

  All above are in fact finite-temperature formulas.

  The obstacle is assumed to be point-like: c << K/2e, and its
  avg size a provides the smallest length scale to the problem.
  The resistance mechanism is assumed to be APB generation,
  so K/2e = 2*X*a, X=0.1-0.5J/m^2. The dislocation is arrested
  (freed) if the extending angle of nearby segments <> K/e, and
  these processes take on the smallest length/time scales.

  When K/2e>1, the obstacles is impenetrable (precipitate ->
  dispersed-particle), but the dislocation can still circumvent
  the obstacle by Orowan loop shedding. The actual process is
  complicated because of dislocation self-interaction, so I
  simply stipulate that upper limit of extending angle is 2*Pi/3.
  
  One can further include thermal activation, which allows the
  arrested [a/b] atoms to collectively overcome the barrier with
  probability rate [vs/b] / [a/b] * exp(G/kT), where [vs/b] is
  the Debye frequency, G = [DK*a] is the free energy barrier as
  DK has been subtracted off the currently applied force.

  Full documentation on the Lagrangian formulation of Nodal Dynamics
  and the abacus discretization model can be found at
  http://www.mit.edu/~liju99/paper_archive.html.

*/

/* Threaded 2D X-Window visualization (SGI) */
#include "V.h"
/** Universal Constants **/
#define PI 3.14159265358979323846
#define AVO 6.0221367E+23
#define BOLZ 1.380658E-23  /* in Joule/K */
#define HBAR 1.054589E-34  /* in Joule.s */
#define EV_IN_J 1.60217733E-19
/* Copper parameters */
#define Cu_AMU (63.546) /* in g/mol */
#define Cu_DENSITY (8920) /* in kg/m^3 */
#define Cu_SHEAR_MODULUS (53.96E9) /* in Pa */
#define Cu_BURGER_VECTOR (pow(Cu_AMU*1E-3/AVO/Cu_DENSITY*4,1./3)/sqrt(2.))
#define Cu_MELTING_POINT (1357.77) /* in K */
/* Fundamental Units in SI */
#define UPROTOTYPE_NAME "Cu"
#define u_scale_law "u(T)=u(0)*(1-T/2/T_melt)"
#define u_scale(T) (1-(T)/2/Cu_MELTING_POINT)
#define ULENGTH_IN_M(T) Cu_BURGER_VECTOR /* assume no thermal expansion */
#define USPEED_IN_MS(T) sqrt(Cu_SHEAR_MODULUS*u_scale(T)/Cu_DENSITY)
#define UTENSION_IN_JM(T) \
(Cu_SHEAR_MODULUS*u_scale(T)/2*Cu_BURGER_VECTOR*Cu_BURGER_VECTOR)
/* Derived Units in SI */
#define UENERGY_IN_J(T) (UTENSION_IN_JM(T)*ULENGTH_IN_M(T))
#define UTIME_IN_S(T)   (ULENGTH_IN_M(T)/USPEED_IN_MS(T))
#define UMASS_IN_KG(T)  (UENERGY_IN_J(T)/USPEED_IN_MS(T)/USPEED_IN_MS(T))
#define USTRESS_IN_PA(T) (UTENSION_IN_JM(T)/ULENGTH_IN_M(T)/ULENGTH_IN_M(T))
#define MY_BOLZ(T) (BOLZ/UENERGY_IN_J(T))

/* Periodic BC */
#define pre(y,i) ((y)[((i)+mesh-1)%mesh])
#define aft(y,i) ((y)[((i)+1)%mesh])
#define prn(y,i,n) ((y)[((i)+mesh-(n))%mesh])
#define afn(y,i,n) ((y)[((i)+(n))%mesh])
#define TINY (1E-8) /* slightly behind the obstacle */
#define STARTY (0.5) /* starting position */
#define STARTW (0.00) /* starting wavy amplitude */
#define STARTT (1./3) /* begin to collect averages */
#define DN (3) /* surrounding nodes to calculate the angle */
#define POWER (2.) /* nth-norm measure for deviations */
#define CHUNK_OF_WORK (0.1) /* notify when much is done */
#define MOBILE_IN_M2 (1E8) /* assumed mobile dislocation density */
   
static int mesh, nob, iseed, *xob, last_print=0;
static double a, c, maxangle, T, B, ss, xmax, ymax, tmax,
    r, r2, xdel, xdel2, tdel, tdel2, t, *y, *ylast, *ynew, *dl, *sl,
    *d, *s, *yt, *ypast, excess, *yob, yavg, ystart, tstart,
    o_free=0, b_free=0, a_free=0, arrested_time=0, waste_time;
static Bool visualization, *arrested, orowan_free, break_free,
    activate_free, final_average=False;
static char buf[128], token[32]={0}, *p, replace[1]={0};
static FILE *con, *save;
static void print_run_conditions (FILE *f);
static void print_predictions (FILE *f);
static void print_statistics (FILE *f);

void print_run_conditions (FILE *f)
{
    fprintf(f,"\n=================== Running Conditions ==================\n");
    fprintf(f,"Assuming a low Peierls barrier metal (%s) with\n",
	    UPROTOTYPE_NAME);
    fprintf(f,"b(T=0)=%.3f A, vs(T=0)=%.1f m/s, u(T=0)=%.1f GPa\n",
	    ULENGTH_IN_M(0)*1E10, USPEED_IN_MS(0),
	    2*UTENSION_IN_JM(0)/ULENGTH_IN_M(0)/ULENGTH_IN_M(0)*1E-9);
    fprintf(f,"and a u-scaling law of %s;\n", u_scale_law);
    fprintf(f,"--------------------------------------------------------\n");
    fprintf(f,"At T = %.2f K, u(T) = %.2f GPa = %.1f%% u(0),\n",
	    T, 2*UTENSION_IN_JM(T)/ULENGTH_IN_M(T)/ULENGTH_IN_M(T)*1E-9,
	    u_scale(T)*100);
    fprintf(f,"b(T) = %.3f A, vs(T) = %.1f m/s, e = %.3f eV/A,\n",
	    ULENGTH_IN_M(T)*1E10, USPEED_IN_MS(T),
	    UTENSION_IN_JM(T)*1E-10/EV_IN_J);
    fprintf(f,"and phonon drag coeff. B = %f [e/b/vs](T);\n", B);
    fprintf(f,"--------------------------------------------------------\n");
    fprintf(f,"Simulation cell xmax = %d [b](T) = %f um,\n",
	    (int)xmax, xmax*ULENGTH_IN_M(T)*1E6);
    fprintf(f,"                ymax = %d [b](T) = %f um;\n",
	    (int)ymax, ymax*ULENGTH_IN_M(T)*1E6);
    fprintf(f,"Obstacle volume concentration c = %.4e = %.4f%%,\n",
	    c, c*100);
    fprintf(f,"         size a = %d [b](T), nob = %d;\n", (int)a, nob);
    fprintf(f,"--------------------------------------------------------\n");    
    fprintf(f,"Obstacle maxangle = K/2e = %.3f, should >> c=%.5f,\n",
	    maxangle, c);
    fprintf(f,"from which we infer that X_APB(T) = %f J/m^2, \n",
	    maxangle/a*UTENSION_IN_JM(T)/ULENGTH_IN_M(T));
    fprintf(f,"      and total energy barrier(T) = %f eV;\n",
	    2*maxangle*a*UENERGY_IN_J(T)/EV_IN_J);
    fprintf(f,"stress = %.6f [e/b^2](T) = %.5f%% u(T) = %.2f MPa;\n",
	    ss, ss/2*100, ss*USTRESS_IN_PA(T)*1E-6);
    fprintf(f,"--------------------------------------------------------\n");    
    fprintf(f,"mesh = %d, xdel = %d [b](T) = %.2f A,\n",
	    mesh, (int)xdel, xdel*ULENGTH_IN_M(T)*1E10);
    fprintf(f,"r = %.3f,  tdel = %.2f [b/vs](T) = %.3f ps,\n",
	    r, tdel, tdel*UTIME_IN_S(T)*1E12);
    fprintf(f,"tmax = %e [b/vs](T) = %d timesteps.\n",
	    tmax, (int)(tmax/tdel));
    fprintf(f,"random number generator seed = %d.\n", iseed);
    fprintf(f,"=========================================================\n\n");
    return;
} /* end print_run_conditions() */


void print_predictions (FILE *f)
{
    double L, s_theory, scale;
    fprintf (f,"\n-------------- Theoretical Predictions ------------------\n");
    fprintf (f,"Full speed = s/B = %f [vs](T) = %.3f m/s;\n",
	     ss/B, ss/B*USPEED_IN_MS(T));
    L = sqrt(xmax*ymax/nob);
    fprintf (f,"L = sqrt(avg_obstacle_area) = %.2f [b](T) = %.2f A,\n",
	     L, L*ULENGTH_IN_M(T)*1E10);
    fprintf (f,"K/2e = %.2f %s, theor. flow stress 2e/b/L*(K/2e)^1.5\n",
	     maxangle, (maxangle>1.)? "> 1" : "<= 1" );
    s_theory = 2. / L * ( (maxangle>1)? 1: pow(maxangle, 1.5) );
    fprintf(f,"= %.6f [e/b^2](T) = %.5f%% u(T) = %.2f MPa,\n",
	    s_theory, s_theory/2*100, s_theory*USTRESS_IN_PA(T)*1E-6 );
    scale = 1 / ((maxangle>1)?1:sqrt(maxangle));
    fprintf(f,"lambda_critical = L/sqrt(K/2e) = %.2f [b](T),\n", scale*L);
    fprintf(f,"---------------------------------------------------------\n\n");
    return;
} /* print_predictions() */


void print_statistics (FILE *f)
{
    int i;
    double cc = o_free+b_free+a_free, dd, ee, ff,
	yadv=yavg-ystart, tadv=t-tstart;
    if ( tadv < TINY ) return;
    fprintf(f,"\n******************** Statistics *********************\n");
    fprintf(f,"At t = %d timesteps = %.2f%% tmax,\n",
	    (int)(t/tdel), t/tmax*100);
    fprintf(f,"the dislocation has advanced by %d [b](T)\n", (int)yavg);
    fprintf(f,"     = %.3f um = %.3f cycles,\n", 
	    yavg*ULENGTH_IN_M(T)*1E6, yavg/ymax, yavg/t/(ss/B) );
    fprintf(f,"----------------------------------------------------\n");
    fprintf(f,"t - tstart = %e [b/vs](T) = %.2f ps,\n",
	    tadv, tadv*UTIME_IN_S(T)*1E12);
    fprintf(f,"speed after tstart = %f [vs](T) = %.3f m/s\n",
	    yadv/tadv, yadv/tadv*USPEED_IN_MS(T));
    fprintf(f,"                   = %.2f%% full speed;\n",
	    yadv/tadv/(ss/B)*100 );
    fprintf(f,"Since tstart, %d obstacles were overcome:\n", (int)cc);
    if ( (int)cc != 0 )
    {
	fprintf(f,"%.2f%% by Orowan loop formation,\n", o_free/cc*100); 
	fprintf(f,"%.2f%% by simple breaking through,\n", b_free/cc*100); 
	fprintf(f,"%.2f%% by thermal activation --\n", a_free/cc*100);
    }
    for (dd=i=0; i<mesh; i++) if (arrested[i]) dd++;
    if (cc+dd>TINY) ee=arrested_time/(cc+dd); else ee=0;
    ff = (cc+dd)/mesh*ee/tadv;
    fprintf(f,"and %d (avg %.2f) mesh points are still arrested;\n",
	    (int)dd, ff*mesh);
    fprintf(f,"----------------------------------------------------\n");
    fprintf(f,"Average jail sentence = %.2f tdel = %.3f ps;\n",
	    ee/tdel, ee*UTIME_IN_S(T)*1E12);
    fprintf(f,"so a mesh-point spends %.3f%% of its time jammed,\n",
	    ff*100);
    fprintf(f,"other %.3f%% of the time free-flying,\n", (1-ff)*100);
    fprintf(f,"and lambda = %.4f L/sqrt(K/2e) = %.2f [b](T);\n",
	    xdel/ff/sqrt(xmax*ymax/nob/maxangle), xdel/ff);
    fprintf(f,"----------------------------------------------------\n");
    fprintf(f,"If we assume mobile dislocation density = %.2e\n",
	    MOBILE_IN_M2);
    fprintf(f,"1/m^2, then predicted strain rate = %.3e 1/s.\n",
	    MOBILE_IN_M2*ULENGTH_IN_M(T)*yadv/tadv*USPEED_IN_MS(T));
    fprintf(f,"*****************************************************\n\n");
    return;
} /* end print_statistics() */


int main (int argc, char *argv[])
{
    int i, j, k;
    double cc, dd, xx[9];
    /* declarations for V */
    volatile double H[2][2];
    VPOINT *dislocation;
    VCIRCLE *obstacles;
    
    for (buf[0]=(char)0,i=0; i<argc; i++)
	strcat(strcat(buf, argv[i]),(i==argc-1)?"":" ");
    printf ("\nProgram \"%s\" starts (pid=%d).\n", buf, getpid());
    watch ("mission", 0);
    
    if ((con=fopen(argv[1],"r"))==NULL)
    {
	printf ("Cannot open input file \"%s\".\n", argv[1]);
	exit(1);
    }
    
    fscanf (con, "PBC width and height [b]:%lf %lf\n",
	    &xmax, &ymax);
    H[0][0]=xmax; H[0][1]=0; H[1][0]=0; H[1][1]=ymax;

    fscanf (con, "\nVolume ratio of obstacles:%lf\n", &c);
    fscanf (con, "\nAverage obstacle size [b]:%lf\n", &a);
    nob = rint(c*xmax*ymax/a/a);
    c = nob*a*a/xmax/ymax;
    
    fscanf (con, "\nK/2e, maximum obstacle angle:%lf\n",
	    &maxangle);
    
    fscanf (con, "\nTemperature [K]:%lf\n", &T);
    B = MY_BOLZ(T)*T;
    
    fscanf (con, "\nExternal stress [e/b^2]:%lf\n", &ss);

    /* one mesh-size is just one obstacle */
    xdel = a;
    mesh = rint(xmax/xdel);
    xdel = xmax/mesh; xdel2 = xdel*xdel;
    
    fscanf (con, "\nr = vs*tdel/xdel:%lf\n", &r);
    r2 = r*r; tdel = r*xdel; tdel2 = tdel*tdel;
    
    fscanf (con, "\nTraversal cycles if in full speed:%lf\n", &cc);
    tmax = cc*ymax/fabs(ss/B);
    
    fscanf (con, "\nRandom number seed (0->time seed):%d\n",
	    &iseed);
    if (iseed==0) iseed = (int)time(0);
    srandom (iseed);
    
    fscanf (con, "\nVisualization:%d\n", &i);
    visualization = (i!=0);
    
    fscanf (con, "\nWaste time:%lf\n", &waste_time);

    fscanf (con, "\nSave results to file (NULL=don't):\n");
    fscanf (con, "%127s %31s %1s", buf, token, replace);
    fclose (con);
    if (!strcasecmp(buf, "NULL")) save = NULL;
    else save = fopen(buf,(replace[0]==(char)0)?"a+":"w+");
    
    print_run_conditions (stdout);
    print_predictions (stdout);
    
    y    = (double *) malloc(mesh*sizeof(double));
    ynew = (double *) malloc(mesh*sizeof(double));
    ylast = (double *) malloc(mesh*sizeof(double));
    dl = (double *) malloc(mesh*sizeof(double));
    sl = (double *) malloc(mesh*sizeof(double));
    d = (double *) malloc(mesh*sizeof(double));
    s = (double *) malloc(mesh*sizeof(double));
    yt = (double *) malloc(mesh*sizeof(double));
    arrested = (Bool *) malloc(mesh*sizeof(Bool));
    xob = (int *) malloc(nob*sizeof(int));
    yob = (double *) malloc(nob*sizeof(double));
    if (visualization)
	ypast = (double *) malloc(mesh*sizeof(double));

    for (i=0; i<nob; i++)
    {
	xob[i] = floor(frandom()*mesh);
	yob[i] = frandom()*ymax;
    }
    
    for (i=0; i<mesh; i++)
    {  /* initial condition */
	y[i] = ymax*(STARTY+STARTW*sin((double)i*2.*PI/mesh));
	ylast[i] = y[i];
	arrested[i] = False;
    }
    
    if (visualization)
    {
	VAddDefaultWindow (H);
	VLock();
	obstacles = (VCIRCLE *)VCreateOG("obstacles",V_CIRCLE,nob);
	for (i=0; i<nob; i++)
	{
	    obstacles[i].s[0] = (double) xob[i]/mesh;
	    obstacles[i].s[1] = yob[i]/ymax;
	    obstacles[i].cid = V_GREEN;
	    obstacles[i].radius = a/2;
	}
	dislocation = (VPOINT *)VCreateOG("dislocation",V_POINT,mesh);
	for (i=0; i<mesh; i++)
	{ 
	    ypast[i] = y[i];
	    dislocation[i].s[0] = (double) i/mesh;
	    dislocation[i].s[1] = y[i]/ymax;
	    dislocation[i].size = 1;
	    dislocation[i].cid = V_WHITE;
	    dislocation[i].lcid = V_BLUE;
	}
	VUnLock();
    }

    /* "abacus" finite-difference scheme for dy/dx ~ O(1) */
    for (t=0; t<tmax; t+=tdel)
    {
	for (i=0; i<mesh; i++)
	{
	    dl[i] = y[i] - pre(y,i);
	    sl[i] = sqrt( dl[i] * dl[i] + xdel2 );
	    d[i] = (aft(y,i) - pre(y,i)) / 2;
	    s[i] = sqrt( d[i] * d[i] + xdel2 );
	    yt[i] = ( y[i] - ylast[i] ) / tdel;
	}
	/* central difference for d^2y/dx^2[i], d^2y/dt^2[i], dy/dt[i] */
	for (i=0; i<mesh; i++)
	{  /* effective stress */
	    cc = ss * s[i] / xdel + s[i] / 4. *
		( aft(d,i) * aft(yt, i) * aft(yt, i) /
		  aft(s,i) / aft(s,i) / aft(s,i) -
		  pre(d,i) * pre(yt, i) * pre(yt, i) /
		  pre(s,i) / pre(s,i) / pre(s,i) );
	    /* effective dissipation */
	    dd = B - d[i] * ( aft(yt,i) - pre(yt,i) ) / 2 / s[i] / s[i];
	    ynew[i] = ( s[i] * ( aft(dl,i) /  aft(sl,i) - 
				 dl[i] / sl[i] ) * r2 +
			tdel2 * cc - ( 1 - dd * tdel / 2 ) * ylast[i] +
			2 * y[i] ) / ( 1 + dd * tdel / 2 );
	}
	/* advance the dislocation */
	for (i=0; i<mesh; i++)
	{
	    ylast[i] = y[i];
	    y[i] = ynew[i];
	}
	/* loop over obstacles */
	for (j=0; j<nob; j++)
	{  /* nearest obstacle image to the dislocation */
	    cc = yob[j];
	    while (cc-ylast[xob[j]] < -ymax/2) cc += ymax;
	    while (cc-ylast[xob[j]] >  ymax/2) cc -= ymax;
	    /* if that obstacle is traversed */
	    if ( copysign (1., ylast[xob[j]] - cc) *
		 copysign (1., y[xob[j]] - cc) == -1. )
	    {
		/* printf ("%f\n", cc); */
		/* flanking force */
		xx[0] = - xdel * DN;
		xx[1] = prn (ylast, xob[j], DN) - ylast[xob[j]];
		xx[2] = sqrt( xx[0] * xx[0] + xx[1] * xx[1] );
		xx[0] /= xx[2]; xx[1] /= xx[2];
		xx[3] = xdel * DN;
		xx[4] = afn (ylast, xob[j], DN) - ylast[xob[j]];
		xx[5] = sqrt( xx[3] * xx[3] + xx[4] * xx[4] );
		xx[3] /= xx[5]; xx[4] /= xx[5];
		xx[6] = xx[0] + xx[3];
		xx[7] = xx[1] + xx[4];
		xx[8] = sqrt( xx[6] * xx[6] + xx[7] * xx[7] );
		/* free by Orowan loop formation, here just by a */
		/* crude estimate of extended angle limit 2*Pi/3 */
		orowan_free = ( xx[8] > sqrt(3.) );
		/* energy barrier */
		excess = a * ( 2 * maxangle - xx[8] );
		break_free = ( !orowan_free ) && ( excess <= 0 );
		/* thermal activation: exponential survival prob. */
		activate_free = ( !orowan_free ) &&
		    ( !break_free ) &&
		    ( frandom() >
		      exp( -tdel / a * exp(-excess/MY_BOLZ(T)/T) ));
		/* event tree */
		if ( orowan_free || break_free || activate_free )
		{
		    if (orowan_free) o_free ++;
		    else if (break_free) b_free ++;
		    else a_free ++;
		    arrested[xob[j]] = False;
		}
		else
		{
		    /* stays slightly behind the obstacle */
		    y[xob[j]] = cc - copysign(TINY, y[xob[j]]-ylast[xob[j]]);
		    /* already arrested */
		    if (arrested[xob[j]])
			arrested_time += tdel;
		    else
		    {  /* obstacle absorb all the impact */
			ylast[xob[j]] = y[xob[j]];
			arrested[xob[j]] = True;
		    }
		}
	    } /* obstacle traversal */
	    else if ( fabs(cc-ylast[xob[j]]) < 2*TINY )
	    {
		/* printf ("%f\n", cc); */
		/* voluntary retreat, very rare */
		arrested[xob[j]] = False;
	    }
	} /* loop over obstacles */
	/* advance since t=0 and deviation from last rendering */
	for (cc=yavg=0,i=0; i<mesh; i++)
	{
	    yavg += y[i];
	    if ( visualization )
		cc += pow(y[i]-ypast[i], POWER);
	}
	yavg = yavg/mesh - STARTY*ymax;
	if ( (t>=STARTT*tmax) && (!final_average) )
	{ /* begin to collect final averages */
	    tstart = t;
	    ystart = yavg;
	    o_free = b_free = a_free = 0;
	    arrested_time = 0;
	    final_average = True;
	}
	if ( visualization && (pow(cc/mesh, 1./POWER) > 0.1*a) )
	{  /* render the motion */
	    VLock();
	    for (i=0; i<mesh; i++)
	    {
		ypast[i] = y[i];
		dislocation[i].s[1] = y[i]/ymax;
		if ( arrested[i] )
		{
		    dislocation[i].cid = V_RED;
		    dislocation[i].size = 3;
		}
		else
		{
		    dislocation[i].cid = V_WHITE;
		    dislocation[i].size = 1;
		}
	    }
	    VUnLock();
	}
	/* notify the user about progress */
	if ( floor(t/tmax/CHUNK_OF_WORK) > last_print )
	{
	    last_print++;
	    if ( !final_average )
		printf("%.2f%% of the simulation done.\n",
		       last_print*CHUNK_OF_WORK*100);
	    else
		print_statistics (stdout);
	}
	/* waste some time here for the eye */
	if ( visualization )
	    for (i=0; i<rint(waste_time*1000); i++,cc+=i);
    } /* integration loop */
    
    if (visualization)
    {
	VCloseAllWindows();
	VDeleteOG("obstacles");
	VDeleteOG("dislocation");
	/* let the threads die out */
	waitfor(0.1); 
    }

    printf("mission terminates successfully.\n");
    printf("time elapsed = %s.\n",
	   earth_time(watch("mission",3)));
    
    /* final results */
    printf ("\n *** Summary ***\n");
    print_run_conditions (stdout);
    print_predictions (stdout);
    print_statistics (stdout);
    
    if ( save != NULL )
    {
	fprintf (save, "\n");
	if ( token[0] != (char)0 )
	{  /* print separator */
	    k = 60;
	    for (i=0; i<k; i++) buf[i] = '%';
	    fprintf (save, "%s\n", buf);
	    for (i=2; i<k-2; i++) buf[i] = ' ';
	    p = buf+(k-strlen(token))/2;
	    strcpy (p, token); p[strlen(token)]=' ';
	    fprintf (save, "%s\n", buf);
	    for (i=0; i<k; i++) buf[i] = '%';
	    fprintf (save, "%s\n", buf);
	}
	print_run_conditions (save);
	print_predictions (save);
	print_statistics (save);
	fprintf (save, "\n");
	fclose(save);
    }
    return(0);
} /* end main() */