#ifndef _INCLUDED_L3_RAD_H
#define _INCLUDED_L3_RAD_H

// *Radial part of wavefunction - antiproton scattering* Copyright (C) Krzysztof Bosak, 2000-01-16...2000-03-31
// All rights reserved.
// kbosak@box43.pl
// http://www.kbosak.prv.pl

#include "l3_mat5.h"
#include "l3_comp.h"
#include "l3_tool.h"
#include <math.h>

double PhaseShift(
	double (*potential)(double),
	double rmax,
	int steps,
	int l,
	double beta,
	double k
	)
{
	assert(potential!=NULL);
	assert(rmax>=0);
	assert(steps>=10);
	assert(k>0);
	assert(l>=0);
	assert(beta!=0);
	const int stepsm1=steps-1;
	const double stepsize=rmax/steps;
	const double off_diag_elem=1.0/(stepsize*stepsize);
	int i;
	basearray<double> diagonal(steps);
	double r=stepsize;
	const double diag_elem=-2.0/(stepsize*stepsize);
	for(i=0; i<stepsm1; i++, r+=stepsize)
	{
		diagonal[i]=-potential(r)+(-l*(l+1)/(r*r)+k*k+diag_elem);
	}
	diagonal[stepsm1]=1;
	basearray<double> under_diagonal(stepsm1);
	for(i=0; i<stepsm1; i++)
	{
		under_diagonal[i]=off_diag_elem/diagonal[i];
		diagonal[i+1]-=under_diagonal[i]*off_diag_elem;
	}
	basearray<double> y(steps);
	y[0]=0;
	for(i=1; i<stepsm1; i++)
	{
		y[i]=-under_diagonal[i-1]*y[i-1];
	}
	const int stepsm2=steps-2;
	y[stepsm1]=-under_diagonal[stepsm2]*y[stepsm2]+beta*off_diag_elem;
	for(i=stepsm2; i>=0; i--)
	{
		diagonal[i]=(y[i]-diagonal[i+1]*off_diag_elem)/diagonal[i];
	}
	const int phase_look_index=steps-8;
	r=phase_look_index*stepsize;
	//const double deriv1=(diagonal[phase_look_index+1]-diagonal[phase_look_index-1])/(2*stepsize);
	const double deriv1=(
		(diagonal[phase_look_index-2]-diagonal[phase_look_index+2])/12+
		(diagonal[phase_look_index+1]-diagonal[phase_look_index-1])*2/3
		)/stepsize;
	return atan(static_cast<double>(diagonal[phase_look_index]*k/deriv1))+l*_PIdiv2-k*r;
}

const ComplexNumber<double> ComplexPhaseShift(
	const ComplexNumber<double> (*potential)(double), 
	double rmax,
	int steps, 
	int l, 
	double beta, 
	double k
	)
{
	assert(potential!=NULL);
	assert(rmax>=0);
	assert(steps>=10);
	assert(k>0);
	assert(l>=0);
	assert(beta!=0);
	const int stepsm1=steps-1;
	const double stepsize=rmax/steps;
	const double off_diag_elem=1.0/(stepsize*stepsize);
	int i;
	basearray<ComplexNumber<double> > diagonal(steps);
	double r=stepsize;
	const double diag_elem=-2.0/(stepsize*stepsize);
	for(i=0; i<stepsm1; i++, r+=stepsize)
	{
		diagonal[i]=-potential(r)+(-l*(l+1)/(r*r)+k*k+diag_elem);
	}
	diagonal[stepsm1]=1;
	basearray<ComplexNumber<double> > under_diagonal(stepsm1);
	for(i=0; i<stepsm1; i++)
	{
		under_diagonal[i]=ComplexNumber<double>(off_diag_elem)/diagonal[i];
		diagonal[i+1]-=under_diagonal[i]*off_diag_elem;
	}
	basearray<ComplexNumber<double> > y(steps);
	y[0]=0;
	for(i=1; i<stepsm1; i++)
	{
		y[i]=-under_diagonal[i-1]*y[i-1];
	}
	const int stepsm2=steps-2;
	y[stepsm1]=-under_diagonal[stepsm2]*y[stepsm2]+beta*off_diag_elem;
	for(i=stepsm2; i>=0; i--)
	{
		diagonal[i]=(y[i]-diagonal[i+1]*off_diag_elem)/diagonal[i];
	}
	const int phase_look_index=steps-8;
	r=phase_look_index*stepsize;
	//const ComplexNumber<double> deriv1=(diagonal[phase_look_index+1]-diagonal[phase_look_index-1])/(2*stepsize);
	const ComplexNumber<double> deriv1=(
		(diagonal[phase_look_index-2]-diagonal[phase_look_index+2])/12+
		(diagonal[phase_look_index+1]-diagonal[phase_look_index-1])*2/3
		)/stepsize;
	const ComplexNumber<double> phaseshift=ComplexArcTan(diagonal[phase_look_index]*k/deriv1)
		+l*static_cast<double>(3.1415926535897932384626433832795/2)-k*r;
	return phaseshift;
}
double ScatteringLength(
	double (*potential)(double), 
	double rmax, 
	int steps, 
	double k1, 
	double k2
	)
{
	const double a1=PhaseShift(potential, rmax, steps, 0, 1, k1)/k1;
	const double a2=PhaseShift(potential, rmax, steps, 0, 1, k2)/k2;
	const double k21=k2/k1;
	return a1-(a2-a1)/(k21*k21-1);
}
const ComplexNumber<double> ComplexScatteringLength(
	const ComplexNumber<double> (*potential)(double), 
	double rmax, 
	int steps, 
	double k1, 
	double k2
	)
{
	const ComplexNumber<double> a1=ComplexPhaseShift(potential, rmax, steps, 0, 1, k1)/k1;
	const ComplexNumber<double> a2=ComplexPhaseShift(potential, rmax, steps, 0, 1, k2)/k2;
	const ComplexNumber<double> k21=k2/k1;
	return a1-(a2-a1)/(k21*k21-1);
}
double EnergyEigenstate(
	double (*func)(double), 
	double x1, 
	double x2, 
	double precision
	)
{
	if(x2<=x1)
	{
		cerr<<"x2 must be >x1\n";
		exit(1);
	}
	const double f1=func(x1);
	const double f2=func(x2);
	if(f2<=f1)
	{
		cerr<<"func(x2) must be >func(x1)\n";
		cerr<<"f1="<<f1<<endl;
		cerr<<"f2="<<f2<<endl;
		exit(1);
	}
	const double level=(f1+f2)*0.5;
	//cout.precision(10);
	do
	{
		double xtest=(x1+x2)*0.5;
		cout<<x1<<'\t'<<xtest<<'\t'<<x2<<endl;
		double ytest=func(xtest);
		if(ytest<level)
			x1=xtest;
		else
			x2=xtest;
	}while(x2-x1>=precision);
	return (x1+x2)*0.5;
}
// EXPERIMENTAL 4TH ORDER METHOD //
const ComplexNumber<double> ComplexPhaseShift5(
	const ComplexNumber<double> (*potential)(double),
	double rmax, 
	int steps, 
	int l, 
	double beta, 
	double k
	)
{
	assert(potential!=NULL);
	assert(rmax>=0);
	assert(steps>=10);
	assert(k>0);
	assert(l>=0);
	assert(beta!=0);

	const double stepsize=rmax/steps;

	assert(beta!=0);
	Mat5Diag<ComplexNumber<double> > mat(steps);

	const int stepsm1=steps-1;
	const int stepsm2=steps-2;
	int i;

	mat.Diag(0)=-2.0/(stepsize*stepsize);
	mat.UDiag(0)=1.0/(stepsize*stepsize);
	mat.U2Diag(0)=0;

	mat.LDiag(0)=+16.0/(12*stepsize*stepsize);
	mat.Diag(1)=-30.0/(12*stepsize*stepsize);
	mat.UDiag(1)=+16.0/(12*stepsize*stepsize);
	mat.U2Diag(1)=-1.0/(12*stepsize*stepsize);

	for(i=2; i<steps-2; i++)
	{
		mat.U2Diag(i)=-1.0/(12*stepsize*stepsize);
		mat.UDiag(i)=+16.0/(12*stepsize*stepsize);
		mat.Diag(i)=-30.0/(12*stepsize*stepsize);
		mat.LDiag(i-1)=+16.0/(12*stepsize*stepsize);
		mat.L2Diag(i-2)=-1.0/(12*stepsize*stepsize);
	}

	mat.L2Diag(stepsm2-2)=0;
	mat.LDiag(stepsm2-1)=1.0/(stepsize*stepsize);
	mat.Diag(stepsm2)=-2.0/(stepsize*stepsize);
	mat.UDiag(stepsm2)=1.0/(stepsize*stepsize);

	mat.L2Diag(stepsm1-2)=0;
	mat.LDiag(stepsm1-1)=0;
	mat.Diag(stepsm1)=1;

	double r=stepsize;
	for(i=0; i<stepsm1; i++, r+=stepsize)
	{
		mat.Diag(i)+=-potential(r)-l*(l+1)/(r*r)+k*k;
	}

	mat.Decompose();
	basearray<ComplexNumber<double> > right_side(steps);
	for(i=0; i<stepsm1; i++)
	{
		right_side[i]=0;
	}
	right_side[stepsm1]=beta;

	basearray<ComplexNumber<double> > result(right_side.size());
	mat.Solve(result, right_side);

	const int phase_look_index=steps-16;
	r=phase_look_index*stepsize;
	const ComplexNumber<double> deriv1=(
		(result[phase_look_index-2]-result[phase_look_index+2])/12+
		(result[phase_look_index+1]-result[phase_look_index-1])*2/3
		)/stepsize;
	const ComplexNumber<double> phaseshift=ComplexArcTan(result[phase_look_index]*k/deriv1)
		+l*static_cast<double>(3.1415926535897932384626433832795l/2)-k*r;
	return phaseshift;
}
const ComplexNumber<double> ComplexScatteringLength5(
	const ComplexNumber<double> (*potential)(double),
	double rmax, 
	int steps, 
	double k1, 
	double k2
	)
{
	const ComplexNumber<double> a1=ComplexPhaseShift5(potential, rmax, steps, 0, 1, k1)/k1;
	const ComplexNumber<double> a2=ComplexPhaseShift5(potential, rmax, steps, 0, 1, k2)/k2;
	const ComplexNumber<double> k21=k2/k1;
	return a1-(a2-a1)/(k21*k21-1);
}

//////////////////////////////////////////////////////////////////////////////
inline double NullPotential(double)
{
	return 0;
}

inline const ComplexNumber<double> ComplexNullPotential(double)
{
	return 0;
}

//////////////////////////////////////////////////////////////////////////////
double R_SQUARE=1;	//real positive
ComplexNumber<double> V_SQUARE(1, 0);	//real positive

inline double SquareWellPotentialREffective()
{
	return R_SQUARE;
}
inline double SquareWellPotential(double r)
{
	if(r<=R_SQUARE)
		return V_SQUARE.re;
	return 0;
}
inline const ComplexNumber<double> ComplexSquareWellPotential(double r)
{
	if(r<=R_SQUARE)
		return V_SQUARE;
	return 0;
}
inline double SquareWellPhaseShift_L0(double K, double Vmax, double Rmax)
{
	const double kappa=sqrt(static_cast<double>(K*K+Vmax));
	return atan(static_cast<double>(K*tan(static_cast<double>(kappa*Rmax))/kappa))-K*Rmax;
}
inline const ComplexNumber<double> ComplexSquareWellPhaseShift_L0(
	double K, 
	const ComplexNumber<double> &Vmax, 
	double Rmax
	)
{
	const ComplexNumber<double> kappa=ComplexSqrt(Vmax+K*K);
	return ComplexArcTan((ComplexTan(kappa*Rmax)*K)/kappa) - K*Rmax;
}
inline double SquareWellScatteringLength_L0(double Vmax, double Rmax)
{
	assert(Rmax>0);
	assert(Vmax>0);
	const double kappa=sqrt(static_cast<double>(Vmax));
	return tan(static_cast<double>(kappa*Rmax))/kappa-Rmax;
}
inline const ComplexNumber<double> ComplexSquareWellScatteringLength_L0(
	ComplexNumber<double> Vmax, 
	double Rmax
	)
{
	assert(Rmax>0);
	assert(Vmax!=0);
	const ComplexNumber<double> kappa=ComplexSqrt(Vmax);
	return ComplexTan(kappa*Rmax)/kappa-Rmax;
}
double SquareWellEnergyEigenState(double v)
{// To search for discontinuities corresponding to eigenstates
	R_SQUARE=1;
	V_SQUARE=v;
	const double rmax=SquareWellPotentialREffective();
	const int steps=200000;//max 2'000'000 for 128MB RAM
	const double k1=1e-2;
	const double k2=1e-7;
	return ComplexScatteringLength(ComplexSquareWellPotential, rmax, steps, k1, k2).re;
}
//////////////////////////////////////////////////////////////////////////////
ComplexNumber<double> G_EXPON(1, 0);	// real + imag, positive or negative
double MJU_EXPON=1;	// real positive

inline double ExponPotentialREffective()
{
	return log(static_cast<double>(G_EXPON.Module()/1e-5));// /MJU_EXPON
}
inline double ExponPotential(double r)
{
	return G_EXPON.re*exp(static_cast<double>(-MJU_EXPON*r));
}
inline const ComplexNumber<double> ComplexExponPotential(double r)
{
	return G_EXPON*static_cast<double>(-MJU_EXPON*r);
}
double ExponEnergyEigenState(double v)
{// To search for discontinuities corresponding to eigenstates
	MJU_EXPON=1;
	G_EXPON=v;
	const double rmax=ExponPotentialREffective();
	const int steps=1000000;//max 2'000'000 for 128MB RAM
	const double k1=1e-3;
	const double k2=1e-4;
	return ComplexScatteringLength(ComplexExponPotential, rmax, steps, k1, k2).re;
}
//////////////////////////////////////////////////////////////////////////////
ComplexNumber<double> G_WOODS(1, 0);	// real + imag, positive or negative
double A_WOODS=1;		// real positive
double R_WOODS=1;		// real positive

inline double WoodsSaxonPotentialREffective()
{
	return R_WOODS*10;//oooverestimated
}
inline double WoodsSaxonPotential(double r)
{
	return G_WOODS.re/(1+exp(static_cast<double>((r-R_WOODS)/A_WOODS)));
}
inline const ComplexNumber<double> ComplexWoodsSaxonPotential(double r)
{
	return G_WOODS/static_cast<double>(1+exp(static_cast<double>((r-R_WOODS)/A_WOODS)));
}

#endif //_INCLUDED_L3_RAD_H