#include <stdio.h>
#include <string.h>
#include <math.h>
#include <GL/glut.h>

GLint width=480, height=240;
GLint mainWindow;

const double TOL = 1E-6;

// Important Constants
const double PI = 3.14159265358979323846;
const double CD = 0.5;
const double RHO = 1.29;
const double g = 9.81, G = 6.6726E-11;

// Timestep, window settings
const double DT = 0.001;
double xmax=1000, ymax=500;
double xmin = xmax/-10.0, ymin = ymax/-10.0;
double dx = (xmax-xmin)/width, dy = (ymax-ymin)/height;

// Projectile and Target settings
double xtarget=xmax*0.9, rtarget = xmax*0.025;
double rball = 0.1;
double aball = PI*rball*rball;
double mball = 1E2;

// Initial position/velocity values
double xp0=0.0, yp0=0.0, v0=100.0, theta0=45;
double xpos=xp0, ypos=yp0, vlast[2]={v0*cos(theta0),v0*sin(theta0)}, vnew[2];

// OpenGL Window functions
void display();
void reshape(int w, int h);
void dispGraph();

// Helper functions for OpenGL Window functions
void init();
void dispString(int x, int y, int z, char *string, void *font);
void dispAxes();

// Math-ish stuff for the calculations
void normalize(double *v);
double abs(double a);

int main(int argc, char **argv)
{
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_RGB);
  glutInitWindowSize(width, height);

  mainWindow = glutCreateWindow("Projectile Fun, Part I");
  init();
  glutDisplayFunc(display);
  glutReshapeFunc(reshape);

  glutMainLoop();

  return 0;
}

void reshape(int w, int h)
{
  // My Favorite Evil Technique (tm) -
  //   if they want to resize the window, too bad.

  glutSetWindow(mainWindow);
  glutReshapeWindow(width, height);

  glFlush();
  glutPostRedisplay();
}

void display()
{
  double force[2], d[2], dmag2, fmag, vmag, vang, i;

  glClear(GL_COLOR_BUFFER_BIT);

  // DRAW AXES
  dispAxes();
  //DRAW TARGET
  glColor3f(0.5, 0.5, 1.0);
  glBegin(GL_QUADS);
    glVertex3f(xtarget-rtarget, -3*dy, 0);
    glVertex3f(xtarget-rtarget, 3*dy, 0);
    glVertex3f(xtarget+rtarget, 3*dy, 0);
    glVertex3f(xtarget+rtarget, -3*dy, 0);
  glEnd();

  // The Yellow Path - gravity plus air resistance
  glColor3f(1.0, 1.0, 0.5);
  glBegin(GL_LINE_STRIP);
    xpos=xp0, ypos=yp0, vnew[0]=v0*cos(theta0), vnew[1]=v0*sin(theta0);
    do {
      vlast[0]=vnew[0], vlast[1]=vnew[1];

      // Update position based on velocity
      if (ypos+vlast[1]*DT < 0)
      {
	xpos += (vlast[0] * abs(ypos/vlast[1]));
	ypos = 0;
      }
      else
	xpos += vlast[0]*DT, ypos += vlast[1]*DT;
      glVertex3f(xpos, ypos, 0);
      
      // GRAVITY - EARTH
      vnew[1] -= g*DT;

      // DRAG FORCE
      vmag = hypot(vlast[0], vlast[1]);
      vang = atan2(vlast[1], vlast[0]);
      fmag =  CD*RHO*aball*(vmag*vmag)/2.0;
      vnew[0] -= (fmag/mball)*DT*cos(vang);
      vnew[1] -= (fmag/mball)*DT*sin(vang);

    } while (ypos > 0);
  glEnd();

  // The Green Path - just gravity
  glColor3f(0.5, 1.0, 0.5);
  glBegin(GL_LINE_STRIP);
    xpos=xp0, ypos=yp0, vnew[0]=v0*cos(theta0), vnew[1]=v0*sin(theta0);
    do {
      vlast[0]=vnew[0], vlast[1]=vnew[1];

      // Update position based on velocity
      if (ypos+vlast[1]*DT < 0)
      {
	xpos += (vlast[0] * abs(ypos/vlast[1]));
	ypos = 0;
      }
      else
	xpos += vlast[0]*DT, ypos += vlast[1]*DT;
      glVertex3f(xpos, ypos, 0);
      
      // GRAVITY - EARTH
      vnew[1] -= g*DT;

    } while (ypos > 0);
  glEnd();

}

void init()
{
  // Window settings
  glClearColor(0.0, 0.0, 0.0, 0.0);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glOrtho(xmin, xmax, ymin, ymax, -1.0, 1.0);
}

void dispString(int x, int y, int z, char *string, void *font)
{
  // Display a string. I'm not sure why I kept this in here,
  // but if I ever decide to stick a string in, I can.
  int len, i;
  glRasterPos3f(x, y, 0);
  len = (int) strlen(string);
  for (i = 0; i < len; i++)
    glutBitmapCharacter(font, string[i]);
}

void dispAxes()
{
  // Draw the axes and tick marks along them.

  double i;
  glColor3f(0.35, 0.35, 0.35);
  glBegin(GL_LINES);
   glVertex3f(xmin, 0.0, 0.0);
   glVertex3f(xmax, 0.0, 0.0);
   glVertex3f(0.0, ymin, 0.0);
   glVertex3f(0.0, ymax, 0.0);
   for (i = 0; i <= xmax+TOL; i += (xmax / 4))
   {
     glVertex3f(i, 0.0, 0.0);
     glVertex3f(i, 4.0*dy, 0.0);
   }
   for (i = 0; i <= ymax+TOL; i += (ymax / 4))
   {
     glVertex3f(0.0, i, 0.0);
     glVertex3f(4.0*dx, i, 0.0);
   }
  glEnd();
}

void normalize(double *v)
{
  // Find a unit vector with the same direction as v.
  double dist = hypot(v[0], v[1]);
  v[0] /= dist;
  v[1] /= dist;
}

double abs(double a)
{
  // the stdlib.h abs() only seems to like ints. ~shrug~
  if (a >= 0) return a;
  return (-1*a);
}

---

	My Supercomp front page.
	This page was created by Gary Sivek for 7th Period Supercomp.