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

#define MAXPTS 20
#define EPS 0.000001

int mousex, mousey;
int sirka=800, vyska=600;
int N,M,D1,D2;
double V[MAXPTS][MAXPTS][3],w[MAXPTS][MAXPTS];
int sel_pt_i=0,sel_pt_j=0,color=0,mode=0;
int is_ctrl,is_alt,is_shift;
double ZOOM=2.0,XOFFSET=-1.0,YOFFSET=0.0,ROTATEX=15.0,ROTATEY=0.0,STEP=0.05;
GLuint wlist1,wlist2,slist1,slist2;
int SHOW_PTS=1;

double ALPHA=0.3, BETA1=0.3, BETA2=0.7;
double MM[4][4];
double A[4][4][3],B[4][4][3],C[4][4][3];
//
//double s[100],t[100],N1[50][50],N2[50][50];

void PisVektorovyText(GLfloat x, GLfloat y, char *format,...) {
  va_list args;
  char buffer[200], *p;

  va_start(args, format); vsprintf(buffer, format, args); va_end(args); 
  
  glPushMatrix(); glLoadIdentity(); glTranslatef(-1.0,0.0,-6.0f); // glRotatef(10,1.0,0.0,0.0);
  glRotatef(0,0.0,1.0,0.0); glScalef(2,2,-2); glTranslatef(x, y, 0); glScalef(0.0007,0.0007,0.0007);
  for (p = buffer; *p; p++) glutStrokeCharacter(GLUT_STROKE_ROMAN, *p);
  glPopMatrix();
}
void SetVertexColor(double u, double v, int I, int J) { 
  double uu,vv; uu=u+I; u/=(N+1); vv=v+J; vv/=(M+1); glColor3d(uu,vv,1.0-uu-vv); 
}

//int ComputeAnalyticPoint(double u, double v, double *x, double *y, double *z){
//}

void DajBod(int i, int j, double *x, double *y, double *z){
  double x1,y1,z1,x2,y2,z2;
  
  if ((i>=0) && (i<=N) && (j>=0) && (j<=M)) { *x=V[i][j][0]; *y=V[i][j][1]; *z=V[i][j][2]; }
  if ((mode==1)||(mode==2)){
    if (fabs(ALPHA)<EPS) { *x=*y=*z=0; return; }
    if ((i==-1) && (j==-1)) {
      DajBod(-1,0,&x1,&y1,&z1); x2=(3-ALPHA)*x1/(2*ALPHA); y2=(3-ALPHA)*y1/(2*ALPHA); z2=(3-ALPHA)*z1/(2*ALPHA);
      DajBod(-1,1,&x1,&y1,&z1); x2+=(2*ALPHA-3)*x1/(2*ALPHA); y2+=(2*ALPHA-3)*y1/(2*ALPHA); z2+=(2*ALPHA-3)*z1/(2*ALPHA);
      DajBod(-1,2,&x1,&y1,&z1); x2+=(ALPHA)*x1/(2*ALPHA); y2+=(ALPHA)*y1/(2*ALPHA); z2+=(ALPHA)*z1/(2*ALPHA);
      *x=x2; *y=y2; *z=z2;
      return;
    }
    if ((i==N+1) && (j==M+1)) {
      DajBod(N,j,&x1,&y1,&z1); x2=(3-ALPHA)*x1/(2*ALPHA); y2=(3-ALPHA)*y1/(2*ALPHA); z2=(3-ALPHA)*z1/(2*ALPHA);
      DajBod(N-1,j,&x1,&y1,&z1); x2+=(2*ALPHA-3)*x1/(2*ALPHA); y2+=(2*ALPHA-3)*y1/(2*ALPHA);z2+=(2*ALPHA-3)*z1/(2*ALPHA);
      DajBod(N-2,j,&x1,&y1,&z1); x2+=(ALPHA)*x1/(2*ALPHA); y2+=(ALPHA)*y1/(2*ALPHA); z2+=(ALPHA)*z1/(2*ALPHA);
      *x=x2; *y=y2; *z=z2;
      return;
    }
    if (i==-1) {
      DajBod(0,j,&x1,&y1,&z1); x2=(3-ALPHA)*x1/(2*ALPHA); y2=(3-ALPHA)*y1/(2*ALPHA); z2=(3-ALPHA)*z1/(2*ALPHA);
      DajBod(1,j,&x1,&y1,&z1); x2+=(2*ALPHA-3)*x1/(2*ALPHA); y2+=(2*ALPHA-3)*y1/(2*ALPHA); z2+=(2*ALPHA-3)*z1/(2*ALPHA);
      DajBod(2,j,&x1,&y1,&z1); x2+=(ALPHA)*x1/(2*ALPHA); y2+=(ALPHA)*y1/(2*ALPHA); z2+=(ALPHA)*z1/(2*ALPHA);
      *x=x2; *y=y2; *z=z2;
      return;
    }
    if (j==-1) {
      DajBod(i,0,&x1,&y1,&z1); x2=(3-ALPHA)*x1/(2*ALPHA); y2=(3-ALPHA)*y1/(2*ALPHA); z2=(3-ALPHA)*z1/(2*ALPHA);
      DajBod(i,1,&x1,&y1,&z1); x2+=(2*ALPHA-3)*x1/(2*ALPHA); y2+=(2*ALPHA-3)*y1/(2*ALPHA); z2+=(2*ALPHA-3)*z1/(2*ALPHA);
      DajBod(i,2,&x1,&y1,&z1); x2+=(ALPHA)*x1/(2*ALPHA); y2+=(ALPHA)*y1/(2*ALPHA); z2+=(ALPHA)*z1/(2*ALPHA);
      *x=x2; *y=y2; *z=z2;
      return;
    }
    if (i==N+1) {
      DajBod(N,j,&x1,&y1,&z1); x2=(3-ALPHA)*x1/(2*ALPHA); y2=(3-ALPHA)*y1/(2*ALPHA); z2=(3-ALPHA)*z1/(2*ALPHA);
      DajBod(N-1,j,&x1,&y1,&z1); x2+=(2*ALPHA-3)*x1/(2*ALPHA); y2+=(2*ALPHA-3)*y1/(2*ALPHA);z2+=(2*ALPHA-3)*z1/(2*ALPHA);
      DajBod(N-2,j,&x1,&y1,&z1); x2+=(ALPHA)*x1/(2*ALPHA); y2+=(ALPHA)*y1/(2*ALPHA); z2+=(ALPHA)*z1/(2*ALPHA);
      *x=x2; *y=y2; *z=z2;
      return;
    }
    if (j==M+1) {
      DajBod(i,M,&x1,&y1,&z1); x2=(3-ALPHA)*x1/(2*ALPHA); y2=(3-ALPHA)*y1/(2*ALPHA); z2=(3-ALPHA)*z1/(2*ALPHA);
      DajBod(i,M-1,&x1,&y1,&z1); x2+=(2*ALPHA-3)*x1/(2*ALPHA); y2+=(2*ALPHA-3)*y1/(2*ALPHA);z2+=(2*ALPHA-3)*z1/(2*ALPHA);
      DajBod(i,M-2,&x1,&y1,&z1); x2+=(ALPHA)*x1/(2*ALPHA); y2+=(ALPHA)*y1/(2*ALPHA); z2+=(ALPHA)*z1/(2*ALPHA);
      *x=x2; *y=y2; *z=z2;
      return;
    }
  }
  if ((mode==3)||(mode==4)){
    if (i==-1) i=0; if (j==-1) j=0;
    if (i==N+1) i=N; if (j==M+1) j=M;
    *x=V[i][j][0]; *y=V[i][j][1]; *z=V[i][j][2];
  }
}

void ComputePoint(double u, double v, double *x, double *y, double *z){
  int i,j,k;
  double P[4][3],Q[3];
  for (i=0;i<4;i++) for (k=0;k<3;k++) P[i][k]=0;
  for (i=0;i<4;i++) for (k=0;k<3;k++) for (j=0;j<4;j++) P[i][k]+=pow(u,3-j)*C[j][i][k];
  for (k=0;k<3;k++) Q[k]=0;
  for (k=0;k<3;k++) for (j=0;j<4;j++) Q[k]+=P[j][k]*pow(v,3-j);
  *x=Q[0]; *y=Q[1]; *z=Q[2];
}

void DrawMeshSegmentWireframe(int I, int J) {
  int i,j,k,l;
  double u,v,x,y,z;

  for (i=-2;i<=1;i++) for (j=-2;j<=1;j++) DajBod(i+I,j+J,&A[i+2][j+2][0],&A[i+2][j+2][1],&A[i+2][j+2][2]);
  for (i=0;i<4;i++) for (j=0;j<4;j++) for (k=0;k<3;k++) B[i][j][k]=0.0;
  for (i=0;i<4;i++) for (j=0;j<4;j++) for (k=0;k<3;k++) C[i][j][k]=0.0;
  for (i=0;i<4;i++) for (j=0;j<4;j++) for (k=0;k<3;k++) for (l=0;l<4;l++) B[i][j][k]+=MM[i][l]*A[l][j][k];
  for (i=0;i<4;i++) for (j=0;j<4;j++) for (k=0;k<3;k++) for (l=0;l<4;l++) C[i][j][k]+=B[i][l][k]*MM[j][l];

  glBegin(GL_LINE_STRIP);
    for (v=0.0;v<=1.0+STEP-EPS;v+=STEP) { ComputePoint(0,v,&x,&y,&z); glVertex3d(x,y,z); }
  glEnd();
  glBegin(GL_LINE_STRIP);
    for (v=0.0;v<=1.0+STEP-EPS;v+=STEP) { ComputePoint(1,v,&x,&y,&z); glVertex3d(x,y,z); }
  glEnd();
  glBegin(GL_LINE_STRIP);
    for (v=0.0;v<=1.0+STEP-EPS;v+=STEP) { ComputePoint(v,0,&x,&y,&z); glVertex3d(x,y,z); }
  glEnd();
  glBegin(GL_LINE_STRIP);
    for (v=0.0;v<=1.0+STEP-EPS;v+=STEP) { ComputePoint(v,1,&x,&y,&z); glVertex3d(x,y,z); }
  glEnd();
}

void DrawMeshSegmentSolid(int I, int J) {
  int i,j,k,l;
  double u,v,x,y,z;

  for (i=-2;i<=1;i++) for (j=-2;j<=1;j++) DajBod(i+I,j+J,&A[i+2][j+2][0],&A[i+2][j+2][1],&A[i+2][j+2][2]);
  for (i=0;i<4;i++) for (j=0;j<4;j++) for (k=0;k<3;k++) B[i][j][k]=0.0;
  for (i=0;i<4;i++) for (j=0;j<4;j++) for (k=0;k<3;k++) C[i][j][k]=0.0;
  for (i=0;i<4;i++) for (j=0;j<4;j++) for (k=0;k<3;k++) for (l=0;l<4;l++) B[i][j][k]+=MM[i][l]*A[l][j][k];
  for (i=0;i<4;i++) for (j=0;j<4;j++) for (k=0;k<3;k++) for (l=0;l<4;l++) C[i][j][k]+=B[i][l][k]*MM[j][l];

  for (u=0.0;u<=1.0-EPS;u+=STEP) {
    glBegin(GL_QUAD_STRIP);
      for (v=0.0;v<=1.0+STEP-EPS;v+=STEP) {
        ComputePoint(u,v,&x,&y,&z);      SetVertexColor(u,v,I,J);      glVertex3d(x,y,z);
        ComputePoint(u+STEP,v,&x,&y,&z); SetVertexColor(u+STEP,v,I,J); glVertex3d(x,y,z);
      }
    glEnd();
  }
}

void GenerateCardinalWireframe(void){
  int i,j;
  
  MM[0][0]=-ALPHA;  MM[0][1]=2-ALPHA; MM[0][2]=ALPHA-2;   MM[0][3]=ALPHA;
  MM[1][0]=2*ALPHA; MM[1][1]=ALPHA-3; MM[1][2]=3-2*ALPHA; MM[1][3]=-ALPHA;
  MM[2][0]=-ALPHA;  MM[2][1]=0;       MM[2][2]=ALPHA;     MM[2][3]=0;
  MM[3][0]=0;       MM[3][1]=1;       MM[3][2]=0;         MM[3][3]=0;

  glNewList(wlist1,GL_COMPILE);
  glColor3d(0.0,0.0,0.7); glLineWidth(3.0);
  for (i=1;i<=N;i++) for (j=1;j<=M;j++) DrawMeshSegmentWireframe(i,j);
  glEndList();
}

void GenerateCardinalSolid(void){
  int i,j;
  
  MM[0][0]=-ALPHA;  MM[0][1]=2-ALPHA; MM[0][2]=ALPHA-2;   MM[0][3]=ALPHA;
  MM[1][0]=2*ALPHA; MM[1][1]=ALPHA-3; MM[1][2]=3-2*ALPHA; MM[1][3]=-ALPHA;
  MM[2][0]=-ALPHA;  MM[2][1]=0;       MM[2][2]=ALPHA;     MM[2][3]=0;
  MM[3][0]=0;       MM[3][1]=1;       MM[3][2]=0;         MM[3][3]=0;
  
  glNewList(slist1,GL_COMPILE);
  for (i=1;i<=N;i++) for (j=1;j<=M;j++) DrawMeshSegmentSolid(i,j);
  glEndList();
}

void GenerateBetaWireframe(void){
  int i,j;
  
  MM[0][0]=-2*pow(BETA1,3);                  MM[0][1]=2*(pow(BETA1,3)+pow(BETA1,2)+BETA1+BETA2);
  MM[0][2]=-2*(pow(BETA1,2)+BETA1+BETA2+1);  MM[0][3]=2;
  MM[1][0]=6*pow(BETA1,3);                   MM[1][1]=-3*(2*pow(BETA1,3)+2*pow(BETA1,2)+BETA2);
  MM[1][2]=3*(2*pow(BETA1,2)+BETA2);         MM[1][3]=0;
  MM[2][0]=-6*pow(BETA1,3);                  MM[2][1]=6*BETA1*(pow(BETA1,2)-1);
  MM[2][2]=6*BETA1;                          MM[2][3]=0;
  MM[3][0]=2*pow(BETA1,3);                   MM[3][1]=4*pow(BETA1,2)+4*BETA1+BETA2;
  MM[3][2]=2;                                MM[3][3]=0;
  for (i=0;i<4;i++) for (j=0;j<4;j++)
    MM[i][j] /= ( 2*pow(BETA1,3) + 4*pow(BETA1,2) + 4*BETA1 + BETA2 + 2 );
  
  glNewList(wlist2,GL_COMPILE);
  glColor3d(0.0,0.0,0.7); glLineWidth(3.0);
  for (i=1;i<=N;i++) for (j=1;j<=M;j++) DrawMeshSegmentWireframe(i,j);
  glEndList();
}

void GenerateBetaSolid(void){
  int i,j;
  
  MM[0][0]=-2*pow(BETA1,3);                  MM[0][1]=2*(pow(BETA1,3)+pow(BETA1,2)+BETA1+BETA2);
  MM[0][2]=-2*(pow(BETA1,2)+BETA1+BETA2+1);  MM[0][3]=2;
  MM[1][0]=6*pow(BETA1,3);                   MM[1][1]=-3*(2*pow(BETA1,3)+2*pow(BETA1,2)+BETA2);
  MM[1][2]=3*(2*pow(BETA1,2)+BETA2);         MM[1][3]=0;
  MM[2][0]=-6*pow(BETA1,3);                  MM[2][1]=6*BETA1*(pow(BETA1,2)-1);
  MM[2][2]=6*BETA1;                          MM[2][3]=0;
  MM[3][0]=2*pow(BETA1,3);                   MM[3][1]=4*pow(BETA1,2)+4*BETA1+BETA2;
  MM[3][2]=2;                                MM[3][3]=0;
  for (i=0;i<4;i++) for (j=0;j<4;j++)
    MM[i][j] /= ( 2*pow(BETA1,3) + 4*pow(BETA1,2) + 4*BETA1 + BETA2 + 2 );
  
  glNewList(slist2,GL_COMPILE);
  for (i=1;i<=N;i++) for (j=1;j<=M;j++) DrawMeshSegmentSolid(i,j);
  glEndList();
}

void display(void) {
  int i,j;

  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
  glTranslatef(XOFFSET,YOFFSET,-6.0f); // move into the screen 6.0
  glRotatef(ROTATEX,1.0,0.0,0.0);
  glRotatef(ROTATEY,0.0,1.0,0.0);
  glScalef(ZOOM,ZOOM,-ZOOM);

  glColor3d(1.0,1.0,1.0); glLineWidth(2.0);
  PisVektorovyText(-1.0,1.0,"N=%d M=%d alpha=%.2f beta1=%.2f beta2=%.2f point=[%d,%d] step=%.2f",
      N,M,ALPHA,BETA1,BETA2,sel_pt_i,sel_pt_j,STEP);

  if (SHOW_PTS) {
    glColor3d(0.5,0.5,0.5); glLineWidth(1.0);
    for (i=0;i<=N;i++) { glBegin(GL_LINE_STRIP); for (j=0;j<=M;j++) glVertex3dv(V[i][j]); glEnd(); }
    for (j=0;j<=M;j++) { glBegin(GL_LINE_STRIP); for (i=0;i<=N;i++) glVertex3dv(V[i][j]); glEnd(); }
  }

  if (mode==1) glCallList(wlist1);
  if (mode==3) glCallList(wlist2);
  if (mode==2) glCallList(slist1);
  if (mode==4) glCallList(slist2);
  
  if (SHOW_PTS) {
    glColor3d(0.5,0.5,0.5);
    glPointSize(5.0);
    glBegin(GL_POINTS); for (i=0;i<=N;i++) for (j=0;j<=M;j++) glVertex3dv(V[i][j]); glEnd();
    glColor3d(0.0,1.0,0.0);
    glPointSize(7.0);
    glBegin(GL_POINTS); glVertex3dv(V[sel_pt_i][sel_pt_j]); glEnd();
  }
  
  glutSwapBuffers();
}

void reshape(int w, int h) {
  sirka=w; vyska=h;
  glViewport(0, 0, w, h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective(45.0f,(GLfloat)sirka/(GLfloat)vyska,0.1f,100.0f);
  glMatrixMode(GL_MODELVIEW);
}

void GenerujDisplayList(void){
  switch (mode) {
    case 1: GenerateCardinalWireframe(); break;
    case 2: GenerateCardinalSolid(); break;
    case 3: GenerateBetaWireframe(); break;
    case 4: GenerateBetaSolid(); break;
  }
}

void keypressed(unsigned char key, int x, int y) { 
//  int i,j,k;

  is_ctrl=glutGetModifiers(); 
  is_alt=is_ctrl&GLUT_ACTIVE_ALT; 
  is_shift=is_ctrl&GLUT_ACTIVE_SHIFT; 
  is_ctrl=is_ctrl&GLUT_ACTIVE_CTRL;
  switch (key) {
    case 27:  exit(0); break;
    case '1': BETA1-=0.05; GenerujDisplayList(); break;
    case '!': BETA1+=0.05; GenerujDisplayList(); break;
    case '2': BETA2-=0.05; GenerujDisplayList(); break;
    case '@': BETA2+=0.05; GenerujDisplayList(); break;
    case 'a': ALPHA-=0.05; GenerujDisplayList(); break;
    case 'A': ALPHA+=0.05; GenerujDisplayList(); break;
    case 'm': if (M-2) { M--; sel_pt_i=sel_pt_j=0; } GenerujDisplayList(); break;
    case 'M': if (!M-MAXPTS+5) M++;                  GenerujDisplayList(); break;
    case 'n': if (N-2) { N--; sel_pt_i=sel_pt_j=0; } GenerujDisplayList(); break;
    case 'N': if (!N-MAXPTS+5) N++;                  GenerujDisplayList(); break;
    case '+': ZOOM+=STEP; break;
    case '-': ZOOM-=STEP; break;
    case 'S': STEP+=0.01;                GenerujDisplayList(); break;
    case 's': if (STEP>0.01) STEP-=0.01; GenerujDisplayList(); break;
    case 'p': SHOW_PTS^=1; break;
  }
  glutPostRedisplay();
}

void keypressed_special(int key, int x, int y) {
  is_ctrl=glutGetModifiers(); 
  is_alt=is_ctrl&GLUT_ACTIVE_ALT; 
  is_shift=is_ctrl&GLUT_ACTIVE_SHIFT; 
  is_ctrl=is_ctrl&GLUT_ACTIVE_CTRL;
 
  if (is_alt) {
    switch (key) {
      case GLUT_KEY_UP:    ROTATEX+=2; break;
      case GLUT_KEY_DOWN:  ROTATEX-=2; break;
      case GLUT_KEY_LEFT:  ROTATEY+=2; break;
      case GLUT_KEY_RIGHT: ROTATEY-=2; break;
    }
  }
  if (is_ctrl) {
    switch (key) {
      case GLUT_KEY_UP:    YOFFSET-=STEP; break;
      case GLUT_KEY_DOWN:  YOFFSET+=STEP; break;
      case GLUT_KEY_LEFT:  XOFFSET+=STEP; break;
      case GLUT_KEY_RIGHT: XOFFSET-=STEP; break;
    }
  }
  if (is_shift) {
    switch (key) {
      case GLUT_KEY_PAGE_UP:    V[sel_pt_i][sel_pt_j][1]+=0.025; break;
      case GLUT_KEY_PAGE_DOWN:  V[sel_pt_i][sel_pt_j][1]-=0.025; break;
      case GLUT_KEY_UP:         V[sel_pt_i][sel_pt_j][2]+=0.025; break;
      case GLUT_KEY_DOWN:       V[sel_pt_i][sel_pt_j][2]-=0.025; break;
      case GLUT_KEY_LEFT:       V[sel_pt_i][sel_pt_j][0]-=0.025; break;
      case GLUT_KEY_RIGHT:      V[sel_pt_i][sel_pt_j][0]+=0.025; break;
    }
    GenerujDisplayList();
  }
  if ((!is_alt)&&(!is_ctrl)&&(!is_shift)) {
    switch (key) {
      case GLUT_KEY_UP:    if (sel_pt_j<M) sel_pt_j++; break;
      case GLUT_KEY_DOWN:  if (sel_pt_j>0) sel_pt_j--; break;
      case GLUT_KEY_LEFT:  if (sel_pt_i>0) sel_pt_i--; break;
      case GLUT_KEY_RIGHT: if (sel_pt_i<N) sel_pt_i++; break;
    }
  }
  glutPostRedisplay();
}

void motion(int x, int y) { mousex = x; mousey = y; }
void mouse(int b, int s, int x, int y) { 
  mousex = x; mousey = y; 
  if ((b==3)||(b==4)) {
    w[sel_pt_i][sel_pt_j]+=-0.05+0.1*(4-b);
    GenerujDisplayList();
    glutPostRedisplay();
  }
  motion(x,y); 
}
void menu(int choice) { 
  mode=choice; 
  GenerujDisplayList();
  glutPostRedisplay(); 
}

int main(int argc, char **argv) {
  int i,j; 
  
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH );
  
  glutInitWindowSize(800, 600); glutInitWindowPosition(0, 0);
  glutCreateWindow("Miso Forisek -- Krivky a plochy v PG");
  glClearColor(0.0, 0.0, 0.0, 0.0); glColor3d(1.0, 0.0, 0.0); glLineWidth(1.0);
 
  glClearDepth(10.0);                           // Enables Clearing Of The Depth Buffer
  glDepthFunc(GL_LESS);                         // The Type Of Depth Test To Do
  glEnable(GL_DEPTH_TEST);                      // Enables Depth Testing
  glShadeModel(GL_SMOOTH);                      // Enables Smooth Color Shading

  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();                             // Reset The Projection Matrix
  gluPerspective(45.0f,(GLfloat)sirka/(GLfloat)vyska,0.1f,100.0f);     // Calculate The Aspect Ratio Of The Window
  glMatrixMode(GL_MODELVIEW);
 
  glutDisplayFunc(display);
  glutReshapeFunc(reshape);
  glutKeyboardFunc(keypressed);
  glutSpecialFunc(keypressed_special);
  glutMotionFunc(motion);
  glutMouseFunc(mouse);
  
  glutCreateMenu(menu); 
  glutAddMenuEntry("Grid only", 0); 
  glutAddMenuEntry("Cardinal mesh (wireframe)", 1); 
  glutAddMenuEntry("Cardinal mesh (solid)", 2); 
  glutAddMenuEntry("Beta-spline mesh (wireframe)", 3); 
  glutAddMenuEntry("Beta-spline mesh (solid)", 4); 
  glutAttachMenu(GLUT_RIGHT_BUTTON);

  M=N=5; D1=D2=3;
  for (i=0;i<MAXPTS;i++) for (j=0;j<MAXPTS;j++) 
    { V[i][j][0]=i/5.0; V[i][j][2]=j/5.0; V[i][j][1]=0.0; w[i][j]=1.0; }
  wlist1=glGenLists(1); 
  wlist2=glGenLists(1); 
  slist1=glGenLists(1); 
  slist2=glGenLists(1); 

  glutMainLoop();
  return 0;
}