#include<GL/glut.h>
#include<iostream.h>
#include<math.h>
#include<stdlib.h>

double v, ang, tdist=600, range, t=0, tstep=0.1, a=32;

const double PI=3.141592653, G=10;

const double objx = 500, objy = 200, objm = 150, pm = 10;

int xsize=300, ysize=300;

void init() {
  glClearColor(0,0,0,0);
  glShadeModel(GL_FLAT);
}

void display() {
  double x, y, vx, vy, v, t=0, ang;

  glClear(GL_COLOR_BUFFER_BIT);
  glLoadIdentity();

  for(int i=0; i<xsize*1.1; i++) 
    for(int j=0; j<ysize*1.1; j++) {
      v=tdist*double(i)/xsize/1.1;
      ang=90.0*double(j)/ysize/1.1;
      vx=v*cos(ang*PI/180);
      vy=v*sin(ang*PI/180);
      x=y=0;  
      glBegin(GL_POINTS);
 
      do {
        t+=tstep;
  
      // Comment the next two lines out to remove the point force
        vx+=G*objm*pm*(objx-x)/pow(pow(objx-x,2)+pow(objy-y,2),3.0/2);
        vy+=G*objm*pm*(objy-y)/pow(pow(objx-x,2)+pow(objy-y,2),3.0/2);
  
        vy-=a*tstep;
   
        y+=vy*tstep;
        x+=vx*tstep;

      } while (y>=0);
      if (fabs(x-tdist) < tdist) {
        glColor3f(1-fabs(x-tdist)/tdist, 1-fabs(x-tdist)/tdist, 0);
        glVertex3f(v, ang, 0);
      }
    glEnd();
    }
}

void reshape(int w, int h) {
  xsize=w; ysize=h;
  glViewport(0, 0, GLsizei(w), GLsizei(h));
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glOrtho(0, tdist, 0, 90, -1, 1);
  glMatrixMode(GL_MODELVIEW);
}

void main(int argc, char* argv[]) {
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
  glutInitWindowSize(xsize, ysize);
  glutInitWindowPosition(100, 100);
  glutCreateWindow("Projectile Plot");
  init();
  glutDisplayFunc(display);
  glutReshapeFunc(reshape);
  glutMainLoop();
}