// Compile with the command line:
// g++ prot.cpp -l SDL
// Type q, or click the upper right x to quit the program
//
// Requires the SDL library included with Linux
// Requires veblib.h for the structures
//

#include <stdlib.h>
#include <SDL/SDL.h>
#include <stdio.h>
#include <math.h>
#include "veblib.h"

void swapv(prot & p1, prot & p2){
  double tmp = p1.vx;
  p1.vx = p2.vx;
  p2.vx = tmp;
  tmp = p1.vy;
  p1.vy = p2.vy;
  p2.vy = tmp;
}

void bounce(prot & p1, prot & p2)
{ 
  double dx = p1.x - p2.x;
  double dy = p1.y - p2.y;
  dy *= dy; 
  dx *= dx;
  double dc = sqrt(dx + dy);          // distance between centers
  double n = (p1.x - p2.x) / dc;      // Normalized movement vectors
  double o = (p1.y - p2.y) / dc;
  double a1 = p1.vt * n;              // distance along new directions
  double a2 = p2.vt * n;
  double b1 = p1.vt * o;
  double b2 = p2.vt * o;

  double P = (2.0 * (a1 - a2)) / (p1.ms + p2.ms); // Change in momentum
  double v1 = p1.vt - P * p2.ms * n;
  double v2 = p2.vt + P * p1.ms * n;


  p1.dx = a1 / v1;                // New movement direction
  p2.dx = a2 / v2;
  p1.dy = b1 / v1;
  p2.dy = b2 / v2;
  p1.vx = a1;
  p2.vx = a2;
  p1.vy = b1;
  p2.vy = b2;
  
  // Code below forces the bounce in case the calculation above does not
  // work. Which, so  far, it does not work.
  switch((int) p1.q){
    case 1:
      if(p1.vx > 0) p1.vx *= -1.0;
      if(p1.vy < 0) p1.vy *= -1.0;
      if(p2.vx < 0) p2.vx *= -1.0;
      if(p2.vy > 0) p2.vy *= -1.0;
      break;
    case 2:
      if(p1.vx < 0) p1.vx *= -1.0;
      if(p1.vy < 0) p1.vy *= -1.0;
      if(p2.vx > 0) p2.vx *= -1.0;
      if(p2.vy > 0) p2.vy *= -1.0;
      break;
    case 3:
      if(p1.vx < 0) p1.vx *= -1.0;
      if(p1.vy > 0) p1.vy *= -1.0;
      if(p2.vx > 0) p2.vx *= -1.0;
      if(p2.vy < 0) p2.vy *= -1.0;
      break;
    case 4:
      if(p1.vx > 0) p1.vx *= -1.0;
      if(p1.vy > 0) p1.vy *= -1.0;
      if(p2.vx < 0) p2.vx *= -1.0;
      if(p2.vy < 0) p2.vy *= -1.0;
      break;
    default:
      p1.vx += 0.1; // Just to see if it ever gets here.
      p2.vx += 0.1;
      p1.vy += 0.1;
      p2.vy += 0.1;
      break;
  }
  
  p1.setvt();
  p2.setvt();
}
void resize(prot & p, double val){
  if(p.z < 1) p.z += val;
  else p.z = .1;
}

int main(int argc, char *argv[])
{
  char c;
  double s4 = .2;
  double s3 = s4 * 42.0;
  double s2 = s4 * 274.0;
  double s1 = s4 * 700.0;
  double x, y;
  int bump;
  int bumptol = 4;
  // Eight protons two of which made into neutrons p1, and p8
  prot p1, p2, p3, p4, p5, p6, p7, p8;
  p1.x = 50;
  p2.x = 250;
  p3.x = 450;
  p4.x = 750;
  p5.x = 650;
  p6.x = 400;   // x and y starting positions
  p7.x = 600;
  p8.x = 650;
  p1.y = 200; 
  p2.y = 550;
  p3.y = 750;
  p4.y = 50;
  p5.y = 150;
  p6.y = 400;
  p7.y = 550;
  p8.y = 340;  // z direction used for size s4 is inside shell four
               // other sizes are calculated via square of the shells rule
  p1.z = s4; p2.z = s4; p3.z = s4;
  p4.z = s4;
  p5.z = s4;
  p6.z = s4;
  p7.z = s4;
  p8.z = s4;
  p1.setmass();
  p1.ms += 2.549920405;  // Extra shell mass of the neutron
  p2.setmass();
  p3.setmass();
  p4.setmass();
  p5.setmass();
  p6.setmass();
  p7.setmass();
  p8.setmass();
  p8.ms += 2.549920405;  // Neutron uses proton struct plus shell one
  p1.vx = 0.12; p1.vy = 6.12;  // Starting movement vectors in cart
  p2.vx = 2.34; p2.vy = 2.34;
  p3.vx = 3.45; p3.vy = 3.45;
  p4.vx = 4.56; p4.vy = 4.56;
  p5.vx = 5.67; p5.vy = 5.67;
  p6.vx = 6.78; p6.vy = 6.78;
  p7.vx = 2.55; p7.vy = 1.40;
  p8.vx = 3.45; p8.vy = 2.87;

  p1.setvt();   // Calculate movement vector from vx and vy
  p2.setvt();
  p3.setvt();
  p4.setvt();
  p5.setvt();
  p6.setvt();
  p7.setvt();
  p8.setvt();

  p1.setvp();   // Calculate vx and vy from movement vector
  p2.setvp();   // Redundant, but checks the math
  p3.setvp();
  p4.setvp();
  p5.setvp();
  p6.setvp();
  p7.setvp();
  p8.setvp();
               // For ending printout to check conservation of momentum
  double tbp = p1.p + p2.p + p3.p + p4.p + p5.p; // beginning momentum
  tbp += p6.p;
  tbp += p7.p;
  tbp += p8.p;
 
  if((SDL_Init(SDL_INIT_VIDEO|SDL_INIT_AUDIO)==-1)) { 
    printf("Could not initialize SDL: %s.\n", SDL_GetError());
    exit(-1);
  }
  init_screen(maxx, maxy);
  ColorFill(color, maxx, maxy);
  while(1){ 
    ColorFill(color, maxx, maxy);
  
    p1.movevtd();
    p2.movevtd();
    p3.movevtd();
    p4.movevtd();
    p5.movevtd();
    p6.movevtd();
    p7.movevtd();
    p8.movevtd();

    neutron(p1.x, p1.y, p1.z);
    neutron(p8.x, p8.y, p8.z);

    if(collide(p1, p2) < bumptol) bounce(p1, p2);
    if(collide(p1, p3) < bumptol) bounce(p1, p3);
    if(collide(p1, p4) < bumptol) bounce(p1, p4);
    if(collide(p1, p5) < bumptol) bounce(p1, p5);
    if(collide(p1, p6) < bumptol) bounce(p1, p6);
    if(collide(p1, p7) < bumptol) bounce(p1, p7);
    if(collide(p1, p8) < bumptol) bounce(p1, p8);
    if(collide(p2, p3) < bumptol) bounce(p2, p3);
    if(collide(p2, p4) < bumptol) bounce(p2, p4);
    if(collide(p2, p5) < bumptol) bounce(p2, p5);
    if(collide(p2, p6) < bumptol) bounce(p2, p6);
    if(collide(p2, p7) < bumptol) bounce(p2, p7);
    if(collide(p2, p8) < bumptol) bounce(p2, p8);
    if(collide(p3, p4) < bumptol) bounce(p3, p4);
    if(collide(p3, p5) < bumptol) bounce(p3, p5);
    if(collide(p3, p6) < bumptol) bounce(p3, p6);
    if(collide(p3, p7) < bumptol) bounce(p3, p7);
    if(collide(p3, p8) < bumptol) bounce(p3, p8);
    if(collide(p4, p5) < bumptol) bounce(p4, p5);
    if(collide(p4, p6) < bumptol) bounce(p4, p6);
    if(collide(p4, p7) < bumptol) bounce(p4, p7);
    if(collide(p4, p8) < bumptol) bounce(p4, p8);
    if(collide(p5, p6) < bumptol) bounce(p5, p6);
    if(collide(p5, p7) < bumptol) bounce(p5, p7);
    if(collide(p5, p8) < bumptol) bounce(p5, p8);
    if(collide(p6, p7) < bumptol) bounce(p6, p7);
    if(collide(p6, p8) < bumptol) bounce(p6, p8);
    if(collide(p7, p8) < bumptol) bounce(p7, p8);
    SDL_UpdateRect(screen, minx, miny, maxx, maxy);
    c = keystroke();
    if(c == 'q') break;
    if(c == 'w') SDL_Delay(5000);
    if(c == 'x'){ p1.vx *= 10; p1.setvt(); }
    if(c == 'f'){
      p1.vt ++;
      p2.vt ++;
      p3.vt ++; 
      p1.setvp();
      p2.setvp();
      p3.setvp();     
    }
    if(c == 's'){
      p1.vt += 10;
      p2.vt -=10;
      p3.vt -=10;
      p1.setvp();
      p2.setvp();
      p3.setvp();
    }
    //SDL_Delay(5);
    
  }   
  SDL_Quit();
  double tep = p1.p + p2.p + p3.p + p4.p + p5.p + p6.p;
  tep += p7.p;
  tep += p8.p;
  printf("\nBeginning Momentum = %5.5f: Ending Momemtun = %5.5f ", tbp, tep);
  printf("\nOn exit p1: vt = %5.2f d = %5.2f p = %5.2f", p1.vt, p1.dx, p1.p);
  printf(" vx = %5.2f vy = %5.2f ", p1.vx, p1.vy);
  printf("\nOn exit p2: vt = %5.2f d = %5.2f p = %5.2f", p2.vt, p2.dx, p2.p);
  printf(" vx = %5.2f vy = %5.2f ", p2.vx, p2.vy);
  printf("\nOn exit p3: vt = %5.2f d = %5.2f p = %5.2f", p3.vt, p3.dx, p3.p);
  printf(" vx = %5.2f vy = %5.2f ", p3.vx, p3.vy);
  printf("\nOn exit p4: vt = %5.2f d = %5.2f p = %5.2f", p4.vt, p4.dx, p4.p);
  printf(" vx = %5.2f vy = %5.2f ", p4.vx, p4.vy);
  printf("\nOn exit p5: vt = %5.2f d = %5.2f p = %5.2f", p5.vt, p5.dx, p5.p);
  printf(" vx = %5.2f vy = %5.2f ", p5.vx, p5.vy);
  printf("\nOn exit p6: vt = %5.2f d = %5.2f p = %5.2f", p6.vt, p6.dx, p6.p);
  printf(" vx = %5.2f vy = %5.2f ", p6.vx, p6.vy);
  printf("\n");
  exit(0);   
}