#include  <stdio.h>
#include  <stdlib.h>
#include  <math.h>

#define   MAX_SIZE  100

struct P3D {
     double  x, y, z;
};

typedef  struct P3D  POINT3D;
typedef  struct P3D  NORM3D;

POINT3D  *r1, *r2;
NORM3D   *n1, *n2;

POINT3D  Row1[MAX_SIZE], Row2[MAX_SIZE];
NORM3D   Norm1[MAX_SIZE], Norm2[MAX_SIZE];

double   u1, u2, v1, v2, u, v;
double   Xmin, Xmax;
double   Ymin, Ymax;
double   Zmin, Zmax;
int      StepU, StepV;
double   DeltaU, DeltaV;

void     Funct(double u, double v, POINT3D *p, NORM3D *Pu, NORM3D *Pv);
void     CrossProduct(NORM3D A, NORM3D B, NORM3D *C);
void     CrossProduct(NORM3D A, NORM3D B, NORM3D *C);
void     OneRow(POINT3D *r, NORM3D *n);
void     SmoothTriangle(POINT3D A, POINT3D B, POINT3D C,
                        NORM3D a, NORM3D b, NORM3D c);
void     GridGen(void);

void  CrossProduct(NORM3D A, NORM3D B, NORM3D *C)
{
     NORM3D  h;
     double   d;

     h.x =   A.y*B.z - B.y*A.z;
     h.y = -(A.x*B.z - B.x*A.z);
     h.z =   A.x*B.y - B.x*A.y;

     d   = sqrt(h.x*h.x + h.y*h.y + h.z*h.z);

     C->x = h.x/d;
     C->y = h.y/d;
     C->z = h.z/d;
}

void  OneRow(POINT3D *r, NORM3D *n)
{
     POINT3D  P;
     NORM3D   Pu, Pv, N;
     double   nn;
     int      j;

     for (j = 0, v = v1; v <= v2; j++, v += DeltaV) {
          Funct(u, v, &P, &Pu, &Pv);
          if (P.x < Xmin) Xmin = P.x;
          if (P.x > Xmax) Xmax = P.x;
          if (P.y < Ymin) Ymin = P.y;
          if (P.y > Ymax) Ymax = P.y;
          if (P.z < Zmin) Zmin = P.z;
          if (P.z > Zmax) Zmax = P.z;
          CrossProduct(Pu, Pv, &N);
          nn = sqrt(N.x*N.x + N.y*N.y + N.z*N.z);
          r[j].x = P.x;
          r[j].y = P.y;
          r[j].z = P.z;
          n[j].x = N.x/nn;
          n[j].y = N.y/nn;
          n[j].z = N.z/nn;
     }
}

void  SmoothTriangle(POINT3D A, POINT3D B, POINT3D C, NORM3D a, NORM3D b, NORM3D c)
{
     printf("          smooth_triangle {\n");
     printf("               < %f, %f, %f >, ",  A.x, A.y, A.z);
     printf(               "< %f, %f, %f >,\n", a.x, a.y, a.z);
     printf("               < %f, %f, %f >, ",  B.x, B.y, B.z);
     printf(               "< %f, %f, %f >,\n", b.x, b.y, b.z);
     printf("               < %f, %f, %f >, " , C.x, C.y, C.z);
     printf(               "< %f, %f, %f >\n" , c.x, c.y, c.z);
     printf("          }\n");
}

#define  SWAP(a, b, c)    { c = a; a = b; b = c; }

void  GridGen(void)
{
     POINT3D *r;
     NORM3D  *n;
     int     j;

     r1 = Row1;
     n1 = Norm1;
     r2 = Row2;
     n2 = Norm2;
     u  = u1;
     OneRow(r1, n1);
     for (u = u1+DeltaU; u <= u2; u += DeltaU) {
          OneRow(r2, n2);
          for (v = v1, j = 0; v < v2; j++, v += DeltaV) {
               SmoothTriangle(r1[j], r2[j], r1[j+1], n1[j], n2[j], n1[j+1]);
               SmoothTriangle(r2[j], r1[j+1], r2[j+1], n2[j], n1[j+1], n2[j+1]);
          }
          SWAP(r1, r2, r);
          SWAP(n1, n2, n);
     }
}

void  main(int argc, char *argv[])
{
     if (argc != 3) {
          printf("Use uv StepU StepV\n");
          exit(1);
     }
     StepU = atoi(argv[1]);
     StepV = atoi(argv[2]);
     if (StepU < 1) {
          printf("*** ERROR: StepU (%d) must be greater than or equal to 2\n",
                 StepU);
          exit(1);
     }
     if (StepU >= MAX_SIZE-1) {
          printf("*** ERROR: The max. of StepU (%d) is %d\n",
                 StepU, MAX_SIZE-1);
          exit(1);
     }
     if (StepV < 1) {
          printf("*** ERROR: StepV (%d) must be greater than or equal to 2\n",
                 StepV);
          exit(1);
     }
     if (StepV >= MAX_SIZE-1) {
          printf("*** ERROR: The max. of StepV (%d) is %d\n",
                 StepV, MAX_SIZE-1);
          exit(1);
     }
     u1 = 0.0;
     u2 = 2*acos(-1.0);
     v1 = 0.0;
     v2 = 2*acos(-1.0);
     DeltaU = (u2 - u1)/StepU;
     DeltaV = (v2 - v1)/StepV;
     printf("#declare  TRIANGULATED_TORUS = \n");
     printf("     union {\n");
     GridGen();
     printf("          clipped_by {\n");
     printf("               box { < %f, %f, %f >,\n",
                                  Xmin, Ymin, Zmin);
     printf("                     < %f, %f, %f >\n",
                                  Xmax, Ymax, Zmax);
     printf("               }\n");
     printf("          }\n");
     printf("          bounded_by {\n");
     printf("               box { < %f, %f, %f >,\n",
                                  Xmin, Ymin, Zmin);
     printf("                     < %f, %f, %f >\n",
                                  Xmax, Ymax, Zmax);
     printf("               }\n");
     printf("          }\n");
     printf("     }\n");
}

void  Funct(double u, double v, POINT3D *p, NORM3D *Pu, NORM3D *Pv)
{
     const double  a = 8.0;
     const double  b = 3.0;
     const double  c = 7.0;

     p->x = (a + b*cos(v))*cos(u); 	/* point on torus           */
     p->y = c*sin(v);
     p->z = (a + b*cos(v))*sin(u);

     Pu->x = (a + b*cos(v))*(-sin(u));  /* partial w.r.t. u         */
     Pu->y = 0.0;
     Pu->z = (a + b*cos(v))*cos(u);

     Pv->x = cos(u)*b*(-sin(v));        /* partial w.r.t. v         */
     Pv->y = c*cos(v);
     Pv->z = sin(u)*b*(-sin(v));
}