#ifndef _SpatialIndex_h
#define _SpatialIndex_h
/*
//#     Filename:       SpatialIndex.h
//#
//#     SpatialIndex is the class for the the sky indexing routines.
//#
//#
//#     Author:         Peter Z. Kunszt, based on A. Szalay s code
//#
//#     Date:           October 15, 1998
//#
//#
//#
//# (c) Copyright The Johns Hopkins University 1998
//# All Rights Reserved
//#
//# The software and information contained herein are proprietary to The
//# Johns Hopkins University, Copyright 1998.  This software is furnished
//# pursuant to a written license agreement and may be used, copied,
//# transmitted, and stored only in accordance with the terms of such
//# license and with the inclusion of the above copyright notice.  This
//# software and information or any other copies thereof may not be
//# provided or otherwise made available to any other person.
//#
//#
*/

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include "SpatialVector.h"

/*
//########################################################################
//
// <GROUP>
// <SUMMARY>Class declarations</SUMMARY>
// 

//########################################################################
//
// <SUMMARY> Spatial Index class </SUMMARY>
//

// The Spatial Index is a quad tree of spherical triangles. The tree
// is built in the following way: Start out with 8 triangles on the
// sphere using the 3 main circles to determine them. Then, every
// triangle can be decomposed into 4 new triangles by drawing main
// circles between midpoints of its edges:
//<srcblock>
//
//                            /\
//                           /  \
//                          /____\
//                         /\    /\
//                        /  \  /  \
//                       /____\/____\
//
//</srcblock>
// This is how the quad tree is built up to a certain level by decomposing
// every triangle again and again.
//
*/

typedef struct  {
  size_t 	level_;		/* layer level */
  size_t 	nVert_;		/* number of vertices in this layer */
  IDTYPE 	nNode_;		/* number of nodes */
  size_t 	nEdge_;		/* number of edges */
  IDTYPE 	firstIndex_;	/* index of first node of this layer */
  size_t 	firstVertex_;	/* index of first vertex of this layer */
} SpatialLayer;

typedef struct  {
  IDTYPE	index_;		/* its own index */
  size_t	v_[3];		/* The three vertex vector indices */
  size_t	w_[3];		/* The three middlepoint vector indices */
  IDTYPE	childID_[4];	/* ids of children */
  IDTYPE	parent_;	/* id of the parent node (needed for sorting)*/
  IDTYPE	id_;		/* numeric id -> name */
} SpatialQuadNode;

typedef struct {
  SpatialQuadNode *vector_;
  size_t 	  length_;
} SpatialQuadNodeVec;

typedef struct {
  SpatialLayer  *vector_;
  size_t 	length_;
} SpatialLayerVec;

typedef struct {
  size_t 		maxlevel_;	/* the depth of the Layer */
  size_t 	        buildlevel_;	/* the depth of the Layer stored */
  IDTYPE		leaves_;	/* number of leaf nodes */
  IDTYPE		storedleaves_;	/* number of stored leaf nodes */
  SpatialQuadNodeVec 	nodes_;		/* the array of nodes */
  SpatialLayerVec	layers_;	/* array of layers */
  SpatialVectorVec	vertices_;	/* array of vertices */
  IDTYPE 		index_;		/* the current index_ of vertices */
} SpatialIndex ;

/* Constructor : give the level of the index
 * and optionally the level to build - i.e. the depth to keep in memory.
 * if maxlevel - buildlevel > 0 , that many levels are generated on the
 * fly each time the index is called.
 */
SpatialIndex * spatialIndexNew( size_t maxlevel, size_t buildlevel );

/* name string conversion to int32/64 */
IDTYPE spatialIndexIdByName( const char * );

/* int32/64 conversion to a string (name of database)
 * WARNING: if name is already allocated, a size of at least 17 is required.
 * The conversion is done by directly calculating the name from a number.
 * To calculate the name of a certain level, the mechanism is that the
 * name is given by (#of nodes in that level) + (id of node).
 * So for example, for the first level, there are 8 nodes, and we get
 * the names from numbers 8 through 15 giving S0,S1,S2,S3,N0,N1,N2,N3.
 * The order is always ascending starting from S0000.. to N3333...
 */
char * spatialIndexNameById( IDTYPE ID, char * name );

/* return leaf number in bitlist for a certain ID. Since the ID here
 * means the number computed from the name, this is simply returning
 * ID -leafCount().
 */
IDTYPE spatialIndexLeafNumberById( const SpatialIndex *si, IDTYPE ID );

/* return leaf id for a certain bitlist index. Same as the function above */
IDTYPE spatialIndexIdByLeafNumber( const SpatialIndex *si, IDTYPE n ) ;

/* return name for a certain leaf index (to be used for name lookup from
 * a bitlist).
 * This function is simply shorthand for nameById(n + leafCount()).
 */
char * spatialIndexNameByLeafNumber( const SpatialIndex *si, 
				     IDTYPE n, char * name );

/* find a node by giving a vector. The ID of the node is returned. */
IDTYPE spatialIndexIdByPoint( const SpatialIndex *si, const 
			      SpatialVector * vector );

/* find a node by giving a ra,dec in degrees. */
IDTYPE spatialIndexIdByPointRaDec( const SpatialIndex *si, 
				   const float64 ra, 
				   const float64 dec );

/* find a node by giving a vector. The ID of the node is returned. */
char* spatialIndexNameByPoint( const SpatialIndex *si, 
			       const SpatialVector * vector );

/* find a node by giving a ra,dec in degrees. */
char* spatialIndexNameByPointRaDec( const SpatialIndex *si, 
				    const float64 ra, const float64 dec);

/* The area in steradians for a given index ID */
float64 spatialIndexAreaID( const SpatialIndex *si, IDTYPE ID );

/* The area in steradians for a given spatial triangle */
float64 spatialIndexArea( const SpatialIndex *si,
			  const SpatialVector * v1, 
			  const SpatialVector * v2, 
			  const SpatialVector * v3 );

/* return the actual vertex vectors */
void spatialIndexNodeVertexLeaf( const SpatialIndex *si,
				 const IDTYPE leaf, 
				 SpatialVector * v1, 
				 SpatialVector * v2, 
				 SpatialVector * v3 ); 

/* return index of vertices for a node */
void spatialIndexNodeVertexID( const SpatialIndex *si,
			       const size_t idx, 
			       size_t * v1, 
			       size_t * v2, 
			       size_t * v3 );

/* print all vertices to output stream */
void spatialIndexShowVertices( const SpatialIndex *si, FILE *out );


/* insert a new node_[] into the list. The vertex indices are given by
 * v1,v2,v3 and the id of the node is set.
 */
size_t spatialIndexNewNode( SpatialIndex *si, size_t v1, size_t v2,size_t v3,
			    IDTYPE id,IDTYPE parent );

/* make new nodes in a new layer. */
void spatialIndexMakeNewLayer( SpatialIndex *si, size_t oldlayer );

/* return the total number of nodes and vertices */
void spatialIndexVMax( SpatialIndex *si, size_t *nodes, size_t *vertices );

/* sort the index so that the leaf nodes are at the beginning */
void spatialIndexSortIndex( SpatialIndex *si );

/* Test whether a vector v is inside a triangle v0,v1,v2. Input
 * triangle has to be sorted in a counter-clockwise direction.
 */
bool spatialIndexIsInside( const SpatialIndex *si, 
			   const SpatialVector * v, 
			   const SpatialVector * v0,
			   const SpatialVector * v1, 
			   const SpatialVector * v2);


#endif