/* $Id: //depot/main/rmanprod/docs/projects/rps/reference/examples/runuquads.c#1 $  (Pixar - RenderMan Division)  $Date: 2009/05/07 $ */
/*
** Copyright (c) 1998 PIXAR.  All rights reserved.  This program or
** documentation contains proprietary confidential information and trade
** secrets of PIXAR.  Reverse engineering of object code is prohibited.
** Use of copyright notice is precautionary and does not imply
** publication.
**
**                      RESTRICTED RIGHTS NOTICE
**
** Use, duplication, or disclosure by the Government is subject to the
** following restrictions:  For civilian agencies, subparagraphs (a) through
** (d) of the Commercial Computer Software--Restricted Rights clause at
** 52.227-19 of the FAR; and, for units of the Department of Defense, DoD
** Supplement to the FAR, clause 52.227-7013 (c)(1)(ii), Rights in
** Technical Data and Computer Software.
**
** Pixar
** 1001 West Cutting Blvd.
** Richmond, CA  94804
*/
#include <stdio.h>
#include <math.h>
#include "runuquads.h"
#ifdef _WINDOWS
#define M_SQRT1_2	0.70710678118654752440
#define M_PI_2		1.57079632679489661923
#define M_PI		3.14159265358979323846
#endif

void
RuNuSphere(RtFloat radius, RtFloat zmin, RtFloat zmax, RtFloat thetamax)
{
   RtFloat umax;
   RtFloat uknots[12];
   static RtFloat uk[12] =
	{ 0.0, 0.0, 0.0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1.0, 1.0, 1.0};
   RtFloat vmin, vmax;
   RtFloat vknots[8];
   static RtFloat vk[8] =
	{ 0.0, 0.0, 0.0, 0.5, 0.5, 1.0, 1.0, 1.0};
   static RtFloat uspline[9][4] =
	{ {1.0, 0.0, 1.0, 1.0},
	  {M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {0.0, 1.0, 1.0, 1.0},
	  {-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {-1.0, 0.0, 1.0, 1.0},
	  {-M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {0.0, -1.0, 1.0, 1.0},
	  {M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {1.0, 0.0, 1.0, 1.0}};
   static RtFloat vspline[5][4] =
	{ {0.0, 0.0, -1.0, 1.0},
	  {M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2},
	  {1.0, 1.0, 0.0, 1.0},
	  {M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {0.0, 0.0, 1.0, 1.0}};
   RtFloat cp[5][9][4];
   RtFloat r[4];
   int i, j, k;

   umax = thetamax / 360.0;
   for (i=0; i<12; i++)
	uknots[i] = uk[i] / umax;

   vmin = zmin / radius;
   vmax = zmax / radius;
   if (vmin < -1.0) vmin = -1.0;
   if (vmax > 1.0) vmax = 1.0;
   vmin = (asin(vmin) + M_PI_2) / M_PI;
   vmax = (asin(vmax) + M_PI_2) / M_PI;
   for (i=0; i<8; i++)
	vknots[i] = (vk[i] - vmin) / (vmax - vmin);

   r[0] = radius;
   r[1] = radius;
   r[2] = radius;
   r[3] = 1.0;
   for (i=0; i<5; i++) {
	for (j=0; j<9; j++) {
	    for (k=0; k<4; k++) {
		cp[i][j][k] = r[k] * uspline[j][k] * vspline[i][k];
	    }
	}
   }

   RiNuPatch((RtInt) 9, (RtInt) 3, uknots, (RtFloat) 0.0, (RtFloat) 1.0,
	     (RtInt) 5, (RtInt) 3, vknots, (RtFloat) 0.0, (RtFloat) 1.0,
	     "Pw", cp, RI_NULL);
}

void
RuNuHyperboloid(RtPoint point1, RtPoint point2, RtFloat thetamax)
{
   RtFloat umax;
   RtFloat uknots[12];
   static RtFloat uk[12] =
	{0.0, 0.0, 0.0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1.0, 1.0, 1.0};
   static RtFloat vknots[4] =
	{0.0, 0.0, 1.0, 1.0};
   static RtFloat xuspline[9][2] =
	{{1.0, 0.0},
	  {M_SQRT1_2, -M_SQRT1_2},
	  {0.0, -1.0},
	  {-M_SQRT1_2, -M_SQRT1_2},
	  {-1.0, 0.0},
	  {-M_SQRT1_2, M_SQRT1_2},
	  {0.0, 1.0},
	  {M_SQRT1_2, M_SQRT1_2},
	  {1.0, 0.0}};
   static RtFloat yuspline[9][2] =
	{{0.0, 1.0},
	  {M_SQRT1_2, M_SQRT1_2},
	  {1.0, 0.0},
	  {M_SQRT1_2, -M_SQRT1_2},
	  {0.0, -1.0},
	  {-M_SQRT1_2, -M_SQRT1_2},
	  {-1.0, 0.0},
	  {-M_SQRT1_2, M_SQRT1_2},
	  {0.0, 1.0}};
   static RtFloat zwuspline[9] =
	 {1.0,
	  M_SQRT1_2,
	  1.0,
	  M_SQRT1_2,
	  1.0,
	  M_SQRT1_2,
	  1.0,
	  M_SQRT1_2,
	  1.0};
   RtFloat cp[2][9][4];
   int i, j;

   umax = thetamax / 360.0;
   for (i=0; i<12; i++)
	uknots[i] = uk[i] / umax;

   for (j=0; j<9; j++) {
	cp[0][j][0] = point1[0] * xuspline[j][0] +
		      point1[1] * xuspline[j][1];
	cp[0][j][1] = point1[0] * yuspline[j][0] +
		      point1[1] * yuspline[j][1];
	cp[0][j][2] = point1[2] * zwuspline[j];
	cp[0][j][3] = zwuspline[j];
	cp[1][j][0] = point2[0] * xuspline[j][0] +
		      point2[1] * xuspline[j][1];
	cp[1][j][1] = point2[0] * yuspline[j][0] +
		      point2[1] * yuspline[j][1];
	cp[1][j][2] = point2[2] * zwuspline[j];
	cp[1][j][3] = zwuspline[j];
    }

   RiNuPatch(9, 3, uknots, 0.0, 1.0,
	     2, 2, vknots, 0.0, 1.0,
	     "Pw", cp, RI_NULL);
}

void
RuNuCone(RtFloat height, RtFloat radius, RtFloat thetamax)
{
    RtPoint point1, point2;

    point1[0] = radius;
    point1[1] = 0.0;
    point1[2] = 0.0;
    point2[0] = 0.0;
    point2[1] = 0.0;
    point2[2] = height;
    RuNuHyperboloid(point1, point2, thetamax);
}

void
RuNuCylinder(RtFloat radius, RtFloat zmin, RtFloat zmax, RtFloat thetamax)
{
    RtPoint point1, point2;

    point1[0] = radius;
    point1[1] = 0.0;
    point1[2] = zmin;
    point2[0] = radius;
    point2[1] = 0.0;
    point2[2] = zmax;
    RuNuHyperboloid(point1, point2, thetamax);
}

void
RuNuDisk(RtFloat height, RtFloat radius, RtFloat thetamax)
{
    RtPoint point1, point2;

    point1[0] = radius;
    point1[1] = 0.0;
    point1[2] = height;
    point2[0] = 0.0;
    point2[1] = 0.0;
    point2[2] = height;
    RuNuHyperboloid(point1, point2, thetamax);
}

void
RuNuTorus(RtFloat majorradius, RtFloat minorradius,
	RtFloat phimin, RtFloat phimax, RtFloat thetamax)
{
   RtFloat umax;
   RtFloat uknots[12];
   RtFloat vmin, vmax;
   RtFloat vknots[12];
   static RtFloat knots[12] =
	{ 0.0, 0.0, 0.0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1.0, 1.0, 1.0};
   static RtFloat uspline[9][4] =
	{ {1.0, 0.0, 1.0, 1.0},
	  {M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {0.0, 1.0, 1.0, 1.0},
	  {-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {-1.0, 0.0, 1.0, 1.0},
	  {-M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {0.0, -1.0, 1.0, 1.0},
	  {M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {1.0, 0.0, 1.0, 1.0}};
   RtFloat vspline[9][4];
   RtFloat cp[9][9][4];
   int i, j, k;

   umax = thetamax / 360.0;
   vmin = phimin / 360.0;
   vmax = phimax / 360.0;

   for (i=0; i<12; i++) {
	uknots[i] = knots[i] / umax;
	vknots[i] = (knots[i] - vmin) / (vmax - vmin);
   }


   vspline[0][0] = majorradius - minorradius;
   vspline[0][1] = majorradius - minorradius;
   vspline[0][2] = 0.0;
   vspline[0][3] = 1.0;
   vspline[1][0] = (majorradius - minorradius) * M_SQRT1_2;
   vspline[1][1] = (majorradius - minorradius) * M_SQRT1_2;
   vspline[1][2] = -minorradius * M_SQRT1_2;
   vspline[1][3] = M_SQRT1_2;
   vspline[2][0] = majorradius;
   vspline[2][1] = majorradius;
   vspline[2][2] = -minorradius;
   vspline[2][3] = 1.0;
   vspline[3][0] = (majorradius + minorradius) * M_SQRT1_2;
   vspline[3][1] = (majorradius + minorradius) * M_SQRT1_2;
   vspline[3][2] = -minorradius * M_SQRT1_2;
   vspline[3][3] = M_SQRT1_2;
   vspline[4][0] = majorradius + minorradius;
   vspline[4][1] = majorradius + minorradius;
   vspline[4][2] = 0.0;
   vspline[4][3] = 1.0;
   vspline[5][0] = (majorradius + minorradius) * M_SQRT1_2;
   vspline[5][1] = (majorradius + minorradius) * M_SQRT1_2;
   vspline[5][2] = minorradius * M_SQRT1_2;
   vspline[5][3] = M_SQRT1_2;
   vspline[6][0] = majorradius;
   vspline[6][1] = majorradius;
   vspline[6][2] = minorradius;
   vspline[6][3] = 1.0;
   vspline[7][0] = (majorradius - minorradius) * M_SQRT1_2;
   vspline[7][1] = (majorradius - minorradius) * M_SQRT1_2;
   vspline[7][2] = minorradius * M_SQRT1_2;
   vspline[7][3] = M_SQRT1_2;
   vspline[8][0] = majorradius - minorradius;
   vspline[8][1] = majorradius - minorradius;
   vspline[8][2] = 0.0;
   vspline[8][3] = 1.0;

   for (i=0; i<9; i++) {
	for (j=0; j<9; j++) {
	    for (k=0; k<4; k++) {
		cp[i][j][k] = uspline[j][k] * vspline[i][k];
	    }
	}
   }

   RiNuPatch(9, 3, uknots, 0.0, 1.0,
	     9, 3, vknots, 0.0, 1.0,
	     "Pw", cp, RI_NULL);
}

void
RuNuParaboloid(RtFloat rmax, RtFloat zmin, RtFloat zmax, RtFloat thetamax)
{
   RtFloat umax;
   RtFloat uknots[12];
   static RtFloat uk[12] =
	{ 0.0, 0.0, 0.0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1.0, 1.0, 1.0};
   RtFloat vmin, vmax;
   RtFloat vknots[6];
   static RtFloat vk[6] =
	{ 0.0, 0.0, 0.0, 1.0, 1.0, 1.0};
   static RtFloat uspline[9][4] =
	{ {1.0, 0.0, 1.0, 1.0},
	  {M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {0.0, 1.0, 1.0, 1.0},
	  {-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {-1.0, 0.0, 1.0, 1.0},
	  {-M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {0.0, -1.0, 1.0, 1.0},
	  {M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2},
	  {1.0, 0.0, 1.0, 1.0}};
   static RtFloat vspline[3][4] =
	{ {0.0, 0.0, 0.0, 1.0},
	  {0.5, 0.5, 0.0, 1.0},
	  {1.0, 1.0, 1.0, 1.0}};
   RtFloat cp[3][9][4];
   RtFloat r[4];
   int i, j, k;

   umax = thetamax / 360.0;
   for (i=0; i<12; i++)
	uknots[i] = uk[i] / umax;

   vmin = zmin / zmax;
   if (vmin < 0.0) vmin = 0.0;
   vmax = 1.0;
   for (i=0; i<6; i++)
	vknots[i] = (vk[i] - vmin) / (vmax - vmin);

   r[0] = rmax;
   r[1] = rmax;
   r[2] = zmax;
   r[3] = 1.0;
   for (i=0; i<3; i++) {
	for (j=0; j<9; j++) {
	    for (k=0; k<4; k++) {
		cp[i][j][k] = r[k] * uspline[j][k] * vspline[i][k];
	    }
	}
   }

   RiNuPatch(9, 3, uknots, 0.0, 1.0,
	     3, 3, vknots, 0.0, 1.0,
	     "Pw", cp, RI_NULL);
}