#include <math.h>
#include <stdio.h>
#include <gl/gl.h>
#include <gl.h>
#include <device.h>
#include <stdlib.h>

#define RDIM 61
#define TDIM 61
#define DR 1.0 / (float) (RDIM - 1)
#define DT (2*3.1415927) / (float) (TDIM - 1)

/* functions are defined: */
float f1 (float x, float y);
float f2 (float x, float y);
float f3 (float x, float y);

float g1 (float x, float y);
float g2 (float x, float y);
float g3 (float x, float y);

float power(float base, float p);

/* global variables: */
typedef long Window;
long xorigin, yorigin, xwidth, ywidth;
float minx, miny, maxx, maxy;
Window ParamWindow;
float c1, c2, Pi;
float point[3000][3];
int k = 1;
int q = 1;
int check, okay;
float red, blue, green;
float coord[2];

/* will be used to call TCL windows */
int system (const char *string); 

/* this gets points from the mouse in the C window and translates them into x and y coordinates */
void getmouseXYcoord()
{
  long xx, yy;
  getorigin(&xorigin,&yorigin);
  getsize(&xwidth,&ywidth);
  xx = getvaluator(MOUSEX);
  yy = getvaluator(MOUSEY);
  c1 = (maxx-minx)*(xx-xorigin-(minx/(minx-maxx))*xwidth)/xwidth;
  c2 = (maxy-miny)*(yy-yorigin-(miny/(miny-maxy))*ywidth)/ywidth;
}

/* creates Bryant's Boy's Surface: */
void Boys()
{ 
  float r, t, x, y;
  int i, j;
   

  printf("(bbox-draw \"world\" off) \n");
  printf("(normalization \"world\" none) \n");
  fflush(stdout);

  printf("(geometry boys { :pretzel }) \n");
  printf("(read geometry { define pretzel \n");
  fflush(stdout);

  printf("{ CMESH \n");
  printf("%d %d \n", RDIM, TDIM);
  for(i = 0, r = 0; i < RDIM; ++i, r += DR) {
    for(j = 0, t = 0; j< TDIM; ++j, t += DT) {
      x = r*fcos(t);
      y = r*fsin(t);
      printf("%f %f %f 1 1 0 1 \n", f1(x, y), f2(x, y), f3(x, y));
    }
  }
 printf("} } ) \n");
 fflush(stdout);
}

/* makes the small spheres which serve as "landmarks" on Boy's
Surface: */
void Landmarks()
{
  printf("(geometry sphere1 { : peach } ) \n");
  printf("(read geometry { define peach \n");
  printf("{ SPHERE \n");
  printf("0.03 \n");
  printf("0 0 -2 \n");
  printf("} } ) \n");
  printf("(merge-ap \"sphere1\" appearance { \n");
  printf("material { \n");
  printf("\t ambient 1 0 0 \n");
  printf("\t *diffuse 1 0 0 \n");
  printf("} } ) \n");
  fflush(stdout);

  printf("(geometry sphere2 { : plum } ) \n");
  printf("(read geometry { define plum \n");
  printf("{ SPHERE \n");
  printf("0.04 \n");
  printf("0 0 0 \n");
  printf("} } ) \n");
  printf("(merge-ap \"sphere2\" appearance { \n");
  printf("material { \n");
  printf("\t ambient 0 0 1 \n");
  printf("\t *diffuse 0 0 1 \n");
  printf("} } ) \n");
  fflush(stdout); 
}

/* makes the line on Boy's Surface: */
void Line()
{
  float x, y;
  int j;
  
  if(check == 0) {
    printf("(geometry line%d { : pastry%d }) \n", q, q);
    printf("(read geometry { define pastry%d \n", q);
    printf("{ VECT \n");
    printf("1 %d 1 \n", k - 1);
    printf("%d \n", k - 1);
    printf("1 \n");
    for(j = 1; j < k; j++) {
      printf("%f %f %f \n", point[j][0], point[j][1], point[j][2]);
    }
    printf("%.3f %.3f %.3f 1 \n", red, green, blue);
    printf("} } ) \n");
    printf("(merge-ap \"line%d\" appearance { \n", q);
    printf("\t *linewidth 4 \n");
    printf("} ) \n");
    fflush(stdout);
  
    q++;
  }
}

/* this is used to prevent the line from jumping: */
float Length () {
  return sqrt((c1 - coord[0])*(c1 - coord[0]) + (c2 - coord[1])*(c2 - coord[1]));
}

void main()
{
  /* variables are defined and initialized */
  float linecoord[2];
  char title[38];
  float blah;
  float x, y;
  int u;
  Device dev;      
  short val, draw, clean, finish, blueline, redline, greenline;
  draw = 0; clean = 0; finish = 0, blueline = 0, redline = 0; greenline = 0;  
  Pi = 4.0*atan(1.0);
  minx = -1.0;
  maxx = 1.0;
  miny = -1.0;
  maxy = 1.0;	

  /* brings up TCL window with instructions */
  system ("instructions2.tcl &");

   /* sets up the C window */
  foreground();
  keepaspect(1, 1);
  sprintf(title,"(%f,%f)  (%f,%f)", minx, miny, maxx, maxy);
  ParamWindow = winopen(title);
  winset(ParamWindow);
  pntsize(4.0);
  gconfig();
  qdevice(LEFTMOUSE);
  qdevice(MIDDLEMOUSE);
  qdevice(RIGHTMOUSE);     
  qdevice(BKEY);
  qdevice(RKEY);
  qdevice(GKEY);
  ortho2(minx, maxx, miny, maxy);
  color(BLACK);
  clear();
  color(WHITE);
  circf(0, 0, 1);
  color(RED);
  circf(0, 0, 0.02);
  color(BLUE);
  circf(0.72555, 0, 0.02);
  circf(-.362775,0.628345, 0.02);
  circf(-.362775,-0.628345, 0.02);
  doublebuffer(); 

  Boys();
  Landmarks();

  /* this loop checks for input from the the mouse, and responds to the input */ 
  while(1) {
    while (qtest()) {
      dev = qread(&val); 
      switch(dev) 
	{
	case LEFTMOUSE:
	  draw = val;
	  break;
	case RIGHTMOUSE:
	  finish = val; 
	  break;
	case MIDDLEMOUSE:
	  clean = val;
	  break;
	case BKEY:
	  blueline = val;
	  break;
	case RKEY:
	  redline = val;
	  break;
	case GKEY:
	  greenline = val;
	  break;
	}
    }
    
    /* exits C program and Geomview if right mouse button is hit */
    if(finish) {
      printf("(exit)");
      fflush(stdout);
    }

    /* deletes lines in C window and on Boy's Surface if middle mouse button is hit */
    if(clean) {
      if(okay == 0) {
	for(u = 1; u < q; u++) {
	  printf("(delete \"line%d\") \n", u);
	  fflush(stdout);
	}
      }
      okay = 1;
      k = 1;
      q = 1;
      color(BLACK);
      clear();
      color(WHITE);
      circf(0, 0, 1);
      color(RED);
      circf(0, 0, 0.02);
      color(BLUE);
      circf(0.72555, 0, 0.02);
      circf(-.362775,0.628345, 0.02);
      circf(-.362775,-0.628345, 0.02);
      swapbuffers();
    }

    /* draws lines in C window*/
    if(draw){
      blah = Length();
      getmouseXYcoord();
      if ((blah < .1 && blah > .001) || k < 3) {
       coord[0] = c1; coord[1] = c2;
       color(BLACK);
       bgnpoint(); 
       v2f(coord); 
       endpoint();
       color(WHITE);     
       swapbuffers();
     
       x = c1; y = c2;
       point[k][0] = f1(x,y);
       point[k][1] = f2(x,y);
       point[k][2] = f3(x,y);
       k++;
       check = 0;
       okay = 0;
      }
    }
/* maps lines onto Boy's Surface in blue, red, or green: */
    if(blueline) {
      red = 0; green = 0; blue = 1;
      Line();
      check = 1;
      k = 1;
    }

    if(redline) {
      red = 1; green = 0; blue = 0;
      Line();
      check = 1;
      k = 1;
    }

    if(greenline) {
      red = 0; green = .545; blue = .1;
      Line();
      check = 1;
      k = 1;
    }

  }
}

/* functions that parametrize the surface: */
float g1(float x, float y)
{
  float g;
  g = (-3*(-(y*(3*power(x,2) - power(y,2))*
          (sqrt(5) + 2*power(x,3) - 6*x*power(y,2))*
          (x - power(x,5) + 10*power(x,3)*power(y,2) - 
            5*x*power(y,4))) - 
       y*(-1 + 5*power(x,4) - 10*power(x,2)*power(y,2) + 
          power(y,4))*(-1 + sqrt(5)*power(x,3) + power(x,6) - 
          3*sqrt(5)*x*power(y,2) - 15*power(x,4)*power(y,2) + 
          15*power(x,2)*power(y,4) - power(y,6))))/
   (2.*(power(y,2)*power(3*power(x,2) - power(y,2),2)*
        power(sqrt(5) + 2*power(x,3) - 6*x*power(y,2),2) + 
       power(-1 + sqrt(5)*power(x,3) + power(x,6) - 
         3*sqrt(5)*x*power(y,2) - 15*power(x,4)*power(y,2) + 
         15*power(x,2)*power(y,4) - power(y,6),2)));
  return g;
}

float g2 (float x, float y)
{
  float g;
  g = (-3*(power(y,2)*(3*power(x,2) - power(y,2))*
        (sqrt(5) + 2*power(x,3) - 6*x*power(y,2))*
        (1 + 5*power(x,4) - 10*power(x,2)*power(y,2) + 
          power(y,4)) + (x + power(x,5) - 
          10*power(x,3)*power(y,2) + 5*x*power(y,4))*
        (-1 + sqrt(5)*power(x,3) + power(x,6) - 
          3*sqrt(5)*x*power(y,2) - 15*power(x,4)*power(y,2) + 
          15*power(x,2)*power(y,4) - power(y,6))))/
   (2.*(power(y,2)*power(3*power(x,2) - power(y,2),2)*
        power(sqrt(5) + 2*power(x,3) - 6*x*power(y,2),2) + 
       power(-1 + sqrt(5)*power(x,3) + power(x,6) - 
         3*sqrt(5)*x*power(y,2) - 15*power(x,4)*power(y,2) + 
         15*power(x,2)*power(y,4) - power(y,6),2)));
  return g;
}

float g3 (float x, float y)
{
  float g;
  g = -0.5 + (-(y*(3*power(x,2) - power(y,2))*
         (sqrt(5) + 2*power(x,3) - 6*x*power(y,2))*
         (1 + power(x,6) - 15*power(x,4)*power(y,2) + 
           15*power(x,2)*power(y,4) - power(y,6))) + 
      2*x*y*(power(x,2) - 3*power(y,2))*
       (3*power(x,2) - power(y,2))*
       (-1 + sqrt(5)*power(x,3) + power(x,6) - 
         3*sqrt(5)*x*power(y,2) - 15*power(x,4)*power(y,2) + 
         15*power(x,2)*power(y,4) - power(y,6)))/
    (power(y,2)*power(3*power(x,2) - power(y,2),2)*
       power(sqrt(5) + 2*power(x,3) - 6*x*power(y,2),2) + 
      power(-1 + sqrt(5)*power(x,3) + power(x,6) - 
        3*sqrt(5)*x*power(y,2) - 15*power(x,4)*power(y,2) + 
        15*power(x,2)*power(y,4) - power(y,6),2));
  return g;
}

float f1 (float x, float y)
{ 
  float j;
  float gg1, gg2, gg3;
  gg1 = g1(x, y);
  gg2 = g2(x, y);
  gg3 = g3(x, y);

  j = gg1/(power(gg1,2) + power(gg2,2) + power(gg3,2));

  return j;
}

float f2 (float x, float y)
{
  float k;
  float gg1, gg2, gg3;
  gg1 = g1(x, y);
  gg2 = g2(x, y);
  gg3 = g3(x, y);

  k = gg2/(power(gg1,2) + power(gg2,2) + power(gg3,2));

  return k;
}

float f3 (float x, float y)
{
  float l;
  float gg1, gg2, gg3;
  gg1 = g1(x,y);
  gg2 = g2(x,y);
  gg3 = g3(x,y);

  l = gg3/(power(gg1,2) + power(gg2,2) + power(gg3,2));

  return l;
}

float power (float base, float p)
{
  int i;
  float ans;
  ans = base;
  for (i = 2; i <= p; i++) {
    ans = ans*base;
  }
  return ans;
}