#include <windows.h>   // use as needed for your system
#include <gl/Gl.h>
#include <gl/glu.h>
#include <gl/glut.h>
#include <math.h>

#define PI 3.14159265

GLfloat light_diffuse[] = {1.0, 1.0, 1.0, 1.0};  /* White diffuse light. */
GLfloat light_position[] = {10.0, 10.0, 10.0, 1.0};  /* Light location. */
GLfloat dome_mat[] = {0.f, .5f, 1.f, 1.f}; /* Material for the dome. */

/* makeMeshNewell is a generic surface drawing routine. Pass in three
   coordinate functions for the parametric representation of the
   surface and a number of steps to use in parameter space. This routine
   uses the Newell method to compute a normal. */
void makeMeshNewell(double x(double,double),double y(double,double),
			  double z(double,double),int steps)
{
  double h,u,v;
  double x0,x1,x2,x3;
  double y0,y1,y2,y3;
  double z0,z1,z2,z3;
  double nx,ny,nz;
  
  h = 1.0/steps; /* Step size */
  
  glBegin(GL_QUADS);
  for(u = 0.0;u < 1.0;u += h)
    for(v = 0.0;v < 1.0;v += h)
    {
	/* Compute the vertices from the parameterization. */
    x0 = x(u,v);
    y0 = y(u,v);
    z0 = z(u,v);
    
    x1 = x(u,v+h);
    y1 = y(u,v+h);
    z1 = z(u,v+h);

    x2 = x(u+h,v+h);
    y2 = y(u+h,v+h);
    z2 = z(u+h,v+h);

    x3 = x(u+h,v);
    y3 = y(u+h,v);
    z3 = z(u+h,v);
     
	/* Use the Newell method to compute the normal. */
	nx = (y0-y1)*(z0+z1)+(y1-y2)*(z1+z2)+(y2-y3)*(z2+z3)+(y3-y0)*(z3+z0);
	ny = (z0-z1)*(x0+x1)+(z1-z2)*(x1+x2)+(z2-z3)*(x2+x3)+(z3-z0)*(x3+x0);    
	nz = (x0-x1)*(y0+y1)+(x1-x2)*(y1+y2)+(x2-x3)*(y2+y3)+(x3-x0)*(y3+y0);
	glNormal3f( nx, ny, nz);
	glVertex3f( x0, y0, z0 );
	glVertex3f( x1, y1, z1 );
	glVertex3f( x2, y2, z2 );
	glVertex3f( x3, y3, z3 );
	}
	glEnd();
}

/* makeMeshExact is a drawing routine for drawing a surface for which you 
   know the normals you want to use for each vertex. */
void makeMeshExact(double x(double,double),double y(double,double),
			  double z(double,double),double nx(double,double),
			  double ny(double,double),double nz(double,double),int steps)
{
  double h,u,v;
  double v_x,v_y,v_z;
  double n_x,n_y,n_z;
  
  h = 1.0/steps; /* Step size */
  
  glBegin(GL_QUADS);
  for(u = 0.0;u < 1.0;u += h)
    for(v = 0.0;v < 1.0;v += h)
    {
	/* Compute the vertices and normals from the parameterization. */
    v_x = x(u,v);
    n_x = nx(u,v);
    v_y = y(u,v);
	n_y = ny(u,v);
    v_z = z(u,v);
	n_z = nz(u,v);
	glNormal3f( n_x, n_y, n_z);
	glVertex3f( v_x, v_y, v_z );
   
    v_x = x(u,v+h);
    n_x = nx(u,v+h);
    v_y = y(u,v+h);
	n_y = ny(u,v+h);
    v_z = z(u,v+h);
	n_z = nz(u,v+h);
	glNormal3f( n_x, n_y, n_z);
	glVertex3f( v_x, v_y, v_z );

    v_x = x(u+h,v+h);
    n_x = nx(u+h,v+h);
    v_y = y(u+h,v+h);
	n_y = ny(u+h,v+h);
    v_z = z(u+h,v+h);
	n_z = nz(u+h,v+h);
	glNormal3f( n_x, n_y, n_z);
	glVertex3f( v_x, v_y, v_z );

    v_x = x(u+h,v);
    n_x = nx(u+h,v);
    v_y = y(u+h,v);
	n_y = ny(u+h,v);
    v_z = z(u+h,v);
	n_z = nz(u+h,v);
	glNormal3f( n_x, n_y, n_z);
	glVertex3f( v_x, v_y, v_z );
	}
	glEnd();
}

/* These three coordinate functions are used to render the dome shape. */
double pointedDomeX(double u,double v) 
{
   return (2*cos(PI*u/4.0)-sqrt(2.0))*cos(2.0*PI*v);
}

double pointedDomeY(double u,double v) 
{
   return 2*sin(PI*u/4.0);
}

double pointedDomeZ(double u,double v) 
{
   return (-2*cos(PI*u/4.0)+sqrt(2.0))*sin(2.0*PI*v);
}

/* These three functions are used to compute the dome shape normals. */
double pointedDomeNX(double u,double v) 
{
   return (2*cos(PI*u/4.0))*cos(2.0*PI*v);
}

double pointedDomeNY(double u,double v) 
{
   return 2*sin(PI*u/4.0);
}

double pointedDomeNZ(double u,double v) 
{
   return (-2*cos(PI*u/4.0))*sin(2.0*PI*v);
}

void pointedDome(int steps,bool newell) 
{
  /* Set material properties for the dome. */
  glMaterialfv(GL_FRONT, GL_DIFFUSE, dome_mat);

  if(newell)
	makeMeshNewell(pointedDomeX,pointedDomeY,pointedDomeZ,steps);
  else
	makeMeshExact(pointedDomeX,pointedDomeY,pointedDomeZ,
		pointedDomeNX,pointedDomeNY,pointedDomeNZ,steps);
}

void
resize(int w,int h)
{
  glViewport(0,0,(GLsizei) w,(GLsizei) h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective( /* field of view in degree */ 40.0,
    /* aspect ratio */ (float) w/(float) h,
    /* Z near */ 5.0, /* Z far */ 50.0);
  glutPostRedisplay();
}

void
display(void)
{
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  /* Be sure to set the matrix mode before drawing. */
  glMatrixMode(GL_MODELVIEW);

  /* Draw four copies of the basic dome shape translated to the four corners of
     a square. */
   glPushMatrix();
   glTranslatef(2.0, 0.0, -2.0);
   pointedDome(20,true);
   glPopMatrix();

   glPushMatrix();
   glTranslatef(-2.0, 0.0, -2.0);
   pointedDome(20,true);
   glPopMatrix();

   glPushMatrix();
   glTranslatef(2.0, 0.0, 2.0);
   pointedDome(20,false);
   glPopMatrix();

   glPushMatrix();
   glTranslatef(-2.0, 0.0, 2.0);
   pointedDome(20,false);
   glPopMatrix(); 

  glutSwapBuffers();
}

void
init(void)
{
  /* Enable a single OpenGL light. */
  glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse);
  glLightfv(GL_LIGHT0, GL_POSITION, light_position);
  glEnable(GL_LIGHT0);
  glEnable(GL_LIGHTING);

  /* Turn on automatic normalization. */
  glEnable(GL_NORMALIZE);

  /* Use depth buffering for hidden surface elimination. */
  glEnable(GL_DEPTH_TEST);

  /* Make sure smooth shading is turned on. */
  glShadeModel(GL_SMOOTH);

  /* Set the clear color to white. */
  glClearColor(1.0, 1.0, 1.0, 1.0);

  /* Setup the view of the cube. */
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective( /* field of view in degree */ 40.0,
    /* aspect ratio */ 1.0,
    /* Z near */ 5.0, /* Z far */ 50.0);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
  gluLookAt(0.0, 2.0, 10.0,  /* eye is at (0,2,10) */
    0.0, 0.0, 0.0,      /* center is at (0,0,0) */
    0.0, 1.0, 0.);      /* up is in positive Y direction */

}

int
main(int argc, char **argv)
{
  glutInit(&argc, argv);
  glutInitWindowSize(512, 512);
  glutInitWindowPosition(100,100);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
  glutCreateWindow("Pointed domes");
  glutDisplayFunc(display);
  glutReshapeFunc(resize);
  
  init();
  glutMainLoop();
  return 0;             /* ANSI C requires main to return int. */
}