using System;
using Spherical;
/*=====================================================================

  File:      Trixel.cs for Spatial Sample
  Summary:   Implements conversions between several representations of a HID
  Date:	     August 16, 2005

---------------------------------------------------------------------
  This spatial search library was developed by Alexander Szalay, Robert Brunner, 
  Peter Kunszt, and Gyorgy Fekete at the Department of Physics and Astronomy, 
  The Johns Hopkins University, in collaboration with Jim Gray, Microsoft Research.
  Details of the algorithms can be found at http://www.sdss.jhu.edu/htm/.
 
  
 * Copyright (C) The Johns Hopkins University.
 * No warranty of the suitability of the code for anything is given
 * This code may be used or modified for any purpose, as long as these
 * credits are included in all future distributions. The authors request any
 * modifications, enhancements to this software to be brought to the attention
 * of mailto:gfekete@pha.jhu.edu
 * 
======================================================= */

namespace Spherical.HTM {
	/// <summary>
	/// Trixel static methods operate on trixels.
	/// 
	/// operations needed to convert between location and HtmID
	/// several representations of a HID, and locations expressed
	/// in either Cartesian or J2000 (ra/dec) space.
	/// </summary>
	public class Trixel {
		/// 
		private enum Startindex {
			_S2 = 0, 
			_N1 = 1,
			_S1 = 2,
			_N2 = 3,
			_S3 = 4,
			_N0 = 5,
			_S0 = 6,
			_N3 = 7
		};
		/// <summary>
		/// Number of bits in an HID
		/// </summary>
		const int eHIDBits = 64;

		/// <summary>
		/// Maximum number of characters in text descriptions
		/// Whenever character arrays are used to keep the text
		/// description of a trixel, they must be exactly eMaxNameSize
		/// characters long
		/// </summary>
		public const int eMaxNameSize = 32;
		private const long IDHIGHBIT  = 0x2000000000000000L;
		private const long IDHIGHBIT2 = 0x1000000000000000L;

		// Private constuctor, all functions are performed with static methods.
		private Trixel() {
		}
		/// <summary>
		/// Convert a Cartesian coordinate to a HID.
		/// 
        /// Given a Cartesian coordinate
		/// and a level number, it returns the 64 bit HID.
		/// <b>WARNING!</b>
		/// x, y, z are assumed to be normalized, so this function
		/// doesn't waste time normalzing.
		/// </summary>
		/// <param name="x">X coordinate of location</param>
		/// <param name="y">Y coordinate of location</param>
		/// <param name="z">Z coordinate of location</param>
		/// <param name="depth">The level of the HID</param>
		/// <returns>64-bit HID</returns>
		public static Int64 CartesianToHid(double x, double y, double z, int depth){
			int len;
			return _cartesianToNameOrHid(x, y, z, depth, null, out len);
		}
        /// <summary>
        /// Find the string representation of the HtmID to a given level of a unit vector
        /// 
        /// The vector <strong>must</strong> be normalized.
        /// </summary>
        /// <param name="x">double</param>
        /// <param name="y">double</param>
        /// <param name="z">double</param>
        /// <param name="depth">int</param>
        /// <returns>text representation of HtmID</returns>
        /// 
        public static String CartesianToString(double x, double y, double z, int depth) {
            char[] name = new char[eMaxNameSize];
			int  len;
			_cartesianToNameOrHid(x, y, z, depth, name, out len);
			return new String(name, 0, len);
		}
        /// <summary>
        /// Convert a HtmID to its textual representation
        /// </summary>
        /// <param name="hid">64-bit HtmID</param>
        /// <returns>The textual representation</returns>
        public static String ToString(Int64 hid) {
            char[] name = new char[eMaxNameSize];
            int len;
            len = ToName(name, hid);
            return new String(name, 0, len);
        }
        /// <summary>
        /// Decide whether or not a node in the tree of a given HtmID is an 
        /// ancestor of another node
        /// </summary>
        /// <param name="grandpa">64-bit HtmID, candidate for ancestry</param>
        /// <param name="hid">64-bit HtmID who is the potential descendant</param>
        /// <returns>true if grandpa is an ancestor of hid</returns>
        public static bool IsAncestor(Int64 grandpa, Int64 hid) {
            int shifts = LevelOfHid(hid) - LevelOfHid(grandpa);
            if(shifts < 0)
                return false;
            Int64 descendant = hid >> (2 * shifts);
            return (descendant == grandpa);
        }
		private static Int64 _cartesianToNameOrHid(double x, double y, double z, int depth, char[] name, out int len) {

			long startID;
			Int64 hid;
			Cartesian v0, v1, v2;
			Cartesian w0, w1, w2;
			Cartesian p = new Cartesian(x, y, z, false);

			// Get the ID of the level0 triangle, and its starting vertices
			//
			//v0 = Cartesian.NaN;
			//v1 = Cartesian.NaN;
			//v2 = Cartesian.NaN;
			//w0 = Cartesian.NaN;
			//w1 = Cartesian.NaN;
			//w2 = Cartesian.NaN;
			w0 = Cartesian.NaN;
			w1 = Cartesian.NaN;
			w2 = Cartesian.NaN;
			startID = _startpane(out v0, out v1, out v2, x, y, z, name);
			// We have two copies of almost identical code, for speed's sake
			//
			len = 2;
			if (name == null) {
				if (startID < 8) {
					hid = 1;
					return hid;
				}
				hid = startID;
				//
				// Start searching for the children
				///
				while (depth-- > 0) {
					hid <<= 2;
					w2.SetMiddlePoint(v0, v1, true);
					w0.SetMiddlePoint(v1, v2, true);
					w1.SetMiddlePoint(v2, v0, true);
					if (Contains(p, v0, w2, w1)) {
						// hid |= 0;

						v1.Set(w2, false);
						v2.Set(w1, false);
					} else if (Contains(p, v1, w0, w2)) {
						hid |= 1;
						v0.Set(v1, false);
						v1.Set(w0, false);
						v2.Set(w2, false);
					} else if (Contains(p, v2, w1, w0)) {
						hid |= 2;
						v0.Set(v2, false);
						v1.Set(w1, false);
						v2.Set(w0, false);
					} else if (Contains(p, w0, w1, w2)) {
						hid |= 3;
						v0.Set(w0, false);
						v1.Set(w1, false);
						v2.Set(w2, false);
					} else {
						// CATASTROPHIC ERROR!!! THROW Something!
						throw new Exception("Panic in Cartesian2hid");
					}
				}
			} else {
				if (startID < 8) {
					name[0] = 'X';
					name[1] = 'X';
					name[2] = '\0';
				}
				hid = startID;
				//
				// Start searching for the children
				///
				
				while (depth-- > 0) {
					w2.SetMiddlePoint(v0, v1, true);
					w0.SetMiddlePoint(v1, v2, true);
					w1.SetMiddlePoint(v2, v0, true);
					if (Contains(p, v0, w2, w1)) {
						name[len++] = '0';
						v1.Set(w2, false);
						v2.Set(w1, false);
					} else if (Contains(p, v1, w0, w2)) {
						name[len++] = '1';
						v0.Set(v1, false);
						v1.Set(w0, false);
						v2.Set(w2, false);
					} else if (Contains(p, v2, w1, w0)) {
						name[len++] = '2';
						v0.Set(v2, false);
						v1.Set(w1, false);
						v2.Set(w0, false);
					} else if (Contains(p, w0, w1, w2)) {
						name[len++] = '3';
						v0.Set(w0, false);
						v1.Set(w1, false);
						v2.Set(w2, false);
					} else {
						// CATASTROPHIC ERROR!!! THROW Something!
						throw new Exception("Panic in Cartesian2hid");
					}
				}
				name[len] = '\0';
			}

			return hid;
		}
		/// <summary>
		/// Convert the location given by (x, y, z) to the symbolic text
		/// name of the HID
		/// <strong>WARNING</strong>:
		/// x, y, z are assumed to be normalized, so this function
		/// doesn't waste time normalizing.
		/// </summary>
		/// <param name="nam">Character array that holds the text for
		/// the trixel's name</param>
		/// <param name="x">X coordinate of location</param>
		/// <param name="y">Y coordinate of location</param>
		/// <param name="z">Z coordinate of location</param>
		/// <param name="depth">The level of the HID</param>
		/// <returns>true if conversion was successful</returns>
		
		/// <summary>
		/// Convert the named trixel to a triangle described
		/// by three vertices.
		/// The vertices are given by three arrays of three doubles.
		/// The coordinates of the triangles are given in
		/// the order (x, y, z) and so that the location
		/// is on the surface of a unit sphere
		/// </summary>
		/// <param name="name">The text description of the trixel</param>
		/// <param name="c0">The X coordinate</param>
		/// <param name="c1">The Y coordinate</param>
		/// <param name="c2">The Z coordinate</param>
		/// <returns>true, if the conversion succeeded, false otherwise</returns>
		public static bool NameToTriangle(char[] name, out Cartesian c0, out Cartesian c1, out Cartesian c2) {
			bool rstat = false;

			Cartesian v0, v1, v2;
			Cartesian w0, w1, w2;


			//
			// Get the top level hemi-demi-semi space
			//
			int k;
			int off0, off1, off2;
			k = (int) name[1] - '0';
			c0 = Cartesian.NaN;
			c1 = Cartesian.NaN;
			c2 = Cartesian.NaN;
			if (k < 0 || k > 3) {// Do not have a valid name
				return rstat;
			}
			if (name[0] != 'N' && name[0] != 'S'){
				return rstat;
			}

			if (name[0] == 'N') k += 4;
			off0 = HtmState.Instance.vertexIndices[k,0];
			off1 = HtmState.Instance.vertexIndices[k,1];
			off2 = HtmState.Instance.vertexIndices[k,2];

			v0 = new Cartesian(
				HtmState.Instance.anchor[off0, 0],
				HtmState.Instance.anchor[off0, 1],
				HtmState.Instance.anchor[off0, 2], false);

			v1 = new Cartesian(
				HtmState.Instance.anchor[off1, 0],
				HtmState.Instance.anchor[off1, 1],
				HtmState.Instance.anchor[off1, 2], false);
			v2 = new Cartesian(
				HtmState.Instance.anchor[off2, 0],
				HtmState.Instance.anchor[off2, 1],
				HtmState.Instance.anchor[off2, 2], false);
			w0 = Cartesian.NaN;
			w1 = Cartesian.NaN;
			w2 = Cartesian.NaN;
			k = 2;
			while(name[k] != '\0'){
				w2.SetMiddlePoint(v0, v1, true);
				w0.SetMiddlePoint(v1, v2, true);
				w1.SetMiddlePoint(v2, v0, true);
				switch(name[k]) {
					case '0':
						v1.Set(w2, false);
						v2.Set(w1, false);
						break;
					case '1':
						v0.Set(v1, false);
						v1.Set(w0, false);
						v2.Set(w2, false);
						break;
					case '2':
						v0.Set(v2, false);
						v1.Set(w1, false);
						v2.Set(w0, false);
						break;
					case '3':
						v0.Set(w0, false);
						v1.Set(w1, false);
						v2.Set(w2, false);
						break;
				}
				k++;
			}
			c0 = new Cartesian(v0, false);
			c1 = new Cartesian(v1, false);
			c2 = new Cartesian(v2, false);
			rstat = true;
			return rstat;
		}

		/// <summary>
		/// 
		/// </summary>
		/// <param name="index">a number between 0 and 7 (inclusive)</param>
		/// <param name="v0">First vertex</param>
		/// <param name="v1">Second vertex</param>
		/// <param name="v2">Third vertex</param>
		public static void TopVertices(int index, out Cartesian v0, out Cartesian v1, out Cartesian v2) {
			//
			// index 0 is for S0, 1 for S1, etc...
			//
			// t = vertex_index[index, j] is vertex j
			// xj = HtmState.Instance.anchor(j,k)
			//
			int t;
			//
			//
			t = HtmState.Instance.vertexIndices[index, 0];
			v0 = new Cartesian(
				HtmState.Instance.anchor[t, 0],
				HtmState.Instance.anchor[t, 1],
				HtmState.Instance.anchor[t, 2], false);

			t = HtmState.Instance.vertexIndices[index, 1];
			v1 = new Cartesian(
				HtmState.Instance.anchor[t, 0],
				HtmState.Instance.anchor[t, 1],
				HtmState.Instance.anchor[t, 2], false);

			t = HtmState.Instance.vertexIndices[index, 2];
			v2 = new Cartesian(
				HtmState.Instance.anchor[t, 0],
				HtmState.Instance.anchor[t, 1],
				HtmState.Instance.anchor[t, 2], false);
			return;
		}
		// /////////////////////////////////////////// private METHODS
		// / startpane, strcpy, midpoint, isinside, copyvec
		// /

		/// <summary>
		/// 
		/// </summary>
		/// <param name="v1"></param>
		/// <param name="v2"></param>
		/// <param name="v3"></param>
		/// <param name="xin"></param>
		/// <param name="yin"></param>
		/// <param name="zin"></param>
		/// <param name="name"></param>
		/// <returns></returns>
		private static Int64 _startpane(
			out Cartesian v1, out Cartesian v2, out Cartesian v3,
			double xin, double yin, double zin,
			char[] name) {
			int t;
			long  baseID;
			int baseindex;

			v1 = Cartesian.NaN;
			v2 = Cartesian.NaN;
			v3 = Cartesian.NaN;
			if ((xin > 0) && (yin >= 0)){
				baseindex = (zin >= 0) ? (int) Startindex._N3 : (int) Startindex._S0;

			} else if ((xin <= 0) && (yin > 0)){
				baseindex = (zin >= 0) ? (int) Startindex._N2 : (int) Startindex._S1;

			} else if ((xin < 0) && (yin <= 0)){
				baseindex = (zin >= 0) ? (int) Startindex._N1 : (int) Startindex._S2;
 
			} else if ((xin >= 0) && (yin < 0)){
				baseindex = (zin >= 0) ? (int) Startindex._N0 : (int) Startindex._S3;

			} else {
				baseindex = (zin >= 0) ? (int) Startindex._N3 : (int) Startindex._S0;
				// was: baseindex = -1;
			}
			if (baseindex < 0)
				return -1;

			baseID = HtmState.Instance.bases[baseindex].HID;
	
			t = HtmState.Instance.bases[baseindex].v1;

			v1 = new Cartesian(
				HtmState.Instance.anchor[t, 0],
				HtmState.Instance.anchor[t, 1],
				HtmState.Instance.anchor[t, 2], false);

			t = HtmState.Instance.bases[baseindex].v2;
			v2 = new Cartesian(
				HtmState.Instance.anchor[t, 0],
				HtmState.Instance.anchor[t, 1],
				HtmState.Instance.anchor[t, 2], false);

			t = HtmState.Instance.bases[baseindex].v3;
			v3 = new Cartesian(
				HtmState.Instance.anchor[t, 0],
				HtmState.Instance.anchor[t, 1],
				HtmState.Instance.anchor[t, 2], false);

			if (name != null){
				name[0] = HtmState.Instance.bases[baseindex].name[0];
				name[1] = HtmState.Instance.bases[baseindex].name[1];
				name[2] = '\0';
			}
			return HtmState.Instance.bases[baseindex].HID;
		}
		/// <summary>
		/// Test if p is inside triangle given by v1, v2, v3
		/// </summary>
		/// <param name="p">point to test</param>
		/// <param name="v1">first vertex of triangle</param>
		/// <param name="v2">second vertex of triangle</param>
		/// <param name="v3">third vertex of triangle</param>
		/// <returns></returns>
		public static bool Contains(Cartesian p, Cartesian v1, Cartesian v2, Cartesian v3) {
			// if (v1 X v2) . p < epsilon then false 
			// same for v2 X v3 and v3 X v1.
			// else return true..
			if (Cartesian.TripleProduct(v1, v2, p) < -Const.Epsilon)
				return false;
			if (Cartesian.TripleProduct(v2, v3, p) < -Const.Epsilon)
				return false;
			if (Cartesian.TripleProduct(v3, v1, p) < -Const.Epsilon)
				return false;
			return true;
		}
		/// <summary>
		/// Create a patch made of the three arcs that the edges of the trixel with the given Htm ID
		/// </summary>
		/// <param name="hid"></param>
		/// <returns></returns>
		public static Patch ToPatch(Int64 hid){
			Cartesian a, b, c;
			ToTriangle(hid, out a, out b, out c);
			Arc[] list = new Arc[3];
			list[0] = new Arc(a, b);
			list[1] = new Arc(b, c);
			list[2] = new Arc(c, a);
			Patch p = new Patch(list);
			return p;
		}
		/// <summary>
		/// 
		/// </summary>
		/// <param name="hid"></param>
		/// <param name="a"></param>
		/// <param name="b"></param>
		/// <param name="c"></param>
		public static void ToTriangle(Int64 hid, out Cartesian a, out Cartesian b, out Cartesian c)
		{
			// check for legal hid value

			char[] name = new char[eMaxNameSize];
			ToName(name, hid);
			NameToTriangle(name, out a, out b, out c);
			// Exceptions
		}
        /// <summary>
        /// Decide if the given 64-bit integer is a valid HtmID
        /// </summary>
        /// <param name="hid">64-bit input number</param>
        /// <returns>true if given number is a valid HtmID, false otherwise</returns>
		public static bool IsValid(Int64 hid) {
			Int32 i;

			// determine index of first set bit
			for (i = 0; i < eHIDBits; i += 2) {
				if (0 != ((hid << i) & IDHIGHBIT))
					break;
				if (0 != ((hid << i) & IDHIGHBIT2))  // invalid id
					return false;
			}
			if (hid < 8)			/* was: ==0, gyuri */
				return false;

			return true;
		}

		/// <summary>
		/// Convert a 64-bit HID to a text desciption of the trixel.
		/// </summary>
		/// <param name="name">An array of HtmTrixel.eMaxNameSize (a const in this class)
		/// character The array is null terminated.</param>
		/// <param name="hid">The HID</param>
		/// <returns>The size, or the length of the text</returns>
		public static int ToName(char[] name, long hid) {
			int size=0, i;
			int c;                                   // a spare character;
			long shifted_bit;
			long shifted_hid;
			// determine index of first set bit, top 2 always assumed 0
			// this is to eliminate the problem of mixing ulongs with longs
			//
			if(hid < 0)
				return -2; // higy bit set

			if(hid < 8)
				return -1; // -1 means Bad hid

			for(i=2; i< eHIDBits; i+=2){
				shifted_hid = hid << (i-2);
				shifted_bit = (shifted_hid) & IDHIGHBIT;
				if (shifted_bit != 0)
					break;
				if ((shifted_hid & IDHIGHBIT2) != 0)
					return -1;
			}
			

			size=(eHIDBits - i) >> 1;
			//
			// fill characters starting with the last one
			//
			for(i = 0; i < size-1; i++) {
				c =  '0' + (int) ((hid >> i*2) & (int) 3);
				name[size-i-1] = (char ) c;
			}
			//
			// put in first character
			//
			shifted_bit = (hid >> (size*2-2)) & 1;
			if(shifted_bit != 0) {
				name[0] = 'N';
			} else {
				name[0] = 'S';
			}
			name[size] = '\0'; // end string
			return size;
		}

		/// <summary>
		/// Convert the trixel from text to 64 bit HID.
		/// </summary>
		/// <param name="sname">The string with the text representation of 
		/// the trixel.</param>
		/// <returns>The 64 bit HID</returns>
		public static Int64 NameToHid(string sname){
			int length = sname.Length;
			char[] name = sname.ToCharArray();
			return NameToHid(name, length);
		}
		/// <summary>
		/// Convert the trixel from text to 64 bit HID
		/// The character array must be of size eMaxNameSize, and
		/// it is not necessary to have null termination. The
		/// number of siginificant characters in the trixel's name
		/// is given as a parameter.
		/// </summary>
		/// <param name="name">The character array with the text 
		/// representation of the trixel.</param>
		/// <param name="effectivelength">The number of significant characters
		/// in the array.</param>
		/// <returns>The 64 bit HID</returns>
		public static Int64 NameToHid(char[] name, int effectivelength){

			long result_hid=0;
			long shifted;
			int i;
			int siz = 0;
			long bits;
			int shift;
			siz = name.Length;
			if(name.Length < 2)
				return 0;	// 0 is an illegal HID

			if(name[0] != 'N' && name[0] != 'S')  // invalid name
				return result_hid;
			if(siz > eMaxNameSize)
				return 0;

			siz = effectivelength;
			for(i = siz-1; i > 0; i--) {// set bits starting from the end
				if(name[i] > '3' || name[i] < '0') {// invalid name
					return 0;
				}
				bits = ( (int)(name[i]-'0'));
				shift = 2*(siz - i -1);
				shifted = bits << shift;
				result_hid +=  shifted;
			}
			bits = 2;                     // set first pair of bits, first bit always set
			if(name[0]=='N') bits++;      // for north set second bit too
			shift = (2*siz - 2);
			shifted = bits << shift;
			result_hid += shifted;
			return result_hid;
		}

		/// <summary>
		/// Truncate the HID to fewer bit
		/// The HID of a trixel implicitly contains the level of the trixel.
		/// When you need a lower level trixel that contains the trixel
		/// of the given HID, this function gives it to you.
		/// If the level of
		/// the given htmid is less than or equal to the desired level,
		/// then there is no change.
		/// </summary>
		/// <param name="htmid"></param>
		/// <param name="level"></param>
		/// <returns>HID of lower level trixel</returns>
		public static Int64 Truncate(Int64 htmid, int level){
			Int64 result = htmid;
			int currentlevel = LevelOfHid(htmid);
			if (level < currentlevel){
				result = (htmid >> 2 * (currentlevel - level));
			}
			return result;
		}
		/// <summary>
		/// Extend given HID to a desired level.
		/// The opposite of truncate. However,
		/// because there are many descendents, the result is a range of
		/// consecutive HIDs. The low and hi values are returned in 
		/// the out variables.
		/// </summary>
		/// <param name="htmid">HID to extend</param>
		/// <param name="level">New level to which to extend</param>
		/// <returns>A pair of 64-bit HtmIDs</returns>
		public static Int64Pair Extend(Int64 htmid, int level){
			Int64 lo, hi;
			int currentlevel;
			int changelevel;
			Int64 dif;
			currentlevel = LevelOfHid(htmid);
			if (level > currentlevel){
				// amount to extend:
				changelevel = level - currentlevel;
				lo = htmid << (2 * changelevel);
				dif = 1L << (2 * changelevel); // Make sure 64 bit stuff works
				dif--;
				hi = lo + dif;
			} else {
				//truncate
				changelevel = currentlevel - level; // could be 0
				lo = htmid >> (2 * changelevel);
				hi = lo;
			}
			return new Int64Pair(lo, hi);
		}

		/// <summary>
        /// Returns the level number of an HID
		/// <param name="htmid">The HID</param>
		/// <returns>The level number</returns>
		public static int LevelOfHid(Int64 htmid){
			Int32 size=0, i;
			if (htmid < 0)
				return -1;
			for(i = 2; i < eHIDBits; i+=2) {
				if ( 0 != ((htmid << (i-2)) & IDHIGHBIT) )
					break;
			}
			size=(eHIDBits-i) / 2;
  			return size-2;
		}

		/// <summary>
		/// Get the internal angles of the trixel described by
		/// a given HID.
		/// </summary>
		/// <param name="hid">The 64-bit HID</param>
		/// <param name="t1">First internal angle</param>
		/// <param name="t2">Second internal angle</param>
		/// <param name="t3">Third internal angle</param>
		/// <returns>true if successful, false if HID is invalid</returns>
		public static bool GetAngles(Int64 hid, out double t1, out double t2, out double t3){
			// formula for area of spherical triangle
			// A = R^2[(t1 + t2 + t3 ) - Pi], but R=1.
			//
			Cartesian n0, n1, n2;
			Cartesian v0, v1, v2;
			char[] name = new char[eMaxNameSize];
			int level = ToName(name, hid);
			if (level < 0) {
				t1 = t2 = t3 = 0.0;
				return false;
			} else {
				NameToTriangle(name, out v0, out v1, out v2);
				//
				// Compute the plane normals for the three GC arcs
				//
				n2 = v1.Cross(v0, true);
				n0 = v2.Cross(v1, true);
				n1 = v0.Cross(v2, true);
				t1 = Math.PI - Math.Acos(n0.Dot(n1));
				t2 = Math.PI - Math.Acos(n1.Dot(n2));
				t3 = Math.PI - Math.Acos(n2.Dot(n0));
			}
			return true	;
		}

		/// <summary>
		/// Provides a rough idea of the length of a triangle side at a given
		/// level. This is based solely on the number of times a great circle
		/// segment is bisected. Actual lengths of trixels sides vary.
		/// </summary>
		/// <param name="level"></param>
		/// <returns>Nominal size of trixel edge in degrees.</returns>
		public static double FeatureSize(int level){
			double angle = Math.PI/2.0; // level 0 is 90 degrees
			if (level == 0)
				return angle;

			for(int i=0; i<level; i++){
				angle /= 2.0;
			}
			return angle;
		}
		/// <summary>
		/// Provides a rough idea of the level of a feature with the given size.
		/// </summary>
		/// <param name="size">degrees expressed as double</param>
		/// <returns>Nominal size of trixel edge in degrees.</returns>
		public static int EstimateLevelFromSize(double size) {
			int level = 0;
			double angle = Math.PI / 2.0; // level 0 is 90 degrees
			
			while(level < 25){
				if (size > angle)
					return level;

				angle /= 2.0;
				level++;
			}
			return level;
		}
	}
}
