// 2 Roads Merging NaSch CA model// By David Clarke ( dave@maths.tcd.ie )// Original source available at www.maths.tcd.ie/~dave/traffic///Compile with cc -o traffic traffic.c#include<stdio.h>#include<stdlib.h>#include<math.h>#define leadin_d 500 #define overlap_d 10#define leadout_d 500#define max_length 2100  // greater then leadin+overlap+leadout#define speed_limit1_d 5#define speed_limit2_d 5#define steps_d 20000 		//Number of Steps updates#define cars1_d 0		// cars to start with (0 is empty road)#define cars2_d 0#define cars0_d 0 // leadout cars#define prob_d 0.5		// probability of stopping#define slow_start_d 0.0#define relax_d 600		// steps ignore before measuring#define flux_time_d 2000		//steps flux measured over#define print_d 0#define split 10000000#define once 0 			// if 1, run for only one density, for scatter plot.#define max_density 1 	// limiter 0->1, can stop program for a lower density typedef struct {  int car;  int vel;  int moved;  int pre_vel;			//previous velocity, used in slow to start  int jtime;  int origin;  int num;} road_t;FILE *out1;FILE *outScatter1, *outScatter2;FILE *out_jtime1, *out_jtime2;//Function definitionsvoid print_road (int);void print_jams ();void print_flow_density();void print_scatter (int c, int tstep);void setup (int);void setup_leadout ();void setspeed (int);void update (int c, int tstep);void measure_density();void variables (int argc, char *argv[]);//global variablesroad_t Road1[max_length],Road2[max_length];int leadin,overlap,leadout, speed_limit[3], steps, cars[3], relax, density_width, flux_time,tot_length,print;double prob, slow_start;int count[3],density_start[3],density_stop[3], time_count, density[3], jtimeMax[3],jtimeMin[3];int main (int argc, char *argv[]){		/* Start up */  int i,j, n, d,pr;  road_t *roadP;  variables (argc, argv);	//take in variables from command line  tot_length=leadin+overlap+leadout;//density measurement parameters.  density_start[1]=density_start[2]=leadin/2;  density_start[0]=leadin+overlap+leadout/2;  density_stop[1]=density_stop[2]=leadin;  density_stop[0]=leadin+overlap+leadout;  out1 = fopen ("density_flow.dat", "w");  outScatter1 = fopen ("trafficYScatter1.dat", "w");  outScatter2 = fopen ("trafficYScatter2.dat", "w");  out_jtime1 = fopen ("jtime1.dat", "w");  out_jtime2 = fopen ("jtime2.dat", "w");  srand (time (0));		// random seedwhile(cars[1]<leadin*max_density&&cars[2]<leadin*max_density){ // Looping for increasing density 		/* Set Up */	n = 0;	setup (1);			// initial positioning of cars	setup (2);			// initial positioning of cars	setup_leadout();	count[1]=count[2]=1;	if(print){	print_road(1);       // initial print		print_road(2);       // initial print	}		/* Relax */		while(n < relax) {n++; //allow system to settle		setspeed(1); 		update(1,-1); 		setspeed(2); 		update(2,-1);}   			//	Give each car a number:	for(i=0;i<leadin+overlap+leadout;i++){		if(Road1[i].car==1){Road1[i].num=count[Road1[i].origin];			count[Road1[i].origin]++;}		if(Road2[i].car==1){Road2[i].num=count[Road2[i].origin];			count[Road2[i].origin]++;}		}	time_count = flux_time;	count[0]=count[1]=count[2] = n = 0;	density[0]=density[1]=density[2]=0;  				/* Run Model */  while (n < steps) {	n++;	time_count--;	setspeed (1);   	update (1,n); 	setspeed (2);   	update (2,n);	if(print){printf("r1"); print_road(1);		 	printf("r2"); print_road(2);	}	if(once){			print_scatter(1,n); // Print Scatter diagram			print_scatter(2,n); 	}	measure_density();	if (time_count == 0)print_flow_density();		//    print_jams ();  } // steps  if(once) break; // do program just once (for scatter)	cars[1]+=5;	cars[2]+=5;	cars[0]+=5;}  fclose (out1);  fclose (outScatter1);  fclose (outScatter2);  fclose (out_jtime1);  fclose (out_jtime2);fprintf(stderr,"\nprogram complete\n");exit (0);}		//end of mainvoid setup (int c){	int i;	int cars_used;	road_t *roadP;	switch(c){		case 1: roadP=&Road1[0];		break;		case 2:  roadP=&Road2[0];		break;		default: fprintf(stderr,"Error! Unknown case in setup\n");}		cars_used=cars[c];	for (i = 0; i < tot_length; i++) {		(roadP+i)->car = 0;		(roadP+i)->vel = 0; }	while (cars_used > 0) {    	for (i = 0; i < leadin; i++) {	// initialise array      		if (((int) (0 + (1.0*leadin/cars_used)*(double)rand()				/(RAND_MAX + 1.0)))==1				&& cars_used > 0 && (roadP+i)->car == 0) {			(roadP+i)->car = 1;			(roadP+i)->vel =(int)(0+(speed_limit[c]+1)*(double)rand()/(RAND_MAX + 1.0));			cars_used--;			(roadP+i)->jtime=-split;			(roadP+i)->origin=c;      		}	}}}void setup_leadout (){	int i;	int cars_used=cars[0];	road_t *roadP;	roadP=&Road1[leadin+overlap];	while (cars_used > 0) {    	for (i = 0; i < leadout; i++) {	// initialise array      		if (((int) (0 + (1.0*leadout/cars_used)*(double)rand()				/(RAND_MAX + 1.0)))==1				&& cars_used > 0 && (roadP+i)->car == 0) {			(roadP+i)->car = 1;			(roadP+i)->vel =(int)(0+(speed_limit[1]+1)*(double)rand()/(RAND_MAX + 1.0));			cars_used--;			(roadP+i)->jtime=-split;			if(cars[1]==0)			(roadP+i)->origin=2;			else if(cars[2]==0)			(roadP+i)->origin=1;			else if((int)((1+(1.0*cars[2])/cars[1])*(double)rand()/(RAND_MAX +1.0))==0)			(roadP+i)->origin=1;			else			(roadP+i)->origin=2;//fprintf(stderr,"rand: %d\n",(int)((1+(1.0*cars[2])/cars[1])*(double)rand()/(RAND_MAX +1.0)));      		}	}}}voidsetspeed (int c){	int i, j, dist, k, d, p,length;	road_t *roadP, *roadO;	switch(c){		case 1: length=leadin+overlap+leadout;		roadP=&Road1[0];		break;		case 2: length=leadin+overlap;		roadP=&Road2[0];		break;		default: fprintf(stderr,"Error! Unknown case in setspeed\n");	}	for (i = 0; i < length; i++) {	// for each road segment    	if ((roadP+i)->car == 1) {	// only look at segments with cars			(roadP+i)->pre_vel = (roadP+i)->vel;			//1. accelerate			if ((roadP+i)->vel < speed_limit[c])(roadP+i)->vel++; 						dist = 0;			// Find distance to next car (up to speed limit)			for(j=1;j<=speed_limit[c];j++){				if (i+j>=length) {					if(c==1){  //Main branch -off the edge					// Find average 						if((roadP+i)->origin==1)roadO=&Road1[0];						else roadO=&Road2[0];						if((roadO+i+j-length)->car==0) dist++;						else { 	break;}					}					else  break;   // Slip lane, stop at end				}			//loop road				else if ((roadP+i+j)->car == 0) dist++;				else break; 			}			//debug			//printf("i:%d j:%d k:%d dist:%d\n",i,j,k,dist);			//printf("vel:%d\n",(roadP+i)->vel);			//2. Deceleration (with infinate breaking power)			if ((roadP+i)->vel > dist) (roadP+i)->vel = dist;							//printf("vel:%d\n",Road[i].vel);						// 3. Random Deceleration      		if((roadP+i)->vel>0&&((int)(0+(1.0/prob)*(double)rand()/				(RAND_MAX+1.0)))==0) 				(roadP+i)->vel--;				}	}}void update (int c,int tstep){	int i, n = 0, d,length,v,k,m;	road_t *roadP, *roadO;	switch(c){		case 1: length=leadin+overlap+leadout;		roadP=&Road1[0];		break;		case 2: length=leadin+overlap;		roadP=&Road2[0];		break;		default: fprintf(stderr,"Error! Unknown case in setspeed\n");	}	roadO=roadP;	for(i=0;i<length;i++)(roadP+i)->moved=0;// Don't move cars twice		for(i=length-1;i>=0;i--){ //move backwards because of merging lane		if ((roadP+i)->car==1&&(roadP+i)->moved==0){ //if car&not moved			k=m=0;			/* ### Merging  ### */			/* ###  Rules   ### *//* Move over at the earliest opportunity, then move up lane as much as possible ie. min(vel,dist to next car) */			if(c==2&&i>=leadin){	//in merging zone				v=Road2[i].vel;				for(k=0;k<=v;k++){  //what k can we move over at?					if(Road1[i+k].car==0&&Road1[i+k-1].car==0){						m++; //succeeded in moving over (and 1 space to spare)						d=k; } // max advancement					if(m>0&&Road1[i+k].car==1) break; //moved over and behind next car				}				if(m==0){ //couldn't move, move in lane at reduced speed					if(v>1){						v--;						if(Road2[i+v].car==1)fprintf(stderr,"Crash while not merging!!");						Road2[i+v].car=1;						Road2[i].car=0;						Road2[i+v].vel=v;						Road2[i+v].jtime=Road2[i].jtime+1;						Road2[i+v].origin=2;						Road2[i+v].num=Road2[i].num;						}				}				else{					if(Road1[i+d].car==1)fprintf(stderr,"Crash while merging!!");					Road1[i+d].car=1;					Road1[i+d].vel=v;					Road1[i+d].jtime=Road2[i].jtime+1;					Road1[i+d].num=Road2[i].num;					Road1[i+d].origin=2;					Road2[i].car=0;				}			}	/* ### Other motion ### */				else if((roadP+i)->vel == 0)(roadP+i)->jtime++; // car doesn't move			else if(i+(roadP+i)->vel>=length){ // car loops.				count[0]++; 				d=(roadP+i)->vel+i-length;//car now loops around				if((roadP+i)->origin==2){ roadO=&Road2[0];}					/* calculate its journey time */					// if time <0 car came from middle,ignore 				if((roadP+i)->jtime>0&&tstep!=-1){					if((roadP+i)->origin==1)						fprintf(out_jtime1,"1\t%d\t%d\n",tstep,(roadP+i)->jtime);					else fprintf(out_jtime2,"2\t%d\t%d\n",tstep,(roadP+i)->jtime);				}				if((roadO+d)->car == 1)fprintf(stderr,"\n *** ### Car crash! ### ***\n");				(roadO+d)->car = 1;				(roadO+d)->vel = (roadP+i)->vel;				(roadO+d)->jtime = 0;				(roadO+d)->num = (roadP+i)->num;				(roadO+d)->origin = (roadP+i)->origin;				(roadO+d)->moved = 1;				(roadP+i)->car = 0;			}				// 4. Update & Move			else {				d = (roadP+i)->vel + i;								if((roadP+d)->car == 1)fprintf(stderr,"\n *** ### Car crash! ### ***\n");				(roadP+d)->car = 1;				(roadP+d)->vel = (roadP+i)->vel;				(roadP+d)->jtime = (roadP+i)->jtime+1;				(roadP+d)->num = (roadP+i)->num;				(roadP+d)->origin = (roadP+i)->origin;				(roadP+d)->moved = 1;				(roadP+i)->car = 0;				if(i<leadin&&d>=leadin){count[c]++;}//car has left the zone			} 					}	}}  //end of update()void measure_density(){	int i,j;	road_t *roadP;	for(j=0;j<3;j++){		switch(j){			case 0: roadP=&Road1[0]; break;			case 1: roadP=&Road1[0]; break;			case 2: roadP=&Road2[0]; break;			default: fprintf(stderr,"error measuring density");}	for (i=density_start[j];i<density_stop[j];i++){		if ((roadP+i)->car == 1) density[j]++;      	}  	//	if(print)  printf("D%d:%d ",j,density[j]);	}}voidprint_road (int c){	int i,length;	road_t *roadP;	switch(c){		case 1: roadP=&Road1[0];		length=leadin+overlap+leadout;		break;		case 2:  roadP=&Road2[0];  		length=leadin+overlap;  		break;		default: fprintf(stderr,"Error! Unknown case in print road\n");}	for (i = 0; i < tot_length; i++) {		if(i==0||i==leadin||i==leadin+overlap||i==leadin+overlap+leadout)	  	printf("|");    	if ((roadP+i)->car == 0) {			printf ("-");    	}	    else {      printf ("%d", (roadP+i)->vel);    // print velocities//      printf ("%d", (roadP+i)->origin); // print origin//      printf ("%d", (roadP+i)->num%10); //print car number mod 10    	} 	}	printf ("\n");}void print_scatter(int c,int tstep){	int i,length;	road_t *roadP;	FILE *Scatter;	switch(c){		case 1: roadP=&Road1[0];		length=leadin+overlap+leadout;		Scatter=outScatter1;		break;	case 2:  roadP=&Road2[0];		length=leadin+overlap;		Scatter=outScatter2;		break;		default:fprintf(stderr,"Error! Unknown case in print scatter\n");}	for (i = 0; i < length; i++) {    	if ((roadP+i)->car == 1&&(roadP+i)->num%5==0) {			if(c==2)fprintf(Scatter,"%d\t%d\n",tstep,i);	/* Separate for plotting different colours based on origin */		else if(Road1[i].origin==1)fprintf(Scatter,"%d\t%d\t-\n",tstep,i);		else if(Road1[i].origin==2)fprintf(Scatter,"%d\t-\t%d\n",tstep,i);		else fprintf(stderr,"Error, origin not 1 or 2");		}	}}void print_flow_density(){int j;	for(j=0;j<3;j++){		fprintf(out1,"%lf\t%lf\t",(density[j]*1.0)/(flux_time*(density_stop[j]-density_start[j])),(count[j]*1.0)/flux_time);		time_count = flux_time;		count[j] = 0;		density[j] = 0;	}	fprintf(out1,"\n");}//Take in variables from command linevoid variables (int argc, char *argv[]){	if (argc == 1) {		leadin = leadin_d;		overlap = overlap_d;		leadout = leadout_d;		speed_limit[1] = speed_limit1_d;		speed_limit[2] = speed_limit2_d;		steps = steps_d;		cars[1] = cars1_d;		cars[2] = cars2_d;		cars[0] = cars0_d;		relax = relax_d;		flux_time = flux_time_d;		prob = prob_d;		slow_start = slow_start_d;		leadin = leadin_d;		print= print_d;	}	else if(argc !=14 ) {		fprintf (stderr,"\nUsage: %s <leadin> <overlap> <leadout> <cars1> <cars2> <cars leadout> <speed1> <speed2> <steps> <relax> <prob> <fluxTime> <print>\n\t \n\n",argv[0]);		fprintf(stderr,"Try: %s  %d %d %d %d %d %d %d %d %d %d %lf %d %d\n\n", argv[0],leadin_d,overlap_d,leadout_d,cars1_d,cars2_d,cars0_d,speed_limit1_d,speed_limit2_d,steps_d,relax_d,prob_d,flux_time_d,print_d);		exit (1);    }	else {	if (sscanf (argv[1], "%d", &leadin) != 1) {		fprintf(stderr,"\narg1: road length not an int\n\n");		exit (1);}	if (sscanf (argv[2], "%d", &overlap) != 1) {		fprintf(stderr,"\narg2: road length not an int\n\n");		exit (1); }	if (sscanf (argv[3], "%d", &leadout) != 1) {		fprintf(stderr,"\narg3: road length not an int\n\n");		exit (1); }	if(leadin+overlap+leadout>max_length){		fprintf(stderr,"\n Total length greater then max length\n\n");		exit(1);}	if (sscanf (argv[4], "%d", &cars[1]) != 1) {		fprintf(stderr,"\narg4:Cars1 not an int\n\n");		exit (1); }	if (sscanf (argv[5], "%d", &cars[2]) != 1) {		fprintf(stderr,"\narg5:Cars2 not an int\n\n");		exit (1); }	if (sscanf (argv[6], "%d", &cars[0]) != 1) {		fprintf(stderr,"\narg5:Cars leadout not an int\n\n");		exit (1);}	if (sscanf (argv[7], "%d", &speed_limit[1]) != 1) {		fprintf(stderr,"\narg7:Speed Limit 1 not an int\n\n");		exit (1); }	if (sscanf (argv[8], "%d", &speed_limit[2]) != 1) {		fprintf(stderr,"\narg8:Speed Limit 2 not an int\n\n");		exit (1);  }	if (sscanf (argv[9], "%d", &steps) != 1) {		fprintf(stderr,"\narg9: Steps not an int\n\n");		exit (1); }	if (sscanf (argv[10], "%d", &relax) != 1) {		fprintf(stderr,"\narg10: relax not an int\n\n");		exit (1);}	if (sscanf (argv[11], "%lf", &prob) != 1 || prob > 1) {		fprintf(stderr,"\narg11:Slowing prob not less then 1\n\n");		exit (1);}	if(cars[1]>leadin||cars[2]>leadin){		fprintf(stderr,"\nToo many cars for road length\n\n");		exit(1);}	if (sscanf (argv[12], "%d", &flux_time) != 1) {		fprintf(stderr,"\narg12:Flux time not an int\n\n");		exit (1);}	if (sscanf (argv[13], "%d", &print) != 1&&(print==0||print==1)) {		fprintf(stderr,"\narg13:Print is not 0 or 1\n\n");		exit (1);}    }}