/******************************************************************************
* CagdSRev.c - Surface of revolution out of a given profile.		      *
*******************************************************************************
* Written by Gershon Elber, Mar. 91.					      *
******************************************************************************/

#include "cagd_loc.h"

static int 
    CircRationalKnotVector[12] = { 0, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 4 },
    CircPolynomialKnotVector[17] = { 0, 0, 0, 0, 1, 1, 1, 2, 2, 2,
							3, 3, 3, 4, 4, 4, 4 };

static CagdRType
    PolyApproxRotAngles[] = {
	0,
	33.523898, /* arcsin(4 (sqrt(2) - 1) / 3) */
	90 - 33.523898,
	90
    };

/*****************************************************************************
* DESCRIPTION:                                                               M
* Constructs a surface of revolution around the Z axis of the given profile  M
* curve. Resulting surface will be a Bspline surface, while input may be     M
* either a Bspline or a Bezier curve.                                        M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To create surface of revolution around Z with.               M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  Surface of revolution.                                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   CagdSurfaceRev, surface of revolution, surface constructors              M
*****************************************************************************/
CagdSrfStruct *CagdSurfaceRev(CagdCrvStruct *Crv)
{
    int i, j, i9,
	Len = Crv -> Length;
    CagdPointType
	PType = Crv -> PType;
    CagdRType **SrfPoints,
	sin45 = sin(M_PI / 4.0),
	**CrvPoints = Crv -> Points;
    CagdSrfStruct
	*Srf = BspPeriodicSrfNew(9, Len, 3, Crv -> Order,
				 FALSE, Crv -> Periodic, CAGD_PT_P3_TYPE);

    for (i = 0; i < 12; i++)
	Srf -> UKnotVector[i] = (CagdRType) CircRationalKnotVector[i];

    switch (Crv -> GType) {
	case CAGD_CBEZIER_TYPE:
	    BspKnotUniformOpen(Len, Crv -> Order, Srf -> VKnotVector);
	    break;
	case CAGD_CBSPLINE_TYPE:
	    CAGD_GEN_COPY(Srf -> VKnotVector, Crv -> KnotVector,
			  sizeof(CagdRType) * (CAGD_CRV_PT_LST_LEN(Crv) +
					                       Crv -> Order));
	    break;
	case CAGD_CPOWER_TYPE:
	    CAGD_FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    return NULL;
	default:
	    CAGD_FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
	    return NULL;
    }

    SrfPoints = Srf -> Points;

    /* For each control point in the original curve - generate 9 points that */
    /* Form a circle perpendicular to the Z axis.			     */
    for (i = i9 = 0; i < Len; i++, i9 += 9) {
	SrfPoints[W][i9] = 1.0;
	switch (PType) {
	    case CAGD_PT_P3_TYPE:
		SrfPoints[W][i9] = CrvPoints[W][i];
	    case CAGD_PT_E3_TYPE:
		SrfPoints[X][i9] = CrvPoints[X][i];
		SrfPoints[Y][i9] = CrvPoints[Y][i];
		SrfPoints[Z][i9] = CrvPoints[Z][i];
		break;
	    default:
		CAGD_FATAL_ERROR(CAGD_ERR_UNSUPPORT_PT);
		break;
	}

	/* Last point is exactly the same as first one in circle - copy it.  */
	for (j = W; j <= Z; j++)
	    SrfPoints[j][i9 + 8] = SrfPoints[j][i9];

	/* The Z components are identical in all circle, while the XY        */
	/* components are rotated 45 degrees at a time:			     */
	for (j = 1; j < 8; j++) {
	    SrfPoints[W][i9 + j] = SrfPoints[W][i9];
	    SrfPoints[X][i9 + j] = SrfPoints[X][i9 + j - 1] * sin45 -
				   SrfPoints[Y][i9 + j - 1] * sin45;
	    SrfPoints[Y][i9 + j] = SrfPoints[X][i9 + j - 1] * sin45 +
				   SrfPoints[Y][i9 + j - 1] * sin45;
	    SrfPoints[Z][i9 + j] = SrfPoints[Z][i9];
	}

	/* And finally we need to compensate for the W's on every second pt. */
	for (j = 1; j < 8; j += 2) {
	    SrfPoints[W][i9 + j] *= sin45;
	    SrfPoints[Z][i9 + j] *= sin45;
	}
    }

    return Srf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Constructs a surface of revolution around the Z axis of the given profile  M
* curve. Resulting surface will be a Bspline surface, while input may be     M
* either a Bspline or a Bezier curve.                                        M
*   Resulting surface will be a polynomial Bspline surface, approximating a  M
* surface of revolution using a polynomial circle approx.                    M
* (See Faux & Pratt "Computational Geometry for Design and Manufacturing").  M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To approximate a surface of revolution around Z with.        M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  Surface of revolution approximation.                   M
*                                                                            *
* KEYWORDS:                                                                  M
*   CagdSurfaceRevPolynomialApprox, surface of revolution, surface           M
*   constructors						             M
*****************************************************************************/
CagdSrfStruct *CagdSurfaceRevPolynomialApprox(CagdCrvStruct *Crv)
{
    int i, j, i13,
	Len = Crv -> Length;
    CagdPointType
	PType = Crv -> PType;
    CagdRType **SrfPoints,
	**CrvPoints = Crv -> Points;
    CagdSrfStruct
	*Srf = BspPeriodicSrfNew(13, Len, 4, Crv -> Order,
				 FALSE, Crv -> Periodic, CAGD_PT_E3_TYPE);

    if (CAGD_IS_RATIONAL_CRV(Crv)) {
	CAGD_FATAL_ERROR(CAGD_ERR_POLYNOMIAL_EXPECTED);
	return NULL;
    }

    for (i = 0; i < 17; i++)
	Srf -> UKnotVector[i] = (CagdRType) CircPolynomialKnotVector[i];

    switch (Crv -> GType) {
	case CAGD_CBEZIER_TYPE:
	    BspKnotUniformOpen(Len, Crv -> Order, Srf -> VKnotVector);
	    break;
	case CAGD_CBSPLINE_TYPE:
	    CAGD_GEN_COPY(Srf -> VKnotVector, Crv -> KnotVector,
			  sizeof(CagdRType) * (CAGD_CRV_PT_LST_LEN(Crv) +
								Crv -> Order));
	    break;
	case CAGD_CPOWER_TYPE:
	    CAGD_FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    return NULL;
	default:
	    CAGD_FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
	    return NULL;
    }

    SrfPoints = Srf -> Points;

    /* For each control point in original curve - generate 13 points that    */
    /* Form a circle 9approximation) perpendicular to the Z axis.	     */
    for (i = i13 = 0; i < Len; i++, i13 += 13) {
	int Quad;

	switch (PType) {
	    case CAGD_PT_E3_TYPE:
		SrfPoints[X][i13] = CrvPoints[X][i];
		SrfPoints[Y][i13] = CrvPoints[Y][i];
		SrfPoints[Z][i13] = CrvPoints[Z][i];
		break;
	    default:
		CAGD_FATAL_ERROR(CAGD_ERR_UNSUPPORT_PT);
		break;
	}

	/* Last point is exactly the same as first one in circle - copy it.  */
	for (j = X; j <= Z; j++)
	    SrfPoints[j][i13 + 12] = SrfPoints[j][i13];

	/* The Z components are identical in all circle, while the XY        */
	/* components are functions of PolyApproxRotAngles:		     */
	for (Quad = 0, j = 1; j < 12; j++) {
	    CagdRType Angle, CosAngle, SinAngle;

	    if (j % 3 == 0)
		Quad++;
	    Angle = Quad * 90 + PolyApproxRotAngles[j % 3];
	    CosAngle = cos(DEG2RAD(Angle));
	    SinAngle = sin(DEG2RAD(Angle));

	    if (ABS(CosAngle) > ABS(SinAngle))
		CosAngle /= ABS(CosAngle);
	    else
		SinAngle /= ABS(SinAngle);
	    SrfPoints[X][i13 + j] = SrfPoints[X][i13] * CosAngle;
	    SrfPoints[Y][i13 + j] = SrfPoints[X][i13] * SinAngle;
	    SrfPoints[Z][i13 + j] = SrfPoints[Z][i13];
	}
    }

    return Srf;
}
