/*************************************/
/* Utility subroutines for Argon 2.0 */
/*************************************/

#include "utils.h"

/* return pseudo random number on (0,1) */
double frandom()
{
    return ((random()+0.5)/pow(2.,31.));
} /* end frandom() */


/* generate n normal variables and store in x[] */
void randnorm (int n, double x[])
{
    int i;
    double s1, s2;
    for (i=0; i<n; i+=2)
    {
	s1 = sqrt(-2.*log(frandom()));
	s2 = 2*acos(-1.)*frandom();
	x[i] = s1 * cos(s2);
	if (i+1 < n) x[i+1] = s1 * sin(s2);
    }
    return;
} /* end randnorm() */


/* see if a file exists, and if so ask */
/* the user whether to overwrite it.   */
int freetowrite (char filename[])
{
    FILE *file;
    char c;
    file = fopen (filename, "r");
    if (file != NULL) 
    {
	fclose (file);
	printf ("file \"%s\" exists, overwrite (y/n)? ",
		filename);
        c = getc(stdin);
	if (c!='\n') while(getc(stdin)!='\n');
        if (!((c=='y')||(c=='Y'))) return(0);
    }
    return(1);
} /* end freetowrite() */


void waitfor (double seconds)
{
    struct timeval time;
    time.tv_sec = floor(seconds);
    time.tv_usec = rint((seconds-time.tv_sec)*1E6);
    select (0, NULL, NULL, NULL, &time);
    return;
} /* end waitfor() */


/*****************************************************/
/* Time processes (real time, in seconds) by 16-byte */
/* tokens. 0: init, 1: check, 2: refresh, 3:destroy. */ 
/*****************************************************/
#define TIME (tp.tv_sec+0.000001*tp.tv_usec)
#define MAX_RECORD 512
#define MAX_TOKEN_SIZE 16
double watch (char token[], int command)
{
    struct list
    {
	char token[MAX_TOKEN_SIZE];
	double init_time;
	struct list *next;
    };
    static struct list mempool[MAX_RECORD]
	= {{{(char)0},0.,(struct list *)NULL}};
    static struct list *p, *q;
    static struct timeval tp;
    double cc;
    
    gettimeofday (&tp, NULL);
    switch (command)
    {
	/* start a new watch */
	case 0:  /* root node is mempool[0] */
	    for ( p = mempool;
		  (p->next != NULL) &&
		      strncmp(p->next->token,
			      token, MAX_TOKEN_SIZE-1);
		  p = p->next );
	if ( p->next == NULL )
	{
	    /* look for free space */
	    for (q=mempool+1; q<mempool+MAX_RECORD; q++)
		if ( *(q->token)==(char)0 )
		{
		    p->next = q;
		    q->next = NULL;
		    strncpy (q->token, token, MAX_TOKEN_SIZE-1);
		    q->init_time = TIME;
		    return(TIME);
		}
	    printf("**error watch(): mempool %d records all used up.\n",
		   MAX_RECORD-1);
	    return(-1.);
	}
	else
	{
	    printf("**error watch(): token name \"%s\" already used.\n",
		   token);
	    return(-1.);
	}
	
	/* return the time elapsed in seconds */
	case 1:
	    for ( p = mempool;
		  (p->next != NULL) &&
		      strncmp(p->next->token,
			      token, MAX_TOKEN_SIZE-1);
		  p = p->next);
	if ( p->next == NULL )
	{
	    printf("**error watch(): token name \"%s\" not found.\n",
		   token);
	    return(-1.);
	}
	else
	    return (TIME - p->next->init_time);

	/* return the time elapsed in seconds */
	/* and refresh the watch's startpoint */
	case 2:
	    for ( p = mempool;
		  (p->next != NULL) &&
		      strncmp(p->next->token,
			      token, MAX_TOKEN_SIZE-1);
		  p = p->next);
	if ( p->next == NULL )
	{
	    printf("**error watch(): token name \"%s\" not found.\n",
		   token);
	    return(-1.);
	}
	else
	{
	    cc = p->next->init_time;
	    p->next->init_time = TIME;
	    return (TIME - cc);
	}

	/* destroy the record and free resources */
	case 3:
	    for ( p = mempool;
		  (p->next != NULL) &&
		      strncmp(p->next->token,
			      token, MAX_TOKEN_SIZE-1);
		  p = p->next);
	if ( p->next == NULL )
	{
	    printf("**error watch(): token name \"%s\" not found.\n",
		   token);
	    return(-1.);
	}
	else
	{
	    cc = p->next->init_time;
	    *(p->next->token) = (char)0;
	    p->next = p->next->next;
	    return (TIME - cc);
	}
	
    default:
	printf("**error watch(): 0:init 1:check 2:refresh 3:destroy.\n");
	return(-1.);
    }
} /* end watch() */
#undef TIME
#undef MAX_RECORD 
#undef MAX_TOKEN_SIZE 


char *earth_time (double seconds)
{
    const struct time_unit
    {
	char *name;
	double scale;
    } earthman[] =
      { {"s",  1},
	{"min", 60},
	{"hr", 60},
	{"day", 24},
	{"month", 30},
	{"year", 365} };
    static char time[128], *p;
    int i, max=sizeof(earthman)/sizeof(struct time_unit);
    double scale; int denom;
    seconds /= earthman[0].scale;
    for (scale=1,i=1; i<max; i++) scale*=earthman[i].scale;
    for (p=time,i=max-1; i>0; i--)
    {
	denom = floor(seconds/scale);
	if (denom>0)
	{
	    sprintf(p, "%d %s, ", denom, earthman[i].name);
	    while((*p)!=(char)0) p++;
	}
	seconds -= denom*scale;
	scale /= earthman[i].scale;
    }
    sprintf(p, "%.3f %s", seconds, earthman[0].name);
    return(time);
} /* earth_time() */


void matran(double A[3][3], double B[3][3])
{
    B[0][0] = A[0][0];
    B[0][1] = A[1][0];
    B[0][2] = A[2][0];
    B[1][0] = A[0][1];
    B[1][1] = A[1][1];
    B[1][2] = A[2][1];
    B[2][0] = A[0][2];
    B[2][1] = A[1][2];
    B[2][2] = A[2][2];
    return;
} /* end matran() */


void matadd(double A[3][3], double B[3][3], double C[3][3])
{
    C[0][0] = A[0][0] + B[0][0];
    C[0][1] = A[0][1] + B[0][1];
    C[0][2] = A[0][2] + B[0][2];
    C[1][0] = A[1][0] + B[1][0];
    C[1][1] = A[1][1] + B[1][1];       
    C[1][2] = A[1][2] + B[1][2];       
    C[2][0] = A[2][0] + B[2][0];           
    C[2][1] = A[2][1] + B[2][1];      
    C[2][2] = A[2][2] + B[2][2];
    return;
} /* end matadd() */


void matsub(double A[3][3], double B[3][3], double C[3][3])
{
    C[0][0] = A[0][0] - B[0][0];
    C[0][1] = A[0][1] - B[0][1];
    C[0][2] = A[0][2] - B[0][2]; 
    C[1][0] = A[1][0] - B[1][0];         
    C[1][1] = A[1][1] - B[1][1];       
    C[1][2] = A[1][2] - B[1][2];       
    C[2][0] = A[2][0] - B[2][0];           
    C[2][1] = A[2][1] - B[2][1];      
    C[2][2] = A[2][2] - B[2][2];
    return;
} /* end matsub() */


void matmul(double A[3][3], double B[3][3], double C[3][3])
{
    C[0][0] = A[0][0]*B[0][0] + A[0][1]*B[1][0] + A[0][2]*B[2][0];
    C[0][1] = A[0][0]*B[0][1] + A[0][1]*B[1][1] + A[0][2]*B[2][1];
    C[0][2] = A[0][0]*B[0][2] + A[0][1]*B[1][2] + A[0][2]*B[2][2];
    C[1][0] = A[1][0]*B[0][0] + A[1][1]*B[1][0] + A[1][2]*B[2][0];
    C[1][1] = A[1][0]*B[0][1] + A[1][1]*B[1][1] + A[1][2]*B[2][1];
    C[1][2] = A[1][0]*B[0][2] + A[1][1]*B[1][2] + A[1][2]*B[2][2];
    C[2][0] = A[2][0]*B[0][0] + A[2][1]*B[1][0] + A[2][2]*B[2][0];
    C[2][1] = A[2][0]*B[0][1] + A[2][1]*B[1][1] + A[2][2]*B[2][1];
    C[2][2] = A[2][0]*B[0][2] + A[2][1]*B[1][2] + A[2][2]*B[2][2];
    return;
} /* end matmul() */


double matinv(double A[3][3], double B[3][3])
{
    double D11,D22,D33,D12,D21,D13,D31,D32,D23,determinant;
    D11=A[1][1]*A[2][2]-A[1][2]*A[2][1];
    D22=A[2][2]*A[0][0]-A[2][0]*A[0][2];
    D33=A[0][0]*A[1][1]-A[0][1]*A[1][0];
    D12=A[1][2]*A[2][0]-A[1][0]*A[2][2];
    D23=A[2][0]*A[0][1]-A[2][1]*A[0][0];
    D31=A[0][1]*A[1][2]-A[0][2]*A[1][1];
    D13=A[1][0]*A[2][1]-A[2][0]*A[1][1];
    D21=A[2][1]*A[0][2]-A[0][1]*A[2][2];
    D32=A[0][2]*A[1][0]-A[1][2]*A[0][0];
    determinant=A[0][0]*D11+A[0][1]*D12+A[0][2]*D13;
    B[0][0]=D11/determinant;
    B[1][1]=D22/determinant;
    B[2][2]=D33/determinant;
    B[0][1]=D21/determinant;
    B[1][2]=D32/determinant;
    B[2][0]=D13/determinant;
    B[1][0]=D12/determinant;
    B[2][1]=D23/determinant;
    B[0][2]=D31/determinant;
    return (determinant);
} /* end matinv() */


/******************************************/
/* Diagonalize 3x3 real symmetric matrix: */
/* A = Q^T * diag(eigval) * Q, eigval[i]  */
/* will be stored in descending order,    */
/* Q[i][] forming a right-handed system.  */
/******************************************/
#define DIAG3_TINY (1E-12)
char *diag3 (double A[3][3], double eigval[3], double Q[3][3])
{
    static char info[32];
    int i;
    double a, b, c, d, e, module, angle;

    info[0] = (char) 0;
    if ( (fabs(A[0][1]-A[1][0]) > DIAG3_TINY) ||
	 (fabs(A[0][2]-A[2][0]) > DIAG3_TINY) ||
	 (fabs(A[1][2]-A[2][1]) > DIAG3_TINY) )
    {
	strcpy (info, "matrix is not symmetric.");
	return (info);
    }
    else
    {
	A[0][1] = (A[0][1]+A[1][0])/2.;
	A[1][0] =  A[0][1];
	A[0][2] = (A[0][2]+A[2][0])/2.;
	A[2][0] =  A[0][2];
	A[1][2] = (A[1][2]+A[2][1])/2.;
	A[2][1] =  A[1][2];
    }
    /* solve x^3 + a*x^2 + b*x + c = 0: */
    a = -A[0][0]-A[1][1]-A[2][2];
    b = A[0][0]*A[1][1]+A[0][0]*A[2][2]+A[1][1]*A[2][2]-
	A[0][1]*A[0][1]-A[0][2]*A[0][2]-A[1][2]*A[1][2];
    c = A[0][0]*A[1][2]*A[1][2]+A[1][1]*A[0][2]*A[0][2]+
	A[2][2]*A[0][1]*A[0][1]-A[0][0]*A[1][1]*A[2][2];
    d = 36*a*b-108*c-8*a*a*a;
    e = 12*b*b*b-3*b*b*a*a-54*a*b*c+81*c*c+12*a*a*a*c;
    if (e > DIAG3_TINY)
    {
	strcpy (info, "imaginary eigenvalues.");
	return (info);
    }
    else
	e = 12*sqrt(fabs(e));
    module = sqrt(d*d+e*e);
    angle = (fabs(module)<DIAG3_TINY)?0:asin(e/module);
    if (d < 0) angle += acos(-1.);
    d = pow(module,1./3.)*cos(angle/3)/6;
    e = pow(module,1./3.)*fabs(sin(angle/3))/6;
    eigval[0] = 2*d-a/3;
    eigval[1] = -d-a/3+sqrt(3.)*e;
    eigval[2] = -d-a/3-sqrt(3.)*e;
    /* make eigval[0] >= eigval[1] >= eigval[2] */
    if (eigval[0]<eigval[1])
    {
	eigval[0] = eigval[1];
	eigval[1] = 2*d-a/3;
	if (eigval[1] < eigval[2])
	{
	    eigval[1] = eigval[2];
	    eigval[2] = 2*d-a/3;
	}
    }
    /* calculate the eigenvectors */
    Q[0][0] = Q[0][1] = Q[0][2] = 0.;
    Q[1][0] = Q[1][1] = Q[1][2] = 0.;
    Q[2][0] = Q[2][1] = Q[2][2] = 0.;
    if ((eigval[0]==eigval[1])&&(eigval[0]==eigval[2]))
    {
	Q[0][0] = Q[1][1] = Q[2][2] = 1.;
	return(info);
    }
    /* first eigenvector */
    if ((A[0][0]==eigval[0])&&(A[1][0]==0.)&&(A[2][0]==0.))
	Q[0][0] = 1.;
    else if ((A[0][1]==0.)&&(A[1][1]==eigval[0])&&(A[2][1]==0.))
	Q[0][1] = 1.;
    else if ((A[0][2]==0.)&&(A[1][2]==0.)&&(A[2][2]==eigval[0]))
	Q[0][2] = 1.;
    else
    {
	Q[0][0] = 1.;
	a = (A[1][1]-eigval[0])*(A[2][2]-eigval[0])
	    -A[1][2]*A[1][2];
	Q[0][1] = (-(A[2][2]-eigval[0])*A[0][1]
		   +A[1][2]*A[0][2])/a;
	Q[0][2] = (A[1][2]*A[0][1]-
		   (A[1][1]-eigval[0])*A[0][2])/a;
	module = sqrt( Q[0][0]*Q[0][0] +
		       Q[0][1]*Q[0][1] +
		       Q[0][2]*Q[0][2] );
	Q[0][0] = Q[0][0] / module;
	Q[0][1] = Q[0][1] / module;
	Q[0][2] = Q[0][2] / module;
    }
    /* third eigenvector */
    if ((A[0][2]==0.)&&(A[1][2]==0.)&&(A[2][2]==eigval[2]))
	Q[2][2] = 1.;    
    else if ((A[0][1]==0.)&&(A[1][1]==eigval[2])&&(A[2][1]==0.))
	Q[2][1] = 1.;
    else if ((A[0][0]==eigval[2])&&(A[1][0]==0.)&&(A[2][0]==0.))
	Q[2][0] = 1.;
    else
    {
	Q[2][0] = 1.;
	a = (A[1][1]-eigval[2])*(A[2][2]-eigval[2])
	    -A[1][2]*A[1][2];
	Q[2][1] = (-(A[2][2]-eigval[2])*A[0][1]
		   +A[1][2]*A[0][2])/a;
	Q[2][2] = (A[1][2]*A[0][1]-
		   (A[1][1]-eigval[2])*A[0][2])/a;
	module = sqrt( Q[2][0]*Q[2][0] +
		       Q[2][1]*Q[2][1] +
		       Q[2][2]*Q[2][2] );
	Q[2][0] = Q[2][0] / module;
	Q[2][1] = Q[2][1] / module;
	Q[2][2] = Q[2][2] / module;
    }
    /* second eigenvector = first x third */
    Q[1][0] = -Q[0][1]*Q[2][2] + Q[2][1]*Q[0][2];
    Q[1][1] = Q[0][0]*Q[2][2] - Q[2][0]*Q[0][2];
    Q[1][2] = -Q[0][0]*Q[2][1] + Q[2][0]*Q[0][1];
    return(info);
} /* end diag3() */


#ifdef DIAG3_TEST
void main()
{
    /* double A[3][3] = {{0.4876,    1.1253,    -0.6633}, */
                      /* {1.1253,    1.0372,    1.2965}, */
                      /* {-0.6633,    1.2965,    1.3789}}; */
    double A[3][3] = {{1.,    0.,    0.},
                      {0.,    1.,    0.},
                      {0.,    0.,    1.}};
    double eigval[3],Q[3][3],QT[3][3],B[3][3],L[3][3],QI[3][3],det;
    /* function driver */
    if (*diag3(A,eigval,Q))
    {
        printf ("Error (diag3): %s\n", diag3(A,eigval,Q));
        exit (1);
    }
    /* verification */
    printf ("\n               %f %f %f",A[0][0],A[0][1],A[0][2]);
    printf ("\n    A       =  %f %f %f",A[1][0],A[1][1],A[1][2]);
    printf ("\n               %f %f %f\n",A[2][0],A[2][1],A[2][2]);
    printf ("\neig[0] = %f, vec[0] = %f %f %f\n",
            eigval[0], Q[0][0], Q[0][1], Q[0][2]);
    printf ("eig[1] = %f, vec[1] = %f %f %f\n",
            eigval[1], Q[1][0], Q[1][1], Q[1][2]);
    printf ("eig[2] = %f, vec[2] = %f %f %f\n",
            eigval[2], Q[2][0], Q[2][1], Q[2][2]);
    printf ("\n               %f %f %f",Q[0][0],Q[0][1],Q[0][2]);
    printf ("\n    Q       =  %f %f %f",Q[1][0],Q[1][1],Q[1][2]);
    printf ("\n               %f %f %f\n",Q[2][0],Q[2][1],Q[2][2]);
    matran (Q,QT);
    matmul (Q,A,B);
    matmul (B,QT,L);
    printf ("\n               %f %f %f",L[0][0],L[0][1],L[0][2]);
    printf ("\nQ * A * Q^T =  %f %f %f",L[1][0],L[1][1],L[1][2]);
    printf ("\n               %f %f %f\n",L[2][0],L[2][1],L[2][2]);
    matmul (Q,QT,L);
    printf ("\n               %f %f %f",L[0][0],L[0][1],L[0][2]);
    printf ("\n  Q * Q^T   =  %f %f %f",L[1][0],L[1][1],L[1][2]);
    printf ("\n               %f %f %f\n",L[2][0],L[2][1],L[2][2]);
    det = matinv (Q,QI);
    printf ("\ndet|Q|   =  %f\n\n",det);    
    return;
} /* end main() */
#endif


/* decompose A to B+C, where B = trace/3 */
double trace_decompose (double A[3][3], double B[3][3], double C[3][3])
{
    double trace;
    trace = (A[0][0] + A[1][1] + A[2][2]) / 3.;
    if (B != NULL)
    {
	B[0][0] = B[1][1] = B[2][2] = trace;
	B[0][1] = B[0][2] = 0.;
	B[1][0] = B[1][2] = 0.;
	B[2][0] = B[2][1] = 0.;
    }
    if (C != NULL)
    {
	C[0][0] = A[0][0] - trace;
	C[1][1] = A[1][1] - trace;
	C[2][2] = A[2][2] - trace;
	C[0][1] = C[0][2] = 0.;
	C[1][0] = C[1][2] = 0.;
	C[2][0] = C[2][1] = 0.;
    }
    return (trace);
} /* end trace_decompose() */


/************************************/
/* Quicksort: arr[] will be altered */
/************************************/
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define SWAPi(a,b) tempi=(a);(a)=(b);(b)=tempi;
#define M 7
#define NSTACK 64
void quicksort (int n, double arr[], int idx[])
{
    int i, ir=n, j, k, l=1;
    int istack[NSTACK], jstack=0, b, tempi;
    double a, temp;
    /* added by Ju Li, Fortran convention */
    arr--; idx--;
    for (i=1; i<=n; i++) idx[i] = i-1;
    /* added by Ju Li, Fortran convention */
    for (;;) {
	if (ir-l < M) {
	    for (j=l+1;j<=ir;j++) {
		a=arr[j];
		b=idx[j];
		for (i=j-1;i>=1;i--) {
		    if (arr[i] <= a) break;
		    arr[i+1]=arr[i];
		    idx[i+1]=idx[i];
		}
		arr[i+1]=a;
		idx[i+1]=b;
	    }
	    if (!jstack) return;
	    ir=istack[jstack];
	    l=istack[jstack-1];
	    jstack -= 2;
	} else {
	    k=(l+ir) >> 1;
	    SWAP(arr[k],arr[l+1])
		SWAPi(idx[k],idx[l+1])
		if (arr[l+1] > arr[ir]) {
		    SWAP(arr[l+1],arr[ir])
			SWAPi(idx[l+1],idx[ir])
			}
	    if (arr[l] > arr[ir]) {
		SWAP(arr[l],arr[ir])
		    SWAPi(idx[l],idx[ir])
		    }
	    if (arr[l+1] > arr[l]) {
		SWAP(arr[l+1],arr[l])
		    SWAPi(idx[l+1],idx[l])
		    }
	    i=l+1;
	    j=ir;
	    a=arr[l];
	    b=idx[l];
	    for (;;) {
		do i++; while (arr[i] < a);
		do j--; while (arr[j] > a);
		if (j < i) break;
		SWAP(arr[i],arr[j])
		    SWAPi(idx[i],idx[j])
		    }
	    arr[l]=arr[j];
	    arr[j]=a;
	    idx[l]=idx[j];
	    idx[j]=b;
	    jstack += 2;
	    if (jstack > NSTACK)
	    {
		printf ("NSTACK too small in sort2.");
		exit(1);
	    }
	    if (ir-i+1 >= j-l) {
		istack[jstack]=ir;
		istack[jstack-1]=i;
		ir=j-1;
	    }
	    else
	    {
		istack[jstack]=j-1;
		istack[jstack-1]=l;
		l=i;
	    }
	}
    }
}
#undef M
#undef NSTACK
#undef SWAP
#undef SWAPi


#ifdef QUICKSORT_TEST
#define K 6
void main()
{
    double a[K] = {0.1, -0.2, 3.6, 1.8, 9.9, -21.}, c[K];
    int idx[K], i;
    for (i=0; i<K; i++) c[i] = a[i];
    quicksort (K, a, idx);
    for (i=0; i<K; i++)
        printf ("%f %f\n", a[i], c[idx[i]]);
}
#undef K
#endif


/*************************************/
/* Piecewise Linear Schedule Manager */
/*************************************/
static struct schedule_list schpool[MAX_SCHEDULES]={0};

/* schedule format: .., [Ttime]/[Ffraction], value, .. */
struct schedule_list *add_schedule
(char token[], char schedule[], double tmin, double tmax)
{
    struct schedule_list *p, *q;
    char *c, *d, buf[MAX_SCHEDULE_INPUT];
    int i, j, k, l, pair_start, o[MAX_SCHEDULE_N];
    double tmp[2][MAX_SCHEDULE_N];
    /* copy the input to a buffer */
    strncpy (buf, schedule, MAX_SCHEDULE_INPUT-1);
    buf[MAX_SCHEDULE_INPUT-1] = '\0';
    /* root node is schpool[0] */
    for ( p = schpool;
	  (p->next != NULL) &&
	      strncmp(p->next->token,
		      token, MAX_SCHEDULE_TOKEN_SIZE-1);
	  p = p->next );
    if ( p->next != NULL )
    {
	printf("add_schedule(): token name \"%s\" already used.\n",
	       token);
	exit(1);
    }
    /* look for free pointer space in schpool[] */
    for (q=schpool+1; q<schpool+MAX_SCHEDULES; q++)
	if (*(q->token) == '\0')
	{
	    p->next = q;
	    q->next = NULL;
	    strncpy (q->token, token, MAX_SCHEDULE_TOKEN_SIZE-1);
	    q->token[MAX_SCHEDULE_TOKEN_SIZE-1] = '\0';
	    j = 0; /* total number of control points */
	    pair_start = 1; /* at the start of a control record */
	    for (c=buf; (*c=='\t')||(*c==' '); c++);
	    for (; (j<MAX_SCHEDULE_N) && (*c!='\0'); c=d)
	    {
		for (d=c; (*d!='\t')&&(*d!=' ')&&(*d!='\0'); d++);
		if (*d != '\0')
		{
		    while ((*d=='\t')||(*d==' ')) d++;
		    *(d-1) = '\0';
		}
		if ( ((*c>='0')&&(*c<='9')) || (*c=='-') || (*c=='+') )
		{
		    if ( pair_start )  /* no control time is given */
			tmp[0][j] = SCHEDULE_NONSENSE; 
		    tmp[1][j] = atof(c); /* control values */
		    pair_start = 1; /* end of this pair */
		    j++;
		}
		else if ( ((*c=='T') || (*c=='t')) && pair_start )
		{ /* explicit control time */
		    tmp[0][j] = atof(c+1);
		    pair_start = 0;  /* wait for value input */
		}
		else if ( ((*c=='F') || (*c=='f')) && pair_start )
		{ /* fractional control time */
		    tmp[0][j] = tmin + (tmax-tmin)*atof(c+1);
		    pair_start = 0;
		}
		else
		{
		    printf ("add_schedule: cannot parse substring "
			    "\"%s\".\n", c);
		    exit(1);
		}
	    }
	    if (j >= MAX_SCHEDULE_N)
	    {
		printf("add_schedule: out of temporary storage"
		       " MAX_SCHEDULE_N = %d\n"
		       "for \"%s\".\n", MAX_SCHEDULE_N, buf);
		exit(1);
	    }
	    q->n = j;
	    q->s = malloc(2*j*sizeof(double));
	    if (tmp[0][0] == SCHEDULE_NONSENSE) tmp[0][0] = tmin;
	    if (tmp[0][j-1] == SCHEDULE_NONSENSE) tmp[0][j-1] = tmax;
	    for (i=1; i<j-1; i++)
		if (tmp[0][i] == SCHEDULE_NONSENSE)
		{ /* equal interpolation between sensible controls */
		    for (k=i-1; tmp[0][k]==SCHEDULE_NONSENSE; k--);
		    for (l=i+1; tmp[0][l]==SCHEDULE_NONSENSE; l++);
		    tmp[0][i] = tmp[0][k] + (tmp[0][l]-tmp[0][k]) *
			(i-k) / (l-k);
		}
	    /* sort control times in ascending order */
	    quicksort (j, &tmp[0][0], o);
	    for (q->nontrivial=i=0; i<j; i++)
	    { /* time already sorted */
		q->s[2*i] = tmp[0][i];  
		q->s[2*i+1] = tmp[1][o[i]];
		q->nontrivial = q->nontrivial ||
		    (fabs(q->s[2*i+1]-q->s[1]) > SCHEDULE_TINY);
	    }
	    return (q);
	}
    printf("add_schedule(): schpool %d pointer slots all used up.\n",
	   MAX_SCHEDULES-1);
    exit(1);
    return(NULL);
} /* return add_schedule() */


double schedule (char token[], double t)
{
    int i;
    struct schedule_list *q;
    for ( q = schpool;
	  (q->next != NULL) &&
	      strncmp(q->next->token,
		      token, MAX_SCHEDULE_TOKEN_SIZE-1);
	  q = q->next );
    if (q->next == NULL)
    {
	printf("schedule: token name \"%s\" not found.\n",
	       token);
	exit (1);
    }
    q = q->next;
    if ( !(q->nontrivial) ) return (q->s[1]);
    if (t <= q->s[0]) return (q->s[1]);
    if (t >= q->s[2*q->n-2]) return (q->s[2*q->n-1]);
    for (i=0; i<q->n-1; i++)
	if ( (t >= q->s[2*i]) && (t < q->s[2*i+2]) ) break;
    if ( fabs(q->s[2*i] - q->s[2*i+2]) < SCHEDULE_TINY )
	return( (q->s[2*i+1] + q->s[2*i+3]) / 2 );
    else return (q->s[2*i+1] + (q->s[2*i+3] - q->s[2*i+1]) *
		 (t - q->s[2*i]) / (q->s[2*i+2] - q->s[2*i]));
}  /* end schedule() */


/* See if the control value remains constant during period */
/* [t1,t2]: if it does, then return the upper and lower    */
/* bounds which contain [t1,t2]; otherwise return NULL.    */
double *calm_period (char token[], double t1, double t2)
{
    int i, j, k;
    struct schedule_list *q;
    double ti, tf;
    static double bound[2];
    for ( q = schpool;
	  (q->next != NULL) &&
	      strncmp(q->next->token,
		      token, MAX_SCHEDULE_TOKEN_SIZE-1);
	  q = q->next );
    if ( q->next == NULL )
    {
	printf("schedule: token name \"%s\" not found.\n", token);
	exit (1);
    }
    q = q->next;
    ti = (t1 < t2)? t1:t2;
    tf = (t1 < t2)? t2:t1;
    for (i=0; i<q->n; i++)
	if (fabs(ti - q->s[2*i]) < SCHEDULE_TINY) break;
	else if (ti < q->s[2*i])
	    if (i==0) break;
	    else {i--; break;}
    if (ti < q->s[2*0]) bound[0] = ti;
    else /* search left for lower bounds */
	for (k=i; k>=0; k--)
	    if (fabs(q->s[2*k+1]-q->s[2*i+1])<SCHEDULE_TINY)
		bound[0] = q->s[2*k];
    for (j=q->n-1; j>=0; j--)
	if (fabs(tf - q->s[2*j]) < SCHEDULE_TINY) break;
	else if (tf > q->s[2*j])
	    if (j==q->n-1) break;
	    else {j++; break;}
    if (tf > q->s[2*(q->n-1)]) bound[1] = tf;
    else /* search right for upper bounds */
	for (k=j; k<q->n; k++)
	    if (fabs(q->s[2*k+1]-q->s[2*j+1])<SCHEDULE_TINY)
		bound[1] = q->s[2*k];
    for (; i<j; i++)
	if ((fabs(q->s[2*i+1]-q->s[2*j+1]) > SCHEDULE_TINY))
	    return(NULL);
    return(bound);
} /* end calm_period() */


void cancel_schedule (char token[])
{
    struct schedule_list *q;
    for ( q = schpool;
	  (q->next != NULL) &&
	      strncmp(q->next->token,
		      token, MAX_SCHEDULE_TOKEN_SIZE-1);
	  q = q->next );
    if ( q->next == NULL )
    {
	printf("schedule: token name \"%s\" not found.\n", token);
	exit (1);
    }
    free (q->next->s);
    q->next->n = 0;
    *(q->next->token) = '\0';
    q->next = q->next->next;
    return;
} /* end cancel_schedule() */


void cancel_all_schedules()
{
    struct schedule_list *q;
    for ( q=schpool; q->next != NULL; )
    {
	free (q->next->s);
	q->next->n = 0;
	*(q->next->token) = '\0';
	q->next = q->next->next;
    }
    return;
} /* end cancel_all_schedules() */


#ifdef SCHEDULE_TEST
void main()
{
    struct schedule_list *a =
	add_schedule ("a", "  F0.1 10 20 t0.9 30  ", 0, 1);
    struct schedule_list *b =
	add_schedule ("b", "  F0.1 10 F0.11 20 t0.9 15  ", 0, 1);
    struct schedule_list *c;
    printf ("%d %f %f %f %f %f %f\n", a->n, a->s[0], a->s[1],
	    a->s[2], a->s[3], a->s[4], a->s[5]);
    printf ("%d %f %f %f %f %f %f\n", b->n, b->s[0], b->s[1],
	    b->s[2], b->s[3], b->s[4], b->s[5]);
    printf ("%f %f %f\n", schedule("a", -0.1),
	    schedule("a", 0.2),
	    schedule("a", 0.7));
    cancel_schedule ("a");
    printf ("%f %f %f %f\n", schedule("b", -0.1),
	    schedule("b", 0.2),
	    schedule("b", 1.7),
	    schedule("b", 0.103)
	    );
    cancel_all_schedules();
    c = add_schedule
	("c", "  F0.1 10   10 10 10 10   F0.9 10 F0.95 20 ",
	 0, 1);
    printf ("c = %d %d %d\n", calm_period("c", -10, 1.5), 
	    calm_period("c", 0.4, 0.5),
	    calm_period("c", 0.4, 0.9));
}
#endif


/* Without changing idx[2*min] and idx[2*max], solve */
/* the memory conflict at idx[2*i+1], idx[2*i+2].    */
void make_space (int idx[], int list[], int i, int min, int max)
{
    int j, k, a, left_shift;
    if (idx[2*i+1] < idx[2*i+2]) return;
    else if (idx[2*i+1] > idx[2*i+2])
    {
	printf ("make_space: detects overflow at i=%d and i+1.\n",
		i);
	exit(1);
    }
    /* look for left free space */
    for (j=i; (j>min)&&(idx[2*j]==idx[2*j-1]); j--);
    /* look for right free space */
    for (k=i+1; (k<max)&&(idx[2*k+1]==idx[2*k+2]); k++);
    if ( (j>min) && (k<max) )
	left_shift =  /* to prevent systematic drift */
	    (idx[2*i+1]-idx[2*j] <  idx[2*k+1]-idx[2*i+2]) ||
	  ( (idx[2*i+1]-idx[2*j] == idx[2*k+1]-idx[2*i+2]) &&
	    (frandom() < 0.5) );
    else if ( (j>min) && (k>=max) ) left_shift = 1;
    else if ( (j<=min) && (k<max) ) left_shift = 0;
    else
    {
	printf ("make_space: min=%d max=%d jammed between "
		"i=%d and i+1.\n", min, max, i);
	exit(1);
    }
    /* minimal action: shift only by 1 unit */
    if ( left_shift )
    {
	for (a=idx[2*j]; a<idx[2*i+1]; a++)
	    list[a-1] = list[a];
	for (a=2*j; a<=2*i+1; a++) idx[a]--;
    }
    else
    {
	for (a=idx[2*k+1]; a>idx[2*i+2]; a--)
	    list[a] = list[a-1];
	for (a=2*k+1; a>=2*i+2; a--) idx[a]++;
    }
    return;
} /* end make_space() */


void safe_append
(int idx[], int list[], int value, int i, int min, int max)
{
    make_space (idx, list, i, min, max);
    list[idx[2*i+1]++] = value;
    return;
} /* end safe_append() */