/* spring3gl.c   spring and dropping mass  command line  s for single buffer */

#include <GL/glut.h>
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265
#endif

static int winWidth = 600;
static int winHeight = 700;
static int xm, ym, bm;       /* mouse position                        */
static int xc = 100;         /* disc/rectangle contact location       */
static int mass_zer = 600;   /* dropping mass reset position          */
static int yc = 600;
static int xs = 130;         /* end of spring location                */
static int spring_top = 380; /* spring zero force                     */
static int spring_bot = 170; /* spring fixed bottom                   */
static int ys = 380; 
static int xp;               /* last mouse location during drag       */
static int yp;
static int tracking = 0;     /* tracking mouse motion                 */
static int loose = 0;        /* loose to move                         */
static double y = 0.0;       /* integrated position falling mass      */
static double vy = 0.0;      /* integrated velocity falling mass      */
static double ay = 0.0;      /* integrating acceleration falling mass */
static double s = 0.0;       /* integrated position spring mass       */
static double vs = 0.0;      /* integrated velocity spring mass       */
static double as = 0.0;      /* integrating acceleration spring mass  */
static int stuck = 0;        /* masses stuck                          */
static int double_k1 = 0;    /* yellow buttons double spring k        */
static int double_mass = 0;  /*                double falling mass    */
static int disc = 0;         /*                disc or flat mass      */
static int do_elastic = 1;   /*                collision type         */
static int mass2 = 0;        /*                mass on spring         */
static int double_mass2 = 0; /*                double spring mass     */
static int change_dt = 0;    /*                change step size       */
static double dt_factor = 1.0;
static int double_buffer = 1;
static double t = 0.0;       /* time since reset                      */
static debug = 0;            /* debug print to stdout, redirect       */
static double scale = 65.0;  /* graphics scale                        */
static GLUquadricObj *obj;

static void reset()
{
  if(debug>2) printf("reset\n");
  t = 0.0;
  xc = 100;      /* graphical coordinates of falling mass */
  yc = mass_zer;
  xs = 130;      /* graphical coordinates of spring */
  ys = spring_top;
  /*   spring_bot */
  y  = 0.0;      /* state of falling mass */
  vy = 0.0;
  ay = 0.0;
  s  = 0.0;      /* state of top of springs mass */
  vs = 0.0;
  as = 0.0;
  stuck = 0;
}

static int in_rect(int x, int y, int w, int h)
{
  return xm>x && xm<x+w && (700-ym)>y && (700-ym)<y+h;
}

static void mouseMotion(int x, int y)
{
  xm = x;
  ym = y;
  bm = 1;  
  if(tracking)
  {
    /* printf("mouse motion x=%d, y=%d \n", x, y); */
    xc = xc + (xm-xp); /* mouse moved (new-old) */
    xp = xm;
    yp = ym;
  }
  glutPostRedisplay();  
}

static void mouseButton(int button, int state, int x, int y)
{
  xm = x;
  ym = y;
  bm = 1;  

  if(state==GLUT_DOWN)
  {
    /* printf("mouse down x=%d, y=%d, track=%d \n", x, y, tracking); */
    if(in_rect(xc-75, yc-5, 150, 10))
    {
      tracking = 1;
      loose = 0;
      xp = xm;
      yp = ym;
    }
    if(in_rect(300,640,280,50))
    {
      double_k1 = 1-double_k1;
      if(debug) if(double_k1) printf("double k1\n");
                else          printf("halve k1\n");
    }
    if(in_rect(300,580,280,50))
    {
      double_mass = 1-double_mass;
      if(debug) if(double_mass) printf("double mass\n");
                else            printf("halve mass\n");
    }
    if(in_rect(300,520,280,50))
    {
      disc = 1-disc;
      if(debug) if(disc) printf("show disc\n");
                else     printf("show bar\n");
    }
    if(in_rect(300,460,280,50))
    {
      do_elastic = 1-do_elastic;
      if(debug) if(do_elastic) printf("do elastic\n");
                else           printf("inelastic\n");
    }
    if(in_rect(300,400,280,50))
    {
      double_mass2 = 1-double_mass2;
      if(debug) if(double_mass2) printf("double mass2\n");
                else             printf("halve mass2\n");
    }
    if(in_rect(300,340,280,50))
    {
      mass2 = 1-mass2;
      if(debug) if(mass2) printf("mass2 on\n");
                else      printf("mass2 off\n");
    }
    if(in_rect(300,280,280,50))
    {
      dt_factor = dt_factor/2.0;
      if(debug) printf("halve dt\n");
    }
    if(in_rect(300,220,280,50))
    {
      dt_factor = dt_factor*2.0;
      if(debug) printf("double dt\n");
    }
  }
  if(state==GLUT_UP)
  {
    /* printf("mouse up x=%d, y=%d, track=%d \n", x, y, tracking); */
    tracking = 0;
    loose = 1;
    reset();
  }
  glutPostRedisplay();
}

static void inelastic(double m1, double v1, double m2, double v2, double * vf)
{
  /* compute final velocity of an inelastic collision */
  /* kinetic energy may not be conserved, masses stick together */
  *vf = (m1*v1 + m2*v2) / (m1+m2);
} /* end inelastic */

static void elastic(double m1, double v1, double m2, double v2,
             double * v1f, double * v2f)
{
  /* computer final velocities of an inelastic collision */
  /* kinetic energy is conserved */
  *v1f = ((m1-m2)/(m1+m2))*v1 + ((2.0*m2)/(m1+m2))*v2;
  *v2f = ((m2-m1)/(m1+m2))*v2 + ((2.0*m1)/(m1+m2))*v1;
} /* end elastic */

static void physicsComputation(void)
{
  double f;
  double dt = 0.02;   /* can be doubled or halved    */
  double mass = 1.0;  /* dropping mass               */
  double dmass;
  double smass = 1.0; /* mass on spring              */
  double dsmass;
  double k = 1.0;     /* spring constant             */
  double kk;
  double g = 1.0;     /* acceleration due to gravity */
  double ddt, vvs, vvy, vf;
  
  if(!loose)
  {
    reset();
    return;
  }

  kk = k;
  if(double_k1) kk = 2.0 * k;
  dmass = mass;
  if(double_mass) dmass = 2.0*mass;
  dsmass = smass;
  if(double_mass2) dsmass = 2.0*smass;
  ddt = dt_factor*dt;
  t = t + ddt;
  
  if(mass2==1 && yc<=ys && stuck==0) /* collision */
  {
    if(do_elastic)
    {
      elastic(dmass, vy, dsmass, vs, &vvy, &vvs);
      if(debug) printf("elastic vy=%f, vynew=%f, vs=%f, vsnew=%f\n",
                       vy, vvy, vs, vvs);
      vy = vvy;
      vs = vvs;
      yc = ys +1; /* fudge to prevent re-collision */
      y  = y + 2.0/scale;
    }
    else /* inelastic */
    {
      stuck = 1;
      inelastic(dmass, vy, dsmass, vs, &vf);
      if(debug) printf("inelastic vy=%f, vs=%f, vf=%f\n",
                       vy, vs, vf);
      vy = vf;
      vs = vf;
    }
  }

  if(stuck) /* use dmass+dsmass coordinates stay together */
  {
    f = kk * (spring_top-ys)/20.0; /* opposite g force */
    as = (-g) + f/(dmass+dsmass);
    vs = vs + as * ddt; /* v = integral a dt */
    s = s + vs * ddt;   /* s = integral v dt */
    ys = (int)(spring_top+(s*scale));
    f = kk * (spring_top-yc)/20.0; /* opposite g force */
    ay = (-g) + f/(dmass+dsmass);
    vy = vy + ay * ddt; /* v = integral a dt */
    y = y + vy * ddt;   /* s = integral v dt */
    yc = (int)(mass_zer+(y*scale));
  }
  else /* not stuck */
  {
    /* mass on spring */
    if(mass2==1 && yc>ys)
    {
      f = kk * (spring_top-ys)/20.0; /* opposite g force */
      as = (-g) + f/dsmass;
      vs = vs + as * ddt; /* v = integral a dt */
      s = s + vs * ddt;   /* s = integral v dt */
      ys = (int)(spring_top+(s*scale));
    }

    if(yc>ys) /* dropping mass above spring */
    {
      ay = -g;            /* f = ma  or  a = f / m  */
      vy = vy + ay * ddt; /* v = integral a dt      */
      y = y + vy * ddt;   /* y = integral v dt      */
      yc = (int)(mass_zer+(y*scale));
    }
    if(mass2==0 && yc<=ys) /* dropping mass on spring */
    {
      f = kk * (spring_top-ys)/20.0; /* opposite g force */
      ay = (-g) + f/dmass;
      vy = vy + ay * ddt;  /* v = integral a dt */
      y = y + vy * ddt;    /* y = integral v dt */
      yc = (int)(mass_zer+(y*scale));
      ys = yc;
      if(ys>spring_top) ys = spring_top; 
    }
    if(mass2==1 && yc<=ys) /* dropping mass on spring */
    {
      f = kk * (spring_top-ys)/20.0; /* opposite g force */
      as = (-g) + f/(dmass+dsmass);
      vs = vs + as * ddt; /* v = integral a dt */
      s = s + vs * ddt;   /* s = integral v dt */
      ys = (int)(spring_top+(s*scale));
      f = kk * (spring_top-yc)/20.0; /* opposite g force */
      ay = (-g) + f/(dmass+dsmass);
      vy = vy + ay * ddt; /* v = integral a dt */
      y = y + vy * ddt;   /* s = integral v dt */
      yc = (int)(mass_zer+(y*scale));
    }
  }
  if(debug>1)printf("t=%f, yc=%d, ys=%d, ay=%f, vy=%f, y=%f, as=%f, vs=%f, s=%f\n",
                    t, yc, ys, ay, vy, y, as, vs, s);
  if(yc<spring_bot) { yc=spring_bot+1; y=0.0; vy=0.0; ay=0.0;}
  if(ys<spring_bot) { ys=spring_bot; s=0.0; vs=0.0; as=0.0;}
  if(yc>mass_zer) { yc=mass_zer; y=0.0; vy=0.0; ay=0.0;}
  if(ys>mass_zer) { ys=mass_zer-1; s=0.0; vs=0.0; as=0.0;}
  glutPostRedisplay();
} /* end physicsComputation */

static void reshape(int w, int h)
{
  glViewport(0, 0, w, h);
  winWidth = w;
  winHeight = h;
  glutPostRedisplay();
} /* end reshape */

static void printstring(int x, int y, char *string)
{
   int len, i;

   glRasterPos2i(x, y);
   len = (int) strlen(string);
   for (i = 0; i < len; i++)
      glutBitmapCharacter(GLUT_BITMAP_HELVETICA_12, string[i]);
} /* end printstring */

static void display(void)
{
  double x=0.0, xx=0.0;
  double y=0.0, yy=0.0;
  double xn,yn;
  double a;
  double b;
  double t=0.0;
  double fs;    /* for drawing spring */
  int i;
  
  a = 1.0/(4.0*M_PI);  /* for drawing spring */
  b = 2.0*a;

  glClear(GL_COLOR_BUFFER_BIT);

  glColor3f(0.0, 0.0, 0.0);
  printstring(20, 680, "drop disc on spring to see motion");
  printstring(20, 140, "Basic physics:"); /* 20, 310 */
  printstring(20, 120, "f = k x  the force from a spring is proportional to distance compressed times spring constant, k");
  printstring(20, 100, "f = m a  the force on a body of mass, m, is proportional to the acceleration, a");
  printstring(20,  80, "acceleration a = dv / dt    velocity v = dx / dt");
  printstring(20,  60, "Momentum p = m v is conserved.");
  printstring(20,  40, "kinetic energy is conserved upon elastic collision");
  printstring(20,  20, "kinetic energy is not conserved upon inelastic collision");
  printstring(340, spring_bot, "any change resets simulation");

  glColor3f(0.0, 0.0, 0.0); /* box */
  glLineWidth(1.0);
  glBegin(GL_LINE_LOOP); 
    glVertex2f(20, spring_bot);
    glVertex2f(20, 660);
    glVertex2f(180,660);
    glVertex2f(180,spring_bot);
  glEnd();

  glColor3f(1.0, 0.0, 0.0); /* falling mass */
  glBegin(GL_LINES); 
    glVertex2f(20, mass_zer);
    glVertex2f(40, mass_zer);
  glEnd();
  if(disc)
  {
    if(double_mass) glLineWidth(3.0);
    else            glLineWidth(1.0);
    glTranslatef(xc, yc+20, 0);
    gluDisk(obj, 20, 20, 20, 20);
    glTranslatef(-xc, -(yc+20), 0);
  }
  if(!disc)
  {
    if(!double_mass)
    {
      glLineWidth(1.0);
      glBegin(GL_LINE_LOOP); 
        glVertex2f(xc-34, yc);
        glVertex2f(xc-34, yc+8);
        glVertex2f(xc+34, yc+8);
        glVertex2f(xc+34, yc);
      glEnd();
    }
    else
    {
      glLineWidth(3.0);
      glBegin(GL_LINE_LOOP); 
        glVertex2f(xc-32, yc);
        glVertex2f(xc-32, yc+6);
        glVertex2f(xc+32, yc+6);
        glVertex2f(xc+32, yc);
      glEnd();
    }
  }

  if(mass2)  /* mass on spring */
  {
    if(!double_mass2)
    {
      glLineWidth(1.0);
      glBegin(GL_LINE_LOOP); 
        glVertex2f(xc-34, ys-8);
        glVertex2f(xc-34, ys);
        glVertex2f(xc+34, ys);
        glVertex2f(xc+34, ys-8);
      glEnd();
    }
    else
    {
      glLineWidth(3.0);
      glBegin(GL_LINE_LOOP); 
        glVertex2f(xc-32, ys-6);
        glVertex2f(xc-32, ys);
        glVertex2f(xc+32, ys);
        glVertex2f(xc+32, ys-6);
      glEnd();
    }
  
  }

  glColor3f(0.0, 0.0, 1.0); /* spring ends */
  glPointSize(5.0);
  glBegin(GL_POINTS); 
    glVertex2f(xc, yc);   /* center of falling mass */
    glVertex2f(130, spring_bot);
    glVertex2f(xs, ys); /* should be moving end */
  glEnd();

  /* draw spring vertical */
  glColor3f(0.0, 0.0, 1.0);
  glBegin(GL_LINES); 
    glVertex2f(20, spring_top);
    glVertex2f(40, spring_top);
  glEnd();
  fs = (double)(ys-spring_bot)/643.0; /* fs<0.5 or fs>1.5 not very visible */
  t = 0.0;    
  y = 0.0;
  x = 0.0;
  if(double_k1) glLineWidth(2.0);
  if(!double_k1) glLineWidth(1.0);
  while(t<32.0*M_PI) /* heuristic, more or less loops */
  {
    t = t+0.1; /* 0.05 flicker, depends on speed of computer */
    xn = -(b - b*cos(t));
    yn = fs*a*t + b*sin(t);
    /* 200 is width scale, 80 is height scale */
    glBegin(GL_LINES); 
      glVertex2f(130+(int)(x *200.0), spring_bot+(int)(y *80.0));
      glVertex2f(130+(int)(xn*200.0), spring_bot+(int)(yn*80.0));
    glEnd();
    x = xn;
    y = yn;
  }    

  glColor3f(0.0, 1.0, 1.0); /* button */
  glBegin(GL_POLYGON);
    glVertex2f(300, 640);
    glVertex2f(580, 640);
    glVertex2f(580, 690);
    glVertex2f(300, 690);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  if(!double_k1) printstring(320, 670, "double spring k");
  if(double_k1)  printstring(320, 670, "halve spring k");
  printstring(320, 650, "note frequency change");

  glColor3f(0.0, 1.0, 1.0); /* button */
  glBegin(GL_POLYGON);
    glVertex2f(300, 580);
    glVertex2f(580, 580);
    glVertex2f(580, 630);
    glVertex2f(300, 630);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  if(!double_mass) printstring(320, 610, "double dropped mass");
  if(double_mass)  printstring(320, 610, "halve dropped mass");
  printstring(320, 590, "note frequency change");

  glColor3f(0.0, 1.0, 1.0); /* button */
  glBegin(GL_POLYGON);
    glVertex2f(300, 520);
    glVertex2f(580, 520);
    glVertex2f(580, 570);
    glVertex2f(300, 570);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  if(disc)  printstring(320, 530, "Change disc to rectangle");
  if(!disc) printstring(320, 530, "Change rectangle to disc");

  glColor3f(0.0, 1.0, 1.0); /* button */
  glBegin(GL_POLYGON);
    glVertex2f(300, 460);
    glVertex2f(580, 460);
    glVertex2f(580, 510);
    glVertex2f(300, 510);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  if(do_elastic)  printstring(320, 470, "Change to inelastic collision ");
  if(!do_elastic) printstring(320, 470, "Change to elastic collision");

  glColor3f(0.0, 1.0, 1.0); /* button */
  glBegin(GL_POLYGON);
    glVertex2f(300, 400);
    glVertex2f(580, 400);
    glVertex2f(580, 450);
    glVertex2f(300, 450);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  if(double_mass2)  printstring(320, 410, "halve the mass on the spring ");
  if(!double_mass2) printstring(320, 410, "double the mass on the spring ");

  glColor3f(0.0, 1.0, 1.0); /* button */
  glBegin(GL_POLYGON);
    glVertex2f(300, 340);
    glVertex2f(580, 340);
    glVertex2f(580, 390);
    glVertex2f(300, 390);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  if(mass2)  printstring(320, 350, "remove the mass from the spring ");
  if(!mass2) printstring(320, 350, "add a mass onto the spring ");
  
  glColor3f(0.0, 1.0, 1.0); /* button */
  glBegin(GL_POLYGON);
    glVertex2f(300, 280);
    glVertex2f(580, 280);
    glVertex2f(580, 330);
    glVertex2f(300, 330);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  printstring(320, 290, "double delta t (slow down)");

  glColor3f(0.0, 1.0, 1.0); /* button */
  glBegin(GL_POLYGON);
    glVertex2f(300, 220);
    glVertex2f(580, 220);
    glVertex2f(580, 270);
    glVertex2f(300, 270);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  printstring(320, 230, "halve delta t  (speed up)");
  
  glFlush();  
  if(double_buffer) glutSwapBuffers(); /* only for double */
} /* end display */

int main(int argc, char *argv[])
{
  glutInit(&argc, argv);
  if(argc>1) if(argv[1][0]=='s') double_buffer = 0;
  if(double_buffer)  glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
  if(!double_buffer) glutInitDisplayMode(GLUT_RGB); /* no double buffer */
  glutInitWindowSize(winWidth, winHeight);
  glutInitWindowPosition(100,100); 
  glutCreateWindow(argv[0]);
  glutReshapeFunc(reshape);
  glutDisplayFunc(display);
  glutMouseFunc(mouseButton);
  glutMotionFunc(mouseMotion);
  glutIdleFunc(physicsComputation);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluOrtho2D(0.0, (GLfloat)winWidth, 0.0, (GLfloat)winHeight);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
  glClearColor(1.0, 1.0, 1.0, 0.0);
  obj = gluNewQuadric(); /* used for disc */
  gluQuadricDrawStyle(obj, GLU_LINE);
  glutMainLoop();
  return 0;
} /* end main */
/* end spring3gl.c */