//--------------------------------------------------------------------
// 1. in order to run DWBABoundState() needs to include Potential.js and Constants.js
// 2. Potentials:
//	  2.1. Do not use Yukawa potential! 
//         The integration method in the level-search-procedure is not adapted for such potential.
//--------------------------------------------------------------------
function DWBABoundState()
{
	this.Const = new Constants();
	this.NName = new NucleusName();
	// nuclear masses and charges
	this.M1 = 1;
	this.Z1 = 0;
	this.M2 =12; 
	this.Z2 = 6; 
	// quantum number for current state
	this.J = 3; // total momentum, doubled
	this.L = 1; // orbital momentum
	this.S = 1; // particle spin, doubled
	this.Ms= 1; // spin projection, doubled
	this.Nr= 0; // main number
	//----------------------------------------------------------------
	// array for all quantum numbers
	this.MaxNumberOfStates = 300;
	this.QState = []; 
	// String
	this.ListOfStates = [];
	this.SymbolState  = [];	
	for(var i = 0; i < this.MaxNumberOfStates; i++)
	{   
		this.QState[i] = [];
		this.SymbolState[i] = "";
		for (var j = 0; j < 5; j++) this.QState[i][j] = 0;
	}
	//----------------------------------------------------------------
	this.State;      // index of current state
	this.NState = -1;// number of state in this potential
	this.Ebound;     // bound energy for current state
	this.kbound;     // real(sqrt(2 mu E/h^2))
	this.mu = 1*this.M1*this.M2/(1*this.M1+1*this.M2); // reduce mass=M1*M2/A
	// parameters of the mesh
	this.NR = 401; // number of mesh points
	this.Rmax = 20.0; // relative distance
	this.dR   = 0.05; // mesh step

	this.type = [];  //new Array(this.nMax);
   	this.psi  = [];  // real part of wave function
	this.EState = [];// energies of the i-th state
	this.EPSILON = 1.e-6; // relative accuracy
		
	this.rA;this.rB;
	this.dLnYin,this.dLnYout;		
	//PotentialDWBABoundState
	this.pot;
	//boolean
	this.LevelsNeedToBeSearched = true;
	this.CheckEnergy = false;
	this.Debugging = true;
	
	this.K_p;// 2*mu_p/hbar^2; mu_p~(A-a)a/A
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//  total spin may be S = 0,1,2 it means 0,1/2h,1h
this.setBoundState = function(z1,m1,z2,m2,nr,j,l,s,potential)
{
	//```````````````````````````````````````````````````````````````````````````	
	this.M1=m1; this.Z1=z1; this.M2=m2; this.Z2=z2; 
	this.mu =1.*this.M1*this.M2/(1*this.M1+1*this.M2);
	this.K_p=2.*this.mu/this.Const.h2_mn;
	//```````````````````````````````````````````````````````````````````````````
	this.LevelsNeedToBeSearched = true;
	this.CheckEnergy = false;
	// Total momentum J, Spin S and it's projection Ms are doubled
	this.Nr = nr; this.J = j; this.L = l; this.S = s;
	this.Ms = this.J - 2*this.L;
	//```````````````````````````````````````````````````````````````````````````
	//                   massivy dlya moih podprogramm
	this.NState = 0;
	this.State  = -1; // before searching levels we can not define it
	//```````````````````````````````````````````````````````````````````````````
	this.pot = potential;
	this.pot.setReaction(this.Z1,this.M1,this.Z2,this.M2);
	this.pot.addCoeff = this.K_p;
	this.pot.addConst = 0.;
	this.pot.setL(this.L,false);
	this.pot.setSpin(this.S,false);
	this.pot.setSpinProjection(this.Ms,false);
	this.pot.calculatePotential(true);
	//```````````````````````````````````````````````````````````````````````````
	this.ListOfStates = null;
	//```````````````````````````````````````````````````````````````````````````
	this.Rmax = this.pot.Rmax; 
	this.dR   = this.pot.dR; 
	this.NR   = this.pot.NR;
	
	for(var i = 0; i < this.NR; i++) this.psi[i] = 0.0;
}
//-------------------------------------------------------------------------------
// в потенциале может быть несколько слагаемых (до 10 штук)
// варьируем глубину 1-го слагаемого, считая, что оно главное
this.adjustDepth = function (e,nr,_2j,l,_2s,eps)
{
	//------------------------------------------------------------------------
	// на всякий случай заведем индекc для варьируемого слагаемого
	iPot = 0;
	//------------------------------------------------------------------------
	// ustanovim nuzhnyi spin
   	this.pot.setSpin(_2s,false);
	//------------------------------------------------------------------------
    iState = -1;
   	nr_max = -1;
   	Scale = 1.25;
     //------------------------------------------------------------------------
     // 1. podgonyaem glubinu chtoby bylo sostoyanie s dannym L i chislom uzlov ne menee nr
     while(nr_max<nr)
     {
        if(nr_max>=0)// masshtabiruem glubinu V0 i VD
        {
			this.pot.setAdjustingDepth(iPot,Scale*this.pot.V0[iPot]);
        }
        if(this.Debugging) console.log("Debugging: V0="+this.pot.V0[iPot]+" Rv="+this.pot.rV[iPot]+" aV="+this.pot.aV[iPot]+" L="+l);
		this.levelSearch(-l);//this.saveLevels();
        if(this.NState>0) nr_max = this.QState[this.NState-1][3]; else nr_max = 0;
     }
     if(this.Debugging) console.log("Debugging: rough adjusting stage is finished.\n");     
     //------------------------------------------------------------------------
     // proverka na nalichie pravel'nogo J
     Stop = false;
     while(true)
     {
     	for(i=0;i<this.NState;i++)
     		if(this.QState[i][0] == _2j && this.QState[i][3] == nr)
     		{
	     		iState = i;
	     		Stop = true;
	     		break;
	     	}
        if(Stop) break;
         
		this.pot.setAdjustingDepth(iPot,Scale*this.pot.V0[iPot]);
		this.levelSearch(-l);
     }
     //------------------------------------------------------------------------
     if(iState < 0) return false;
     //------------------------------------------------------------------------
     // podgonka Ebound
     var
     	Iteration = 0,
     	e0 = this.EState[iState] - e, 
		em,ep,
        a,b,c,D,
        x1,x2,
        Vnew,Vm,Vp,
        dV = e0*0.25;

     while(Math.abs(e0)>eps)
     {
	    Iteration++;
	    if(Iteration>1000)
	    {
			if(this.Debugging) console.log("Debugging: Failed to search bound state!\n"+"\tE(current) = "+e0);
		    return false;
	    }
		//---------------------------------------------------------------------
     	Vm = this.pot.V0[iPot];
     	Vp = 1.*Vm+1.*dV;
        //---------------------------------------------------------------------
		this.pot.setAdjustingDepth(iPot,Vp);
	    this.levelSearch(-l);
		ep = this.EState[iState] - e;
		//---------------------------------------------------------------------
		// linear approximation  e = a + b*V0
		//vp = pot.V0; vm = pot.V0 - dV;
		//a = (vp*e0 - vm*ep)/dV;
		//b = (ep-e0)/dV;
		//dv = -(vp*e0 - vm*ep)/(ep-e0);
		Vnew = Vp - dV*ep/(ep-e0);
		this.pot.setAdjustingDepth(iPot,Vnew);
        //---------------------------------------------------------------------
      	this.levelSearch(-l);
        e0 = this.EState[iState] - e;

        if(this.Debugging) console.log("Debugging: current V0="+Vnew+" Ebound = "+this.EState[iState]+" E0="+e+" delta1="+(Vnew-Vp)+" delta2="+(Vnew-Vm));
        //---------------------------------------------------------------------
        // est' sluchai ne dostizheniya tochnosti! Togda proga zacyklivaetsya
        if((Math.abs(Vnew - Vp) < this.Const.small_add) || (Math.abs(Vnew - Vm) < this.Const.small_add)) 
        {
	        if(this.Debugging) console.log("Debugging: required accuracy "+eps+" was not achieved. Current one is "+e0);
	        break;
        }
		dV = 0.5*Math.abs(Vnew - Vp);
		if(dV>0.5*Math.abs(this.EState[iState])) dV = 0.5*Math.abs(this.EState[iState]);
		if(e0>0) dV = -dV;
	}
    //------------------------------------------------------------------------
    this.Nr = nr; this.J = _2j; this.L = l; this.S = _2s;
	if(this.Debugging) console.log("Debugging: adjusting done: nr="+this.Nr+"| 2j="+this.J+" l="+this.L+" 2s="+this.S);
    //------------------------------------------------------------------------
    return true;
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//   Na osnove poluchennogo spektra schitaem w.f.
//   n < 0 - schitaem sostoyanie s tekushimi kvantovymi chislami
this.calculateWF = function (n)
{
	//-------------------------------------------------------------------------
	if(this.LevelsNeedToBeSearched) this.levelSearch(1);
	//-------------------------------------------------------------------------
	if(this.NState<=0) return false;
	//-------------------------------------------------------------------------
	if(n>=this.NState) return false;
	if(n>=0 && n!=this.State) this.setCurrentState(n);
	//-------------------------------------------------------------------------
	this.pot.setL(this.L,false);
	//-------------------------------------------------------------------------
	var
		iMatch,
		_rA,_rB,
		rMatch,
		rMatchMin=1.0;
	this.Ebound = this.EState[this.State]; 
	this.kbound = Math.sqrt(2.*this.mu*Math.abs(this.Ebound)/this.Const.h2_mn); 
	//-------------------------------------------------------------------------
	// Mathching point
	if(!this.turningPoints(this.Ebound)) return false;
	if(this.rB < rMatchMin) this.rB = rMatchMin;
	iMatch = Math.floor(this.rB/this.pot.dR); 
	rMatch = iMatch*this.pot.dR;
	if(rMatch < this.rA || rMatch > this.rB) return false;
	if(iMatch >= this.pot.NR) return false;
	//-------------------------------------------------------------------------
	if(this.getLogDerivatives(this.Ebound,iMatch,true))
	{
        if(Math.abs((this.dLnYin-this.dLnYout)/(1.*Math.abs(this.dLnYin) + 1.*Math.abs(this.dLnYout))) > this.EPSILON)
        return false;
	}
	else return false;
	//-------------------------------------------------------------------------
	return true;
}
this.doDebugging = function (doIt) 
{
	this.Debugging = doIt;
	this.pot.doDebugging(doIt);
}
this.doEnergyCheck = function (DoIt)
{
	this.CheckEnergy = DoIt;
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
this.getAproximateWF = function (r)
{
	var
		ir = Math.floor(r/this.dR+0.5),
		p,re,im;
	//---------------------------------------------------------------------------
	if(ir<0) return 0;
	if(ir>=this.NR-1) return this.getAsymptotics(r);
	//---------------------------------------------------------------------------
	p=(r-ir*this.dR)/this.dR;
	if(ir==0 || ir==1 || ir==this.NR-2 || ir==this.NR-3)
	{
		re = (1.-p)*this.psi[ir] + p*this.psi[1*ir+1];
	}
	else
	{
		re = 0.5*p*(p-1)*this.psi[ir-1] + (1-p*p)*this.psi[ir] + 0.5*p*(1*p+1)*this.psi[1*ir+1];
	}
	//---------------------------------------------------------------------------
	return re;
}
this.getAsymptotics = function (r)
{
	var
		rmax=(this.NR-1)*this.dR,
	   	eKr = 1.0/Math.exp(-rmax*this.kbound),
	   	Cre = this.psi[NR-1]*eKr,
	   	re;
	eKr = Math.exp(-r*this.kbound);
	re = eKr*Cre;
	return re;
}
this.getLogDerivatives = function (Eb,iMatch,doNorm)
{
	this.pot.calculatePotential(true); // vrode by esli chtoto izmenili => sootvetstvuushii flag izmenen
	
	if(iMatch>=this.pot.NR) { this.dLnYin=1.; this.dLnYout=-1.; return false;}
	
	var
		k = Math.sqrt(2.*this.mu*Math.abs(Eb)/this.Const.h2_mn),
		k2 = -k*k, 
		kappa,
		h2_12 = this.pot.dR*this.pot.dR/12.,
		C=1./(120.*this.pot.dR), 
		D, sum, dY, F0,F1,F2, i;
	//-------------------------------------------------------------------------
	// Integration from the origin (Numerov)
	this.psi[0]=0.0; this.psi[1]=1.e-100;
	for(i=2;i<=iMatch+2;i++)
	{  
		F0 = 1. +     h2_12*(1*k2-this.pot.V[i-2]);
        F1 = 2. - 10.*h2_12*(1*k2-this.pot.V[i-1]);
        F2 = 1. +     h2_12*(1*k2-this.pot.V[i]);

        this.psi[i] = (1*this.psi[i-1]*F1 - 1*this.psi[i-2]*F0)/F2;
	}
	this.dLnYin = 
		C*(  -4.*Math.log(Math.abs(this.psi[1*iMatch-3]))+30.*Math.log(Math.abs(this.psi[1*iMatch-2]))
		   -120.*Math.log(Math.abs(this.psi[1*iMatch-1]))+40.*Math.log(Math.abs(this.psi[  iMatch]))
		    +60.*Math.log(Math.abs(this.psi[1*iMatch+1]))- 6.*Math.log(Math.abs(this.psi[1*iMatch+2])));
	D = this.psi[1*iMatch-3]; // dlya posleduushei normirovki
	//-------------------------------------------------------------------------
	// Integration from infinity (Numerov)
	kappa = Math.sqrt(-1.*k2+1.*this.pot.V[this.pot.NR-1]);
	this.psi[this.pot.NR-1]= 1.e-100;
	this.psi[this.pot.NR-2]= this.psi[this.pot.NR-1] * Math.exp(kappa*this.pot.dR);
	for(i=this.pot.NR-3;i>=iMatch-3;i--)
	{  
		F0 = 1. +     h2_12*(1*k2-this.pot.V[i]);
		F1 = 2. - 10.*h2_12*(1*k2-this.pot.V[1*i+1]);
        F2 = 1. +     h2_12*(1*k2-this.pot.V[1*i+2]);

        this.psi[i] =  (1.*this.psi[1*i+1]*F1 - 1.*this.psi[1*i+2]*F2)/F0;
	}
	this.dLnYout = 
		C*(  -4.*Math.log(Math.abs(this.psi[1*iMatch-3]))+30.*Math.log(Math.abs(this.psi[1*iMatch-2]))
		   -120.*Math.log(Math.abs(this.psi[1*iMatch-1]))+40.*Math.log(Math.abs(this.psi[  iMatch]))
		    +60.*Math.log(Math.abs(this.psi[1*iMatch+1]))- 6.*Math.log(Math.abs(this.psi[1*iMatch+2])));
	//-------------------------------------------------------------------------
	//  normiruem w.f.
	if(doNorm)
	{
		D = this.psi[iMatch-3]/D; for(i=0;i<iMatch-3;i++) this.psi[i] = D*this.psi[i];
		sum=0.; F1=0.; 
		for(i=1;i<this.pot.NR;i++) 
		{
			F2  = this.psi[i]*this.psi[i]; 
			sum = 1.*sum+1.*F1+1.*F2; 
			F1  = F2;
		}
		C = 1./Math.sqrt(0.5*sum*this.pot.dR);
		for(i=0;i<this.pot.NR;i++) this.psi[i] = C*this.psi[i];
	}
	//-------------------------------------------------------------------------
	return true;
}
this.getState = function (n)
{
	iL = '';
	switch (this.QState[n][1])
	{   case 0: iL='s'; break; case 1: iL='p'; break; case 2: iL='d'; break;
		case 3: iL='f'; break; case 4: iL='g'; break; case 5: iL='h'; break;
		case 6: iL='i'; break; case 7: iL='j'; break; case 8: iL='k'; break;
	}
	if(this.QState[n][0]%2==0) j = ""+this.QState[n][0]/2; else j = this.QState[n][0]+"/2";
	buf = ""+this.QState[n][3]+""+iL+""+j;
	return buf;
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//                   Searching the energy levels
this.levelSearch = function (_Lmax)
{
	var //double
		Eb,Emax,Emin,dE,
		F1,F2,F, E1,E2,
		rMatch, rMatchMin=1.0,         
		rP,
		acc;
	var //boolean
		Stop;
	var //int
		s,l,j,nr, i,k,
		NStateOld, NodeNumber,
		iMatch, Niter,
		Lmin,Lmax,
		WF_Sign;
	//````````````````````````````````````````````````````````````````````````
	// otcenim primernoe znachenie maximal'nogo L_bound
	// L(L+1) ~ L^2 = 2mu/h^2 *V0*Rv^2/2
	if(_Lmax>0)
	{
		rP = this.pot.getMaximalRadius();
		F1 = Math.abs(this.pot.getPocket(true,false));
		Lmax = Math.floor(1.0 + 1.*rP*Math.sqrt(this.mu*F1/this.Const.h2_mn));
		if(Lmax<_Lmax) Lmax=_Lmax;
		Lmin = 0;
	}
	else {Lmin = Math.abs(_Lmax); Lmax = Lmin;}
	if(this.Debugging) console.log("Debugging: Lmax="+Lmax+"| Lmin="+Lmin);
	//````````````````````````````````````````````````````````````````````````
	//                  Nachinaem poisk sostoyanii
	this.NState = 0; NStateOld=-1;
	//````````````````````````````````````````````````````````````````````````
	for(l=Lmin;l<=Lmax;l++)
	{
		if(NStateOld == this.NState) break;
		NStateOld = this.NState;
		//=====================================================================
		this.pot.setL(l,true);
		//=====================================================================
		// spin-orbit coupling
		if(this.S>0)
		{
			if(this.Debugging) console.log("Debugging: Spin-Orbit coupling!");
			for(s=-1*this.S;s<=1*this.S;s=1*s+2)
			{
				if(l==0 && s>-this.S) continue;
				
				Stop = false;
				this.pot.setSpinProjection(s,true);
				if(this.Debugging) console.log("Debugging: L="+this.pot.L+" S="+this.S+" MS="+s+" LS="+this.pot.LS);
				
				Emin = this.pot.getPocket(true,true); 
				rP   = this.pot.getPocket(false,true);
				
				if(this.Debugging) console.log("Debugging: Emin="+Emin+" rP="+rP);

				if(Emin>0 || rP >= this.Rmax-this.dR) continue;// no pocket or no bound states
				Emax = 0;
				
				NodeNumber=1;
				dE = (1.*Emax-1.*Emin)/50.;
				Eb = 1.*Emin - 0.5*dE;
				WF_Sign = 0;
				while(Eb<0)
				{
					F1 = 1; F2 = 1; Niter=0; dE = (1.*Emax-1.*Eb)/50.0;
					iMatch = Math.floor(this.pot.getMaximalRadius()/this.pot.dR); 
					//`````````````````````````````````````````````````````````````
					// rough search
					while(F1*F2 > 0)
					{
						//```````````````````````````````````````````````````````````
						// ulovka chtoby ne propustit' urovni u poverhnosti
						if(Eb>-10 && this.NState==0 && dE > 0.5) dE = 0.5;
						if(Eb>-2 && dE > 0.05) dE = 0.05;
						//```````````````````````````````````````````````````````````
						Eb = 1.*Eb + 1.*dE; 
						if(Eb > -0.25*dE) {Stop=true; break;}
						if(!this.turningPoints(Eb)) continue;//{Stop=true; break;}
						// knowing rB value we choose a matching radius
						if(this.rB<rMatchMin) this.rB = rMatchMin;
						iMatch = Math.floor(this.rB/this.pot.dR); 
						rMatch = iMatch*this.pot.dR;
						if((rMatch < this.rA || rMatch > this.rB) || (iMatch >= this.pot.NR))
						{
							Stop=true;
							break;
						}
						//---------------------------------------------------------------
						F1 = F2;
						if(!this.getLogDerivatives(Eb,iMatch,false))
						{
							Stop=true; break;
						}
						F2 = this.dLnYin - this.dLnYout;
						if(Niter==0)
						{
							Niter=1;
							F1=F2;
						}
						//-----------------------------------------------------------
						if(WF_Sign!=0)
						{
							F1 = F2;
							if(WF_Sign*this.psi[iMatch-4]>0) continue;
							else WF_Sign=0;
						}
						//-----------------------------------------------------------
						/*/ srazu posle naidennogo sobstvennogo znacheniya w.f.(in)
						// nachinaet "zagibat'sya" vse kruche i kruche i v kakoito
						// moment {d(ln Psi)/dr = 1/Psi dPsi/dr} rezko menyaet znak
						// (naklon staryi, a Psi stala otritcatel'noi). Chto by eto
						// ne interpretirovat' kak novoe s.z., sledim za ee znakom
						if(dLnPsiIn * dLnPsiIn_OLD < 0 && In_Out_OLD >= 0) F1=-F1;
						dLnPsiIn_OLD = dLnPsiIn;
						In_Out_OLD = dLnPsiIn*dLnPsiOut;*/
					}
					if(Stop) break;
					//``````````````````````````````````````````````````````````````
					// accurate search
					E1=1.*Eb-1.*dE; E2=Eb; acc = 1;
					while(acc > this.EPSILON)
					{
						Eb = 0.5*(1.*E1+1.*E2);
						if(!this.getLogDerivatives(Eb,iMatch,false))
						{
							Stop=true; break;
						}
						acc = Math.abs((this.dLnYin - this.dLnYout)/(1.*Math.abs(this.dLnYin) + 1.*Math.abs(this.dLnYout)));
						F = 1.*this.dLnYin - 1.*this.dLnYout;
						if(F*F1 > 0.) {F1=F; E1=Eb;} else E2=Eb;
						Niter=1*Niter+1; if(Niter>1000) {Stop=true; break;}
					}
					if(Stop) break;
					//``````````````````````````````````````````````````````````````
					if(this.Debugging) console.log("Debugging: Ebound="+Eb+"|Nodes="+NodeNumber+"|NState="+this.NState);
					WF_Sign = (this.psi[iMatch-4]<0)?(-1):1;
					this.EState[this.NState] = Eb; 
					this.QState[this.NState][3] = NodeNumber;
					this.QState[this.NState][0] = Math.abs(2*l+1*s); 
					this.QState[this.NState][1] = l; 
					this.QState[this.NState][2] = this.S;
					this.NState++; NodeNumber++;
					if(this.NState == this.MaxNumberOfStates) break;
					//``````````````````````````````````````````````````````````````
				}// end of while(Eb)
			}
		}
		//=====================================================================
		// no spin-orbit coupling
		else
		{
			if(this.Debugging) console.log("Debugging: no Spin-Orbit coupling!");
			Stop = false;
			this.pot.setSpinProjection(0,true);	//this.pot.calculatePotential(true);
			if(this.Debugging) console.log("Debugging: L="+this.pot.L+" S="+this.S+" LS="+this.pot.LS);
			
			Emin = this.pot.getPocket(true,true);
			rP = this.pot.getPocket(false,true);
			
			if(this.Debugging) console.log("Debugging: Emin="+Emin+" rP="+rP);
			
			if(Emin>=0 || rP >= 1.*this.Rmax-1.*this.dR) continue;// no pocket or no bound states
			Emax = 0;
			
			dE = (1.*Emax-1.*Emin)/50.; if(dE>1) dE=1;
			Eb = 1.*Emin - 0.5*dE;
			NodeNumber = 1;
			WF_Sign = 0;
			while(Eb<0)
			{
				F1 = 1; F2 = 1; Niter=0; dE = (1.*Emax-1.*Eb)/50.0; if(dE>1) dE=1;
				iMatch = Math.floor(this.pot.getMaximalRadius()/this.pot.dR); 
				//`````````````````````````````````````````````````````````````````
				// rough search
				while(F1*F2 > 0)
				{
					//```````````````````````````````````````````````````````````````
					// ulovka chtoby ne propustit' urovni u poverhnosti
					if(Eb>-10 && this.NState==0 && dE > 0.5) dE = 0.5;
					if(Eb>-2 && dE > 0.05) dE = 0.05;
					//```````````````````````````````````````````````````````````````
					Eb = 1.*Eb + 1.*dE; 
					if(Eb > -0.25*dE) {Stop=true; break;} //if(Eb>0) {Stop=true; break;}					
					if(!this.turningPoints(Eb)) continue;//{Stop=true; break;}
					// on the basis of rB value we choose a matching radius.
					// It could not be to large, because in this case
					// we need a very small energy step in order to do not miss
					// that we are looking for
					if(this.rB<rMatchMin) this.rB = rMatchMin;
					iMatch = Math.floor(this.rB/this.pot.dR);
					rMatch = iMatch*this.pot.dR;
					if((rMatch < this.rA || rMatch > this.rB) || (iMatch >= this.pot.NR))
					{
						Stop = true;
						break;
					}
					//---------------------------------------------------------------
					F1 = F2;
					if(!this.getLogDerivatives(Eb,iMatch,false))
					{
						Stop=true; break;
					}
					F2 = 1.*this.dLnYin - 1.*this.dLnYout;
					if(Niter==0)
					{
						Niter=1;
						F1 = F2;
					}
					//---------------------------------------------------------------
					if(WF_Sign!=0)
					{
						F1 = F2;
						if(WF_Sign*this.psi[iMatch-4]>0) continue;
						else WF_Sign=0;
					}
					//---------------------------------------------------------------
					/*/ srazu posle naidennogo sobstvennogo znacheniya w.f.(in)
					// nachinaet "zagibat'sya" vse kruche i kruche i v kakoito
					// moment {d(ln Psi)/dr = 1/Psi dPsi/dr} rezko menyaet znak
					// (naklon staryi, a Psi stala otritcatel'noi). Chto by eto
					// ne interpretirovat' kak novoe s.z., sledim za ee znakom
					if(dLnPsiIn * dLnPsiIn_OLD < 0 && In_Out_OLD >= 0) F1=-F1;
					dLnPsiIn_OLD = dLnPsiIn;
					In_Out_OLD = dLnPsiIn*dLnPsiOut;*/
				}
				if(Stop) break;
				//`````````````````````````````````````````````````````````````````
				// accurate search
				E1=1.*Eb-1.*dE; E2=Eb; acc = 1;
				while(acc > this.EPSILON)
				{
					Eb = 0.5*(1.*E1+1.*E2);
					if(!this.getLogDerivatives(Eb,iMatch,false))
					{
						Stop=true; break;
					}
					acc = Math.abs((this.dLnYin - this.dLnYout)/(Math.abs(this.dLnYin) + Math.abs(this.dLnYout)));
					F = 1.*this.dLnYin - 1.*this.dLnYout;
					if(F*F1 > 0.) {F1=F; E1=Eb;} else E2=Eb;
					Niter=1*Niter+1; if(Niter>1000) {Stop=true; break;}
				}
				if(Stop) break;
				//`````````````````````````````````````````````````````````````````
				if(this.Debugging) console.log("Debugging: Ebound="+Eb+"|Nodes="+NodeNumber+"|NState="+this.NState);
				WF_Sign = (this.psi[iMatch-4]<0)?(-1):1;
				this.EState[this.NState] = Eb; 
				this.QState[this.NState][3]=NodeNumber;
				this.QState[this.NState][0]=2*l; 
				this.QState[this.NState][1]=l; 
				this.QState[this.NState][2]=this.S;
				this.NState++; NodeNumber++;
				if(this.NState==this.MaxNumberOfStates) break;
				//`````````````````````````````````````````````````````````````````
			}// end of while(Eb)
		}// end of else - no spin coupling
		if(l==Lmax-1 && this.NState!=NStateOld) Lmax=1*Lmax+1;
		if(this.NState==this.MaxNumberOfStates)
		{
			console.log("Number of levels in the potential is more then "+this.MaxNumberOfStates+
				".\nOnly "+this.MaxNumberOfStates+" ones are were recognized.\n");
			break;
		}
	}// end of l-loop
	//````````````````````````````````````````````````````````````````````````
	//  Sort of the states
	for (i=0;i<this.NState-1;i++)
		for (k=1*i+1;k<this.NState;k++)
			if (this.EState[k]<this.EState[i])
			{
				F = this.EState[k];    this.EState[k]   =this.EState[i];    this.EState[i]=F;
				j = this.QState[k][0]; this.QState[k][0]=this.QState[i][0]; this.QState[i][0]=j;
				l = this.QState[k][1]; this.QState[k][1]=this.QState[i][1]; this.QState[i][1]=l;
				s = this.QState[k][2]; this.QState[k][2]=this.QState[i][2]; this.QState[i][2]=s;
				nr= this.QState[k][3]; this.QState[k][3]=this.QState[i][3]; this.QState[i][3]=nr;
			}
	OrbMom = '';
	for(k=0;k<this.NState;k++)
	{
		switch (this.QState[k][1])
		{   
			case 0: OrbMom='S'; break; case 1: OrbMom='P'; break;
			case 2: OrbMom='D'; break; case 3: OrbMom='F'; break;
			case 4: OrbMom='G'; break; case 5: OrbMom='H'; break;
			case 6: OrbMom='I'; break; case 7: OrbMom='J'; break;
			case 8: OrbMom='K'; break; default: OrbMom='?';
		}
		buf = ""+(this.QState[k][3])+""+OrbMom;
		if(this.QState[k][0]%2 != 0) buf +=""+this.QState[k][0]+"/2";
		else                         buf +=""+(1*this.QState[k][0]/2);
		this.SymbolState[k] = buf;
	}
	//````````````````````````````````````````````````````````````````````````
	this.LevelsNeedToBeSearched=false;
	//````````````````````````````````````````````````````````````````````````
	//   vostanovim tekushie znacheniya orb. momemnta i proektcii spina
	this.setCurrentState(-1);
	this.pot.setL(this.L,false);
	this.pot.setSpinProjection(this.Ms,true);
	//````````````````````````````````````````````````````````````````````````
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
this.saveLevels = function ()
{
	buf  = "Level scheme for the system:";
	buf += ""+this.NName.getNucleusName(this.M1,this.Z1)+" + ";
	buf += ""+this.NName.getNucleusName(this.M2,this.Z2)+" (";
	buf += ""+((this.NState>0)?this.NState:0)+" levels have been found) \n\n";

	if(this.NState<=0) buf += "No levels!";
	else
	{
		Ntot = [];//this.ListOfStates = [];
		//---------------------------------------------------------------------
		// Opredelim chislo sostoyanii na dannom urovne
		// N = Sum(2J_i+1)
		for(k=0;k<this.NState;k++)
		{
			if(this.QState[k][0]%2 != 0) Ntot[k] = (1*this.QState[k][0]+1);
			else					     Ntot[k] = (2*this.QState[k][0]+1);
			//if(k>0) Ntot[k] += Ntot[k-1];
		}
		//---------------------------------------------------------------------
		buf += "   E(MeV)  J(h)  L(h) S(h)  symbol    Nr   Ntot\n";
		for(k=this.NState-1;k>=0;k--)
		{
			buf += this.EState[k].toFixed(5)+"  ";
			if(this.QState[k][0]%2 != 0) 
				buf += ""+this.QState[k][0]+"/2  "+this.QState[k][1]+"  "+this.QState[k][2]+"/2";
			else 
				buf += ""+(this.QState[k][0]/2)+"  "+this.QState[k][1]+"  "+(this.QState[k][2]/2);
				   
			buf += "  "+this.SymbolState[k]+"  "+this.QState[k][3]+"  "+Ntot[k]+"\n";
		}
	}
	if(this.Debugging) console.log(new String(buf));
	return buf;
}

this.saveWaveFunction = function ()
{
	buf = "wave function for the state: " + this.getState(this.State)+" E = "+
		  this.EState[this.State].toFixed(4)+" MeV\n";
	for (i=0;i<this.NR;i++) 
		buf += Math.floor(this.pot.dR*i*1000.)/1000.+"  "+ this.psi[i].toString()+"\n";
	return buf;
}
this.setCurrentEnergy = function (n)
{
	this.Ebound = this.EState[n];
	this.kbound=Math.sqrt(2.*this.mu*Math.abs(this.Ebound)/this.Const.h2_mn);
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
this.setCurrentState = function (n)
{
	if(n>=0 && n<this.NState)
	{
		this.setCurrentEnergy(n);
		this.J = this.QState[n][0];
		this.L = this.QState[n][1];
		this.S = this.QState[n][2];
		this.Nr= this.QState[n][3];
		this.Ms = this.J - 2*this.L;
		this.State = n;
	}
	else if(this.NState>0)
	{
		// ishem sostoyanie blizkoe po energii v pervuu ochered', a potom uzhe ostal'noe
		if(this.CheckEnergy)
		{
			this.CheckEnergy = false;
			for(i=0;i<this.NState;i++)
				if(Math.abs(this.Ebound-this.EState[i])<0.05)
					if(this.J == this.QState[i][0] && this.L == this.QState[i][1] && this.S == this.QState[i][2])
					{
						this.setCurrentEnergy(i);
						this.Ms = this.J - 2*this.L; 
						this.State = i; 
						return;
					}
		}
		// poisk polnogo sootvetstviya
		for(i=0;i<this.NState;i++)
			if(this.J == this.QState[i][0] && this.L == this.QState[i][1] && this.S == this.QState[i][2] && this.Nr == this.QState[i][3])
			{
				this.setCurrentEnergy(i);
				this.Ms = this.J - 2*this.L; 
				this.State = i; 
				return;
			}
		// poisk chastichnogo sootvetstviya, aktual'no posle supersimetrizatcii
		for(i=0;i<this.NState;i++)
			if(this.J == this.QState[i][0] && this.L == this.QState[i][1] && this.S == this.QState[i][2])
			{
				this.setCurrentEnergy(i);
				this.Ms = this.J - 2*this.L; 
				this.State = i; 
				return;
			}
		this.State = 0;
		this.setCurrentEnergy(0);
		this.J = this.QState[0][0];
		this.L = this.QState[0][1];
		this.S = this.QState[0][2];
		this.Nr= this.QState[0][3];
		this.Ms = this.J - 2*this.L;
	}
}

this.setM1 = function (val)
{ 
	this.M1 = val; 
	this.mu = 1.*this.M1*this.M2/(1*this.M1+1*this.M2);
	this.K_p = 2.*this.mu/_Constants.h2_mn; 
	this.pot.setReaction(this.Z1,this.M1,this.Z2,this.M2);	
	this.pot.addCoeff = K_p;
	this.LevelsNeedToBeSearched=true;
}
this.setM2 = function (val)
{ 
	this.M2 = val; 
	this.mu = 1.*this.M1*this.M2/(1*this.M1+1*this.M2);
	this.K_p = 2.*this.mu/_Constants.h2_mn; 
	this.pot.setReaction(this.Z1,this.M1,this.Z2,this.M2);
	this.pot.addCoeff = this.K_p;     
	this.LevelsNeedToBeSearched=true;
}
this.setRmax = function (rmax,dr,recalculate)
{
	this.pot.setDefaultParameters(rmax,dr);
	this.Rmax = this.pot.Rmax; 
	this.dR = this.pot.dR; 
	this.NR = this.pot.NR;
	
	this.psi= [this.NR];
}
this.setZ1 = function (val)
{ 
	this.Z1=val; 
	this.pot.setReaction(this.Z1,this.M1,this.Z2,this.M2);
	this.LevelsNeedToBeSearched=true;
}
this.setZ2 = function (val)
{ 
	this.Z2=val; 
	this.pot.setReaction(this.Z1,this.M1,this.Z2,this.M2);
	this.LevelsNeedToBeSearched=true;
}
this.turningPoints = function (Eb)
{
	eps = 1.e-5;
	//-------------------------------------------------------------------------
	xmax = 5.*this.pot.getMaximalRadius(); 
	if(xmax>20) xmax = 20; 
	if(xmax< 5) xmax = 5; 
	dx=this.pot.getMaximalRadius()/50.;  
	if(dx>0.25) dx=0.25;
	//-------------------------------------------------------------------------
	this.rA=-1.; this.rB=-1.;
	x = 0.; v = this.pot.getPotential(x,true);
	//-------------------------------------------------------------------------
	// left turning point
	if(Eb > v) this.rA=0.;
	else
	{
		vb = 1.*Eb - 1.*v; xa = xb = x; va=vb;
		// poisk intervala
		while(vb*va > 0.)
		{ 
			xa=xb; va=vb;
			xb= 1.*xb + 1.*dx; 
			vb = 1.*Eb - 1.*this.pot.getPotential(xb,true);
			if(xb>xmax) return false;
		}
		// metod deleniya otrezkov popolam
		v=1.; Nstep=0;
		while(Math.abs(v)>eps)
		{ 
			x=0.5*(1.*xa+1.*xb); 
			v = 1.*Eb - 1.*this.pot.getPotential(x,true);
			if(v*va > 0.) {va=v; xa=x;}
			else xb=x;
			Nstep++; if(Nstep > 100) break;//kakie-to problemy !?!
		}
		this.rA=x;
	}
	//-------------------------------------------------------------------------
	// right turning point
	x  = 1.*x + 0.5*dx; 
	vb = 1.*Eb - 1.*this.pot.getPotential(x,true); 
	xa = xb = x; va=vb;
	// poisk intervala
	while(vb*va > 0.)
	{
		xa=xb; va=vb;
		xb=1.*xb+1.*dx; vb = 1.*Eb - 1.*this.pot.getPotential(xb,true);
		if(xb>xmax)
		{
			// est' tolko odna tochka povorta i ona vneshnya
			this.rB  = this.rA; 
			this.rA  = 1.*this.rA + 1.*this.pot.ComFun.small_add;
			return true;
		}
	}
	// metod deleniya otrezkov popolam
	v=1.; Nstep=0;
	while(Math.abs(v)>eps)
	{
		x=0.5*(1.*xa+1.*xb); v = 1.*Eb - 1.*this.pot.getPotential(x,true);
		if(v*va > 0.) {va=v; xa=x;}
		else xb=x;
		Nstep++; if(Nstep > 100) break;//kakie-to problemy !?!
	}
	this.rB=x;
	return true;
}

this.getRMS = function () 
{
	if(this.NState<=0 || this.LevelsNeedToBeSearched) return -1.0;
	norm = 0.0;
	rms  = 0.0;
	
	if(this.calculateWF(this.State))
		for(ir=0;ir<this.NR;ir++)
		{
			r = ir*this.dR;
			norm = 1.*norm + 1.*this.psi[ir]*this.psi[ir];
			rms  = 1.*rms  + 1.*this.psi[ir]*this.psi[ir]*r*r;
		}
		
	return Math.sqrt(rms/norm);
}
}