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; } } }