#include "utils.h"

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() */


/* 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() */


/*****************************************************/
/* 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.);
    }
}
#undef TIME
#undef MAX_RECORD 512
#undef MAX_TOKEN_SIZE 16
/* watch() */

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() */


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() */


/************************************/
/* 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 test_quicksort
#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\n", c[idx[i]]);
}
#undef K
#endif