/******************************************************************************
* Cagd_Arc.c - Curve representation of arcs and circles			      *
*******************************************************************************
* Written by Gershon Elber, Jun. 90.					      *
******************************************************************************/

#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include "cagd_loc.h"

#define UNIT_CIRCLE_ORDER 3		        /* Quadratic rational curve. */
#define UNIT_CIRCLE_LENGTH 9	  /* Nine control points in the unit circle. */
#define UNIT_PCIRCLE_ORDER 4			  /* Cubic polynomial curve. */
#define UNIT_PCIRCLE_LENGTH 13	       /* control points in the unit circle. */

static int
    UnitCircleKnots[UNIT_CIRCLE_ORDER + UNIT_CIRCLE_LENGTH] =
					{ 0, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 4 };
static int
    UnitCircleX[UNIT_CIRCLE_LENGTH] = { 1, 1, 0, -1, -1, -1, 0, 1, 1 };
static int
    UnitCircleY[UNIT_CIRCLE_LENGTH] = { 0, 1, 1, 1, 0, -1, -1, -1, 0 };

static int
    UnitPCircleKnots[UNIT_PCIRCLE_ORDER + UNIT_PCIRCLE_LENGTH] =
			 { 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
* Creates an arc at the specified position as a rational quadratic Bezier    M
* curve.								     M
*   The arc is assumed to be less than 180 degrees from Start to End in the  M
* shorter path as arc where Center as arc center.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Start:       Point of beginning of arc.                                  M
*   Center:      Point of arc.                                               M
*   End:         Point of end of arc. 	                                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *: A rational quadratic Bezier curve representing the arc. M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvCreateArc, circle, arc                                             M
*****************************************************************************/
CagdCrvStruct *BzrCrvCreateArc(CagdPtStruct *Start,
			       CagdPtStruct *Center,
			       CagdPtStruct *End)
{
    int i;
    CagdCrvStruct
	*Arc = BzrCrvNew(3, CAGD_PT_P3_TYPE);
    CagdRType Len, CosAlpha, Radius,
	**Points = Arc -> Points;
    CagdVecStruct V1, V2, V;

    /* Copy first point. */
    Points[X][0] = Start -> Pt[0];
    Points[Y][0] = Start -> Pt[1];
    Points[Z][0] = Start -> Pt[2];
    Points[W][0] = 1.0;

    /* Copy last point. */
    Points[X][2] = End -> Pt[0];
    Points[Y][2] = End -> Pt[1];
    Points[Z][2] = End -> Pt[2];
    Points[W][2] = 1.0;

    /* Compute position of middle point. */
    Len = 0.0;
    for (i = 0; i < 3; i++) {
	V1.Vec[i] = Start -> Pt[i] - Center -> Pt[i];
	V2.Vec[i] = End -> Pt[i] - Center -> Pt[i];
	V.Vec[i] = V1.Vec[i] + V2.Vec[i];
	Len += SQR(V.Vec[i]);
    }

    if (APX_EQ(Len, 0.0)) {
	CagdCrvFree(Arc);
	CAGD_FATAL_ERROR(CAGD_ERR_180_ARC);
	return NULL;
    }
    else
	Len = sqrt(Len);

    for (i = 0; i < 3; i++)
	V.Vec[i] /= Len;

    /* Compute cosine alpha (where alpha is the angle between V and V1. */
    Radius = sqrt(DOT_PROD(V1.Vec, V1.Vec));
    CosAlpha = DOT_PROD(V1.Vec, V.Vec) / Radius;

    CAGD_DIV_VECTOR(V, CosAlpha);
    CAGD_MULT_VECTOR(V, Radius);

    /* And finally fill in the middle point with CosAlpha as the Weight. */
    Points[X][1] = (Center -> Pt[0] + V.Vec[0]) * CosAlpha;
    Points[Y][1] = (Center -> Pt[1] + V.Vec[1]) * CosAlpha;
    Points[Z][1] = (Center -> Pt[2] + V.Vec[2]) * CosAlpha;
    Points[W][1] = CosAlpha;

    return Arc;
}

/******************************************************************************
* DESCRIPTION:                                                               M
* Creates a circle at the specified position as a rational quadratic Bspline M
* curve. 								     M
*   Constructs a unit circle as 4 90 degrees arcs of rational quadratic      M
* Bezier segments using a predefined constants.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   None                                                                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *: A rational quadratic bsplinecurve representing a unit   M
*                    circle.                                                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvCreateUnitCircle, circle                                           M
*****************************************************************************/
CagdCrvStruct *BspCrvCreateUnitCircle(void)
{
    int i;
    CagdRType Weight,
        W45 = sin(M_PI / 4.0);
    CagdCrvStruct
	*Circle = BspCrvNew(UNIT_CIRCLE_LENGTH, UNIT_CIRCLE_ORDER,
							      CAGD_PT_P3_TYPE);
    CagdRType
	**Points = Circle -> Points;

    for (i = 0; i < UNIT_CIRCLE_LENGTH + UNIT_CIRCLE_ORDER; i++)
	Circle -> KnotVector[i] = UnitCircleKnots[i];

    for (i = 0; i < UNIT_CIRCLE_LENGTH; i++) {
	Weight = Points[W][i] = i % 2 ? W45: 1.0;
	Points[X][i] = UnitCircleX[i] * Weight;
	Points[Y][i] = UnitCircleY[i] * Weight;
	Points[Z][i] = 0.0;
    }

    return Circle;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Creates a circle at the specified position as a rational quadratic Bspline M
* curve. Circle is always paralell to the XY plane.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Center:   Of circle to be created.                                       M
*   Radius:   Of circle to be created.                                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A circle centered at Center and radius Radius that is M
*                      parallel to the XY plane represented as a rational    M
*                      quadratic Bspline curve.                              M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvCreateCircle, circle                                               M
*****************************************************************************/
CagdCrvStruct *BspCrvCreateCircle(CagdPtStruct *Center, CagdRType Radius)
{
    CagdPtStruct OriginPt;
    CagdCrvStruct
	*Circle = BspCrvCreateUnitCircle();

    /* Do it in two stages: 1. scale, 2. translate */
    OriginPt.Pt[0] = OriginPt.Pt[1] = OriginPt.Pt[2] = 0.0;
    CagdCrvTransform(Circle, OriginPt.Pt, Radius);
    CagdCrvTransform(Circle, Center -> Pt, 1.0);

    return Circle;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Approximates a unit circle as a cubic polynomial Bspline curve.	     M
*   Construct a circle as four 90 degrees arcs of polynomial cubic Bezier    M
* segments using predefined constants.					     M
*   See Faux & Pratt "Computational Geometry for Design and Manufacturing"   M
* for a polynomial approximation to a circle.                                M
*                                                                            *
* PARAMETERS:                                                                M
*   None                                                                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *: A cubic polynomial Bspline curve approximating a unit   M
*                    circle                                                  M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvCreateUnitPCircle, circle                                          M
*****************************************************************************/
CagdCrvStruct *BspCrvCreateUnitPCircle(void)
{
    int i, j, Quad;
    CagdCrvStruct
	*PCircle = BspCrvNew(UNIT_PCIRCLE_LENGTH, UNIT_PCIRCLE_ORDER,
							      CAGD_PT_E3_TYPE);
    CagdRType
	**Points = PCircle -> Points;

    for (i = 0; i < UNIT_PCIRCLE_LENGTH + UNIT_PCIRCLE_ORDER; i++)
	PCircle -> KnotVector[i] = UnitPCircleKnots[i];

    Points[X][0] = Points[X][12] = 1.0;
    Points[Y][0] = Points[Y][12] = 0.0;
    Points[Z][0] = Points[Z][12] = 0.0;

    /* 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);
	Points[X][j] = Points[X][0] * CosAngle;
	Points[Y][j] = Points[X][0] * SinAngle;
	Points[Z][j] = 0.0;
    }

    return PCircle;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Approximates a circle as a cubic polynomial Bspline curve at the specified M
* position and radius.							     M
*   Construct the circle as four 90 degrees arcs of polynomial cubic Bezier  M
* segments using predefined constants.					     M
*   See Faux & Pratt "Computational Geometry for Design and Manufacturing"   M
* for a polynomial approximation to a circle.                                M
*                                                                            *
* PARAMETERS:                                                                M
*   Center:   Of circle to be created.                                       M
*   Radius:   Of circle to be created.                                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A circle approximation centered at Center and radius  M
*                      Radius that is parallel to the XY plane represented   M
*                      as a polynomial cubic Bspline curve.                  M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvCreatePCircle, circle                                              M
*****************************************************************************/
CagdCrvStruct *BspCrvCreatePCircle(CagdPtStruct *Center, CagdRType Radius)
{
    CagdPtStruct OriginPt;
    CagdCrvStruct
	*Circle = BspCrvCreateUnitPCircle();

    /* Do it in two stages: 1. scale, 2. translate */
    OriginPt.Pt[0] = OriginPt.Pt[1] = OriginPt.Pt[2] = 0.0;
    CagdCrvTransform(Circle, OriginPt.Pt, Radius);
    CagdCrvTransform(Circle, Center -> Pt, 1.0);

    return Circle;
}
