#include <iostream>
#include <math.h>
#include <fstream>
#include <iomanip.h>

//function declarations


//General:

//converts barycentric coordinates to SG? array index
unsigned long int bary_to_index(unsigned long int bary[]);

//converts SG? array index to barycentric coordinates
unsigned long int index_to_bary(unsigned long int bary[], unsigned long int index);

//(1) finds indices of neighbors of a given level for a given index
//(2) finds smallest triangle on level of index containing index 
//option controlled by harmonicx
//true for (1)
//false for (2)
void neighbors(unsigned long int index, unsigned long int neighbors[],bool& harmonicx, int level);


//Spherical harmonic maps:

//computes euclidean energy of a harmonic map in a given coordinate
double e(int level, double SG[]);

//computes euclidean energy of harmonic map
double energy(int level);  

//computes euclidean energy with respect to neighbors
double penergy(int level, long int index); 

//finds euclidean distance squared between points index1 and index2
double distsq(unsigned long int index1, unsigned long int index2);
//computes spherical energy of harmonic map on sphere
double senergy(int level);

//takes the x and y values of index and returns the associated
//positive z-value.
double zcoord(long unsigned int index);

//finds harmonic map to sphere
void sharmonic(bool);

//finds spherical distance between point at index1 and index2
double sdist(unsigned long int index1, unsigned long int index2); 

//takes index as point as north pole (or south),
//assigns neighbors of index geodesic coordinates,
//and computes center of mass of the neighbors in geodesic 
//coordinates. Shows if index is at the center of mass
//of its neighbors.
void cmass(long int index,double &mx, double &my, long int level);

//computes barycentric coordinates of index in terms of the 
//points defining the smallest triangle on the level of index containing
//index. Should be close to (2/5, 2/5, 1/5) as in the euclidean
//harmonic extension algorithm
double baryw(long int index, double weights[], long int level );

//assigns solution of Lagrange Multipliers Problem to gasket
void assigns(long int index, double tempx, double tempy, double tempz, double lambda);

//writes data for harmonic map to a file in columnar text format, suitable
//for importing to Maple with the "read" command
void swriteout(char s[]);



//Hyperbolic (-x^2 - y^2 + z^2 = 1) harmonic maps:

//computes Minkowski energy with respect to neighbors
double mpenergy(int level, long int index);

//computes Minkowski energy hyperbolic harmonic map
double menergy(int level);

//finds Minkowski distance squared between points index1 and index2
double mdistsq(unsigned long int index1, unsigned long int index2);

//computes hyperbolic distance between index1 and index2
double hdist(long int index1, long int index2); 

//computes hyperbolic energy
double henergy(int level);

//computes hyperbolic harmonic map with given boundary image
//and possibly initial values
void hharmonic(bool);

//given x and y coordinates of a point on the (positive branch of the) hyperboloid, returns 
//z coordinate
double zhcoord(long unsigned int index);

//gives x coordinate of point in Poincare disc
double poincx(long int index);

//gives y coordinate of point in Poincare disc
double poincy(long int index);

double hcmass(long int index, double coord[]);

//writes data for harmonic map to a file in columnar text format, suitable
//for importing to Maple with the "read" command. Coordinates on the 
//Poincare disc given
void hwriteout(char s[]);


//Computational functions:

//returns absolute value of arg
double absolute(double arg);

//returns square of arg
double sq(double arg); 

//returns square of arg
int sq(int arg);

//returns arg^power. (Note that in C++ "^" denotes the bitwise XOR. 
//It does not denote exponentiation. Therefore, do not use expressions
//of the form x^y in your program unless you really want x XOR y). 
//Will not give reliable answers for large values of power. Though
//the for loop with j in this function allows j to run from 1 to 10,
//presumbably allowing power to be as high as 2^11-1=2047,
//don't good results with powers even close to this high. Overflow is liking to occur
//well before that. If you really need to use larger exponents, consider using the 
//function double pow provided in "math.h"
long int pwr(int arg, int power);

const int m=8;
//const int SIZE =(pwr(3,m+1)-3)/2 + 3;
const int SIZE = 88575;
const double sqrt32 = sqrt(3)/2;
const double tol = 5e-7;
const double inf = 1e-16;
const double pi = 3.14159265359;

const bool INITIAL_VALUES = true;
const bool NO_INITIAL_VALUES = false;

double SGx[SIZE];
double SGy[SIZE];
double SGz[SIZE];

//for the function pwr, powers of 2
const int powers[31] = 
{ 
	1,
	2,
	4,
	8,
	16,
	32,
	64,
	128,
	256,
	512,
	1024,
	2048,
	4096, 
	8192,
	16384,
	32768,
	65536,
	131072,
	262144,
	524288,
	1048576,
	2097152,
	4194304,
	8388608,
	16777216,
	33554432,
	67108864,
	134217728,
	268435456,
	536870912,
	1073741824
};	


void main()
{

//define gasket

		int i;
	
		for(i=1; i<=8; i++)
		  {
		    long int mag = pwr(10,i);
		    cout << mag << endl;
		    cout << pwr(10,i) << endl;

		    SGx[0]=mag*cos(0.);
		    SGy[0]=mag*sin(0.);
		    SGz[0]=zhcoord(0);

		    SGx[1]=mag*cos(2*pi/3);
		    SGy[1]=mag*sin(2*pi/3);
		    SGz[1]=zhcoord(1);


		    SGx[2]=mag*cos(4*pi/3);
		    SGy[2]=mag*cos(4*pi/3);
		    SGz[2]=zhcoord(2);

		/*
		SGx[3]=1;
		SGy[3]=0;
		SGz[3]=0;

		SGx[4]=1;
		SGy[4]=0;
		SGz[4]=0;

		SGx[5]=1;
		SGy[5]=0;
		SGz[5]=0;
		*/
	
 		hharmonic(NO_INITIAL_VALUES);
		
	/*	double c[2];
		long int index;
		for(index=(pwr(3,m)-3)/2 + 3;index<=(pwr(3,m+1)-3)/2+2;index++)
		  {
		    hcmass(index, c);
		    cout << index << " " << c[0] << " " << c[1] << endl;

		  }
*/

		char s[6] = "hyp";
		s[3] = 48+(i/10);
		s[4] = 48+(i-10*(i/10));

		hwriteout(s);
		}
}


//General Functions:

unsigned long int bary_to_index(unsigned long int bary[])
{	
	int twogcd=100; 
	int i,j;
	
//find greatest power of two (twogcd) that divides all three coordinates
//of the barycentric coordinates


	for(i=0; i<=2; i++)
	{
		j=0;

		if(bary[i]!=0)
		{
			while( int(bary[i]/powers[j]) == double(bary[i])/double(powers[j]))
			{
				j++;
			}
			if(j-1 < twogcd)
				twogcd=j-1;
		}
	}	

	unsigned long int lbary[3];
	unsigned long int levelindex=0;	

	unsigned long int p = powers[twogcd];
	
	lbary[0]=bary[0]/p;
	lbary[1]=bary[1]/p;
	lbary[2]=bary[2]/p;

//Takes care of boundary points
	if( twogcd==m )
	{
	for(i=0; i<=2; i++)
		if(lbary[i]==1)
		return i;	
	}	
 

// m-twogcd-1 is the highest power of two that will appear 
// as a bit in all of the scaled barycentric coordinates,
// excluding the case for the boundary points
	for(i=m-twogcd-1; i>=1; i--)
	{
		j=0;

		for(j=0; j<=2; j++)
		{
			if( (powers[i] & lbary[j]) != powers[i]  )
				;
			else
			{
				levelindex = 3 * levelindex + j;
				break;
			}
		}
		

	}
	switch( (1 & lbary[0]) + 2*(1 & lbary[1]) + 4*(1 & lbary[2]) )
	{
		case 3:
			levelindex = 3*levelindex;
			break;
		case 5:
			levelindex = 3*levelindex + 1;
			break;
		case 6:
			levelindex = 3*levelindex + 2;
			break;
		default:
			cout << "Shouldn't be here!" << endl;
			break;
	}
//Note on indexing.

//The value of twogcd indicates the first level n=m-twogcd on which
//the point represented by the barycentric coordinates appears.
//On each level n (excluding level 0), there are 3^n vertices.
//levelindex assigns a number between 0 and 3^n for a point on level n.
//levelindex is then added to (3 + 3^2 + 3^3 + . . . 3^(n-1)) + 3
//=(3^n - 3)/2 + 3 to find the index.  
//The boundary points are assigned indices of
//0, 1, and 2. This explains the "+ 3" in the above formula.

	return levelindex + (pwr(3,m-twogcd) - 3)/2 + 3;

}

unsigned long int index_to_bary(unsigned long int bary[], unsigned long
int index)
{
unsigned long int twogcd;
		
//find twogcd from index
//if loop goes all the way to twogcd=m,
//we know we have a boundary point

//boundary point case
	if(index <= 2)
	{
		bary[0]=0;
		bary[1]=0;
		bary[2]=0;
		bary[index] = powers[m];
		return 0;
	}


	for(twogcd=0; twogcd <= m-1; twogcd++)
		if( (pwr(3,m-twogcd) - 3)/2 + 3 <= index )
			break;

	index -= (pwr(3,m-twogcd) - 3)/2 + 3;
	
	switch(index % 3)
	{
		case 0:
			bary[0]=1;
			bary[1]=1;
			bary[2]=0;
			break;

		case 1:
			bary[0]=1;
			bary[1]=0;
			bary[2]=1;
			break;

		case 2:
			bary[0]=0;
			bary[1]=1;
			bary[2]=1;
			break;
	}
						
	index = (index - (index % 3))/3;
	
	int i;
	unsigned long int k;

//decode bits
	for(i=1; i<=m-twogcd-1; i++)
	{
		k=index % 3;
		bary[k] = bary[k] | powers[i];

		index=(index - k)/3;
	} 	

//multiply by to bring to proper level
	unsigned long int p=powers[twogcd];
	for(i=0; i<=2; i++)
		bary[i] *= p;

	return p;
		
}

void neighbors(unsigned long int index, unsigned long int neighbors[],
bool& harmonicx, int level)
{
//when harmonicx == true, the three indices of the vertices of the
//triangle used in the harmonic extension are returned. Those
//in neighbors[0] and neighbors[1] are the near neighbors on level
//k, where k is the level whereon the vertex was added. The 
//one in neighbors[2] is the far neighbor. 

//when harmonicx == false, the four neighbors of the vertex
//represented by index are returned in neighbors[0], . . . ,neighbors[3].
//Neighbors are ordered by the counterclockwise angle made with the
//positive x-axis. The variable level controls at which level the
//neighbors are computed. 


//function should never be called for boundary points
	if(index <= 2 || index >= (pwr(3,m+1)-3)/2 + 3 )
	{
		cout << "Invalid index in function neighbors" << endl;
		cout << "Called with index = " << index << endl;
	}

	unsigned long int bary[3];
	unsigned long int p=index_to_bary(bary, index);
	unsigned long int k=powers[m-level];

//for all points besides boundary points
//neighbor structure is determined by index mod 3

//(0,-1,1) ->    0 deg
//(1,-1,0) ->   60 deg
//(1,0,-1) ->  120 deg
//(0,1,-1) ->  180 deg 
//(-1,1,0) ->  240 deg 
//(-1,0,1) ->  300 deg

//case    -----neighbors---- | far extension
	
//case 0: 0  , 60 , 240, 300 | (-1, -1,  2)
//case 1: 120, 180, 240, 300 | (-1,  2, -1)
//case 2: 0  , 60 , 120, 180 | ( 2, -1, -1)

//triangle flag determines whether to compute the triangular completion 

	unsigned long int tempbary[3];

	switch(harmonicx)

	{
	case true:
	
		switch(index % 3)
		{
			case 0:
//60
				tempbary[0]=bary[0] + p;
				tempbary[1]=bary[1] - p;
				tempbary[2]=bary[2];
				neighbors[0]=bary_to_index(tempbary);
//240
				tempbary[0]=bary[0] - p;
				tempbary[1]=bary[1] + p;
				tempbary[2]=bary[2];
				neighbors[1]=bary_to_index(tempbary);
//far extension vertex
				tempbary[0]=bary[0] - p;
				tempbary[1]=bary[1] - p;
				tempbary[2]=bary[2] + 2*p;
				neighbors[2]=bary_to_index(tempbary);

				break;
	
			case 1:			
//120
				tempbary[0]=bary[0] + p;
				tempbary[1]=bary[1];
				tempbary[2]=bary[2] - p;
				neighbors[0]=bary_to_index(tempbary);
//300
				tempbary[0]=bary[0] - p;
				tempbary[1]=bary[1];
				tempbary[2]=bary[2] + p;
				neighbors[1]=bary_to_index(tempbary);
//far extension vertex
				tempbary[0]=bary[0] - p;
				tempbary[1]=bary[1] + 2*p;
				tempbary[2]=bary[2] - p;
				neighbors[2]=bary_to_index(tempbary);

				break;

			case 2:
//0
				tempbary[0]=bary[0];
				tempbary[1]=bary[1] - p;
				tempbary[2]=bary[2] + p;
				neighbors[0]=bary_to_index(tempbary);
//180
				tempbary[0]=bary[0];
				tempbary[1]=bary[1] + p;
				tempbary[2]=bary[2] - p;
				neighbors[1]=bary_to_index(tempbary);
//far extension vertex
				tempbary[0]=bary[0] + 2*p;
				tempbary[1]=bary[1] - p;
				tempbary[2]=bary[2] - p;
				neighbors[2]=bary_to_index(tempbary);

				break;
		}	
	
		break;


	case false:
		
		switch(index % 3)
		{
			case 0:
//0
				tempbary[0]=bary[0];
				tempbary[1]=bary[1] - k;
				tempbary[2]=bary[2] + k;
				neighbors[0]=bary_to_index(tempbary);
//60
				tempbary[0]=bary[0] + k;
				tempbary[1]=bary[1] - k;
        			tempbary[2]=bary[2];
				neighbors[1]=bary_to_index(tempbary);
//240
				tempbary[0]=bary[0] - k;
				tempbary[1]=bary[1] + k;
				tempbary[2]=bary[2];
				neighbors[2]=bary_to_index(tempbary);
//300
				tempbary[0]=bary[0] - k;
				tempbary[1]=bary[1];
				tempbary[2]=bary[2] + k;
				neighbors[3]=bary_to_index(tempbary);

				break;
	
			case 1:			
//120
				tempbary[0]=bary[0] + k;
				tempbary[1]=bary[1];
				tempbary[2]=bary[2] - k;
				neighbors[0]=bary_to_index(tempbary);
//180
				tempbary[0]=bary[0];
				tempbary[1]=bary[1] + k;
				tempbary[2]=bary[2] - k;
				neighbors[1]=bary_to_index(tempbary);
//240
				tempbary[0]=bary[0] - k;
				tempbary[1]=bary[1] + k;
				tempbary[2]=bary[2];
				neighbors[2]=bary_to_index(tempbary);
//300
				tempbary[0]=bary[0] - k;
				tempbary[1]=bary[1];
				tempbary[2]=bary[2] + k;
				neighbors[3]=bary_to_index(tempbary);

				break;

			case 2:

//0
				tempbary[0]=bary[0];
				tempbary[1]=bary[1] - k;
				tempbary[2]=bary[2] + k;
				neighbors[0]=bary_to_index(tempbary);
//60
				tempbary[0]=bary[0] + k;
				tempbary[1]=bary[1] - k;
				tempbary[2]=bary[2];
				neighbors[1]=bary_to_index(tempbary);
//120
				tempbary[0]=bary[0] + k;
				tempbary[1]=bary[1];
				tempbary[2]=bary[2] - k;
				neighbors[2]=bary_to_index(tempbary);
//180
				tempbary[0]=bary[0];
				tempbary[1]=bary[1] + k;
				tempbary[2]=bary[2] - k;
				neighbors[3]=bary_to_index(tempbary);

				break;
		}	

		break;

	}

	return;				
}

//Spherical Functions:

double e(int level, double SG[]) 
{
		double eng=0;
		long unsigned int index;
		long unsigned int n[4];
		bool harmonicx = false;

	for(index = (pwr(3,level)-3)/2 + 3;
	 index <=  (pwr(3,level+1)-3)/2 + 2; index++)

	{
		neighbors(index, n, harmonicx, level);			

 		switch( index % 3)
		{
			case 0:
				eng += 
				pow(double(5)/double(3),level)*(

				sq(SG[index]-SG[n[0]])+
				sq(SG[index]-SG[n[1]])+
				sq(SG[index]-SG[n[2]])
				);

				break;

			case 1:
				eng += 
				pow(double(5)/double(3),level)*(

				sq(SG[index]-SG[n[0]])+
				sq(SG[index]-SG[n[2]])+
				sq(SG[index]-SG[n[3]])
				);

				break;

			case 2:

				eng += 
				pow(double(5)/double(3),level)*(

				sq(SG[index]-SG[n[0]])+
				sq(SG[index]-SG[n[2]])+
				sq(SG[index]-SG[n[3]])

				);

				break;
		}

	}
return eng;

}

double energy(int level)
{

//compute total energy 

	double eng=0;
	unsigned long int index;
	bool harmonicx = false;
	unsigned long int n[4];

	if(level == 0)
	{
		return 	distsq(0, 1) +
			distsq(1, 2) +
			distsq(2, 0);
	}
	

	for(index =  (pwr(3,level)-3)/2 + 3; index <=  (pwr(3,level+1)-3)/2 + 2; index++)
	{
			
		neighbors(index, n, harmonicx, level);			


 		switch( index % 3)
		{
			case 0:
				eng += 
				pow(double(5)/double(3),level)*

				(
				distsq(index, n[0])+
				distsq(index, n[1])+
				distsq(index, n[2])
				);

				break;

			case 1:
				eng += 
				pow(double(5)/double(3),level)*

				(
				distsq(index, n[0])+
				distsq(index, n[2])+
				distsq(index, n[3])
				);

				break;

			case 2:

				eng += 
				pow(double(5)/double(3),level)*

				(
				distsq(index, n[0])+
				distsq(index, n[2])+
				distsq(index, n[3])
				);

				break;
		}
			
	}

return eng;

}


double penergy(int level, long int index)
{
//compute euclidean energy of point at index index with respect to its neighbors on 
//level level

	bool harmonicx = false;
	unsigned long int n[4];

	if(level == 0)
	{
		cout << "Improper call to function penergy: called for boundary point"
		 << endl;
	}
	

			
		neighbors(index, n, harmonicx, level);			


 		switch( index % 3)
		{
			case 0:
				return
				pow(double(5)/double(3),level)*
				(
				distsq(index, n[0])+
				distsq(index, n[1])+
				distsq(index, n[2])+
				distsq(index, n[3])
				);

				break;

			case 1:
				return 
				pow(double(5)/double(3),level)*
				(
				distsq(index, n[0])+
				distsq(index, n[1])+
				distsq(index, n[2])+
				distsq(index, n[3])
				);

				break;

			case 2:

				return
				pow(double(5)/double(3),level)*
				(
				distsq(index, n[0])+
				distsq(index, n[1])+
				distsq(index, n[2])+
				distsq(index, n[3])
				);

				break;
			
		}
}

double distsq( unsigned long int index1, unsigned long int index2)
{

	return 	sq( SGx[index1] - SGx[index2] ) + 
		sq( SGy[index1] - SGy[index2] ) + 
		sq( SGz[index1] - SGz[index2] ) ;
			


}		

double senergy(int level)
{ 


//compute total spherical energy 

	double eng=0;
	unsigned long int index;
	bool harmonicx = false;
	unsigned long int n[4];

	if(level == 0)
		return 	  sq(sdist(0, 1)) 
			+ sq(sdist(1, 2))
			+ sq(sdist(2, 0));

	for(index =  (pwr(3,level)-3)/2 + 3; index <=  (pwr(3,level+1)-3)/2 + 2; index++)
	{
			
		neighbors(index, n, harmonicx, level);			

 		switch( index % 3)
		{
			case 0:
				eng += 
				pow(double(5)/double(3),level)*
				(
					sq(sdist(index, n[0])) +
					sq(sdist(index, n[1])) +
					sq(sdist(index, n[2]))
				);

				break;

			case 1:
				eng += 
				pow(double(5)/double(3),level)*
				(
					sq(sdist(index, n[0])) +
					sq(sdist(index, n[2])) +
					sq(sdist(index, n[3]))

				);

				break;

			case 2:

				eng += 
				pow(double(5)/double(3),level)*
				(
					sq(sdist(index, n[0])) +
					sq(sdist(index, n[2])) +
					sq(sdist(index, n[3]))

				);

				break;
		}
			
	}

return eng;

}


double zcoord(long unsigned int index)
{
	return sqrt(1 - sq(SGx[index]) - sq(SGy[index]) );
}

void sharmonic(bool values)
{
//the below algorithm first makes the normalized 2/5, 2/5, 1/5 
//guess for the harmonic extension, then optimizes the energy of
//the guess

	unsigned long int n[4];
	long unsigned int index;
	double length;
	bool harmonicx;	
	int level;
	double oldsgx;
	double oldsgy;	
	int j;

	
	for(level=1; level<=m; level++)
	{

		if(values!=true || level!=1 )
		{
			for(index =  (pwr(3,level)-3)/2 + 3; index <=  (pwr(3,level+1)-3)/2 + 2; index++)

			{
				harmonicx=true;	
				//0 is a dummy argument		
				neighbors(index, n, harmonicx,0);
	
				SGx[index] = (2 * SGx[n[0]] + 2 * SGx[n[1]] + SGx[n[2]])/5;
				SGy[index] = (2 * SGy[n[0]] + 2 * SGy[n[1]] + SGy[n[2]])/5;
				SGz[index] = (2 * SGz[n[0]] + 2 * SGz[n[1]] + SGz[n[2]])/5;

				length =  sqrt( sq(SGx[index]) + sq(SGy[index]) + sq(SGz[index]) );
			
				SGx[index] = SGx[index] / length;
				SGy[index] = SGy[index] / length;
				SGz[index] = SGz[index] / length;
			}

		}
	
	harmonicx = false;	

//oldsg{x,y} are set to 2 so that the program doesn't mistakenly finish
//early

	oldsgx=2;
	oldsgy=2;

//while loop repeats until SG{x,y}[4] settles down 
	while	( 
		absolute(oldsgx-SGx[4]) + 
		absolute(oldsgy-SGy[4]) >=
		tol
		)
	{

//optimize energy one pass

		oldsgx = SGx[4];
		oldsgy = SGy[4];

		for(index = 3;
		 index <=  (pwr(3,level+1)-3)/2 + 2; index++)
		{		
			neighbors(index, n, harmonicx, level);

			double tempx = SGx[n[0]] + SGx[n[1]] + SGx[n[2]] + SGx[n[3]];
			double tempy = SGy[n[0]] + SGy[n[1]] + SGy[n[2]] + SGy[n[3]];
			double tempz = SGz[n[0]] + SGz[n[1]] + SGz[n[2]] + SGz[n[3]];

			double temp  = sqrt(sq(tempx) + sq(tempy) + sq(tempz));

//two possible parameters to Lagrange multipliers problem
			double lambda1 = 4 + temp;
			double lambda2 = 4 - temp;

			double energies[4];
	
//compute local energy (energy with neighbors) for both possibilities of lambda
//and both signs of SGz[index]

			assigns(index, tempx, tempy, tempz, lambda1);
			energies[0]= penergy(level,index);
			
			assigns(index, tempx, tempy, tempz, lambda2);
			energies[1] = penergy(level,index);
			
			
		double min=1000;
		int minindex;

//find minimum 
		for(j=0; j<=1; j++)
		{
			if(energies[j] <= min)
			{
				minindex = j;
				min = energies[j];
			}	


		}	

		
//set point on gasket to minimum
		int sign;
		switch(minindex)
		{
			case 0:
			       
			  assigns(index, tempx, tempy, tempz, lambda1);
			  
			  if(SGz[index] < 0)
			    sign = -1;
			  else 
			    sign=1;

			  SGz[index] = sign*zcoord(index);
			  break;

			case 1:
			  assigns(index, tempx, tempy, tempz, lambda2);
			  
			  if(SGz[index] < 0)
			    sign = -1;
			  else 
			    sign=1;

			  SGz[index] = sign*zcoord(index);
			  break;
				

		}		



		}//end for


	
			
	}//end while


	cout << 
	"Level: " << level <<  
	"\t Spherical energy: " << senergy(level) <<
	"\t Euclidean energy: " << energy(level) << endl;

	}



}

double sdist(unsigned long int index1, unsigned long int index2)
{

//The vector represented by index1 and index2 have magnitude 1. 
//Thus the their dot product is the cosine of the angle between them.
//We compute the spherical distance (arc length) between them.
		
	return acos	
			(
				SGx[index1]*SGx[index2] + 
				SGy[index1]*SGy[index2] + 
				SGz[index1]*SGz[index2] 
			);


}		

void cmass(long int index, double &mx, double &my, long int level)
{
	unsigned long int n[4];
	bool harmonicx = false;
	int i;
	mx=0;
	my=0;
	double phii;
	double z;
	double zi;
	
	z = zcoord(index);

	neighbors(index, n, harmonicx, level);

	for(i=0; i<=3; i++)
	{
		
		zi = zcoord(n[i]);
		phii = sdist(n[i], index);		

		mx += .25 * phii / sin(phii)	
			*(SGx[n[i]] - SGx[index] * (zi + cos(phii)) 
		/ (1 + z));

		my += .25 * phii / sin(phii)	
			*(SGy[n[i]] - SGy[index] * (zi + cos(phii)) 
		/ (1 + z));

	
	}
}


double baryw(long int index, double weights[], long int level )
{
	unsigned long int n[3];
	bool harmonicx = true;
	int i;
	double phii;
	double z;
	double zi;
	double ax[3];
	double ay[3];

	z = zcoord(index);

	neighbors(index, n, harmonicx, level);

	for(i=0; i<=2; i++)
	{
		
		zi = zcoord(n[i]);
		phii = sdist(n[i], index);		

		ax[i] = .25 * phii / sin(phii)	
			*(SGx[n[i]] - SGx[index] * (zi + cos(phii)) 
		/ (1 + z));

		ay[i] = .25 * phii / sin(phii)	
			*(SGy[n[i]] - SGy[index] * (zi + cos(phii)) 
		/ (1 + z));
	}
	double temp = -ax[0]*ay[1]+ax[0]*ay[2]+ay[0]*ax[1] -
			ay[0]*ax[2]-ax[1]*ay[2] + ax[2]*ay[1];

	weights[0] = (-ax[1]*ay[2] + ax[2]*ay[1]) / temp;
	weights[1] = (-ay[0]*ax[2] + ax[0]*ay[2]) / temp;
	weights[2] = (-ax[0]*ay[1] + ay[0]*ax[1]) / temp;	

	return sq((double)(weights[0] - .4)) 
		+ sq((double)(weights[1] - .4)) 
		+ sq((double)(weights[2] - .2));
}


void assigns(long int index, double tempx, double tempy, double tempz, double lambda)
{

  SGx[index] = tempx/(lambda - 4);
  SGy[index] = tempy/(lambda - 4);
  SGz[index] = tempz/(lambda - 4);
}

void swriteout(char file[])
{
	long unsigned int index=0;
	ofstream out(file);
	out << setprecision(8);
	for(index=0; index <= (pwr(3,m+1)-3)/2 + 2; index++)
	  out << SGx[index] << " " << SGy[index] << " "
		<< SGz[index] << endl;

}

//Hyperbolic functions:

double mpenergy(int level, long int index)
{
//compute euclidean energy of point at index index with respect to its neighbors on 
//level level

	bool harmonicx = false;
	unsigned long int n[4];

	if(level == 0)
	{
		cout << "Improper call to function mpenergy: called for boundary point"
		 << endl;
	}
	

			
		neighbors(index, n, harmonicx, level);			


 		switch( index % 3)
		{
			case 0:
				return
				pow(double(5)/double(3),level)*
				(
				mdistsq(index, n[0])+
				mdistsq(index, n[1])+
				mdistsq(index, n[2])+
				mdistsq(index, n[3])
				);

				break;

			case 1:
				return 
				pow(double(5)/double(3),level)*
				(
				mdistsq(index, n[0])+
				mdistsq(index, n[1])+
				mdistsq(index, n[2])+
				mdistsq(index, n[3])
				);

				break;

			case 2:

				return
				pow(double(5)/double(3),level)*
				(
				mdistsq(index, n[0])+
				mdistsq(index, n[1])+
				mdistsq(index, n[2])+
				mdistsq(index, n[3])
				);

				break;
			
		}
}

//computes Minkowski energy hyperbolic harmonic map
double menergy(int level)
{ 


//compute Minkowski energy on hyperboloid 

double eng=0;
	unsigned long int index;
	bool harmonicx = false;
	unsigned long int n[4];

	if(level == 0)
		return 	  mdistsq(0, 1) 
			+ mdistsq(1, 2)
			+ mdistsq(2, 0);

	for(index =  (pwr(3,level)-3)/2 + 3; index <=  (pwr(3,level+1)-3)/2 + 2; index++)
	{
			
		neighbors(index, n, harmonicx, level);			

 		switch( index % 3)
		{
			case 0:
				eng += 
				pow(double(5)/double(3),level)*
				(
				 mdistsq(index, n[0]) +
				 mdistsq(index, n[1]) +
		                 mdistsq(index, n[2])
				);

				break;

			case 1:
				eng += 
				pow(double(5)/double(3),level)*
				(
					mdistsq(index, n[0]) +
					mdistsq(index, n[2]) +
					mdistsq(index, n[3])

				);

				break;

			case 2:

				eng += 
				pow(double(5)/double(3),level)*
				(
					mdistsq(index, n[0]) +
					mdistsq(index, n[2]) +
					mdistsq(index, n[3])

				);

				break;
		}
			
	}

return eng;

}


//finds Minkowski distance squared between points index1 and index2
double mdistsq(unsigned long int index1, unsigned long int index2)
{

  return 
    + sq(SGx[index1]-SGx[index2]) 
    + sq(SGy[index1]-SGy[index2]) 
    - sq(SGz[index1]-SGz[index2]); 
}		


double hdist(long int index1, long int index2)
{
	return	
	  2*asinh(
		  sqrt(
		         sq(SGx[index1]-SGx[index2])
		       + sq(SGy[index1]-SGy[index2])
		       - sq(SGz[index1]-SGz[index2])
		      )/2
		);
		
} 

double henergy(int level)
{
//compute total hyperbolic energy 

	double eng=0;
	unsigned long int index;
	bool harmonicx = false;
	unsigned long int n[4];

	if(level == 0)
		return 	  sq(hdist(0, 1)) 
			+ sq(hdist(1, 2))
			+ sq(hdist(2, 0));

	for(index =  (pwr(3,level)-3)/2 + 3; index <=  (pwr(3,level+1)-3)/2 + 2; index++)
	{
			
		neighbors(index, n, harmonicx, level);			

 		switch( index % 3)
		{
			case 0:
				eng += 
				pow(double(5)/double(3),level)*
				(
					sq(hdist(index, n[0])) +
					sq(hdist(index, n[1])) +
					sq(hdist(index, n[2]))
				);

				break;

			case 1:
				eng += 
				pow(double(5)/double(3),level)*
				(
					sq(hdist(index, n[0])) +
					sq(hdist(index, n[2])) +
					sq(hdist(index, n[3]))

				);

				break;

			case 2:

				eng += 
				pow(double(5)/double(3),level)*
				(
					sq(hdist(index, n[0])) +
					sq(hdist(index, n[2])) +
					sq(hdist(index, n[3]))

				);

				break;
		}
			
	}

return eng;

}

void hharmonic(bool values)
{
//the below algorithm first makes the normalized 2/5, 2/5, 1/5 
//guess for the harmonic extension, then optimizes the energy of
//the guess


 	unsigned long int n[4];
	long unsigned int index;
	double length;
	bool harmonicx;	
	int level;
	double oldsgx;
	double oldsgy;	
	int j;

	
	for(level=1; level<=m; level++)
	{

		if(values!=true || level!=1 )
		{
			for(index =  (pwr(3,level)-3)/2 + 3; index <=  (pwr(3,level+1)-3)/2 + 2; index++)

			{
				harmonicx=true;	
				//0 is a dummy argument		
				neighbors(index, n, harmonicx,0);
	
				SGx[index] = (2 * SGx[n[0]] + 2 * SGx[n[1]] + SGx[n[2]])/5;
				SGy[index] = (2 * SGy[n[0]] + 2 * SGy[n[1]] + SGy[n[2]])/5;
				SGz[index] = (2 * SGz[n[0]] + 2 * SGz[n[1]] + SGz[n[2]])/5;

				length = sqrt( 
					      - sq(SGx[index])
					      - sq(SGy[index])
					      + sq(SGz[index])
					     );
			
				SGx[index] = SGx[index] / length;
				SGy[index] = SGy[index] / length;
				SGz[index] = SGz[index] / length;
			}

		}
	
	harmonicx = false;	

//oldsg{x,y} are set to 2 so that the program doesn't mistakenly finish
//early

	oldsgx=2;
	oldsgy=2;

//while loop repeats until SG{x,y}[4] settles down 
	while	( 
		absolute(oldsgx-SGx[4]) + 
		absolute(oldsgy-SGy[4]) >=
		tol
		)
	{

//optimize energy one pass

		oldsgx = SGx[4];
		oldsgy = SGy[4];

		for(index = 3;
		 index <=  (pwr(3,level+1)-3)/2 + 2; index++)
		{		
			neighbors(index, n, harmonicx, level);

			double tempx = SGx[n[0]] + SGx[n[1]] + SGx[n[2]] + SGx[n[3]];
			double tempy = SGy[n[0]] + SGy[n[1]] + SGy[n[2]] + SGy[n[3]];
			double tempz = SGz[n[0]] + SGz[n[1]] + SGz[n[2]] + SGz[n[3]];
	
			double temp  = sqrt(-sq(tempx) - sq(tempy) + sq(tempz) );

			double lambda[2];
			  
			lambda[0] = 4 + temp;
			lambda[1] = 4 - temp;
			
//test points on gasket		

			double energies[2]; 
			int minindex=0;

			int i;
			for(i=0; i<=1; i++)
			  {
			    //this code differs from the analogous code in sharmonic because we
			    //are restricting our view to only one branch of the hyperboloid 
			    SGx[index] = tempx/(lambda[i]-4);
			    SGy[index] = tempy/(lambda[i]-4);		
			    SGz[index] = zhcoord(index);			    
			    
			    energies[i]=mpenergy(level, index);
			    
			    if(energies[i]<energies[minindex])
			      minindex=i;
			  }

			SGx[index] = tempx/(lambda[minindex]-4);
			SGy[index] = tempy/(lambda[minindex]-4);		
			SGz[index] = zhcoord(index);

	       }//end for


	
			
	}//end while

	cout <<
	"Level: " << level <<  
	"\t Hyperbolic energy: " << henergy(level) <<
	"\t Minkowski: " << menergy(level) << endl;

	}



}

double zhcoord(long unsigned int index)
{

  return sqrt(sq(SGx[index]) + sq(SGy[index]) + 1);

}

//gives x coordinate of point in Poincare disc
double poincx(long int index)
{

  return SGx[index]/SGz[index];

}


//gives y coordinate of point in Poincare disc
double poincy(long int index)
{
  return SGy[index]/SGz[index];
}

double hcmass(long int index, double coord[])
{
  long unsigned int n[4];
  double a;
  double b;
  double c;
  double d;
  bool harmonicx = false;

  neighbors(index, n, harmonicx, m);

  int i;
  double rpart=0;
  double ipart=0;
  double real;
  double imag;
  double mag;

  for(i=0; i<=3; i++)
    {
      a = SGx[n[i]] - SGx[index];
      b = SGy[n[i]] - SGy[index];
      c = 1 - SGx[n[i]]*SGx[index] + SGy[n[i]]*SGy[index];
      d = SGx[n[i]]*SGy[index] - SGx[index]*SGy[n[i]];

      real = (a*c + b*d)/(sq(c) + sq(d));
      imag = (b*c + a*d)/(sq(c) + sq(d));

      mag = sqrt( sq(real) + sq(imag));
      
      coord[0] += real*atanh(mag)/mag;
      coord[1] += imag*atanh(mag)/mag;

    }
  
  coord[0] = coord[0] / 4;
  coord[1] = coord[1] / 4;
  
}

void hwriteout(char file[])
{
	long unsigned int index=0;
	ofstream out(file);
	out << setprecision(8);
	for(index=0; index <= (pwr(3,m+1)-3)/2 + 2; index++)
	  out << poincx(index) << " " << poincy(index) << endl;

}


//Computational Functions:

double absolute(double arg) 
{
if(arg <= 0)
	return -arg;
else 
	return arg;
}

double sq(double arg)
{
	return arg*arg;
}


int sq(int arg)
{
	return arg*arg;
}



long int pwr(int arg, int power)
{
	
	int curargpower = arg;
	int result = 1;
	int j;
	
	for(j=0; j<=10; j++)
	{
	
		if((power & powers[j])==powers[j])
				result=result*curargpower;

		curargpower=sq(curargpower);	
			

	}

return result;	
	
}