using System;
using ODE = Science.Mathematics.OrdinaryDifferentialEquation;

namespace Science.Physics.QuantumMechanics
{
	/// <summary>
	/// SchroedingerEquation for One Dimensional Space
	/// In units of hbar=m=1.
	/// </summary>

	public class SchroedingerEquation1D
	{

		public SchroedingerEquation1D()
		{
		}
        private int numberOfVariables = 5;
        private int numberOfParameters = 2;
        private double start = -10.0;
        private double end = 0.0;
		private double left = -9.0, right = 9.0;
		public double LeftBoundary
		{
			set{left = value;}
		}
		public double RightBoundary
		{
			set{right = value;}
		}
        private Science.Mathematics.Function.ToLastType<double,double> potential;
        public Science.Mathematics.Function.ToLastType<double, double> PotentialEnergy
		{
			set{potential=value;}
		}
		private double energyguess = 1.0;
		public double EnergyGuess
		{
			set{energyguess=value;}
		}
		private double slope = 1.0E-15;
		public double SlopeGuess
		{
			set{slope=value;}
		}
		private double[] waveFunctionX;
		public double[] WaveFunctionXCoordinate
		{
			get{return waveFunctionX;}
		}
        private double[] waveFunction;
        public double[] WaveFunction
        {
            get { return waveFunction; }
        }
        private double[] waveFunctionD;
        public double[] WaveFunctionDerivative
        {
            get { return waveFunctionD; }
        }
		private double energy;
		public double Energy
		{
			get{return energy;}
		}
        private bool ok;
        public bool SolvedQ
        {
            get { return ok; }
        }
		public void Solve()
		{
			start = (start < left ? start : left);
			start = (right < -start ? start : -right);
 
			Science.Mathematics.Function.ToLastType<double,double[],double[]> load
                = new Science.Mathematics.Function.ToLastType<double, double[], double[]>(ValuesAtStart);
            Science.Mathematics.Function.ToLastType<double, double[], double[], double[]> score
                = new Science.Mathematics.Function.ToLastType<double, double[], double[], double[]>(RestrictionOfContinium);
			ODE.BoundaryCondition bc 
				= new ODE.BoundaryCondition(start,load,end,score);

            Science.Mathematics.Function.ToLastType<double, double[], double[]> func
                = new Science.Mathematics.Function.ToLastType<double, double[], double[]>(Derivs);
			
			ODE.EquationWithBoundaryCondition eq 
				= new ODE.EquationWithBoundaryCondition();
            eq.FirstDerivativeFunction = func;
            eq.Condition = bc;
            eq.HowMany = 20;

			double[] initialParameter = new Double[numberOfParameters];
			initialParameter[0] = energyguess;   
			initialParameter[1] = 0.0;    
			eq.InitialGuess = initialParameter;
            
		    eq.Solve();

            if (Math.Abs((Convert.ToDouble(eq.Solution[1][eq.IndependentVariable.Count - 1]) + slope) / (Convert.ToDouble(eq.Solution[3][eq.IndependentVariable.Count - 1]) + slope) - 1.0) > 0.001 ||
                Math.Abs((Convert.ToDouble(eq.Solution[2][eq.IndependentVariable.Count - 1]) + slope) / (Convert.ToDouble(eq.Solution[4][eq.IndependentVariable.Count - 1]) + slope) - 1.0) > 0.001)
            {
                ok = false;
                try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Bad initial guess. Try again with modification.");
                }
            }
            else ok = true;

			waveFunctionX = new Double[eq.IndependentVariable.Count * 2];
            waveFunction = new Double[eq.IndependentVariable.Count * 2];
            waveFunctionD = new Double[eq.IndependentVariable.Count * 2];


			for(int k = 0; k < eq.IndependentVariable.Count; k++)
			{
				waveFunctionX[k] = Convert.ToDouble(eq.IndependentVariable[k]);
                waveFunction[k] = Convert.ToDouble(eq.Solution[1][k]);
                waveFunctionD[k] = Convert.ToDouble(eq.Solution[2][k]);
			}
			for(int k = 0; k < eq.IndependentVariable.Count; k++)
			{
				waveFunctionX[eq.IndependentVariable.Count + k] = -1.0*Convert.ToDouble(eq.IndependentVariable[eq.IndependentVariable.Count-1-k]);
                waveFunction[eq.IndependentVariable.Count + k] = Convert.ToDouble(eq.Solution[3][eq.IndependentVariable.Count - 1 - k]);
                waveFunctionD[eq.IndependentVariable.Count + k] = Convert.ToDouble(eq.Solution[4][eq.IndependentVariable.Count - 1 - k]);
			}
            energy = Convert.ToDouble(eq.Solution[0][0]);
		}		
	
		private double[] ValuesAtStart(double start, double[] p)
		{
			double[] y = new Double[numberOfVariables];
			y[0] = p[0];    
			y[1] = 0.0;	
			y[2] = slope;  
			y[3] = 0.0;
			y[4] = y[2]*p[1];
			return y;
		}
		private double[] RestrictionOfContinium(double end, double[] v, double[] p)
		{
			double[] f = new Double[numberOfParameters];
			f[0] = v[1]-v[3]; 
			f[1] = v[2]-v[4];
			return f;
		}
		private double[] Derivs(double x, double[] y)  // Schroedinger Equation in First Order
		{
			double[] dydx = new Double[y.Length];
			dydx[0]=0.0;
			dydx[1]=y[2];
			dydx[2]=2.0*(potential(x)*y[1]-y[0]*y[1]);
			dydx[3]=-y[4];
			dydx[4]=-2.0*(potential(-x)*y[3]-y[0]*y[3]);
			return dydx;
		} 	
	}  
}