/********************************************************
Graphics functions and structures for the neutron model
*********************************************************/
Uint32 color = 0x55ddff;
Uint32 Red = 0xff0000, Blue = 0x0000ff;
Uint32 Yellow =  0xffff00, Greenish = 0x00ffff;
#define sign(a) (((a)<0) ? -1 : (a)>0 ? 1 : 0)
#define PI 3.141592653589793238462643383279502884197169399375105
int minx = 0, miny = 0, maxx = 800, maxy = 600;
void proton(double, double, double);
void putpixel(unsigned int, unsigned int, unsigned int);
void line(int, int, double, double, int);
void linexy(int, int, int, int, int);
int ColorFill(int color, int maxx, int maxy);
struct vector
{
  double x;
  double y;
  void normalize(void){
    double vt = sqrt(x * x + y * y);
    x = x / vt;
    y = y / vt;
  }
  void dot(vector & a, double & dotp){
    dotp = x * a.x + y * a.y;
  }
};

struct abc
{            // Write to screen based upon 
  int curx;  // turn and go notation
  int cury;
  int scale;
  int font;
  int color;
  void init(int x, int y, int sc){
    curx = x;
    cury = y;
    scale = sc;
    font = 1;
  }
};

struct control
{
  int view;
  int trigger;
  int lines;
  int vec;
  int twovec;
  void reset(){
    view = 0;
    lines = 0;
    vec = 0;
    trigger = 0;
    twovec = 0;
  }
  void setview(){
    view = 1;
  }
  void setlines(){
    lines = 1;
  }
  void setvec(){
    vec = 1;
  }
  void settwovec(){
    twovec = 1;
  }
};

struct nprot
{             // Proton based upon polar coordinates
  double x;   // x position
  double y;   // y position
  double z;   // z position as size
  double vx;  // x component movement
  double vy;  // y component movement
  double vz;  // z component movement
  double vt;  // Movement vector magnitude
  double rad; // Movement vector direction
  double rs1; // Polarity vector of shell one
  double rs2; // Polarity vector of shell two
  double rs3; // Polarity vector of shell three
  double rse; // Electron polarity vector
  double ms1;
  double ms2;
  double ms3;
  double ms4;
  double ms;  // Total Mass
  double p;   // Total Momemtum
  int half;   // Time to live when it is a neutron
  int lines;  // Special effects flag
  int electron;
  int id;
  int bumper;
  void addelec(){
    electron = 1;
  }
  void makeneut(){
    ms1 = 2.549920405;
    ms = ms1 + ms2 + ms3 + ms4;
    half = rand() % 1000;  // About 6 minutes time to live
    p = vt * ms;
  }
  void setmass(){  // In units of electron masses
    ms1 = 2.549920405;
    ms2 = ms1 * ms1;
    ms3 = ms2 * ms2;
    ms4 = ms3 * ms3;
    int twopi = 2 * PI * 1000;
    vt = (float) (rand() % 1000) / 100;
    x = rand() % maxx;
    y = rand() % maxy;
    z = .10;
    rad = (float) (rand() % twopi) / 1000;
    rs1 = (float) (rand() % twopi) / 1000;
    rs2 = (float) (rand() % twopi) / 1000;
    rs3 = (float) (rand() % twopi) / 1000;
    rse = (float) (rand() % twopi) / 1000;
    lines = 0;
    electron = 0;
    makeneut();
  }
  void makeprot(){
    ms1 = 0;
    vt += 5;
    ms = ms2 + ms3 + ms4;
    p = ms * vt;
  }
  void rotate(){  // Rotate the amplitude around the shell
    rs1 += .310;  // Speed of rotation not to scale
    rs2 -= .610;  // Middle shell rotating counter clockwise
    rs3 += .810;  // Direction of rotation is SWAG
    rse -= .180;
    if(rse < 0)      rse += 2 * PI;
    if(rs1 > 2 * PI) rs1 -= 2 * PI; // Keep it positive between
    if(rs2 < 0)      rs2 += 2 * PI; // zero and 2 PI
    if(rs3 > 2 * PI) rs3 -= 2 * PI;
  }
  void ellipse(){ // Draw the proton neutron and electron
    double c = 0; // Shell sizes in units of shell 4
    double sr4  = z;       if(sr4 < 1) sr4 = 1;
    double sr3  = z * 42;  if(sr3 < 1) sr3 = 1;
    double sr2  = z * 274; if(sr2 < 1) sr2 = 1;
    double sr1  = z * 701; if(sr1 < 1) sr1 = 1;
    double elec = z * 3497;
    double inc = PI / sr1; // Dot density around circle
    double prot = ms2 + ms3 + ms4;
    int minamp, maxamp;
    double xx, yy, xrs, yrs;
    int ofst = 0;
    int pos = 0xff0000;   // Red for positive field
    int neg = 0x0000ff;   // Blue for negative field
    rotate();
    while(c < 2 * PI){
      c += inc;
      rs1 += inc;
      rs2 += inc;
      rs3 += inc;
      rse += inc;
      xx = x + cos(c) * sr4;
      yy = y + sin(c) * sr4;
      putpixel(xx, yy, pos);          // One dot for shell four
      minamp = sin(rs3) * z * 20;   
      maxamp = sin(rs3) * z * 20;   
      if(minamp > 0) minamp *= -1.0;
      if(maxamp < 0) maxamp *= -1.0;               
      for(ofst = minamp; ofst < maxamp; ofst ++){ 
        xx = x + cos(c) * (sr3 + ofst);  
        yy = y + sin(c) * (sr3 + ofst);
        if(ofst < 0){                 // Shell three hash marks         
          putpixel(xx, yy, pos);
        }
        else{        
          putpixel(xx, yy, neg);
        }
      }
      minamp = sin(rs2) * z * 60;
      maxamp = sin(rs2) * z * 60;
      if(minamp > 0) minamp *= -1.0;
      if(maxamp < 0) maxamp *= -1.0;
      for(ofst = minamp; ofst < maxamp; ofst ++){
        xx = x + cos(c) * (sr2 + ofst);
        yy = y + sin(c) * (sr2 + ofst);
        if(ofst < 0){               // Shell two hash marks  
          putpixel(xx, yy, neg);
        }
        else{
          putpixel(xx, yy, pos);
        }
      }
      if(ms > prot){                   // It is a neutron
        minamp = sin(rs1) * z * 200;   // so draw shell one
        maxamp = sin(rs1) * z * 200;   // maxamp can be a different value
        if(minamp > 0) minamp *= -1.0; // minamp will be less than zero
        if(maxamp < 0) maxamp *= -1.0; // maxamp will be more than zero
        for(ofst = minamp; ofst < maxamp; ofst ++){
          xx = x + cos(c) * (sr1 + ofst);
          yy = y + sin(c) * (sr1 + ofst);
          if(ofst < 0){               // Shell one hash marks  
            putpixel(xx, yy, pos);
          }
          else{
            putpixel(xx, yy, neg);
          }
        }
        
      }
      if(electron){
        inc = PI / elec;
        minamp = sin(rse) * z * 300;
        maxamp = sin(rse) * z * 300;
        if(minamp > 0) minamp *= -1;
        if(maxamp < 0) maxamp *= -1;
        for(ofst = minamp; ofst < maxamp; ofst ++){
          xx = x + cos(c) * (elec + ofst);
          yy = y + sin(c) * (elec + ofst);
          if(ofst < 0){               // Shell one hash marks  
            putpixel(xx, yy, pos);
          }
          else{
            putpixel(xx, yy, neg);
          }
        }
      }
    }  
  }
  void move(){
    int b = 0;
    double r = z * 274;
    if(x + r > maxx){ rad = PI - rad; x = maxx - r; b = 1; }
    if(x - r < minx){ rad = PI - rad; x = minx + r; b = 1; }
    if(y + r > maxy){ rad = 2 * PI - rad; y = maxy - r; b = 1; }
    if(y - r < miny){ rad = 2 * PI - rad; y = miny + r; b = 1; }
    if(b == 1 && vt < 20){ vt *= 2.50; }
    vx = cos(rad) * vt;
    vy = sin(rad) * vt;
    x += vx;
    y += vy;
    ellipse();
    if(half > 0) half --;
    if(half == 1) makeprot();
    if(lines > 0) lines --;
    p = ms * vt;
  }
  void getrad(void){  // Get radian from movement vectors
    double d = sqrt(vx * vx + vy * vy);
    double v = vx / d;
    rad = acos(v / vx * vy);
    if(v > 0){ rad = 2 * PI - rad; }
    rad += PI / 2;
  }
  void getvt(void){
    vt = sqrt(vx * vx + vy * vy);
  }
  void bump(nprot & p1, control & c1){
    if(c1.lines == 1) linexy(x, y, p1.x, p1.y, 0x000000);
    if(z < p1.z || z > p1.z) return;              // Wrong plane exit
    double dx = x - p1.x;
    double dy = y - p1.y;
    double sr = z * 274 + p1.z * 274;             // Collision distance
    if(dx * dx + dy * dy > sr * sr){ 
      if(p1.id == bumper){ 
        bumper = -1; 
        p1.bumper = -1;
      }      
      return;                                     // No collision return
    }
    if(p1.id == bumper){ return; }                // Same bumper return
    double dc = sqrt(dx * dx + dy * dy);          // Distance between centers
    double nx = dx / dc;
    double ny = dy / dc;
    double a1, a2;
    /*********************************************/
    if(dc < sr - 10){  // Makeshift binding should be modified to rely on
      bumper = -1;     // electrical forces at surface of shells once the
      p1.bumper = -1;  // outside shell is penetrated by the speed of the
      vt = 1;          // particles.
      p1.vt = vt; p1.rad = rad;
      p1.x = x + z * 100;
      p1.y = y + z * 100;
      return; 
    }
    /*********************************************/
    a1 = vx * nx + vy * ny;
    a2 = p1.vx * nx + p1.vy * ny;
    double op = (2 * (a1 - a2)) / (ms + p1.ms);
    vx = vx - op * p1.ms * nx;
    vy = vy - op * p1.ms * ny;
    p1.vx = p1.vx + op * ms * nx;
    p1.vy = p1.vy + op * ms * ny;
    getvt();
    p1.getvt();
    getrad();                                     // Get new radians
    p1.getrad();
    if(bumper <= 0){ bumper = p1.id; }
    if(p1.bumper <= 0){ p1.bumper = id; }
    if(c1.view == 1) c1.trigger = 1;
  }
  void showvec(){                                 // Draw a line along movement vector
    line(x, y, rad, vt * 10, 0xffff00);
  }
};

SDL_Surface *screen;
int init_screen(int srt, int sdn)
{
  screen = SDL_SetVideoMode(srt, sdn, 32, SDL_SWSURFACE);
  if( screen == NULL ){
    fprintf(stderr, "Couldn't set video mode: %s\n",
    SDL_GetError());
    exit(1);
  }
}

void putpixel(unsigned int x, unsigned int y, unsigned int color)
{
  if(x >= minx && x < maxx && y >= miny && y < maxy){
     unsigned int *ptr = (unsigned int*)screen->pixels;
     int lineoffset = y * (screen->pitch / 4);
     ptr[lineoffset + x] = color;
  }
}
Uint32 getpixel(SDL_Surface *surface, int x, int y)
{
  int bpp = surface->format->BytesPerPixel;
  // Here p is the address to the pixel we want to retrieve 
  if(x >= minx && x < maxx && y >= miny && y < maxy){
     Uint8 *p = (Uint8 *)surface->pixels + y * surface->pitch + x * bpp;
     switch(bpp) {
        case 1:
          return *p;
        case 2:
          return *(Uint16 *)p;
        case 3:
          if(SDL_BYTEORDER == SDL_BIG_ENDIAN)
              return p[0] << 16 | p[1] << 8 | p[2];
          else
              return p[0] | p[1] << 8 | p[2] << 16;
        case 4:
          return *(Uint32 *)p;
        default:
          return 0;       /* shouldn't happen, but avoids warnings */
     }
  }
  return 0; 
}

int ColorFill(int color, int maxx, int maxy)
{
    int x, y;
    for(y = 0; y < maxy; y++){
       for(x = 0; x < maxx; x++){
          putpixel(x, y, color);
       }
    }
}

void linexy(int x, int y, int x2, int y2, int hue)
{  // Draw a line using Cartesian coordinates
  double a = x2, b = y2, xx;
  double dx = x - x2, dy = y - y2;
  double vt = sqrt(dx * dx + dy * dy);
  double vx = dx / vt;
  double vy = dy / vt;
  for(xx = 0; xx < vt; xx ++){
    putpixel(a, b, hue);
    a += vx;
    b += vy;
  }
}

void line(int x, int y, double rad, double vt, int hue)
{ // Draw a line using polar coordinates
  double a = x, b = y, xx = 0, vx = cos(rad), vy = sin(rad);
  while(xx++ < vt){
    putpixel(a, b, hue);
    a += vx;
    b += vy;
  }
}

void MakeCircle(float xx, float yy, float radius, Uint32 color)
{                // Larger circle gets more dots per degree
  double inc = 1 / (radius * 2);
  double circle = 0;
  while(circle < 2 * PI){
    putpixel(xx + cos(circle) * radius, yy + sin(circle) * radius, color);
    circle += inc;
  }
}

char keystroke(void)
{
  int count = 0;
  char keydat = 0;
  SDL_Event event;
  /* Poll for events. SDL_PollEvent() returns 0 when there are no  */
  /* more events on the event queue, our while loop will exit when */
  /* that occurs.                                                  */
  while( SDL_PollEvent( &event ) ){
    switch( event.type ){
      case SDL_KEYDOWN:
        keydat = event.key.keysym.sym;
        //printf( "Key press detected\n" );
        count ++;
        break;
      case SDL_KEYUP:
        //printf( "Key release detected\n" );
        count ++;
        break;
      case SDL_QUIT:
        keydat = 'q';
        break;
      case SDLK_KP_ENTER:
        keydat = 0x0a;
        break;
      case SDLK_RETURN:
        keydat = 0x0a;
        break;
      case SDLK_KP0:
        keydat = '0';
        break;
      case SDLK_KP1:
        keydat = '1';
        break;
      case SDLK_KP2:
        keydat = '2';
        break;
      case SDLK_KP3:
        keydat = '3';
        break;
      default:
        break;
    }
  }
  return keydat;
}