/* modified version */ 	else if (VERSION == 2) {		/* probability that the first trip would be a pause */		double prob_pause = PAUSE / (EXP_1_V*EXP_R + PAUSE);			/* the first trip is a pause */		if (prob_pause > uniform()) {			/* constant pause */			if (PAUSETYPE == 1) {				time_arrival = TIME + PAUSE;				// constant pause			}			/* uniform pause */			else if (PAUSETYPE == 2) {				time_arrival = TIME + 2*PAUSE*uniform();	// uniform pause [0, 2*PAUSE]			}							first_trip = 0;						// indicating the first trip is a pause 		}		/* the first trip is a move based on the steady-state pdf */		else {			time_arrival = TIME;			first_trip = 1;						// indicating the first trip is a move 		}	}	time_update = TIME;	time_transition = 0.0;	position.X = position.Y = position.Z = 0.0;	destination.X = destination.Y = destination.Z = 0.0;	direction.X = direction.Y = direction.Z = 0.0;	speed = 0.0;	RandomPosition();	fprintf(stdout, NODE_FORMAT3, index, 'X', position.X);	fprintf(stdout, NODE_FORMAT3, index, 'Y', position.Y);	fprintf(stdout, NODE_FORMAT3, index, 'Z', position.Z);	neighbor = new Neighbor[NODES];	if(neighbor == 0) {		perror("new");		exit(1);	}	for(i = 0; i < NODES; i++) {		neighbor[i].index = i;		neighbor[i].reachable = (index == i) ? 1 : 0;		neighbor[i].time_transition = 0.0;	}}voidNode::RandomPosition(){    position.X = uniform() * MAXX;    position.Y = uniform() * MAXY;	position.Z = 0.0;}voidNode::RandomDestination(){   	destination.X = uniform() * MAXX;   	destination.Y = uniform() * MAXY;	destination.Z = 0.0;	assert(destination != position);}/******************************************************************************************  * Speeds are chosen based on the given type and distribution ******************************************************************************************/voidNode::RandomSpeed(){	/* original version */	if (VERSION == 1) {       	speed = uniform() * MAXSPEED;		assert(speed != 0.0);	}	/* modified version */	else if (VERSION == 2) {		/* uniform speed */		if (SPEEDTYPE == 1) {			/* using steady-state distribution for the first trip */			if (first_trip == 1) {				/* pick a speed by rejection technique */				double temp_v, temp_fv;					do {	        	 	temp_v = uniform() * (MAXSPEED - MINSPEED) + MINSPEED;					temp_fv = uniform() * PDFMAX;				} while (temp_fv > 1/temp_v*EXP_R / (EXP_1_V*EXP_R + PAUSE) / (MAXSPEED-MINSPEED));	 				speed = temp_v;				first_trip = 0;		// reset first_trip flag 				}			/* using the original distribution from the second trip on */			else {	        	speed = uniform() * (MAXSPEED - MINSPEED) + MINSPEED;				assert(speed != 0.0);			}		}		/* normal speed */		else if (SPEEDTYPE == 2) {			/* using steady-state distribution for the first trip */			if (first_trip == 1) {				double temp_v, temp_fv, square, fv;					/* rejection technique */				do {	        	 	temp_v = uniform() * (MAXSPEED - MINSPEED) + MINSPEED;					temp_fv = uniform() * PDFMAX;					square = (temp_v - MEAN)*(temp_v - MEAN)/SIGMA/SIGMA;					fv = 1/KAPPA/sqrt(2.0*PI*SIGMA*SIGMA) * exp(-0.5*square);				} while (temp_fv > 1.0/temp_v*fv*EXP_R / (EXP_1_V*EXP_R + PAUSE));	 				speed = temp_v;				first_trip = 0;			}			/* using the original distribution from the second trip on */			else {				double temp_v, temp_fv, square;				double max_normal = 1.0/KAPPA/sqrt(2.0*PI*SIGMA*SIGMA);		// max of normal distribution					/* rejection technique */				do {	         		temp_v = uniform() * (MAXSPEED - MINSPEED) + MINSPEED;					temp_fv = uniform() * max_normal;					square = (temp_v - MEAN)*(temp_v - MEAN)/SIGMA/SIGMA;				} while (temp_fv > max_normal * exp(-0.5*square));	 				speed = temp_v;				assert(speed != 0.0);			}		}		/* other types of speed for future use */		else			;	}}voidNode::Update(){	position += (speed * (TIME - time_update)) * direction;	if(TIME == time_arrival) {		vector v;		if(speed == 0.0 || PAUSE == 0.0) {           	RandomDestination();			RandomSpeed();			v = destination - position;			direction = v / v.length();			time_arrival = TIME + v.length() / speed;		}		else {			destination = position;			speed = 0.0;			/* original version */			if (VERSION == 1) {				time_arrival = TIME + PAUSE;			}			/* modified version */			else if (VERSION == 2) {				/* constant pause */				if (PAUSETYPE == 1) {					time_arrival = TIME + PAUSE;				}				/* uniform pause */				else if (PAUSETYPE == 2) {					time_arrival = TIME + 2*PAUSE*uniform();				}			}		}		fprintf(stdout, NODE_FORMAT,			TIME, index, destination.X, destination.Y, speed);		}	time_update = TIME;	time_transition = 0.0;}voidNode::UpdateNeighbors(){	static Node *n2;	static Neighbor *m1, *m2;	static vector D, B, v1, v2;	static double a, b, c, t1, t2, Q;	static u_int32_t i, reachable;	v1 = speed * direction;	/*	 *  Only need to go from INDEX --> N for each one since links	 *  are symmetric.	 */	for(i = index+1; i < NODES; i++) {		m1 = &neighbor[i];		n2 = &NodeList[i];		m2 = &n2->neighbor[index];		assert(i == m1->index);		assert(m1->index == n2->index);		assert(index == m2->index);                assert(m1->reachable == m2->reachable);                reachable = m1->reachable;		/* ==================================================		   Determine Reachability		   ================================================== */		{	vector d = position - n2->position;			if(d.length() < RANGE) {#ifdef SANITY_CHECKS				if(TIME > 0.0 && m1->reachable == 0)					assert(RANGE - d.length() < ROUND_ERROR);#endif				m1->reachable = m2->reachable = 1;			}			// Boundary condition handled below.			else {#ifdef SANITY_CHECKS				if(TIME > 0.0 && m1->reachable == 1)					assert(d.length() - RANGE < ROUND_ERROR);#endif				m1->reachable = m2->reachable = 0;			}#ifdef DEBUG                        fprintf(stdout, "# %.6f (%d, %d) %.2fm\n",                                TIME, index, m1->index, d.length());#endif		}		/* ==================================================		   Determine Next Event Time		   ================================================== */		v2 = n2->speed * n2->direction;		D = v2 - v1;		B = n2->position - position;		a = (D.X * D.X) + (D.Y * D.Y) + (D.Z * D.Z);		b = 2 * ((D.X * B.X) + (D.Y * B.Y) + (D.Z * B.Z));		c = (B.X * B.X) + (B.Y * B.Y) + (B.Z * B.Z) - (RANGE * RANGE);		if(a == 0.0) {			/*			 *  No Finite Solution			 */			m1->time_transition= 0.0;			m2->time_transition= 0.0;			goto  next;		}		Q = b * b - 4 * a * c;		if(Q < 0.0) {			/*			 *  No real roots.			 */			m1->time_transition = 0.0;			m2->time_transition = 0.0;			goto next;		}		Q = sqrt(Q);		t1 = (-b + Q) / (2 * a);		t2 = (-b - Q) / (2 * a);		// Stupid Rounding/Boundary Cases		if(t1 > 0.0 && t1 < ROUND_ERROR) t1 = 0.0;		if(t1 < 0.0 && -t1 < ROUND_ERROR) t1 = 0.0;		if(t2 > 0.0 && t2 < ROUND_ERROR) t2 = 0.0;		if(t2 < 0.0 && -t2 < ROUND_ERROR) t2 = 0.0;		if(t1 < 0.0 && t2 < 0.0) {			/*			 *  No "future" time solution.			 */			m1->time_transition = 0.0;			m2->time_transition = 0.0;			goto next;		}		/*		 * Boundary conditions.		 */		if((t1 == 0.0 && t2 > 0.0) || (t2 == 0.0 && t1 > 0.0)) {			m1->reachable = m2->reachable = 1;			m1->time_transition = m2->time_transition = TIME + max(t1, t2);		}		else if((t1 == 0.0 && t2 < 0.0) || (t2 == 0.0 && t1 < 0.0)) {			m1->reachable = m2->reachable = 0;			m1->time_transition = m2->time_transition = 0.0;		}		/*		 * Non-boundary conditions.		 */		else if(t1 > 0.0 && t2 > 0.0) {			m1->time_transition = TIME + min(t1, t2);			m2->time_transition = TIME + min(t1, t2);		}		else if(t1 > 0.0) {			m1->time_transition = TIME + t1;			m2->time_transition = TIME + t1;		}		else {			m1->time_transition = TIME + t2;			m2->time_transition = TIME + t2;		}		/* ==================================================		   Update the transition times for both NODEs.		   ================================================== */		if(time_transition == 0.0 || (m1->time_transition &&		   time_transition > m1->time_transition)) {			time_transition = m1->time_transition;		}		if(n2->time_transition == 0.0 || (m2->time_transition &&		   n2->time_transition > m2->time_transition)) {			n2->time_transition = m2->time_transition;		}        next:                if(reachable != m1->reachable && TIME > 0.0) {                        LinkChangeCount++;                        link_changes++;                        n2->link_changes++;                }	}}voidNode::Dump(){	Neighbor *m;	u_int32_t i;	fprintf(stdout,		"Node: %d\tpos: (%.2f, %.2f, %.2f) dst: (%.2f, %.2f, %.2f)\n",		index, position.X, position.Y, position.Z,		destination.X, destination.Y, destination.Z);	fprintf(stdout, "\tdir: (%.2f, %.2f, %.2f) speed: %.2f\n",		direction.X, direction.Y, direction.Z, speed);	fprintf(stdout, "\tArrival: %.2f, Update: %.2f, Transition: %.2f\n",		time_arrival, time_update, time_transition);	for(i = 0; i < NODES; i++) {		m = &neighbor[i];		fprintf(stdout, "\tNeighbor: %d (%x), Reachable: %d, Transition Time: %.2f\n",			m->index, (int) m, m->reachable, m->time_transition);	}}/* ======================================================================   Dijkstra's Shortest Path Algoritm   ====================================================================== */void dumpall(){	u_int32_t i;	fprintf(stdout, "\nTime: %.2f\n", TIME);	for(i = 0; i < NODES; i++) {		NodeList[i].Dump();	}}voidComputeW(){	u_int32_t i, j;	u_int32_t *W = D2;	memset(W, '\xff', sizeof(int) * NODES * NODES);	for(i = 0; i < NODES; i++) {		for(j = i; j < NODES; j++) {			Neighbor *m = &NodeList[i].neighbor[j];			if(i == j)				W[i*NODES + j] = W[j*NODES + i] = 0;			else				W[i*NODES + j] = W[j*NODES + i] = m->reachable ? 1 : RTG_INFINITY;			}	}}voidfloyd_warshall(){	u_int32_t i, j, k;	ComputeW();	// the connectivity matrix	for(i = 0; i < NODES; i++) {		for(j = 0; j < NODES; j++) {			for(k = 0; k < NODES; k++) {				D2[j*NODES + k] = min(D2[j*NODES + k], D2[j*NODES + i] + D2[i*NODES + k]);			}		}	}#ifdef SANITY_CHECKS	for(i = 0; i < NODES; i++)		for(j = 0; j < NODES; j++) {			assert(D2[i*NODES + j] == D2[j*NODES + i]);			assert(D2[i*NODES + j] <= RTG_INFINITY);		}#endif}/* *  Write the actual GOD entries to a TCL script. */voidshow_diffs(){	u_int32_t i, j;	for(i = 0; i < NODES; i++) {		for(j = i + 1; j < NODES; j++) {			if(D1[i*NODES + j] != D2[i*NODES + j]) {				if(D2[i*NODES + j] == RTG_INFINITY)					DestUnreachableCount++;                                if(TIME > 0.0) {                                        RouteChangeCount++;                                        NodeList[i].route_changes++;                                        NodeList[j].route_changes++;                                }				if(TIME == 0.0) {					fprintf(stdout, GOD_FORMAT2,						i, j, D2[i*NODES + j]);#ifdef SHOW_SYMMETRIC_PAIRS					fprintf(stdout, GOD_FORMAT2,						j, i, D2[j*NODES + i]);#endif				}				else {					fprintf(stdout, GOD_FORMAT, 						TIME, i, j, D2[i*NODES + j]);#ifdef SHOW_SYMMETRIC_PAIRS					fprintf(stdout, GOD_FORMAT, 						TIME, j, i, D2[j*NODES + i]);#endif				}			}		}	}	memcpy(D1, D2, sizeof(int) * NODES * NODES);}voidshow_routes(){	u_int32_t i, j;	fprintf(stdout, "#\n# TIME: %.12f\n#\n", TIME);	for(i = 0; i < NODES; i++) {		fprintf(stdout, "# %2d) ", i);		for(j = 0; j < NODES; j++)			fprintf(stdout, "%3d ", D2[i*NODES + j] & 0xff);		fprintf(stdout, "\n");	}	fprintf(stdout, "#\n");}voidshow_counters(){	u_int32_t i;	fprintf(stdout, "#\n# Destination Unreachables: %d\n#\n",		DestUnreachableCount);	fprintf(stdout, "# Route Changes: %d\n#\n", RouteChangeCount);        fprintf(stdout, "# Link Changes: %d\n#\n", LinkChangeCount);        fprintf(stdout, "# Node | Route Changes | Link Changes\n");	for(i = 0; i < NODES; i++)		fprintf(stdout, "# %4d |          %4d |         %4d\n",                        i, NodeList[i].route_changes,                        NodeList[i].link_changes);	fprintf(stdout, "#\n");}