/*
//#     Filename:       SpatialIndex.c
//#
//#     The SpatialIndex functions are here.
//#
//#     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.
//#
//#
//#     Modification History:
//#
*/
#define HTM_LIB
#include <stdlib.h>//added 8/15/02
#include <string.h>//added 8/15/02
#include "SpatialIndex.h"
#include "SpatialEdge.h"
/*
// ===========================================================================
//
// Macro definitions for readability
//
// ===========================================================================
*/
#define N(x) si->nodes_.vector_[(x)]
#define V(x) si->vertices_.vector_[si->nodes_.vector_[index].v_[(x)]]
#define IV(x) si->nodes_.vector_[index].v_[(x)]
#define W(x) si->vertices_.vector_[si->nodes_.vector_[index].w_[(x)]]
#define IW(x) si->nodes_.vector_[index].w_[(x)]
#define ICHILD(x) si->nodes_.vector_[index].childID_[(x)]

#define IV_(x) si->nodes_.vector_[si->index_].v_[(x)]
#define IW_(x) si->nodes_.vector_[si->index_].w_[(x)]
#define ICHILD_(x) si->nodes_.vector_[si->index_].childID_[(x)]
#define IOFFSET 9
/*
// ===========================================================================
//
// Member functions for class SpatialIndex
//
// ===========================================================================

/////////////CONSTRUCTOR//////////////////////////////////
*/
SpatialIndex * spatialIndexNew( size_t maxlevel, size_t buildlevel  )
{
  SpatialIndex *si;
  SpatialEdge *edge;
  size_t nodes,vertices;
  size_t pl=0;
  size_t level;
  int i;
  float64 v[6][3] = {     /* set the first 6 vertices */
    {0.0L,  0.0L,  1.0L}, /* 0 */
    {1.0L,  0.0L,  0.0L}, /* 1 */
    {0.0L,  1.0L,  0.0L}, /* 2 */
   {-1.0L,  0.0L,  0.0L}, /* 3 */
    {0.0L, -1.0L,  0.0L}, /* 4 */
    {0.0L,  0.0L, -1.0L}  /* 5 */
  };

  /* get space for spatial index */
  si = (SpatialIndex *)malloc(sizeof(SpatialIndex));
  si->maxlevel_ = maxlevel;
  si->buildlevel_ = (buildlevel == 0 || buildlevel > maxlevel) ? maxlevel 
	                                                  : buildlevel;

  /* allocate enough space for structures now */
  si->layers_.vector_ = (SpatialLayer *)malloc((si->buildlevel_+1) * 
					       sizeof(SpatialLayer));
  spatialIndexVMax( si, &nodes, &vertices );

  si->nodes_.vector_ = (SpatialQuadNode *)malloc((nodes+1) * 
						 sizeof(SpatialQuadNode));
  si->nodes_.length_ = nodes+1;
  si->vertices_.vector_ = (SpatialVector *)malloc((vertices+1) *
						  sizeof(SpatialVector));
  si->vertices_.length_   = vertices+1;
  si->vertices_.capacity_ = vertices+1;

  /* initialize invalid node */
  N(0).index_ = 0;

  /* initialize first layer */
  si->layers_.vector_[0].level_ = 0;
  si->layers_.vector_[0].nVert_ = 6;
  si->layers_.vector_[0].nNode_ = 8;
  si->layers_.vector_[0].nEdge_ = 12;
  si->layers_.vector_[0].firstIndex_ = 1;
  si->layers_.vector_[0].firstVertex_ = 0;
  si->layers_.length_++;

  for(i = 0; i < 6; i++)
    spatialVectorSet( &(si->vertices_.vector_[i]), v[i][0], v[i][1], v[i][2]);

  /* create the first 8 nodes - index 1 through 8 */
  si->index_ = 1;
  spatialIndexNewNode(si,1,5,2,8,0);  /* S0 */
  spatialIndexNewNode(si,2,5,3,9,0);  /* S1 */
  spatialIndexNewNode(si,3,5,4,10,0); /* S2 */
  spatialIndexNewNode(si,4,5,1,11,0); /* S3 */
  spatialIndexNewNode(si,1,0,4,12,0); /* N0 */
  spatialIndexNewNode(si,4,0,3,13,0); /* N1 */
  spatialIndexNewNode(si,3,0,2,14,0); /* N2 */
  spatialIndexNewNode(si,2,0,1,15,0); /* N3 */

  /*  loop through buildlevel steps, and build the nodes for each layer */
  level = si->buildlevel_;
  while(level-- > 0) {
    edge = spatialEdgeNew( si, pl );
    spatialEdgeMakeMidPoints( edge );
    spatialIndexMakeNewLayer( si, pl );
    ++pl;
  }
  spatialIndexSortIndex( si );

  return si;
}

/* /////////////SHOWVERTICES/////////////////////////////////
// showVertices: print every vertex to the output stream
*/
void spatialIndexShowVertices( const SpatialIndex *si, FILE *out )
{
  size_t i;
  for(i = 0; i < si->vertices_.length_ - 1; i++) {
    spatialVectorWrite( &(si->vertices_.vector_[i]), out );
    fprintf(out,"\n");
  }
}

/* /////////////NODEVERTEX///////////////////////////////////
// nodeVertex: return index of vertices for a node
*/
void spatialIndexNodeVertexID( const SpatialIndex *si,
			       const size_t index, 
			       size_t * v1, 
			       size_t * v2, 
			       size_t * v3 ) {
  *v1 = IV(0);
  *v2 = IV(1);
  *v3 = IV(2);
}
/* /////////////NODEVERTEX///////////////////////////////////
// nodeVertex: return the vectors of the vertices, based on a bitlist leaf
// index! */
void spatialIndexNodeVertexLeaf( const SpatialIndex *si,
				 const IDTYPE leaf, 
				 SpatialVector * v0, 
				 SpatialVector * v1, 
				 SpatialVector * v2 )
{
  IDTYPE index, id, sid, i;
  char name[HTMNAMEMAX];
  SpatialVector *w0, *w1, *w2;

  if(si->buildlevel_ == si->maxlevel_) {
    index = leaf + IOFFSET;
    spatialVectorAssign( v0, &(V(0)) );
    spatialVectorAssign( v1, &(V(1)) );
    spatialVectorAssign( v2, &(V(2)) );
    return;
  }

  /* buildlevel < maxlevel
   * get the id of the stored leaf that we are in
   * and get the vertices of the node we want */
  id = spatialIndexIdByLeafNumber( si, leaf );
  sid = id >> ((si->maxlevel_ - si->buildlevel_)*2);
  index = sid - si->storedleaves_ + IOFFSET;
  spatialVectorAssign( v0, &(V(0)) );
  spatialVectorAssign( v1, &(V(1)) );
  spatialVectorAssign( v2, &(V(2)) );

  /* get name of leaf node */
  spatialIndexNameByLeafNumber( si, leaf, name );

  /* loop through additional name characters and
   * pick the correct triangle accordingly, storing the
   * vertices in v1,v2,v3 */
  w0 = spatialVectorNew();
  w1 = spatialVectorNew();
  w2 = spatialVectorNew();
  for( i = si->buildlevel_ + 2; i < si->maxlevel_ + 2; i++) {
    spatialVectorAddNorm( v1, v2, w0 ) ;
    spatialVectorAddNorm( v0, v2, w1 ) ;
    spatialVectorAddNorm( v1, v0, w2 ) ;

    switch(name[i]) {
    case '0':
      spatialVectorAssign( v1, w2 );
      spatialVectorAssign( v2, w1 );
      break;
    case '1':
      spatialVectorAssign( v0, v1 );
      spatialVectorAssign( v1, w0 );
      spatialVectorAssign( v2, w2 );
      break;
    case '2':
      spatialVectorAssign( v0, v2 );
      spatialVectorAssign( v1, w1 );
      spatialVectorAssign( v2, w0 );
      break;
    case '3':
      spatialVectorAssign( v0, w0 );
      spatialVectorAssign( v1, w1 );
      spatialVectorAssign( v2, w2 );
      break;
    }
  }
  spatialVectorDelete( w0 );
  spatialVectorDelete( w1 );
  spatialVectorDelete( w2 );

}

/* /////////////MAKENEWLAYER/////////////////////////////////
// makeNewLayer: generate a new layer and the nodes in it
*/
void spatialIndexMakeNewLayer( SpatialIndex *si, size_t oldlayer )
{
  IDTYPE index, id, ioffset;
  size_t newlayer = oldlayer + 1;

  si->layers_.vector_[newlayer].level_  = si->layers_.vector_[oldlayer].level_+1;
  si->layers_.vector_[newlayer].nVert_  = si->layers_.vector_[oldlayer].nVert_ + 
                                  si->layers_.vector_[oldlayer].nEdge_;
  si->layers_.vector_[newlayer].nNode_  = 4 * si->layers_.vector_[oldlayer].nNode_;
  si->layers_.vector_[newlayer].nEdge_  = si->layers_.vector_[newlayer].nNode_ + 
                                  si->layers_.vector_[newlayer].nVert_ - 2;
  si->layers_.vector_[newlayer].firstIndex_ = si->index_;
  si->layers_.vector_[newlayer].firstVertex_ = si->layers_.vector_[oldlayer].firstVertex_ + 
                                       si->layers_.vector_[oldlayer].nVert_;
  si->layers_.length_++;

  ioffset = si->layers_.vector_[oldlayer].firstIndex_ ; 
  for(index = ioffset;
      index < ioffset + si->layers_.vector_[oldlayer].nNode_; index++){
    id = N(index).id_ << 2;
    ICHILD(0) = spatialIndexNewNode(si,IV(0),IW(2),IW(1),id++,index);
    ICHILD(1) = spatialIndexNewNode(si,IV(1),IW(0),IW(2),id++,index);
    ICHILD(2) = spatialIndexNewNode(si,IV(2),IW(1),IW(0),id++,index);
    ICHILD(3) = spatialIndexNewNode(si,IW(0),IW(1),IW(2),id,index);
  }
}

/* /////////////NEWNODE//////////////////////////////////////
// newNode: make a new node
*/
size_t spatialIndexNewNode( SpatialIndex *si, size_t v1, size_t v2,size_t v3,
			    IDTYPE id,IDTYPE parent )
{
  IV_(0) = v1;		/* vertex indices */
  IV_(1) = v2;
  IV_(2) = v3;
  IW_(0) = 0;		/* middle point indices */
  IW_(1) = 0;
  IW_(2) = 0;
  ICHILD_(0) = 0;	/* child indices */
  ICHILD_(1) = 0;	/* index 0 is invalid node. */
  ICHILD_(2) = 0;
  ICHILD_(3) = 0;

  N(si->index_).id_ = id;		/* set the id */
  N(si->index_).index_ = si->index_;	/* set the index */
  N(si->index_).parent_ = parent;	/* set the parent */
  return si->index_++;
}


float64 spatialIndexAreaID( const SpatialIndex *si, IDTYPE ID )
{
  size_t leaf = spatialIndexLeafNumberById( si, ID );
  float64 area;

  SpatialVector n0;
  SpatialVector n1;
  SpatialVector n2;

  spatialIndexNodeVertexLeaf( si, leaf, &n0, &n1, &n2);
  return spatialIndexArea(si,&n0,&n1,&n2);
}


/* /////////////AREA////////////////////////////////////////
// area: routine to precompute the area of a node using
//
//   AREA = 4*arctan sqrt(tan(s/2)tan((s-a)/2)tan((s-b)/2)tan((s-c)/2))
//
//   with s = (a+b+c)/2
//
// (with many thanks to Eduard Masana, emasana@pchpc10.am.ub.es )
*/
float64 spatialIndexArea( const SpatialIndex *si,
			  const SpatialVector * v0, 
			  const SpatialVector * v1, 
			  const SpatialVector * v2 ) {

  float64 a = acos( spatialVectorDot( v0 , v1) );
  float64 b = acos( spatialVectorDot( v1 , v2) );
  float64 c = acos( spatialVectorDot( v2 , v0) );

  float64 s = (a + b + c)/2.0;

  float64 area = 4.0*atan(sqrt(tan(s/2.0)*
			       tan((s-a)/2.0)*
			       tan((s-b)/2.0)*
			       tan((s-c)/2.0)));        
  return area;
}

/* /////////////VMAX/////////////////////////////////////////
// vMax: compute the maximum number of vertices for the
//       polyhedron after buildlevel of subdivisions and
//       the total number of nodes that we store
//       also, calculate the number of leaf nodes that we eventually have.
*/
void spatialIndexVMax( SpatialIndex *si, size_t *nodes, size_t *vertices )
{
  IDTYPE nv = 6;    /* initial values */
  IDTYPE ne = 12;
  IDTYPE nf = 8;
  int32 i  = si->buildlevel_;
  *nodes = nf;

  while(i-->0){
    nv += ne;
    nf *= 4;
    ne  = nf + nv -2;
    *nodes += nf;
  }
  *vertices = nv;
  si->storedleaves_ = nf;

  /* calculate number of leaves */
  i = si->maxlevel_ - si->buildlevel_;
  while(i-- > 0)
    nf *= 4;
  si->leaves_ = nf;
}

/* /////////////SORTINDEX////////////////////////////////////
// sortIndex: sort the index so that the first node is the invalid node
//            (index 0), the next 8 nodes are the root nodes
//            and then we put all the leaf nodes in the following block
//            in ascending id-order.
//           All the rest of the nodes is at the end.
*/
void spatialIndexSortIndex( SpatialIndex *si ) {
  size_t index;
  size_t nonleaf;
  size_t leaf, i;

  /* create a copy of the node list */
  SpatialQuadNodeVec oldnodes;
  oldnodes.vector_ = (SpatialQuadNode *)malloc (si->nodes_.length_ *
						sizeof(SpatialQuadNode));
  memcpy(oldnodes.vector_, si->nodes_.vector_, 
	 si->nodes_.length_ * sizeof(SpatialQuadNode));
  oldnodes.length_ = si->nodes_.length_;

#define ON(x) oldnodes.vector_[(x)]

  /* now refill the nodes_ list according to our sorting. */
  for( index=IOFFSET, leaf=IOFFSET, nonleaf=si->nodes_.length_-1; 
       index < si->nodes_.length_; index++) {

    if( ON(index).childID_[0] == 0 ) { /* childnode */
      /* set leaf into list */
      N(leaf) = ON(index);
      /* set parent's pointer to this leaf */
      for ( i = 0; i < 4; i++) {
	if(N(N(leaf).parent_).childID_[i] == index) {
	  N(N(leaf).parent_).childID_[i] = leaf;
	  break;
	}
      }
      leaf++;
    } else {
      /* set nonleaf into list from the end
       * set parent of the children already to this
       * index, they come later in the list. */
      N(nonleaf) = ON(index);
      ON(N(nonleaf).childID_[0]).parent_ = nonleaf;
      ON(N(nonleaf).childID_[1]).parent_ = nonleaf;
      ON(N(nonleaf).childID_[2]).parent_ = nonleaf;
      ON(N(nonleaf).childID_[3]).parent_ = nonleaf;
      /* set parent's pointer to this leaf */
      for ( i = 0; i < 4; i++) {
	if(N(N(nonleaf).parent_).childID_[i] == index) {
	  N(N(nonleaf).parent_).childID_[i] = nonleaf;
	  break;
	}
      }
      nonleaf--;
    }
  }
  free(oldnodes.vector_);
}
/*
//////////////////IDBYNAME/////////////////////////////////////////////////
// Translate ascii leaf name to a uint32
//
// The following encoding is used:
//
// The string leaf name has the always the same structure, it begins with
// an N or S, indicating north or south cap and then numbers 0-3 follow 
// indicating which child to descend into. So for a depth-5-index we have
// strings like
//                 N012023  S000222  N102302  etc
//
// Each of the numbers correspond to 2 bits of code (00 01 10 11) in the
// uint32. The first two bits are 10 for S and 11 for N. For example
//
//                 N 0 1 2 0 2 3
//                 11000110001011  =  12683 (dec)
//
// The leading bits are always 0.
//
// --- WARNING: This works only up to 15 levels. 
//              (we probably never need more than 7)
//
*/
IDTYPE spatialIndexIdByName( const char *name )
{
  IDTYPE out=0, i;
  uint32 size = 0;

  if(name == 0) return out;  /* null pointer-name */
  if(name[0] != 'N' && name[0] != 'S') return out; /* invalid name */

  size = strlen(name);       /* determine string length */
  /* at least size-2 required, don't exceed max */
  if(size < 2 || size > HTMNAMEMAX) return out;

  for(i = size-1; i > 0; i--) {/* set bits starting from the end */
    if(name[i] > '3' || name[i] < '0')return 0; /* invalid name */
    out += ((IDTYPE)(name[i]-'0')) << 2*(size - i -1);
  }

  i = 2;                     /* set first pair of bits, first bit always set */
  if(name[0]=='N') i++;      /* for north set second bit too */
  out += (i << (2*size - 2) );

  /************************
  // This code may be used later for hashing !
  if(size==2)out -= 8;
  else {
    size -= 2;
    uint32 offset = 0, level4 = 8;
    for(i = size; i > 0; i--) { // calculate 4 ^ (level-1), level = size-2
      offset += level4;
      level4 *= 4;
    }
    out -= level4 - offset;
  }
  **************************/
  return out;
}

/*
//////////////////NAMEBYID/////////////////////////////////////////////////
// Translate uint32 to an ascii leaf name
//
// The encoding described above may be decoded again using the following
// procedure:
//
//  * Traverse the uint32 from left to right.
//  * Find the first 'true' bit.
//  * The first pair gives N (11) or S (10).
//  * The subsequent bit-pairs give the numbers 0-3.
//
*/

char * spatialIndexNameById( IDTYPE id, char * name )
{
  uint32 size=0, i;
  bool ok=true;


  /*************
  // This code might be useful for hashing later !!

  // calculate the level (i.e. 8*4^level) and add it to the id:
  uint32 level=0, level4=8, offset=8;
  while(id >= offset) {
    if(++level > 13) { ok = false; offset = 0; break; }// level too deep
    level4 *= 4;
    offset += level4;
  }
  id += 2 * level4 - offset;
  **************/

  /* determine index of first set bit */
  for(i = 0; i < IDSIZE; i+=2) {
    if ( (id << i) & IDHIGHBIT ) break;
    if ( (id << i) & IDHIGHBIT2 ) {  /* invalid id */
      ok = false;
      break;
    }
  }
  if(!ok || id == 0) {
    if (!name)
      name = (char *)malloc(sizeof(char)); /* single character, zero. */
    *name = 0;
    return name;
  }

  size=(IDSIZE-i) >> 1;
  /* allocate characters */
  if(!name)
    name = (char *)malloc(sizeof(char) *(size+1));

  /* fill characters starting with the last one */
  for(i = 0; i < size-1; i++)
    name[size-i-1] = '0' + (char)( (id >> i*2) & 3);

  /* put in first character */
  if( (id >> (size*2-2)) & 1 ) {
    name[0] = 'N';
  } else {
    name[0] = 'S';
  }
  name[size] = 0; /* end string */

  return name;
}
/*
//////////////////IDBYPOINT////////////////////////////////////////////////
// Find a leaf node where a vector points to
*/

IDTYPE spatialIndexIdByPoint( const SpatialIndex *si, 
			      const SpatialVector * v ) {
  size_t index, oldindex, i, len, level;
  SpatialVector t;
  SpatialVector v0, v1, v2;
  SpatialVector w0, w1, w2;
  char name[HTMNAMEMAX];

  /* start with the 8 root triangles, find the one which v points to */
  for(index=1; index <=8; index++) {
    spatialVectorCross ( &(V(0)), &(V(1)), &t );
    if( spatialVectorDot( &t, v) < -gEpsilon) continue;
    spatialVectorCross ( &(V(1)), &(V(2)), &t );
    if( spatialVectorDot( &t, v) < -gEpsilon) continue;
    spatialVectorCross ( &(V(2)), &(V(0)), &t );
    if( spatialVectorDot( &t, v) < -gEpsilon) continue;
    break;
  }
  /* loop through matching child until leaves are reached */
  while(ICHILD(0)!=0) {
    oldindex = index;
    for(i = 0; i < 4; i++) {
      index = si->nodes_.vector_[oldindex].childID_[i];
      spatialVectorCross ( &(V(0)), &(V(1)), &t );
      if( spatialVectorDot( &t, v) < -gEpsilon) continue;
      spatialVectorCross ( &(V(1)), &(V(2)), &t );
      if( spatialVectorDot( &t, v) < -gEpsilon) continue;
      spatialVectorCross ( &(V(2)), &(V(0)), &t );
      if( spatialVectorDot( &t, v) < -gEpsilon) continue;
      break;
    }
  }
  /* return if we have reached maxlevel */
  if(si->maxlevel_ == si->buildlevel_)return N(index).id_;

  /* from now on, continue to build name dynamically.
   * until maxlevel_ levels depth, continue to append the
   * correct index, build the index on the fly. */
  spatialIndexNameById(N(index).id_,name);
  len = strlen(name);
  spatialVectorAssign( &v0 , &(V(0)) );
  spatialVectorAssign( &v1 , &(V(1)) );
  spatialVectorAssign( &v2 , &(V(2)) );

  level = si->maxlevel_ - si->buildlevel_;
  while(level--) {
    spatialVectorAddNorm( &v1, &v2, &w0 );
    spatialVectorAddNorm( &v0, &v2, &w1 );
    spatialVectorAddNorm( &v1, &v0, &w2 );

    if(spatialIndexIsInside(si, v, &v0, &w2, &w1)) {
      name[len++] = '0';
      spatialVectorAssign( &v1, &w2 );
      spatialVectorAssign( &v2, &w1 );
      continue;
    } else if(spatialIndexIsInside(si, v, &v1, &w0, &w2)) {
      name[len++] = '1';
      spatialVectorAssign( &v0, &v1 );
      spatialVectorAssign( &v1, &w0 );
      spatialVectorAssign( &v2, &w2 );
      continue;
    } else if(spatialIndexIsInside(si, v, &v2, &w1, &w0)) {
      name[len++] = '2';
      spatialVectorAssign( &v0, &v2 );
      spatialVectorAssign( &v1, &w1 );
      spatialVectorAssign( &v2, &w0 );
      continue;
    } else if(spatialIndexIsInside(si, v, &w0, &w1, &w2)) {
      name[len++] = '3';
      spatialVectorAssign( &v0, &w0 );
      spatialVectorAssign( &v1, &w1 );
      spatialVectorAssign( &v2, &w2 );
      continue;
    }
  }
  name[len] = '\0';
  return spatialIndexIdByName(name);
}
/*
//////////////////ISINSIDE/////////////////////////////////////////////////
// Test whether a vector is inside a triangle. 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)
{
  SpatialVector t;
  spatialVectorCross ( v0, v1, &t );
  if( spatialVectorDot( &t, v ) < -gEpsilon) return false;
  spatialVectorCross ( v1, v2, &t );
  if( spatialVectorDot( &t, v ) < -gEpsilon) return false;
  spatialVectorCross ( v2, v0, &t );
  if( spatialVectorDot( &t, v ) < -gEpsilon) return false;
  return true;
}


/* ////////////////LEAFNUMBERBYID///////////////////////////////////////////
 */
IDTYPE spatialIndexLeafNumberById( const SpatialIndex *si, IDTYPE id ) {
  return (id - si->leaves_);
}

/* //////////////////IDBYLEAFNUMBER///////////////////////////////////////////
 */
IDTYPE spatialIndexIdByLeafNumber( const SpatialIndex *si, IDTYPE n ) {
  return (si->leaves_ + n);
}

/* //////////////////IDBYLEAFNUMBER///////////////////////////////////////////
 */
char * spatialIndexNameByLeafNumber( const SpatialIndex *si, IDTYPE n, char * name )
{
  return spatialIndexNameById( spatialIndexIdByLeafNumber(si,n), name);
}

/* //////////////////IDBYPOINT////////////////////////////////////////////////
// Find a leaf node where a ra/dec points to
*/
IDTYPE spatialIndexIdByPointRaDec( const SpatialIndex *si, 
				   const float64 ra, 
				   const float64 dec ) {
  SpatialVector v;
  spatialVectorSetRaDec(&v,ra,dec);
  return spatialIndexIdByPoint(si,&v);
}

/* //////////////////NAMEBYPOINT//////////////////////////////////////////////
// Find a leaf node where a ra/dec points to, return its name
*/
char* spatialIndexNameByPointRaDec( const SpatialIndex *si,
				    const float64 ra, const float64 dec) {
  return spatialIndexNameById(spatialIndexIdByPointRaDec(si,ra,dec),NULL);
}

/* //////////////////NAMEBYPOINT//////////////////////////////////////////////
// Find a leaf node where v points to, return its name
*/

char* spatialIndexNameByPoint( const SpatialIndex *si, 
			       const SpatialVector * vector ){
  return spatialIndexNameById(spatialIndexIdByPoint(si,vector),NULL);
}