#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <getopt.h>
#include <sys/types.h>
#include <unistd.h>

#define PI 3.14159265

char *program_name;

void rdmob ( int nid, double theta, double speed, double etime );

typedef enum {false, true} bool;

void print_usage (FILE *stream, int exit_code);

int nodes;
double minspeed;
double maxspeed;
double x;
double y;
double simtime;

double pause_time;
double theta;  // angle
double etime;   // travel time (epoch time)
double speed;  // speed of travel 
double avgspeed;
double avgdist;

struct npos {
   double x;
   double y;
   double z;
   double t;
};

struct npos *positions;

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

   bool n_given = false;   
   bool m_given = false;
   bool M_given = false;
   bool x_given = false;
   bool y_given = false;
   bool p_given = false;
   bool t_given = false;

   double ptime;

   char *serror = "";

   srand48(time(NULL)/7*(long)getpid());
   program_name = argv[0];

   int next_option;

  // A string listing valid short options letters.
  const char* const short_options = "hn:x:y:m:M:p:t:v";

  // An array describing valid long options.
  const struct option long_options[] = {
    { "help",     0, NULL, 'h' },
    { "nodes",    1, NULL, 'n' },
    { "x",        1, NULL, 'x' },
    { "y",        1, NULL, 'y' },
    { "minspeed", 1, NULL, 'm' },
    { "maxspeed", 1, NULL, 'M' },
    { "pause",    1, NULL, 'p' },
    { "simtime",  1, NULL, 't' },
    { NULL,       0, NULL, 0   }   /* Required at end of array.  */
  };

  do {
     next_option = getopt_long (argc, argv, short_options, long_options, NULL);

     switch (next_option)
     {
       case 'h':
          print_usage (stdout, 0);

       case 'n':   // -n or --nodes
          nodes = atoi(optarg);
          n_given = true;
          break;

       case 'x':   // -x
          x = atof(optarg);
          x_given = true;
          break;

       case 'y':   // -y
          y = atof(optarg);
          y_given = true;
          break;

       case 'm':   // -m or --minspeed
          minspeed = atof(optarg);
          m_given = true;
          break;

       case 'M':   // -M or --maxspeed
          maxspeed = atof(optarg);
          M_given = true;
          break;

       case 'p':   // -p or --pause
          pause_time = atof(optarg);
          p_given = true;
          break;

       case 't':   // -p or --pause
          simtime = atof(optarg);
          t_given = true;
          break;

       case '?':   // The user specified an invalid option.
          print_usage (stderr, 1);

       case -1:    // Done with options.
          break;

       default:    // Something else: unexpected.
          abort ();
     }
   } while ( next_option != -1 );

   if ( !n_given ) {
      serror = "error: number of nodes missing";
   } else if ( nodes <= 0 ) {
      serror = "error: number of nodes not positive";
   } else if ( !x_given ) {
      serror = "error: width (x) missing";
   } else if ( x <= 0 ) {
      serror = "error: x not positive";
   } else if ( !y_given ) {
      serror = "error: length (y) missing";
   } else if ( y <= 0 ) {
      serror = "error: y not positive";
   } else if ( !p_given ) {
      serror = "error: pause time missing";
   } else if ( pause <= 0 ) {
      serror = "error: pause time not positive";
   } else if ( !m_given ) {
      serror = "error: minspeed missing";
   } else if ( minspeed <= 0 ) {
      serror = "error: minspeed not positive";
   } else if ( !M_given ) {
      serror = "error: maxspeed missing";
   } else if ( minspeed > maxspeed ) {
      serror = "error: minspeed larger than maxspeed";
   } else if ( !t_given ) {
      serror = "error: simtime missing";
   } else if ( simtime <= 0 ) {
      serror = "error: simtime non positive";
   }

   if ( strlen(serror) > 0 ) {
      fprintf(stderr, "%s\n", serror);
      exit(1);
   }

   positions = (struct npos *) malloc(nodes * sizeof(struct npos));
   assert(positions != NULL);

   avgspeed = (maxspeed + minspeed) / 2;
   avgdist  = 0.5 * sqrt(x*x+y*y) / sqrt(2);

   // set initial position of nodes randomly
   for ( i = 0; i < nodes; i++ ) {
      positions[i].x = x * drand48();
      positions[i].y = y * drand48();
      positions[i].z = 0.0; 
      positions[i].t = 0.0; 

      fprintf(stdout, "$node_(%d) set X_ %.1f \n", i, positions[i].x);
      fprintf(stdout, "$node_(%d) set Y_ %.1f \n", i, positions[i].y);
      fprintf(stdout, "$node_(%d) set Z_ %.1f \n", i, positions[i].z);
      fprintf(stdout, "\n");
   }

   // generate movements for each node according to the random direction mobility model 
   for ( i = 0; i < nodes; i++ ) {
      while ( positions[i].t < simtime ) {
         theta = 2 * PI * drand48(); // randomly choose an angle
         speed = minspeed + (maxspeed - minspeed) * drand48();
         etime  = -log(drand48())*(avgdist/avgspeed); // exponentially generate epoch time with average avgdist/avgspeed
         ptime  = 2 * pause_time * drand48();

         //fprintf(stdout, "theta: %g, speed %g, time %g\n", theta, speed, etime);
         rdmob(i, theta, speed, etime);
         positions[i].t += ptime;
      }
   }
}

void rdmob ( int nid, double theta, double speed, double etime )
{
   struct npos *p = &positions[nid];
   double dist = etime * speed;
   double nt;
   double cx, cy; // current x, y coordinates
   double nx, ny; // new x, y coordinates

   cx = p->x;
   cy = p->y;

   assert(cx>0);
   assert(cx<x);
   assert(cy>0);
   assert(cy<y);
   assert(etime>0);

   nx = cx + dist * cos(theta); // new x coordinate assuming no border restrictions
   ny = cy + dist * sin(theta); // new y coordinate assuming no border restrictions

   //fprintf(stdout, "node %d: x: %g, y: %g, dist: %g, time: %g, speed: %g, theta: %g, nx: %g, ny: %g\n", nid, cx, cy, dist, etime, speed, theta*180/PI, nx, ny);

   //fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, p->x, p->y, speed);

   bool yfirst = false;

   if ( theta < PI/2 ) { // right/up 
      if ( nx > x - 1 && ny > y - 1 ) {
         if ( tan(theta) > (y-1 - cy)/(x-1 - cx) ) {
            yfirst = true;
         }
      }

      if ( !yfirst && nx > x - 1 ) { // the trajectory exceeds the width, make it bounce
         nx = x - 1;
         ny = cy + tan(theta) * (nx - cx);
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);

         nt = sqrt((nx-cx)*(nx-cx) + (ny-cy)*(ny-cy)) / speed; // new elapsed time till the bouncing
         p->t += nt;
         p->x = nx; 
         p->y = ny; 

         // Recursively call rdmob with the new bounce angle and remaining time. Speed does not change, 
         // but the last (time) argument keeps decreasing (so recursion should stop at some point). Same below.
         rdmob(nid, PI - theta, speed, etime - nt);
      } else if ( ny > y - 1 ) {
         ny = y - 1;
         nx = cx + (ny - cy)/tan(theta); 
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);

         nt = sqrt((nx-cx)*(nx-cx) + (ny-cy)*(ny-cy)) / speed; // new elapsed time just before the bouncing
         p->t += nt;
         p->x = nx; 
         p->y = ny; 
         rdmob(nid, 2*PI - theta, speed, etime - nt);
      } else { // destination is within borders
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);
         p->t += etime;
         p->x = nx; 
         p->y = ny; 
      }
   } else if ( theta < PI ) { // left/up
      if ( nx < 1 && ny > y - 1 ) {
         if ( tan(theta) < (y-1 - cy)/(1 - cx) ) {
            yfirst = true;
         }
      }

      if ( !yfirst && nx < 1 ) {
         nx = 1;
         ny = cy + tan(theta)*(nx-cx);
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);

         nt = sqrt((nx-cx)*(nx-cx) + (ny-cy)*(ny-cy)) / speed; // new elapsed time till the bouncing
         p->t += nt;
         p->x = nx; 
         p->y = ny; 

         rdmob(nid, PI - theta, speed, etime - nt);
      } else if ( ny > y - 1 ) {
         ny = y - 1;
         nx = cx + (ny - cy)/tan(theta); 
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);

         nt = sqrt((nx-cx)*(nx-cx) + (ny-cy)*(ny-cy)) / speed; // new elapsed time just before the bouncing
         p->t += nt;
         p->x = nx; 
         p->y = ny; 
         rdmob(nid, 2*PI - theta, speed, etime - nt);
      } else {
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);
         p->t += etime;
         p->x = nx; 
         p->y = ny; 
      }
   } else if ( theta < 1.5 * PI  ) { // left/down
      if ( nx < 1 && ny < 1 ) {
         if ( tan(theta) > (1 - cy)/(1 - cx) ) {
            yfirst = true;
         }
      }

      if ( !yfirst && nx < 1 ) {
         nx = 1;
         ny = cy + tan(theta)*(nx-cx);
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);

         nt = sqrt((nx-cx)*(nx-cx) + (ny-cy)*(ny-cy)) / speed; // new elapsed time till the bouncing
         p->t += nt;
         p->x = nx;
         p->y = ny;

         if ( theta == PI ) {
            rdmob(nid, 0, speed, etime - nt);
         } else {
            rdmob(nid, 3*PI - theta, speed, etime - nt);
         }
      } else if ( ny < 1 ) {
         ny = 1;
         nx = cx + (ny - cy)/tan(theta); 
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);

         nt = sqrt((nx-cx)*(nx-cx) + (ny-cy)*(ny-cy)) / speed; // new elapsed time just before the bouncing
         p->t += nt;
         p->x = nx; 
         p->y = ny; 
         rdmob(nid, 2*PI - theta, speed, etime - nt);
      } else {
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);
         p->t += etime;
         p->x = nx; 
         p->y = ny; 
      }
   } else { // down/right
      if ( nx > x - 1 && ny < 1 ) {
         if ( tan(theta) < (1 - cy)/(x-1 - cx) ) {
            yfirst = true;
         }
      }

      if ( !yfirst && nx > x-1 ) {
         nx = x - 1;
         ny = cy + tan(theta)*(nx-cx);
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);

         nt = sqrt((nx-cx)*(nx-cx) + (ny-cy)*(ny-cy)) / speed; // new elapsed time till the bouncing
         p->t += nt;
         p->x = nx;
         p->y = ny;

         rdmob(nid, 3*PI - theta, speed, etime - nt);
      } else if ( ny < 1 ) {
         ny = 1;
         nx = cx + (ny - cy)/tan(theta); 
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);

         nt = sqrt((nx-cx)*(nx-cx) + (ny-cy)*(ny-cy)) / speed; // new elapsed time just before the bouncing
         p->t += nt;
         p->x = nx; 
         p->y = ny; 
         rdmob(nid, 2*PI - theta, speed, etime - nt);
      } else {
         fprintf(stdout, "$ns_ at %.1f  \"$node_(%d) setdest %f %f %f\" \n", p->t, nid, nx, ny, speed);
         p->t += etime;
         p->x = nx; 
         p->y = ny; 
      }
   }
}

void print_usage (FILE *stream, int exit_code)
{
   fprintf (stream, "Usage:  %s -n num_nodes -x x -y y -m minspeed -M maxspeed -p pause-time\n", program_name);
   fprintf (stream,
           "  -h  --help             Display this usage information.\n"
           "  -x  --x width          Width of rectangular area (Required).\n"
           "  -y  --y length         Length of rectangular area (Required).\n"
           "  -n  --nodes num_nodes  Number of nodes (Required, must be larger than zero).\n"
           "  -m  --minspeed speed   Minimum speed (Required, must be larger than zero).\n"
           "  -M  --maxspeed speed   Maximum speed (Required, must be larger than minspeed).\n"
           "  -t  --time simtime     Simulation time(Required, must be larger 0).\n"
           "  -p  --pause pause-time Pause time between epochs (Required, must be non-negative).\n");
   exit (exit_code);
}