#include "harmonic.h"
#include <vector>
using namespace std;

double harmonic::hermite(int k, double x) {
	double h0 = 1;
	if (k==0) return h0;
	double h1= 2*x;
	if (k==1) return h1;
	double h2=-2*h0+2*x*h1;
	if (k==2) return h2;
	for (int i=2; i<k; i++) {
		h0=h1;
		h1=h2;
		h2=-2*i*h0+2*x*h1;
	}
	return h2;
}


double harmonic::wav1d(int k, double x) {
	double fac=1; // 2^n * n!
	for (int i=1; i<(k+0.5); i++) {
		fac *= 2*i;
	}
	double norm = atan(1)*4; // (pi)
	double param = sqrt(mass*string);
	norm = sqrt (param/norm);
	norm = sqrt(norm/fac);
	x *= sqrt(param);
	norm *= hermite(k,x);
	norm *= exp (-x*x/2);
	return norm;
}
	
double harmonic::wav3d(int k, int l, double x) { 
	// 
	double param=sqrt(mass*string);
	double omega=sqrt(string/mass);
	double a0=1;
	double a1=2*param; // guess
	double res=a0;
	if (x<1e-20) x=1e-20;
	for (int i=0; i<k+0.5; i++) {
		if (i>1) {
			double a= (-2*mass*(i-0.5+l)*omega+param*(2*i+2*l-1))/i/(i+2*l-1)*a0;
			res+=exp((i+l)*log(x))*a;
			a0=a1;
			a1=a;
		}
		else if (i==1) res+=a1*exp((l+1)*log(x));
		else for (int m=1; m<l+0.5; m++) res*=x;
	}
	res*=wav1d(0,x);
	return res;
}

double harmonic::secdiv(int k, double x) {
	double E = (k+0.5) * sqrt(string/mass);
        double V=string/2*x*x;
	double secdiv= -2*mass*(E-V);
	return secdiv;
}