#include <iostream>
#include <cmath>
#include <plotter.h>
#include <vector>

using namespace std;
void GetStablePoint(double mu, int Nskip, vector<double>& xn, double x0 = 0.01)
{////////////////////////////////////////////////////////////////////
 // simulation of the bug population
 // Input:
 //         mu     --- the growth parameter
 //         Nskip  --- number of generations before steady state
 //         x0     --- starting number of bugs
 // Output: 
 //         xn     --- generations after steady state is reached
 ///////////////////////////////////////////////////////////////////
  double x = x0;
  for (int i=0; i<Nskip; i++)  x = mu*x*(1-x); // before the stedy state
  // we have stedy state now
  for (int i=0; i<xn.size(); i++){  // store next few generations
    x = mu*x*(1-x);
    xn[i] = x;
  }
}

bool InitializePlotter(XPlotter& plotter)
{
  //************** Initialization of plotter **************************/
  if (plotter.openpl () < 0){ // open Plotter      
    cerr<<"Couldn't open Plotter"<<endl;;
    return false;
  }

  plotter.fspace (0, -0.2, 5, 1.2);  // specify user coordinate system
  plotter.filltype (0);              // objects will be filled
  plotter.bgcolorname ("white");     // background color for the window
  
  plotter.fline(1,0,4,0); // coordinate system
  plotter.fline(1,0,1,1); // is plotted
  
  plotter.flinewidth (0.01);         // line thickness in user coordinates 
  plotter.colorname("black");
  plotter.ffontsize(0.1);
  plotter.ffontname("HersheySerif-Bold"); // ticks labels
  plotter.fmove(1.,-0.07);
  plotter.alabel('c', 'c', "1.0");
  plotter.fmove(4.,-0.07);
  plotter.alabel('c', 'c', "4.0");
  plotter.fmove(0.85, 0.0);
  plotter.alabel('c', 'c', "0.0");
  plotter.fmove(0.85, 1.0);
  plotter.alabel('c', 'c', "1.0");
  
  plotter.ffontsize(0.12);
  plotter.fmove(0.8, 0.5);
  plotter.alabel('c', 'c', "x\\sbn");    // axis labels
  plotter.fmove(2.5,-0.1);
  plotter.ffontname("Symbol");
  plotter.alabel('c', 'c', "m");
  
  plotter.colorname("red");
  return true;
}

int main(int argc, char *argv[], char *env[])
{
  int N=2000;      // Number of points for parameter mu
  int Nskip = 400; // Number of bug generations before reaching stedy state
  int Nplot = 100; // Number of generations plotted (after Nskip)
  int Nstart = 10; // Number of different generations of bugs (with different starting value)
  
  int i=0;
  while (++i<argc){
    std::string str(argv[i]);
    if (str=="-N" && i<argc-1)      N       = atoi(argv[++i]);
    if (str=="-Nskip" && i<argc-1)  Nskip   = atoi(argv[++i]);
    if (str=="-Nplot" && i<argc-1)  Nplot   = atoi(argv[++i]);
    if (str=="-Nstart" && i<argc-1) Nstart  = atoi(argv[++i]);
    if (str=="-h" || str=="--help"){
      std::clog<<"************************* Bug population simulation *********\n";
      std::clog<<"**            Copyright KH, January 2007                   **\n";
      std::clog<<"*************************************************************\n";
      std::clog<<"Usage: bugs [Options]\n";
      std::clog<<"Options:   -N         Number of mu points ("<<N<<")\n";
      std::clog<<"           -Nskip     Number of bug generations before reaching stedy state  ("<<Nskip<<")\n";
      std::clog<<"           -Nplot     Number of generations plotted (after Nskip)  ("<<Nplot<<")\n";
      std::clog<<"           -Nstart    Number of different initial conditions simulated ("<<Nstart<<")\n";
      return 0;
    }
  }

  //************** Initialization of plotter **************************/
  PlotterParams params; // set a Plotter parameter
  params.setplparam ("BITMAPSIZE", (char *)"1000x280");
  XPlotter plotter(cin, cout, cerr, params); // declare Plotter
  //  GIFPlotter plotter(cin, cout, cerr, params); // declare Plotter
  if (!InitializePlotter(plotter)) return 1;
  //************** Initialization of plotter **************************/
  
  vector<double> xn(Nplot);  // this vector will contain bug population after the stedy state is reached
  srand48(time(0));          // initialization of random number generator
  
  for (int l=0; l<Nstart; l++){ // Nstart number of different initial conditions are simulated
    double x0 = drand48();      // initial number of bugs is choosen randomly
    
    for (int i=0; i<N; i++){
    
      double mu = 1 + 3.*i/(N-1.);       // growth parameter is between [1-4]
      GetStablePoint(mu, Nskip, xn, x0); // few (xn.size()) generations after stedy state is reached are computed
    
      for (int j=0; j<xn.size(); j++) plotter.fpoint (mu, xn[j]); // plotting of bug populations
    }
    
    plotter.flushpl();// to flush output, otherwise display will be empty untill we are finished
  }
  
  /************** Plotter Done**************************/
  if (plotter.closepl () < 0) cerr<<"Couldn't close Plotter"<<endl;
  /*****************************************************/
  return 0;
}