/******************************************************************************
* Sph_Pts.c - distribute uniformly points on a sphere.			      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, May 1996.					      *
******************************************************************************/

#include <math.h>
#include <stdio.h>
#include "irit_sm.h"
#include "geomat3d.h"
#include "iritprsr.h"
#include "allocate.h"

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes a honey comb distribution of points on a sphere. Result is      M
* an approximation.						             M
*                                                                            *
* PARAMETERS:                                                                M
*   HoneyCombSize:      Size of honey comb.                                  M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:    A point list object.                                M
*                                                                            *
* KEYWORDS:                                                                  M
*   CGPointCoverOfUnitHemiSphere                                             M
*****************************************************************************/
IPObjectStruct *CGPointCoverOfUnitHemiSphere(RealType HoneyCombSize)
{
    RealType
	Cos30 = cos(M_PI / 6),
	Alpha = 2 * asin(HoneyCombSize / 2);
    int x, y,
	d = (int) (20 * M_PI / Alpha);
    IPPolygonStruct *Poly;
    IPVertexStruct *V,
	*VHead = NULL;
    IPObjectStruct *PObj;

    for (y = -d; y <= d; y++) {
	for (x = -d; x <= d; x++) {
	    RealType CoordZ = 0, DAlpha,
		CoordX = x + (y & 0x01 ? 0.5 : 0.0),
		CoordY = Cos30 * y,
		Len = sqrt(SQR(CoordX) + SQR(CoordY));

	    if (APX_EQ(Len, 0.0))
	        Len = IRIT_EPS;

	    if ((DAlpha = Len * Alpha) < M_PI / 2.0) {
		RealType
		    NewLen = tan(DAlpha);

		/* Angular span fits the southern hemisphere. */
		CoordX *= NewLen / Len;
		CoordY *= NewLen / Len;

		Len = sqrt(SQR(CoordX) + SQR(CoordY) + 1.0);
		CoordX /= Len;
		CoordY /= Len;
		CoordZ = 1.0 / Len;

		V = IPAllocVertex(0, NULL, VHead);
		VHead = V;
		V -> Coord[0] = CoordX;
		V -> Coord[1] = CoordY;
		V -> Coord[2] = CoordZ;
	    }
	}
    }

    Poly = IPAllocPolygon(0, VHead, NULL);

    PObj = GenPOLYObject(Poly);
    IP_SET_POINTLIST_OBJ(PObj);
    return PObj;
}