using System;

public class gaussian{
	public double A;
	public gaussian(double a){this.A=a;}
	}

public class nbody{

static Random rnd = new Random();

public static double rexp(double x){
	double r = rnd.NextDouble();
	return -r.log()/x;
	}

static double hbarc = 137.035999;
static double pi = 3.1415927;
static double mass = 1;

public static double L;

public static double[] matrix_elements(gaussian g1,gaussian g2){
	var B=g1.A+g2.A;
	var detB = B;
	var overlap = (pi/detB).pow(1.5);
	var kin = 6*(g1.A*L*g2.A/B)*overlap;
	var beta=B;
	var pot = -2*beta.sqrt()/pi.sqrt()*overlap;
	double[] result = {overlap,kin+pot};
	return result;
	}

public static gaussian[] basis;
public static matrix N,H;
public static GEVD z;

public static double energy(vector x){
	for(int i=0;i<x.size;i++) basis[i].A=x[i].abs();
	for(int i=0;i<x.size;i++)
	for(int j=0;j<x.size;j++)
		{
		double[] r = matrix_elements(basis[i],basis[j]);
		N[i,j]=r[0];
		H[i,j]=r[1];
		}
	z=new GEVD(H,N);
	double e=z.e[0];
	return e;
	}	

public static int Main(){
	amoeba.function master = energy;
	L = 0.5/mass;
	//var x = new vector("0.2 1.0");
	//var x = new vector("0.2");
	//var x = new vector("0.2 1.0 5");
	var x = new vector("0.2 0.4 2 11");
	//var x = new vector("0.1 0.2 2 11 15");
	int K=x.size;
	var ss = new vector(K);
	for(int i=0;i<K;i++)ss[i]=0.05;
	basis = new gaussian[K];
	for(int i=0;i<K;i++)basis[i]=new gaussian(1);
	N = new matrix(K,K);
	H = new matrix(K,K);
	amoeba A = new amoeba(master,x,ss);
	double acc=0.001;
	A.downhill(acc);
	x=A.p[A.lo];
	x.print("x=","{0,9:g4}");
	energy(x);
	vector c=z.V[0];
	var f = new Func<double,double>( r => {
		double sum=0;
		for(int i=0;i<c.size;i++)sum+=c[i]*Math.Exp(-x[i]*r*r);
		return sum;
		});
	var fo = new Func<double,double> (r => {return Math.Exp(-r)*f(1)/Math.Exp(-1.0);});
	for(double r=0;r<=7;r+=0.05)Console.Error.WriteLine("{0} {1} {2}",r,f(r),fo(r));
	return 0;
	}

}