开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > Windows编程 > 通讯编程 > 查看源码
setdest.cc
资源名称：ns-allinone-2.33.tar.gz　[点击查看]
上传用户：rrhhcc
上传日期：2015-12-11
资源大小：54129k
文件大小：25k
源码类别：

通讯编程

开发平台：

Visual C++

setdest.cc：源码内容
  1. /*
  2.  *
  3.  * The original code of setdest was included in ns-2.1b8a.
  4.  * This file is the modified version by J. Yoon <jkyoon@eecs.umich.edu>,
  5.  * Department of EECS, University of Michigan, Ann Arbor. 
  6.  *
  7.  * (1) Input parameters
  8.  * <Original version>
  9.  * => -M maximum speed (minimum speed is zero as a default)
  10.  * => -p pause time (constant) 
  11.  * => -n number of nodes
  12.  * => -x x dimension of space
  13.  * => -y y dimension of space
  14.  *
  15.  * <Modified version>
  16.  * => -s speed type (uniform, normal)
  17.  * => -m minimum speed > 0 
  18.  * => -M maximum speed
  19.  * => -P pause type (constant, uniform)
  20.  * => -p pause time (a median if uniform is chosen)
  21.  * => -n number of nodes
  22.  * => -x x dimension of space
  23.  * => -y y dimension of space
  24.  *
  25.  * (2) In case of modified version, the steady-state speed distribution is applied to 
  26.  * the first trip to eliminate any speed decay. If pause is not zero, the first 
  27.  * trip could be either a move or a pause depending on the probabilty that the 
  28.  * first trip is a pause. After the first trip regardless of whether it is 
  29.  * a move or a pause, all subsequent speeds are determined from the given speed 
  30.  * distribution (e.g., uniform or normal).
  31.  *
  32.  * (3) Refer to and use scenario-generating scripts (make-scen.csh for original version, 
  33.  * make-scen-steadystate.csh for modified version).
  34.  *
  35.  *
  36.  */
  40. extern "C" {
  41. #include <assert.h>
  42. #include <fcntl.h>
  43. #include <math.h>
  44. #include <stdio.h>
  45. #include <stdlib.h>
  46. #include <string.h>
  47. #include <sys/time.h>
  48. #include <sys/types.h>
  49. #include <sys/uio.h>
  50. #include <unistd.h>
  51. #if !defined(sun)
  52. #include <err.h>
  53. #endif
  54. };
  55. #include "../../../tools/rng.h"
  57. #include "setdest.h"
  60. // #define DEBUG
  61. #define SANITY_CHECKS
  62. //#define SHOW_SYMMETRIC_PAIRS
  65. #define GOD_FORMAT "$ns_ at %.12f "$god_ set-dist %d %d %d"n"
  66. #define GOD_FORMAT2 "$god_ set-dist %d %d %dn"
  67. #define NODE_FORMAT "$ns_ at %.12f "$node_(%d) setdest %.12f %.12f %.12f"n"
  68. #define NODE_FORMAT2 "$node_(%d) setdest %.12f %.12f %.12fn"
  69. #define NODE_FORMAT3 "$node_(%d) set %c_ %.12fn"
  71. #define RTG_INFINITY 0x00ffffff
  72. #ifndef min
  73. #define min(x,y) ((x) < (y) ? (x) : (y))
  74. #define max(x,y) ((x) > (y) ? (x) : (y))
  75. #endif
  76. #define ROUND_ERROR 1e-9
  77. #ifndef PI
  78. #define  PI   3.1415926
  79. #endif
  82. static int count = 0;
  84. /* ======================================================================
  85.    Function Prototypes
  86.    ====================================================================== */
  87. void usage(char**);
  88. void init(void);
  89. double uniform(void);
  91. void dumpall(void);
  92. void ComputeW(void);
  93. void floyd_warshall(void);
  95. void show_diffs(void);
  96. void show_routes(void);
  97. void show_counters(void);
  100. /* ======================================================================
  101.    Global Variables
  102.    ====================================================================== */
  103. const double RANGE = 250.0; // transmitter range in meters
  104. double TIME = 0.0; // my clock;
  105. double MAXTIME = 0.0; // duration of simulation
  107. double MAXX = 0.0; // width of space
  108. double MAXY = 0.0; // height of space
  109. double PAUSE = 0.0; // pause time
  110. double MAXSPEED = 0.0; // max speed
  111. double MINSPEED = 0.0; // min speed 
  112. double SS_AVGSPEED = 0.0; // steady-state avg speed 
  113. double  KAPPA = 0.0; // normalizing constant 
  114. double  MEAN = 0.0; // mean for normal speed
  115. double  SIGMA = 0.0; // std for normal speed
  116. double  EXP_1_V = 0.0; // expactation of 1/V
  117. double  EXP_R = 0.0; // expectation of travel distance R
  118. double  PDFMAX = 0.0; // max of pdf for rejection technique
  119. u_int32_t SPEEDTYPE = 1; // speed type (default = uniform)
  120. u_int32_t PAUSETYPE = 1; // pause type (default = constant)
  121. u_int32_t VERSION = 1; // setdest version (default = original by CMU) 
  122. u_int32_t NODES = 0; // number of nodes
  123. u_int32_t RouteChangeCount = 0;
  124. u_int32_t LinkChangeCount = 0;
  125. u_int32_t DestUnreachableCount = 0;
  128. Node *NodeList = 0;
  129. u_int32_t *D1 = 0;
  130. u_int32_t *D2 = 0;
  133. /* ======================================================================
  134.    Random Number Generation
  135.    ====================================================================== */
  136. #define M 2147483647L
  137. #define INVERSE_M ((double)4.656612875e-10)
  139. char random_state[32];
  140. RNG *rng;
  142. double
  143. uniform()
  144. {
  145.         count++;
  146.         return rng->uniform_double() ;
  147. }
  150. /* ======================================================================
  151.    Misc Functions...
  152.    ====================================================================== */
  155. /* compute the expectation of travel distance E[R] in a rectangle */
  156. void
  157. compute_EXP_R()
  158. {
  159. #define csc(x)  (1.0/sin(x)) // csc function
  160. #define sec(x) (1.0/cos(x)) // sec function
  161. #define sin2(x) (sin(x)*sin(x)) // sin^2
  162. #define sin3(x) (sin2(x)*sin(x)) // sin^3
  163. #define cos2(x) (cos(x)*cos(x)) // cos^2
  164. #define cos3(x) (cos2(x)*cos(x)) // cos^3
  166. double x = MAXX, y = MAXY; // max x and max y
  167. double x2 = x*x, x3 = x*x*x; // x^2 and x^3
  168. double y2 = y*y, y3 = y*y*y; // y^2 and y^3
  170. double term1 = sin(atan2(y,x)) / 2.0 / cos2(atan2(y,x));
  171. double term2 = 0.5 * log( sec(atan2(y,x)) + y/x );
  172. double term3 = -1.0 * x3 / y2 / 60.0 / cos3(atan2(y,x)) + 1.0/60.0 * x3 / y2;
  173. double term4 = (term1 + term2) * x2 / 12.0 / y + term3;
  175. double term5 = -1.0 * cos(atan2(y,x)) / 2.0 / sin2(atan2(y,x));
  176. double term6 = 0.5 * log( csc(atan2(y,x)) - x/y );
  177. double term7 = -1.0 * y3 / x2 / 60.0 / sin3(atan2(y,x)) + 1.0/60.0 * y3 / x2;
  178. double term8 = -1.0 * (term5 + term6) * y2 / 12.0 / x + term7;
  180. EXP_R = (4 * (term4 + term8));  // E[R]
  181. }
  183. void
  184. usage(char **argv)
  185. {
  186. fprintf(stderr, "nusage:n");
  187. fprintf(stderr,
  188. "n<original 1999 CMU version (version 1)>n %st-v <1> -n <nodes> -p <pause time> -M <max speed>n",
  189. argv[0]);
  190. fprintf(stderr,
  191. "tt-t <simulation time> -x <max X> -y <max Y>n");
  192. fprintf(stderr,
  193. "nORn<modified 2003 U.Michigan version (version 2)>n %st-v <2> -n <nodes> -s <speed type> -m <min speed> -M <max speed>n",
  194. argv[0]);
  195. fprintf(stderr,
  196. "tt-t <simulation time> -P <pause type> -p <pause time> -x <max X> -y <max Y>n");
  197. fprintf(stderr,
  198. "tt(Refer to the script files make-scen.csh and make-scen-steadystate.csh for detail.) nn");
  199. }
  201. void
  202. init()
  203. {
  204. /*
  205.  * Initialized the Random Number Generation
  206.  */
  207. /* 
  208. This part of init() is commented out and is replaced by more
  209. portable RNG (random number generator class of ns) functions.
  211. struct timeval tp;
  212. int fd, seed, bytes;
  214. if((fd = open("/dev/random", O_RDONLY)) < 0) {
  215. perror("open");
  216. exit(1);
  217. }
  218. if((bytes = read(fd, random_state, sizeof(random_state))) < 0) {
  219. perror("read");
  220. exit(1);
  221. }
  222. close(fd);
  224. fprintf(stderr, "*** read %d bytes from /dev/randomn", bytes);
  226. if(bytes != sizeof(random_state)) {
  227.   fprintf(stderr,"Not enough randomness. Reading `.rand_state'n");
  228.   if((fd = open(".rand_state", O_RDONLY)) < 0) {
  229.     perror("open .rand_state");
  230.     exit(1);
  231.   }
  232.   if((bytes = read(fd, random_state, sizeof(random_state))) < 0) {
  233.     perror("reading .rand_state");
  234.     exit(1);
  235.   }
  236.   close(fd);
  237. }
  239.          if(gettimeofday(&tp, 0) < 0) {
  240. perror("gettimeofday");
  241. exit(1);
  242. }
  243. seed = (tp.tv_sec  >> 12 ) ^ tp.tv_usec;
  244.         (void) initstate(seed, random_state, bytes & 0xf8);*/
  246. /*
  247.  * Allocate memory for globals
  248.  */
  249. NodeList = new Node[NODES];
  250. if(NodeList == 0) {
  251. perror("new");
  252. exit(1);
  253. }
  255. D1 = new u_int32_t[NODES * NODES];
  256. if(D1 == 0) {
  257. perror("new");
  258. exit(1);
  259. }
  260. memset(D1, 'xff', sizeof(u_int32_t) * NODES * NODES);
  262. D2 = new u_int32_t[NODES * NODES];
  263. if(D2 == 0) {
  264. perror("new");
  265. exit(1);
  266. }
  267. memset(D2, 'xff', sizeof(u_int32_t) * NODES * NODES);
  268. }
  270. extern "C" char *optarg;
  272. int
  273. main(int argc, char **argv)
  274. {
  275. char ch;
  277. while ((ch = getopt(argc, argv, "v:n:s:m:M:t:P:p:x:y:i:o:")) != EOF) {       
  279. switch (ch) { 
  281. case 'v':
  282.   VERSION = atoi(optarg);
  283.   break;
  285. case 'n':
  286.   NODES = atoi(optarg);
  287.   break;
  289. case 's':
  290. SPEEDTYPE = atoi(optarg);
  291. break;
  293. case 'm':
  294. MINSPEED = atof(optarg);
  295. break;
  297. case 'M':
  298. MAXSPEED = atof(optarg);
  299. break;
  301. case 't':
  302. MAXTIME = atof(optarg);
  303. break;
  305. case 'P':
  306. PAUSETYPE = atoi(optarg);
  307. break;
  309. case 'p':
  310. PAUSE = atof(optarg);
  311. break;
  313. case 'x':
  314. MAXX = atof(optarg);
  315. break;
  317. case 'y':
  318. MAXY = atof(optarg);
  319. break;
  321. default:
  322. usage(argv);
  323. exit(1);
  324. }
  325. }
  327. if(MAXX == 0.0 || MAXY == 0.0 || NODES == 0 || MAXTIME == 0.0) {
  328. usage(argv);
  329. exit(1);
  330. }
  332. /* specify the version */
  333. if (VERSION != 1 && VERSION != 2) {
  334.   printf("Please specify the setdest version you want to use. For original 1999 CMU version use 1; For modified 2003 U.Michigan version use 2n");
  335.   exit(1);
  336. }
  338. if (VERSION == 2 && MINSPEED <= 0) {
  339.   usage(argv);
  340.   exit(1);
  341. } else if (VERSION == 1 && MINSPEED > 0) {
  342.   usage(argv);
  343.   exit(1);
  344. }
  347. // The more portable solution for random number generation
  348. rng = new RNG;
  349. rng->set_seed(RNG::HEURISTIC_SEED_SOURCE); 
  353. /****************************************************************************************
  354.    * Steady-state avg speed and distribution depending on the initial distirbutions given
  355.  ****************************************************************************************/
  357. /* original setdest */
  358. if (VERSION == 1) {
  359. fprintf(stdout, "#n# nodes: %d, pause: %.2f, max speed: %.2f, max x: %.2f, max y: %.2fn#n",
  360. NODES, PAUSE, MAXSPEED, MAXX, MAXY);
  361. }
  363. /* modified version */
  364. else if (VERSION == 2) {
  365. /* compute the expectation of travel distance in a rectangle */
  366. compute_EXP_R();
  368. /* uniform speed from min to max */
  369. if (SPEEDTYPE == 1) {
  370. EXP_1_V = log(MAXSPEED/MINSPEED) / (MAXSPEED - MINSPEED); // E[1/V]
  371. SS_AVGSPEED = EXP_R / (EXP_1_V*EXP_R + PAUSE); // steady-state average speed
  372. PDFMAX = 1/MINSPEED*EXP_R / (EXP_1_V*EXP_R + PAUSE) / (MAXSPEED-MINSPEED); // max of pdf for rejection technique
  373. }
  375. /* normal speed clipped from min to max */
  376. else if (SPEEDTYPE == 2) {
  377. int bin_no = 10000; // the number of bins for summation
  378. double delta = (MAXSPEED - MINSPEED)/bin_no;  // width of each bin 
  379. int i;
  380. double acc_k, acc_e, square, temp_v;
  382. MEAN = (MAXSPEED + MINSPEED)/2.0; // means for normal dist.
  383. SIGMA = (MAXSPEED - MINSPEED)/4.0; // std for normal dist.
  384. /* computing a normalizing constant KAPPA, E[1/V], and pdf max */
  385. KAPPA = 0.0;
  386. EXP_1_V = 0.0;
  387. PDFMAX = 0.0;
  389. /* numerical integrals */
  390. for (i=0; i<bin_no; ++i) {
  391. temp_v = MINSPEED + i*delta; // ith v from min speed
  392. square = (temp_v - MEAN)*(temp_v - MEAN)/SIGMA/SIGMA;
  394. acc_k = 1.0/sqrt(2.0*PI*SIGMA*SIGMA)*exp(-0.5*square);
  395. KAPPA += (acc_k*delta); // summing up the area of rectangle
  397. acc_e = 1.0/temp_v/sqrt(2.0*PI*SIGMA*SIGMA)*exp(-0.5*square);
  398. EXP_1_V += (acc_e*delta); // summing up for the denominator of pdf
  400. /* find a max of pdf */
  401. if (PDFMAX < acc_e) PDFMAX = acc_e;
  402. }
  403. EXP_1_V /= KAPPA; // normalizing
  404. SS_AVGSPEED = EXP_R / (EXP_1_V*EXP_R + PAUSE); // steady-state average speed
  405. PDFMAX = EXP_R*PDFMAX/KAPPA / (EXP_1_V*EXP_R + PAUSE); // max of pdf for rejection technique
  406. }
  407. /* other types of speed for future use */
  408. else
  409. ;
  411. fprintf(stdout, "#n# nodes: %d, speed type: %d, min speed: %.2f, max speed: %.2fn# avg speed: %.2f, pause type: %d, pause: %.2f, max x: %.2f, max y: %.2fn#n",
  412. NODES , SPEEDTYPE, MINSPEED, MAXSPEED, SS_AVGSPEED, PAUSETYPE, PAUSE, MAXX, MAXY);
  413. } 
  416. init();
  418. while(TIME <= MAXTIME) {
  419. double nexttime = 0.0;
  420. u_int32_t i;
  422. for(i = 0; i < NODES; i++) {
  423. NodeList[i].Update();
  424. }
  426. for(i = 0; i < NODES; i++) {
  427. NodeList[i].UpdateNeighbors();
  428. }
  430. for(i = 0; i < NODES; i++) {
  431. Node *n = &NodeList[i];
  433. if(n->time_transition > 0.0) {
  434. if(nexttime == 0.0)
  435. nexttime = n->time_transition;
  436. else
  437. nexttime = min(nexttime, n->time_transition);
  438. }
  440. if(n->time_arrival > 0.0) {
  441. if(nexttime == 0.0)
  442. nexttime = n->time_arrival;
  443. else
  444. nexttime = min(nexttime, n->time_arrival);
  445. }
  447. }
  449. floyd_warshall();
  451. #ifdef DEBUG
  452. show_routes();
  453. #endif
  455.   show_diffs();
  457. #ifdef DEBUG
  458. dumpall();
  459. #endif
  461. assert(nexttime > TIME + ROUND_ERROR);
  462. TIME = nexttime;
  463. }
  465. show_counters();
  467. int of;
  468. if ((of = open(".rand_state",O_WRONLY | O_TRUNC | O_CREAT, 0777)) < 0) {
  469.   fprintf(stderr, "open rand staten");
  470.   exit(-1);
  471.   }
  472. for (unsigned int i = 0; i < sizeof(random_state); i++)
  473.           random_state[i] = 0xff & (int) (uniform() * 256);
  474. if (write(of,random_state, sizeof(random_state)) < 0) {
  475.   fprintf(stderr, "writing rand staten");
  476.   exit(-1);
  477.   }
  478. close(of);
  481. }
  484. /* ======================================================================
  485.    Node Class Functions
  486.    ====================================================================== */
  487. u_int32_t Node::NodeIndex = 0;
  489. Node::Node()
  490. {
  491. u_int32_t i;
  493. index = NodeIndex++;
  495. //if(index == 0)
  496. // return;
  498. route_changes = 0;
  499.     link_changes = 0;
  503.     /*******************************************************************************
  504.  * Determine if the first trip is a pause or a move with the steady-state pdf
  505.      *******************************************************************************/
  507. /* original version */
  508. if (VERSION == 1) {
  509. time_arrival = TIME + PAUSE; // constant pause
  510. }
  512. /* modified version */ 
  513. else if (VERSION == 2) {
  514. /* probability that the first trip would be a pause */
  515. double prob_pause = PAUSE / (EXP_1_V*EXP_R + PAUSE);
  517. /* the first trip is a pause */
  518. if (prob_pause > uniform()) {
  519. /* constant pause */
  520. if (PAUSETYPE == 1) {
  521. time_arrival = TIME + PAUSE; // constant pause
  522. }
  523. /* uniform pause */
  524. else if (PAUSETYPE == 2) {
  525. time_arrival = TIME + 2*PAUSE*uniform(); // uniform pause [0, 2*PAUSE]
  526. }
  528. first_trip = 0; // indicating the first trip is a pause 
  529. }
  530. /* the first trip is a move based on the steady-state pdf */
  531. else {
  532. time_arrival = TIME;
  533. first_trip = 1; // indicating the first trip is a move 
  534. }
  535. }
  538. time_update = TIME;
  539. time_transition = 0.0;
  541. position.X = position.Y = position.Z = 0.0;
  542. destination.X = destination.Y = destination.Z = 0.0;
  543. direction.X = direction.Y = direction.Z = 0.0;
  544. speed = 0.0;
  546. RandomPosition();
  548. fprintf(stdout, NODE_FORMAT3, index, 'X', position.X);
  549. fprintf(stdout, NODE_FORMAT3, index, 'Y', position.Y);
  550. fprintf(stdout, NODE_FORMAT3, index, 'Z', position.Z);
  552. neighbor = new Neighbor[NODES];
  553. if(neighbor == 0) {
  554. perror("new");
  555. exit(1);
  556. }
  558. for(i = 0; i < NODES; i++) {
  559. neighbor[i].index = i;
  560. neighbor[i].reachable = (index == i) ? 1 : 0;
  561. neighbor[i].time_transition = 0.0;
  562. }
  563. }
  567. void
  568. Node::RandomPosition()
  569. {
  570.     position.X = uniform() * MAXX;
  571.     position.Y = uniform() * MAXY;
  572. position.Z = 0.0;
  573. }
  576. void
  577. Node::RandomDestination()
  578. {
  579.     destination.X = uniform() * MAXX;
  580.     destination.Y = uniform() * MAXY;
  581. destination.Z = 0.0;
  582. assert(destination != position);
  583. }
  586. /****************************************************************************************** 
  587.  * Speeds are chosen based on the given type and distribution
  588.  ******************************************************************************************/
  589. void
  590. Node::RandomSpeed()
  591. {
  592. /* original version */
  593. if (VERSION == 1) {
  594.         speed = uniform() * MAXSPEED;
  595. assert(speed != 0.0);
  596. }
  598. /* modified version */
  599. else if (VERSION == 2) {
  600. /* uniform speed */
  601. if (SPEEDTYPE == 1) {
  602. /* using steady-state distribution for the first trip */
  603. if (first_trip == 1) {
  604. /* pick a speed by rejection technique */
  605. double temp_v, temp_fv;
  607. do {
  608.            temp_v = uniform() * (MAXSPEED - MINSPEED) + MINSPEED;
  609. temp_fv = uniform() * PDFMAX;
  610. } while (temp_fv > 1/temp_v*EXP_R / (EXP_1_V*EXP_R + PAUSE) / (MAXSPEED-MINSPEED));
  611.  
  612. speed = temp_v;
  613. first_trip = 0; // reset first_trip flag 
  614. }
  615. /* using the original distribution from the second trip on */
  616. else {
  617.          speed = uniform() * (MAXSPEED - MINSPEED) + MINSPEED;
  618. assert(speed != 0.0);
  619. }
  620. }
  621. /* normal speed */
  622. else if (SPEEDTYPE == 2) {
  623. /* using steady-state distribution for the first trip */
  624. if (first_trip == 1) {
  625. double temp_v, temp_fv, square, fv;
  627. /* rejection technique */
  628. do {
  629.            temp_v = uniform() * (MAXSPEED - MINSPEED) + MINSPEED;
  630. temp_fv = uniform() * PDFMAX;
  631. square = (temp_v - MEAN)*(temp_v - MEAN)/SIGMA/SIGMA;
  632. fv = 1/KAPPA/sqrt(2.0*PI*SIGMA*SIGMA) * exp(-0.5*square);
  633. } while (temp_fv > 1.0/temp_v*fv*EXP_R / (EXP_1_V*EXP_R + PAUSE));
  634.  
  635. speed = temp_v;
  636. first_trip = 0;
  637. }
  638. /* using the original distribution from the second trip on */
  639. else {
  640. double temp_v, temp_fv, square;
  641. double max_normal = 1.0/KAPPA/sqrt(2.0*PI*SIGMA*SIGMA); // max of normal distribution
  643. /* rejection technique */
  644. do {
  645.           temp_v = uniform() * (MAXSPEED - MINSPEED) + MINSPEED;
  646. temp_fv = uniform() * max_normal;
  647. square = (temp_v - MEAN)*(temp_v - MEAN)/SIGMA/SIGMA;
  648. } while (temp_fv > max_normal * exp(-0.5*square));
  649.  
  650. speed = temp_v;
  651. assert(speed != 0.0);
  652. }
  653. }
  654. /* other types of speed for future use */
  655. else
  656. ;
  657. }
  658. }
  661. void
  662. Node::Update()
  663. {
  664. position += (speed * (TIME - time_update)) * direction;
  666. if(TIME == time_arrival) {
  667. vector v;
  669. if(speed == 0.0 || PAUSE == 0.0) {
  671.             RandomDestination();
  672. RandomSpeed();
  674. v = destination - position;
  675. direction = v / v.length();
  676. time_arrival = TIME + v.length() / speed;
  677. }
  678. else {
  680. destination = position;
  681. speed = 0.0;
  683. /* original version */
  684. if (VERSION == 1) {
  685. time_arrival = TIME + PAUSE;
  686. }
  687. /* modified version */
  688. else if (VERSION == 2) {
  689. /* constant pause */
  690. if (PAUSETYPE == 1) {
  691. time_arrival = TIME + PAUSE;
  692. }
  693. /* uniform pause */
  694. else if (PAUSETYPE == 2) {
  695. time_arrival = TIME + 2*PAUSE*uniform();
  696. }
  697. }
  698. }
  700. fprintf(stdout, NODE_FORMAT,
  701. TIME, index, destination.X, destination.Y, speed);
  703. }
  705. time_update = TIME;
  706. time_transition = 0.0;
  707. }
  710. void
  711. Node::UpdateNeighbors()
  712. {
  713. static Node *n2;
  714. static Neighbor *m1, *m2;
  715. static vector D, B, v1, v2;
  716. static double a, b, c, t1, t2, Q;
  717. static u_int32_t i, reachable;
  719. v1 = speed * direction;
  721. /*
  722.  *  Only need to go from INDEX --> N for each one since links
  723.  *  are symmetric.
  724.  */
  725. for(i = index+1; i < NODES; i++) {
  727. m1 = &neighbor[i];
  728. n2 = &NodeList[i];
  729. m2 = &n2->neighbor[index];
  731. assert(i == m1->index);
  732. assert(m1->index == n2->index);
  733. assert(index == m2->index);
  734.                 assert(m1->reachable == m2->reachable);
  736.                 reachable = m1->reachable;
  738. /* ==================================================
  739.    Determine Reachability
  740.    ================================================== */
  741. { vector d = position - n2->position;
  743. if(d.length() < RANGE) {
  744. #ifdef SANITY_CHECKS
  745. if(TIME > 0.0 && m1->reachable == 0)
  746. assert(RANGE - d.length() < ROUND_ERROR);
  747. #endif
  748. m1->reachable = m2->reachable = 1;
  749. }
  750. // Boundary condition handled below.
  751. else {
  752. #ifdef SANITY_CHECKS
  753. if(TIME > 0.0 && m1->reachable == 1)
  754. assert(d.length() - RANGE < ROUND_ERROR);
  755. #endif
  756. m1->reachable = m2->reachable = 0;
  757. }
  758. #ifdef DEBUG
  759.                         fprintf(stdout, "# %.6f (%d, %d) %.2fmn",
  760.                                 TIME, index, m1->index, d.length());
  761. #endif
  762. }
  764. /* ==================================================
  765.    Determine Next Event Time
  766.    ================================================== */
  767. v2 = n2->speed * n2->direction;
  769. D = v2 - v1;
  770. B = n2->position - position;
  772. a = (D.X * D.X) + (D.Y * D.Y) + (D.Z * D.Z);
  773. b = 2 * ((D.X * B.X) + (D.Y * B.Y) + (D.Z * B.Z));
  774. c = (B.X * B.X) + (B.Y * B.Y) + (B.Z * B.Z) - (RANGE * RANGE);
  776. if(a == 0.0) {
  777. /*
  778.  *  No Finite Solution
  779.  */
  780. m1->time_transition= 0.0;
  781. m2->time_transition= 0.0;
  782. goto  next;
  783. }
  785. Q = b * b - 4 * a * c;
  786. if(Q < 0.0) {
  787. /*
  788.  *  No real roots.
  789.  */
  790. m1->time_transition = 0.0;
  791. m2->time_transition = 0.0;
  792. goto next;
  793. }
  794. Q = sqrt(Q);
  796. t1 = (-b + Q) / (2 * a);
  797. t2 = (-b - Q) / (2 * a);
  799. // Stupid Rounding/Boundary Cases
  800. if(t1 > 0.0 && t1 < ROUND_ERROR) t1 = 0.0;
  801. if(t1 < 0.0 && -t1 < ROUND_ERROR) t1 = 0.0;
  802. if(t2 > 0.0 && t2 < ROUND_ERROR) t2 = 0.0;
  803. if(t2 < 0.0 && -t2 < ROUND_ERROR) t2 = 0.0;
  805. if(t1 < 0.0 && t2 < 0.0) {
  806. /*
  807.  *  No "future" time solution.
  808.  */
  809. m1->time_transition = 0.0;
  810. m2->time_transition = 0.0;
  811. goto next;
  812. }
  814. /*
  815.  * Boundary conditions.
  816.  */
  817. if((t1 == 0.0 && t2 > 0.0) || (t2 == 0.0 && t1 > 0.0)) {
  818. m1->reachable = m2->reachable = 1;
  819. m1->time_transition = m2->time_transition = TIME + max(t1, t2);
  820. }
  821. else if((t1 == 0.0 && t2 < 0.0) || (t2 == 0.0 && t1 < 0.0)) {
  822. m1->reachable = m2->reachable = 0;
  823. m1->time_transition = m2->time_transition = 0.0;
  824. }
  826. /*
  827.  * Non-boundary conditions.
  828.  */
  829. else if(t1 > 0.0 && t2 > 0.0) {
  830. m1->time_transition = TIME + min(t1, t2);
  831. m2->time_transition = TIME + min(t1, t2);
  832. }
  833. else if(t1 > 0.0) {
  834. m1->time_transition = TIME + t1;
  835. m2->time_transition = TIME + t1;
  836. }
  837. else {
  838. m1->time_transition = TIME + t2;
  839. m2->time_transition = TIME + t2;
  840. }
  842. /* ==================================================
  843.    Update the transition times for both NODEs.
  844.    ================================================== */
  845. if(time_transition == 0.0 || (m1->time_transition &&
  846.    time_transition > m1->time_transition)) {
  847. time_transition = m1->time_transition;
  848. }
  849. if(n2->time_transition == 0.0 || (m2->time_transition &&
  850.    n2->time_transition > m2->time_transition)) {
  851. n2->time_transition = m2->time_transition;
  852. }
  853.         next:
  854.                 if(reachable != m1->reachable && TIME > 0.0) {
  855.                         LinkChangeCount++;
  856.                         link_changes++;
  857.                         n2->link_changes++;
  858.                 }
  859. }
  860. }
  862. void
  863. Node::Dump()
  864. {
  865. Neighbor *m;
  866. u_int32_t i;
  868. fprintf(stdout,
  869. "Node: %dtpos: (%.2f, %.2f, %.2f) dst: (%.2f, %.2f, %.2f)n",
  870. index, position.X, position.Y, position.Z,
  871. destination.X, destination.Y, destination.Z);
  872. fprintf(stdout, "tdir: (%.2f, %.2f, %.2f) speed: %.2fn",
  873. direction.X, direction.Y, direction.Z, speed);
  874. fprintf(stdout, "tArrival: %.2f, Update: %.2f, Transition: %.2fn",
  875. time_arrival, time_update, time_transition);
  877. for(i = 0; i < NODES; i++) {
  878. m = &neighbor[i];
  879. fprintf(stdout, "tNeighbor: %d (%lx), Reachable: %d, Transition Time: %.2fn",
  880. m->index, (long) m, m->reachable, m->time_transition);
  881. }
  882. }
  885. /* ======================================================================
  886.    Dijkstra's Shortest Path Algoritm
  887.    ====================================================================== */
  888. void dumpall()
  889. {
  890. u_int32_t i;
  892. fprintf(stdout, "nTime: %.2fn", TIME);
  894. for(i = 0; i < NODES; i++) {
  895. NodeList[i].Dump();
  896. }
  897. }
  899. void
  900. ComputeW()
  901. {
  902. u_int32_t i, j;
  903. u_int32_t *W = D2;
  905. memset(W, 'xff', sizeof(int) * NODES * NODES);
  907. for(i = 0; i < NODES; i++) {
  908. for(j = i; j < NODES; j++) {
  909. Neighbor *m = &NodeList[i].neighbor[j];
  910. if(i == j)
  911. W[i*NODES + j] = W[j*NODES + i] = 0;
  912. else
  913. W[i*NODES + j] = W[j*NODES + i] = m->reachable ? 1 : RTG_INFINITY;
  914. }
  915. }
  916. }
  918. void
  919. floyd_warshall()
  920. {
  921. u_int32_t i, j, k;
  923. ComputeW(); // the connectivity matrix
  925. for(i = 0; i < NODES; i++) {
  926. for(j = 0; j < NODES; j++) {
  927. for(k = 0; k < NODES; k++) {
  928. D2[j*NODES + k] = min(D2[j*NODES + k], D2[j*NODES + i] + D2[i*NODES + k]);
  929. }
  930. }
  931. }
  933. #ifdef SANITY_CHECKS
  934. for(i = 0; i < NODES; i++)
  935. for(j = 0; j < NODES; j++) {
  936. assert(D2[i*NODES + j] == D2[j*NODES + i]);
  937. assert(D2[i*NODES + j] <= RTG_INFINITY);
  938. }
  939. #endif
  940. }
  943. /*
  944.  *  Write the actual GOD entries to a TCL script.
  945.  */
  946. void
  947. show_diffs()
  948. {
  949. u_int32_t i, j;
  951. for(i = 0; i < NODES; i++) {
  952. for(j = i + 1; j < NODES; j++) {
  953. if(D1[i*NODES + j] != D2[i*NODES + j]) {
  955. if(D2[i*NODES + j] == RTG_INFINITY)
  956. DestUnreachableCount++;
  958.                                 if(TIME > 0.0) {
  959.                                         RouteChangeCount++;
  960.                                         NodeList[i].route_changes++;
  961.                                         NodeList[j].route_changes++;
  962.                                 }
  964. if(TIME == 0.0) {
  965. fprintf(stdout, GOD_FORMAT2,
  966. i, j, D2[i*NODES + j]);
  967. #ifdef SHOW_SYMMETRIC_PAIRS
  968. fprintf(stdout, GOD_FORMAT2,
  969. j, i, D2[j*NODES + i]);
  970. #endif
  971. }
  972. else {
  973. fprintf(stdout, GOD_FORMAT, 
  974. TIME, i, j, D2[i*NODES + j]);
  975. #ifdef SHOW_SYMMETRIC_PAIRS
  976. fprintf(stdout, GOD_FORMAT, 
  977. TIME, j, i, D2[j*NODES + i]);
  978. #endif
  979. }
  980. }
  981. }
  982. }
  984. memcpy(D1, D2, sizeof(int) * NODES * NODES);
  985. }
  988. void
  989. show_routes()
  990. {
  991. u_int32_t i, j;
  993. fprintf(stdout, "#n# TIME: %.12fn#n", TIME);
  994. for(i = 0; i < NODES; i++) {
  995. fprintf(stdout, "# %2d) ", i);
  996. for(j = 0; j < NODES; j++)
  997. fprintf(stdout, "%3d ", D2[i*NODES + j] & 0xff);
  998. fprintf(stdout, "n");
  999. }
  1000. fprintf(stdout, "#n");
  1001. }
  1003. void
  1004. show_counters()
  1005. {
  1006. u_int32_t i;
  1008. fprintf(stdout, "#n# Destination Unreachables: %dn#n",
  1009. DestUnreachableCount);
  1010. fprintf(stdout, "# Route Changes: %dn#n", RouteChangeCount);
  1011.         fprintf(stdout, "# Link Changes: %dn#n", LinkChangeCount);
  1012.         fprintf(stdout, "# Node | Route Changes | Link Changesn");
  1013. for(i = 0; i < NODES; i++)
  1014. fprintf(stdout, "# %4d |          %4d |         %4dn",
  1015.                         i, NodeList[i].route_changes,
  1016.                         NodeList[i].link_changes);
  1017. fprintf(stdout, "#n");
  1018. }