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

#ifdef __MSDOS__
#include <stdlib.h>
#endif /* __MSDOS__ */

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

#define UNIT_CIRCLE_ORDER 3
#define UNIT_CIRCLE_LENGTH 9	  /* Nine 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 };

/******************************************************************************
* Creates an arc at the specified position as a rational Bezier curve.	      *
*   The arc is assumed to be less than 180 degrees from Start to End in the   *
* shorter path as arc where Center as arc center.			      *
******************************************************************************/
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);
	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;
}

/******************************************************************************
* Creates a circle at the specified position as a rational Bspline curve.     *
*   Construct a circle as 4 90 degrees arcs of rationa bezier segments using  *
* the constants defined above.						      *
******************************************************************************/
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;
}


/******************************************************************************
* Creates a circle at the specified position as a rational bezier curve.      *
******************************************************************************/
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;
}