#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <GLUT/glut.h>

// display functions
void mouse(int button, int state, int x, int y) {
	exit(0);
}
void idle(void) {
}

// linear feedback shift register (LFSR)
// ori: input value
int LFSR(int ori){
  // copy ori to val
  int val = ori;
  
  // specify number of bits in LFSR and integer value of largest bit
  int M = 24;
  int top = 1 << (M-1);   // bitwise leftshift operator, to generate powers of 2
  // build taps array
	int Ntaps = 4;
	int taps_array[Ntaps];
	taps_array[0] = 1;
	taps_array[1] = 2;
	taps_array[2] = 7;
	taps_array[3] = 24;
	
  // iterate taps array M times
  int i, j, counter;
  for(i=0; i<M; i++){
	  // pass ori through taps array
    counter = 0;
	  for(j=0; j<Ntaps; j++){
      // bitwise rightshift and mod operators, to check if jth position 0 or 1
	  	if((val >> (M - taps_array[j])) % 2){
	  		counter += 1;
	  	}
	  }
    
	  // determine output
	  int ter;
	  if(counter % 2){
      // set top bit 1 and shift all bits to the right
	  	ter = top + (val >> 1);
	  }
	  else{
      // set top bit 0 and shift all bits to the right
	  	ter = (val >> 1);
	  }
    
    // update val
    val = ter;
  }
  
  // output
  return val;
}

// define coarse-grained molecule
struct molecule{
  // id = -1 defines a root of linked list
  // id = 0 defines a water (just one hydrophilic)
  // id = 1 defines a lipid (one hydrophilic and two hydrophobic particles)
  int id;
  // one hydrophilic particle
  double x_phil;
  double y_phil;
  double vx_phil;
  double vy_phil;
  double Fx_phil;
  double Fy_phil;
  // two hydrophobic particles
  double x_phob1;
  double y_phob1;
  double vx_phob1;
  double vy_phob1;
  double Fx_phob1;
  double Fy_phob1;
  double x_phob2;
  double y_phob2;
  double vx_phob2;
  double vy_phob2;
  double Fx_phob2;
  double Fy_phob2;
};

// define element of linked list of molecules
struct molEl{
  struct molecule m;              // each element includes one molecule ...
  struct molEl *left;             // ... and a pointer to the previous molecule in the cell ...
  struct molEl *right;            // ... and a pointer to the next molecule in the cell
};

// define function to compute force between two particles
struct forces {
  double Fx, Fy;
};
struct forces compute_force(struct molEl *p1, struct molEl *p2){
  // define output
  struct forces f;
  // define short parameter names
  int id1, id2;
  id1 = p1->m.id;
  id2 = p2->m.id;
  double x1, y1, vx1, vy1, x2, y2, vx2, vy2;
  double a_AA = 25;
  double a_AB = 50;
  double a_AW = 26;
  double a_BB = 25;
  double a_BW = 50;
  double a_WW = 25;
  // force between two waters
  if(id1 == 0 && id2 == 0){
    
  }
  // force between water and lipid
  else if(id1 == 0 && id2 == 1){
    
  }
  // force between lipid and water
  else if(id1 == 1 && id2 == 0){
    
  }
  // force between two lipids
  else if(id1 == 1 && id2 == 1){
    
  }
  
  // output structure
  return f;
}

// functions to add and remove nodes from linked lists
void insert_molEl_right(struct molEl *i, struct molEl *p){
  struct molEl *f;
  f = i->right;
  
  i->right = p;
  p->right = f;
  f->left = p;
  p->left = i;
}
void remove_molEl(struct molEl *p){
  struct molEl *l;
  struct molEl *r;
  l = p->left;
  r = p->right;
  
  l->right = r;
  r->left = l;
}

void display(void) {
  // display information
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_MODELVIEW);
  
  // initialize parameters
  int Nwater = 200;               // number of water molecules
  int Nlipid = 100;               // number of lipid molecules
  int Lbox = 10;                  // length of box
  int vi_max = 2;                 // maximum initial velocity
  int Rcutoff = 1;                // cutoff radius
  double Rparticle = 0.2;         // particle radius
  int Lcell = Rcutoff;            // set cell length to cutoff radius
  int Ncells = Lbox/Rcutoff;      // define number of cells along each axis
  int Nsteps = 100;               // number of integration steps
  double deltaT = 0.01;           // time step
  double lambda = 0.65;           // integration parameter
  
  int M = 24;                     // specify number of bits in LFSR
  int top = 1 << M;               // calculate largest integer value
  int val = 1;                    // initialize LFSR
  double r, theta;                // initialize random numbers
  
  // initialize cells
  struct molEl *cells_left[Ncells][Ncells];    // each element points to left root in cell
  struct molEl *cells_right[Ncells][Ncells];   // each element points to right root in cell
  int m,n;
  for(m=0; m<Ncells; m++){
    for(n=0; n<Ncells; n++){
      // allocate memory
      cells_left[m][n] = malloc(sizeof(struct molEl));
      cells_right[m][n] = malloc(sizeof(struct molEl));
      // declare left and right roots of linked list
      cells_left[m][n]->m.id = -1;
      cells_left[m][n]->left = 0;
      cells_left[m][n]->right = cells_right[m][n];
      cells_right[m][n]->m.id = -1;
      cells_right[m][n]->left = cells_left[m][n];
      cells_right[m][n]->right = 0;
    }
  }
  
  // generate molecules and assign to cells
  int i;
  // create and assign water molecules
  struct molEl *waters[Nwater];
  for(i=0; i<Nwater; i++){
    // allocate and create water molecule
    waters[i] = malloc(sizeof(struct molEl));
    waters[i]->m.id = 0;
    val = LFSR(val);
    r = ((double)val)/top;
    waters[i]->m.x_phil = Lbox*r;
    val = LFSR(val);
    r = ((double)val)/top;
    waters[i]->m.y_phil = Lbox*r;
    val = LFSR(val);
    r = ((double)val)/top;
    waters[i]->m.vx_phil = vi_max*r;
    val = LFSR(val);
    r = ((double)val)/top;
    waters[i]->m.vy_phil = vi_max*r;
    
    // add molecule to appropriate linked list
    m = floor(waters[i]->m.x_phil);
    n = floor(waters[i]->m.y_phil);
    insert_molEl_right(cells_left[m][n], waters[i]);
    
    // display water molecule, translate and color particles
	  glPushMatrix();
	  glTranslatef(waters[i]->m.x_phil, waters[i]->m.y_phil, 0);
	  glColor3f(0.0f,0.0f,0.5f);
	  glutSolidSphere(Rparticle, 40, 40);
	  glPopMatrix();
  }
  // create lipid molecules
  struct molEl *lipids[Nlipid];
  for(i=0; i<Nlipid; i++){
    // allocate and create lipid molecule
    lipids[i] = malloc(sizeof(struct molEl));
    // hydrophilic head
    lipids[i]->m.id = 1;
    val = LFSR(val);
    r = ((double)val)/top;
    lipids[i]->m.x_phil = Lbox*r;
    val = LFSR(val);
    r = ((double)val)/top;
    lipids[i]->m.y_phil = Lbox*r;
    val = LFSR(val);
    r = ((double)val)/top;
    lipids[i]->m.vx_phil = vi_max*r;
    val = LFSR(val);
    r = ((double)val)/top;
    lipids[i]->m.vy_phil = vi_max*r;
    // hydrophobic thorax
    val = LFSR(val);
    r = ((double)val)/top;
    theta = 2*M_PI*r;
    lipids[i]->m.x_phob1 = lipids[i]->m.x_phil + Rparticle*cos(theta);
    lipids[i]->m.y_phob1 = lipids[i]->m.y_phil + Rparticle*sin(theta);
    val = LFSR(val);
    r = ((double)val)/top;
    lipids[i]->m.vx_phob1 = vi_max*r;
    val = LFSR(val);
    r = ((double)val)/top;
    lipids[i]->m.vy_phob1 = vi_max*r;
    // hydrophobic abdomen
    val = LFSR(val);
    r = ((double)val)/top;
    theta = 2*M_PI*r;
    lipids[i]->m.x_phob2 = lipids[i]->m.x_phob1 + Rparticle*cos(theta);
    lipids[i]->m.y_phob2 = lipids[i]->m.y_phob1 + Rparticle*sin(theta);
    val = LFSR(val);
    r = ((double)val)/top;
    lipids[i]->m.vx_phob2 = vi_max*r;
    val = LFSR(val);
    r = ((double)val)/top;
    lipids[i]->m.vy_phob2 = vi_max*r;
    
    // add molecule to appropriate linked list
    m = floor(lipids[i]->m.x_phil);
    n = floor(lipids[i]->m.y_phil);
    insert_molEl_right(cells_left[m][n], lipids[i]);
    
    // display lipid molecule, translate and color particles
	  glPushMatrix();
	  glTranslatef(lipids[i]->m.x_phil, lipids[i]->m.y_phil, 0);
	  glColor3f(0.5f,0.0f,0.0f);
	  glutSolidSphere(0.2, 40, 40);
	  glPopMatrix();
	  glPushMatrix();
	  glTranslatef(lipids[i]->m.x_phob1, lipids[i]->m.y_phob1, 0);
	  glColor3f(0.5f,0.5f,0.0f);
	  glutSolidSphere(0.2, 40, 40);
	  glPopMatrix();
	  glPushMatrix();
	  glTranslatef(lipids[i]->m.x_phob2, lipids[i]->m.y_phob2, 0);
	  glColor3f(0.5f,0.5f,0.0f);
	  glutSolidSphere(0.2, 40, 40);
	  glPopMatrix();
  }
  // create conductor pointer, to navigate linked lists
  struct molEl *conductor_to;
  struct molEl *conductor_from;
  
  // integrate system
  int dx, dy, neighbor_x, neighbor_y;
  double trial_vx, trial_vy, upd_Fx, upd_Fy;
  int t;
  for(t=0; t<Nsteps; t++){
    /*                             */
    /* compute forces on particles */
    /*                             */
    // loop through cells
    for(m=0; m<Ncells; m++){
      for(n=0; n<Ncells; n++){
        // loop through molecules within cell (compute force ON these particles)
        conductor_to = cells_left[m][n];
        while(conductor_to != NULL){
          // do not process roots of linked lists
          if(conductor_to->m.id != -1){
            // if lipid, compute force on self
            if(conductor_to->m.id == 1){
              
            }
            // loop over neighboring cells
            for(dx=-1; dx<2; dx++){
              for(dy=-1; dy<2; dy++){
                neighbor_x = (m + dx) % Ncells;     // periodic boundary conditions
                if(neighbor_x < 0){
                  neighbor_x += Ncells;
                }
                neighbor_y = (n + dy) % Ncells;     // periodic boundary conditions
                if(neighbor_y < 0){
                  neighbor_y += Ncells;
                }
                // loop through neighboring molecules (compute force FROM these particles)
                conductor_from = cells_left[neighbor_x][neighbor_y];
                while(conductor_from != NULL){
                  
                  conductor_to->m.Fx_phil += 
                }
              }
            }
          }
          
          // advance to next molecule in cell
          conductor = conductor->right;
        }
      }
    }
    
    /*                                                                              */
    /* update molecule positions and velocities, reset forces to zero for next loop */
    /*                                                                              */
    for(m=0; m<Ncells; m++){
      for(n=0; n<Ncells; n++){
        // loop through molecules within cell
        conductor_to = cells_left[m][n];
        while(conductor_to != NULL){
          // do not process roots of linked lists
          if(conductor_to->m.id != -1){
            /*                                      */
            /* modified velocity Verlet algorithm   */
            /*                                      */
            // water molecule
            if(conductor_to->m.id == 0){
              // update position
              conductor_to->m.x_phil += deltaT*conductor_to->m.vx_phil + 0.5*pow(deltaT,2)*conductor_to->m.Fx_phil;
              conductor_to->m.y_phil += deltaT*conductor_to->m.vy_phil + 0.5*pow(deltaT,2)*conductor_to->m.Fy_phil;
              // build trial velocity
              trial_vx = conductor_to->m.vx_phil + lambda*deltaT*conductor_to->m.Fx_phil;
              trial_vy = conductor_to->m.vy_phil + lambda*deltaT*conductor_to->m.Fy_phil;
              // update force
              upd_Fx = ;
              upd_Fy = ;
              // update velocity
              conductor_to->m.vx_phil += 
            }
            // lipid molecule
            if(conductor_to->m.id == 1){
              
            }
          }
          
          // advance to next molecule in cell
          conductor = conductor->right;
        }
      }
    }
    
    /*                                                                    */
    /* assign molecules to new cells, impose periodic boundary conditions */
    /*                                                                    */
    
  }
  
  // render graphics
  glFlush();
}

int main(int argc, char **argv) {
	GLfloat light_ambient[] = {0.5, 0.5, 0.5, 1.0};
	GLfloat light_diffuse[] = {0.8, 0.8, 0.8, 1.0};
	GLfloat light_specular[] = {0.3, 0.3, 0.3, 1.0};
	GLfloat light_position[] = {0, 0, 1, 0.0};
  
	glutInit(&argc,argv);
	glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH);
	glutInitWindowSize(500, 500);
	glutCreateWindow("MD simulation");
	glutDisplayFunc(display);
	glutMouseFunc(mouse);
	glutIdleFunc(idle);
	glEnable(GL_LIGHTING);
	glEnable(GL_COLOR_MATERIAL);
	glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, 1.0);
	glLightModelfv(GL_LIGHT_MODEL_AMBIENT, light_ambient);
	glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
	glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse);
	glLightfv(GL_LIGHT0, GL_SPECULAR, light_specular);
	glLightfv(GL_LIGHT0, GL_POSITION, light_position);
	glEnable(GL_LIGHT0);
	glEnable(GL_DEPTH_TEST);
	glClearColor(1.0, 1.0, 1.0, 1.0);
	glMatrixMode(GL_PROJECTION);
  glOrtho(0, 10, 0, 10, -0.3, 0.3);
	glutMainLoop();
  
  return 0;
}