/* 

streetcar.c -- a program to simulate a disordered streetcar route

by Jeffrey S. Rosenthal, 2007

*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <sys/types.h>

#define DATACHOICE 0
#define LEADTIME 24.0

#define NUMCARSPER (28)  /* per unit route length */
#define ROUTELENGTH 3
#define NUMCARS (NUMCARSPER * ROUTELENGTH)
#define NUMREPS 10
#define NUMSIMS 1000

#define SIGMA 0.01
#define DELTA 0.0005
#define DIAMETER (0.05/NUMCARSPER)
#define EPSILON (1.0/60/2)
/* #define EPSILON DIAMETER */
#define MAXTIMES 1000

double dsfactor;

#define DISPLAYLEVEL 2

/* #define STREETCARFLAG true */
#define STREETCARFLAG true

/* #define RANDVEL true */
#define RANDVEL true

#define PI 3.1415926536

#define SCREENWIDTH 80

main(argc, argv)
int argc;
char *argv[];
{

  double pos[NUMCARS], vel[NUMCARS];
  double upper, nextvel;
  double min(), max(), normal(), reflection(), drand48(), sq();
  double time[MAXTIMES], timeuntil, timesince, prevtime;
  int i, j, itnum, done, timenum, counting, spacenum, carcount, PRINTLENGTH;
  int simcount, simnum, numepss, numepsu, numepseps;
  int sumpass, sumtimes;
  double sumsince, sumuntil, sssince, ssuntil, sumprod, thetime, y;
  double meansince, meanuntil, meanprod, cov;
  double meansqsince, meansquntil, sdsince, sduntil, corr;
  double fracepsu, fracepseps;
  double begintime, endtime;
  struct timeval tmptv;

  /* Initialisations. */
  seedrnd();
  dsfactor = sqrt(DELTA) * SIGMA;
  simcount = 0;
  sumsince = sumuntil = sumprod = sssince = ssuntil = 0.0;
  sumpass = sumtimes = numepss = numepsu = numepseps = 0;
  PRINTLENGTH = 10;
  if (LEADTIME*NUMREPS/DELTA >= 1000000)
    PRINTLENGTH = 1;

  /* Start up. */
  printf("\nBeginning %d reps of", NUMREPS);
  if (DATACHOICE==0) {
    printf(" %d sims of %.1f hrs of %d", NUMSIMS, LEADTIME+1, NUMCARS);
#ifdef STREETCARFLAG
  printf(" streetcars (diam=%f)", DIAMETER);
#else
  printf(" buses");
#endif
  } else {
    printf(" data set #%d.", DATACHOICE);
  }
  printf("\n");
  gettimeofday (&tmptv, (struct timezone *)NULL);
  begintime = 1.0 * tmptv.tv_sec + 0.000001 * tmptv.tv_usec;

  /* Do the repeated simulations. */
  for (j=0; j<NUMREPS; j++) {

if (DATACHOICE==0) {
/* Do simulations. */

    /* Initialise the positions, velocities, etc. */
    done = itnum = timenum = counting = 0;
    prevtime = 0.0;
    for (i=0; i<NUMCARS; i++) {
	pos[i] = ((double)i) / NUMCARSPER - ROUTELENGTH;
	vel[i] = 1.0;
    }

    /* Do the iterations. */
    while (done == 0) {

	/* Try moving the cars, in reverse order. */
	for (i=NUMCARS-1; i>=0; i--) {

#ifdef RANDVEL
	    vel[i] = reflection( vel[i] + dsfactor * normal() );
#endif
	    pos[i] = pos[i] + DELTA * vel[i];

#ifdef STREETCARFLAG
	    /* Compute upper, etc. */
	    if (i==NUMCARS-1) {
		upper = pos[0];
		nextvel = vel[0];
	    } else {
		upper = pos[i+1];
		nextvel = vel[i+1];
	    }
	    if (upper < pos[i])
		upper = upper + ROUTELENGTH;
	    upper = upper - DIAMETER;

	    /* Slow down if reached upper. */
	    if (pos[i] > upper) {
#ifdef VERBOSE
	      printf("REDUCING pos[%d] from %f to %f!\n", i, pos[i], upper);
#endif
	      pos[i] = upper;
	      vel[i] = min( vel[i], nextvel );
	    }
#endif

	    /* Check for cars "over the edge". */
	    if (pos[i] >= 0.0) {
		pos[i] = pos[i] - ROUTELENGTH;
		sumpass++;
		if (counting==0) {
		    prevtime = DELTA * itnum;
		} else {
		    time[timenum] = itnum * DELTA - LEADTIME;
		    timenum++;
		    sumtimes++;
		    if (itnum * DELTA >= LEADTIME + 1.0) {
			done = 1;
		    }
		}
            }

	} /* End of "i" for loop. */

	if( (DISPLAYLEVEL>=4) ||
	    ( (DISPLAYLEVEL>=2) && (divisible(itnum,10000)) ) ||
		( (DISPLAYLEVEL>=1) && (done==1) ) ) {
	  for (spacenum=0; spacenum<SCREENWIDTH; spacenum++) {
	    carcount = 0;
	    for (i=0; i<NUMCARS; i++) {
	      if ( ifloor( (pos[i]+1) * SCREENWIDTH) == spacenum )
	        carcount = carcount + 1;
	    }
	    /* printf("%d", carcount); */
	    printthechar(carcount);
	  }
	  printf("\n");
	  if (DISPLAYLEVEL>=4)
	    usleep(20000);
	}

	if ( (counting==0) && (itnum * DELTA >= LEADTIME) ) {
	    /* Start counting arrival times. */
	    counting = 1;
	    time[timenum] = prevtime - LEADTIME;
	    timenum++;
	}

        itnum++;
    } /* End of "while" loop. */

    if (DISPLAYLEVEL >= 2) {
      printf("\nEND OF ITERATION %d.   %d ARRIVAL TIMES.\n", j, timenum);
      if (DISPLAYLEVEL >= 3) {
	for (i=0; i<NUMCARS; i++) {
	    printf("pos[%d] = %f\n", i, pos[i]);
        }
        printf("\n");
	for (i=0; i<timenum; i++) {
	    printf("time[%d] = %f\n", i, time[i]);
        }
      }
      printf("\n");
    }

} else if (DATACHOICE==1) {
/* Use data from January 25, 2007. */

timenum = 30;

time[0] = -0.003888889;
time[1] = 0.018888889;
time[2] = 0.029444444;
time[3] = 0.040277778;
time[4] = 0.080000000;
time[5] = 0.204166667;
time[6] = 0.222222222;
time[7] = 0.257777778;
time[8] = 0.267222222;
time[9] = 0.273888889;
time[10] = 0.293055556;
time[11] = 0.372777778;
time[12] = 0.405555556;
time[13] = 0.418611111;
time[14] = 0.439166667;
time[15] = 0.552777778;
time[16] = 0.587777778;
time[17] = 0.618611111;
time[18] = 0.644166667;
time[19] = 0.665555556;
time[20] = 0.695555556;
time[21] = 0.702777778;
time[22] = 0.715277778;
time[23] = 0.742222222;
time[24] = 0.795833333;
time[25] = 0.818888889;
time[26] = 0.843055556;
time[27] = 0.874166667;
time[28] = 0.923055556;
time[29] = 1.004444444;

} else if (DATACHOICE==2) {
/* Use data from January 29, 2007. */

timenum = 31;

time[0] = -0.0002777778;
time[1] = 0.0202777778;
time[2] = 0.0283333333;
time[3] = 0.0847222222;
time[4] = 0.1202777778;
time[5] = 0.1441666667;
time[6] = 0.1650000000;
time[7] = 0.1886111111;
time[8] = 0.2202777778;
time[9] = 0.2713888889;
time[10] = 0.2911111111;
time[11] = 0.3022222222;
time[12] = 0.3669444444;
time[13] = 0.4250000000;
time[14] = 0.4550000000;
time[15] = 0.4677777778;
time[16] = 0.5083333333;
time[17] = 0.5161111111;
time[18] = 0.5522222222;
time[19] = 0.5952777778;
time[20] = 0.6419444444;
time[21] = 0.6538888889;
time[22] = 0.6630555556;
time[23] = 0.6711111111;
time[24] = 0.8152777778;
time[25] = 0.8538888889;
time[26] = 0.8819444444;
time[27] = 0.9047222222;
time[28] = 0.9161111111;
time[29] = 0.9780555556;
time[30] = 1.0208333333;

} else {

/* DATACHOICE value is out of range. */
fprintf(stderr, "\nError -- invalid data choice.\n\n");
exit(1);

}

    for (simnum=0; simnum<NUMSIMS; simnum++) {

	/* Simulate a rider arriving at position "y". */
	y = drand48();
	timesince = timeuntil = 10000.0;
	for (i=0; i<timenum; i++) {
	    thetime = time[i] - y;
	    if ( (thetime >= 0) && (thetime < timeuntil) )
		timeuntil = thetime;
	    if ( (-thetime >= 0) && (-thetime < timesince) )
		timesince = -thetime;
	}
	simcount++;
	sumsince = sumsince + timesince;
	sumuntil = sumuntil + timeuntil;
	sssince = sssince + sq(timesince);
	ssuntil = ssuntil + sq(timeuntil);
	sumprod = sumprod + timesince * timeuntil;
	if ( timeuntil < EPSILON )
          numepsu++;
	if ( timesince < EPSILON ) {
          numepss++;
	  if ( timeuntil < EPSILON )
            numepseps++;
        }
    }

    /* Output this trial's results to stdout. */
    if ( (DISPLAYLEVEL >= 1) || (divisible(j, PRINTLENGTH)) ) {
      printf("%4d / %d:  sumuntil = %f,  sumsince = %f\n",
			j+1, NUMREPS, sumuntil, sumsince);
    }

  } /* End of "j" for loop. */

  /* Wrap up. */
  printf("\nDone: %d reps of", NUMREPS);
  if (DATACHOICE==0) {
    printf(" %d sims of %.1f hrs of %d", NUMSIMS, LEADTIME+1, NUMCARS);
#ifdef STREETCARFLAG
  printf(" streetcars (diam=%f)", DIAMETER);
#else
  printf(" buses");
#endif
  } else {
    printf(" data set #%d.", DATACHOICE);
  }
  printf("\n");
  printf("NUMCARSPER=%d, SIGMA=%f, DELTA=%f\n", NUMCARSPER, SIGMA, DELTA);
  gettimeofday (&tmptv, (struct timezone *)NULL);
  endtime = 1.0 * tmptv.tv_sec + 0.000001 * tmptv.tv_usec;

  /* Final output. */
  printf("meannumpass=%.1f, meannumtimes=%.1f, ellapsed time = %.2f secs.\n",
    ((double)sumpass)/NUMREPS, ((double)sumtimes)/NUMREPS, endtime-begintime);
  meansince = sumsince/simcount;
  meanuntil = sumuntil/simcount;
  meanprod = sumprod/simcount;
  meansqsince = sssince/simcount;
  meansquntil = ssuntil/simcount;
  sdsince = sqrt(meansqsince - sq(meansince));
  sduntil = sqrt(meansquntil - sq(meanuntil));
  cov = meanprod - meansince*meanuntil;
  corr = cov / sdsince / sduntil;
  printf("meansince = %f,  meanuntil = %f, meanprod = %f\n",
  					meansince, meanuntil, meanprod);
  printf("meansqsince = %f,  meansquntil = %f,\nsdsince = %f, sduntil = %f\n",
  		meansqsince, meansquntil, sdsince, sduntil);
  printf("cov = %f,  corr = %f\n", cov, corr);
  fracepsu = ((double)numepsu)/simcount;
  fracepseps = ((double)numepseps)/numepss;
  printf("fracepsu=%f, fracepseps=%f, perc=%f (eps=%f)\n",
	fracepsu, fracepseps, 100.0 * fracepseps / fracepsu, EPSILON );

  printf("\n");
  return(0);

}


/* REFLECTION -- reflect aa to be within [0,1] */
double reflection(double aa)
{
    double reflection();

    if (aa > 1.0)
	return( reflection(2.0 - aa) );
    else if (aa < 0.0)
	return( reflection(- aa) );
    else
	return( aa );
}


/* PRINTTHECHAR */
printthechar(int nn)
{
    if (nn==0)
      printf(" ");
    else if (nn==1)
      printf("|");
    else if (nn==2)
      printf("+");
    else if (nn==3)
      printf("X");
    else
      printf("*");
}


/* NORMAL:  return a standard normal random number. */
double normal()
{
    double R, theta, drand48();

    R = - log(drand48());
    theta = 2 * PI * drand48();

    return( sqrt(2*R) * cos(theta));
}


/* SEEDRND. */
seedrnd()
{
    int seed;
    struct timeval tmptv;
    gettimeofday (&tmptv, (struct timezone *)NULL);
    seed = (int) (tmptv.tv_usec - 1000000 * ifloor ( ( (double)
tmptv.tv_usec )
                                  / ( (double) 1000000.0 ) ) );
    srand48(seed);
}


/* IFLOOR */
int ifloor(double x)  /* returns floor(x) as an integer */
{
    return((int)floor(x));
}


/* SQ */
double sq(double xxx)
{
    return(xxx * xxx);
}


/* MIN */
double min(double aa, double bb)
{
    if (aa < bb)
	return(aa);
    return(bb);
}


/* DIVISIBLE */
int divisible (int n1, int n2)
/* Check if n1 is divisible by n2. */
{
  return (n1 == n2 * (n1/n2));
}