using System;
using ODE = Science.Mathematics.OrdinaryDifferentialEquation;
namespace Science.Physics.QuantumMechanics
{
///
/// SchroedingerEquation for Two Dimensional Space
/// In units of hbar=m=1.
/// System must have rotational symmetry.
///
public class SchroedingerEquation2D
{
public SchroedingerEquation2D()
{
}
int numberOfVariables = 5;
int numberOfParameters = 2;
double start = 0.0;
double end = 5.0;
public double RadialBoundary
{
set{end = value/2.0;}
}
int qn = 0;
public int AzimuthalQuantumNumber
{
set{qn=value;}
}
double coeff = 0.0;
public double CoefficientOfSingularAtOrigin
{
set{coeff = value;}
}
private Science.Mathematics.Function.ToLastType potential;
public Science.Mathematics.Function.ToLastType 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 atzero = 1.0;
public double ValueAtOriginGuess
{
set{atzero=value;}
}
private bool ok;
public bool SolvedQ
{
get{return ok;}
}
private double[] waveFunctionR;
public double[] WaveFunctionRadialCoordinate
{
get { return waveFunctionR; }
}
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;}
}
public void Solve()
{
Science.Mathematics.Function.ToLastType load
= new Science.Mathematics.Function.ToLastType(ValuesAtStart);
Science.Mathematics.Function.ToLastType score
= new Science.Mathematics.Function.ToLastType(RestrictionOfContinium);
ODE.BoundaryCondition bc
= new ODE.BoundaryCondition(start,load,end,score);
Science.Mathematics.Function.ToLastType func
= new Science.Mathematics.Function.ToLastType(Derivs);
ODE.EquationWithBoundaryCondition eq
= new ODE.EquationWithBoundaryCondition();
eq.FirstDerivativeFunction = func;
eq.Condition = bc;
eq.HowMany = 200;
double[] initialParameter = new Double[numberOfParameters];
initialParameter[0] = energyguess;
initialParameter[1] = atzero;
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;
waveFunctionR = 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++)
{
waveFunctionR[k] = Convert.ToDouble(eq.IndependentVariable[k]);
waveFunction[k] = Math.Pow(waveFunctionR[k],Math.Abs(qn))*Convert.ToDouble(eq.Solution[1][k]);
if (qn == 0) waveFunctionD[k] = Convert.ToDouble(eq.Solution[2][k]);
else waveFunctionD[k] = Math.Abs(qn)*Math.Pow(waveFunctionR[k],Math.Abs(qn)-1)
* Convert.ToDouble(eq.Solution[1][k])
+ Math.Pow(waveFunctionR[k], Math.Abs(qn)) * Convert.ToDouble(eq.Solution[2][k]);
}
for(int k = 0; k < eq.IndependentVariable.Count; k++)
{
waveFunctionR[eq.IndependentVariable.Count+k] = 2.0*end-1.0*Convert.ToDouble(eq.IndependentVariable[eq.IndependentVariable.Count-1-k]);
waveFunction[eq.IndependentVariable.Count + k] = Math.Pow(waveFunctionR[eq.IndependentVariable.Count + k], Math.Abs(qn)) * Convert.ToDouble(eq.Solution[3][eq.IndependentVariable.Count - 1 - k]);
if (qn == 0) waveFunctionD[eq.IndependentVariable.Count + k] = Convert.ToDouble(eq.Solution[4][eq.IndependentVariable.Count - 1 - k]);
else waveFunctionD[eq.IndependentVariable.Count+k] = Math.Abs(qn)
*Math.Pow(waveFunctionR[eq.IndependentVariable.Count+k],Math.Abs(qn)-1)
* Convert.ToDouble(eq.Solution[3][eq.IndependentVariable.Count - 1 - k])
+ Math.Pow(waveFunctionR[eq.IndependentVariable.Count + k], Math.Abs(qn)) * 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] = p[1];
y[2] = coeff*p[1]/(Math.Abs(qn)+0.5);
y[3] = 0.0;
y[4] = slope;
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];
if(Math.Abs(x)>0.000000000001)
dydx[2]=2.0*(potential(x)*y[1]-y[0]*y[1])-(2.0*Math.Abs(qn)+1.0)/x*y[2];
else
dydx[2]=2.0*(potential(0.000000000001)*y[1]-y[0]*y[1])
-(2.0*Math.Abs(qn)+1.0)/0.000000000001*y[2];
dydx[3]=-y[4];
dydx[4]=-2.0*(potential(2.0*end-x)*y[3]-y[0]*y[3])
+(2.0*Math.Abs(qn)+1.0)/(2.0*end-x)*y[4];
return dydx;
}
}
}