PART I : SOLVING THE EIGEN ENERGIES FOR DIFFERENT n-STATE

Solving the Radial function and find the corresponding eigen energy.

In ths section, we solve the radial function by guessing the Eigen energy. If the solution for corresponding energy satisfy the boundary condition, then that energy is the Eigen energy we need.

code for radial function and energy

(***********************************************************************************)
(*To close all notification*)
Off[General::luc]
Off[General::stop]
Off[DSolve::bvsing]
Off[InterpolatingFunction::dmval]
(***********************************************************************************)


(*Initialisation*)
(***********************************************************************************)
Package ="2";(*"1-->let the l= 0, and scan all the intersection point along x-axis","2-->scan the 1st intersection point for different value of l"*)
n=5;(*principle quantum number*)(*the energy will be collected from ground state to correspond (Subscript[E, n])*)
iota=0.00000001;(*add in to l to avoid some singularity*)
rmax=100;




(*to initialise the condition for package*)
(***************************************************************)
If[Package == "1", 
nlastpackage=0;
nini=0.0;
nlast=2000000.0;
]; 
If[Package == "2",
nlastpackage=n-1;
delta=30000.0;
nini=0.0;
nlast=delta/4.0;
 ];
(***********************************************************************************)
(*End ofnitialisation*)

(*visualisation for the progress of the program*)
(***************************************************************)
Print[ProgressIndicator[Dynamic[nn],{nini,nlast}]];
(***************************************************************)



(*main part of program*)
(***************************************************************)
For[npackage=0,npackage≤nlastpackage,npackage++

(*initialise the condition for package for the loop run nlastpackage more than 1*)
(***************************************************************)
If[Package == "1", 
l=0.0+iota;
ex=(0.5/((l+1)^2)+0.5/(10*(l+1)^2));
nini=0.0;
nlast=2000000.0;
epsilon=ex/(nlast-nini);
]; 
If[Package == "2",
l=(npackage-1.0)+iota;
ex=(0.5/((l+1)^2)+0.5/(10*(l+1)^2));
delta=30000.0;
nini=0.0;
nlast=delta/4.0;
epsilon=ex/delta;
(*nl=(n-1.0);*)
 ];
(***************************************************************)


(*radial function*)
(***************************************************************)
einitial=-(ex);(*the initial guess for energy*)
ee[a_]:=einitial+a*epsilon;(*keep changing the guess energy by the step of epsilon*)
sol[n_]:=u[r]/.DSolve[{-r*u'[r]-((r^2)/2)u''[r]-(r-0.5*l*(l+1))u[r]==(r^2)* ee[nn]*u[r],u[rmax]==rmax*Exp[-rmax],u[rmax-epsilon]==(rmax-epsilon)*Exp[-(rmax-epsilon)]},u[r],r][[1]];(*solve for radial function*)
(***************************************************************)


(*Loop to collect the data and plot the U[0] versus n*)
(***************************************************************)
For[nn=nini,nn≤nlast,nn+=500,
soll=sol[nn];
f[n_,rdum_]:=soll/.{r->rdum};
function[nn]=f[nn,0.0000001](*U(0)*)
];
a=Table[{nn,Re[(function[nn])]/. C[1]->1},{nn,nini,nlast,500}];(*store all data in list a*)(*list a for some n will contain the constant and complex number,C[1]->1 is use to replace it with a constant and Re[] execute to consider the real part of number*)
(*ListPlot[a(*,AxesOrigin->{0.0,0.0}*)]*)
xx=Interpolation[a];(*to smooth the function to a line instead of a series of points*)
Print[Plot[xx[nn],{nn,nini,nlast},PlotLabel->"U[0] versus n",AxesLabel->{"n","U[0]"}]];(*plot the U[0] versus n*)
(***************************************************************)

(*Bisection method*)
(***************************************************************)
ft[y_]:=xx[y];(*function gives by interpolation*)
xF=nlast;
x1=nini;
x2=nini+1000;
xm=N[(x1+x2)/2.0];
tol=0.0001;
root=1;
runlimit=50;(*set to limit the loop,because in higher n, the value will never fall below the tolarent*)
nl=0.0;

For[x1=nini,x1<xF,{x1=x2,
       x2=x1+1000;},
run=1;
{While[ft[x1]*ft[x2]<0,

xm=N[(x1+x2)/2.0];

{While[
Abs[ft[xm]]>tol,
xm=N[(x1+x2)/2.0];

If[ft[x1]*ft[xm]<0,
x2=xm;If[run==runlimit,Break[]];
,(*else*)
x1=xm;If[run==runlimit,Break[]];
];(*end of if*)
run++;
];(*end of while*)
Print["The ",
Subscript["E",If[Package =="1",root,npackage]]
," is ",ee[xm]];
If[Package =="1",no=root,no=npackage];
filename="data"<>ToString[no]<>".dat";
Export[filename,ee[xm],"Data"];
root++;nl++;
};
{       x1=x2,
       x2=x1+1000;}
]};(*end of while*)
If[nl==n,Break[]];
    ](*end of for*)
(***************************************************************)
(*end of bisection*)

Clear[ex,a,f,function];(*clear some data which will e*)

]
(*end of the radial part*)



TimeUsed[]
EmitSound[Sound[SoundNote[]]]

PART II : SOLVING THE ANGULAR AND RADIAL FUNCTION

In this part, we are going to solve all the angular (Φ,Θ) and radial functions(R). Once these functions were obtained, we need to normalise them by dividing the functions with (Sqrt[area]).

Source Code for Wavefunction

(******************)
(* PHI function *)

phiupperB=2Pi;(* upper boundary of ϕ *)
philowerB=0.0;(* lower boundary of ϕ *)

phiBC=1;
phiupperBC=phiBC; (* upper boundary condition of ϕ *)
philowerBC=phiBC; (* lower boundary condition of ϕ *)

squarephifunction[ml_,ϕ_]:=(phifunction[ml,ϕ])^2;(* to square the unnormalised ϕ function *)
normalisedphifunction[ml_,ϕ_]:=(1/Sqrt[phiarea])*phifunction[ml,ϕ];(* to normalised the ϕ function *)

(* create a delayed function to calculate NDSolve[ ], area *)
phicalculation:=
(
phisolution[ml_]:=NDSolve[{phi''[ϕ]==-(ml^2)phi[ϕ],
phi[philowerB]==phiupperBC,
phi[phiupperB]==philowerBC},
phi,{ϕ,philowerB,phiupperB},PrecisionGoal->10];
phifunction[ml_,ϕ_]:=Evaluate[phi[ϕ]/.(phisolution[ml])];
phiarea=NIntegrate[squarephifunction[ml,ϕ],{ϕ,philowerB,phiupperB}];
);

PHI[ml_,ϕ_]:=normalisedphifunction[ml,ϕ];(* finally normalised ϕ function *)

(******************)


(******************)
(* THETA function *)
(*+++++++++++++++++++++++++++++++*)

value=N[0.9999999,10];(* to define a value which is very close to 1.0 *)

Clear[thetaBC];

thetaupperB=value;(* to set the upper interval of x *)
thetalowerB=-value;(* to set the lower interval of x *)
(************************)
(* in the following three steps, i m going to count how many decimal place of 'value' *)
realdigit=RealDigits[value,10];
list=TakeWhile[realdigit[[1]],#==9&];
decimalplace=Length[list];
(************************)
(* we need to define the interpolating range to be slightly out of the interval [-value,value] *)
thetaupperinterpolate=N[(thetaupperB+Last[list]*10^((-decimalplace)+(-1))),10];
thetalowerinterpolate=N[(thetalowerB-Last[list]*10^((-decimalplace)+(-1))),10];
(* we need to define the plot range for Θ function to be slightly inside of the interval [-value,value], 
so that no extrapolation will be needed in plotting the Θ function *)
upperplotrange=N[(thetaupperB-Last[list]*10^((-decimalplace)+(-1))),10];
lowerplotrange=N[(thetalowerB+Last[list]*10^((-decimalplace)+(-1))),10];


(****************************)
squarethetafunction[l_,ml_,g_]:=(thetafunction[l,ml,g])^2;(* to square up the Θ(x) function *)
normalisedthetafunction[l_,ml_,g_]:=(1/Sqrt[(thetaarea)])*thetafunction[l,ml,g];(* to divide the Θ(x) with the area of square under the curve *)
(****************************)


(* Case 1 : l is odd, ml is zero *)
(* the Θ(x) is being changed to W(x) *)
thetaconditionoddlzeroml:=
(
thetaBC=1;(* define boundary conditions *)
thetasolution[l_,ml_]:=NDSolve[{(1-g^2)theta''[g]-(2/g)*theta'[g]+((2/(g^2))+l(l+1)-(ml^2)/(1-g^2))theta[g]==0,
theta[thetalowerB]==thetaBC,
theta[thetaupperB]==thetaBC},
theta,{g,thetalowerinterpolate,thetaupperinterpolate},PrecisionGoal->10];(* solve the 2nd order ODE *)
thetafunction[l_,ml_,g_]:=Evaluate[((theta[g])/g)/.(thetasolution[l,ml])];(* obtain the Θ(x) function instead of W(x) *)
thetaarea=NIntegrate[squarethetafunction[l,ml,g],{g,thetalowerinterpolate,thetaupperinterpolate}];(* find the area *)
THETA[l_,ml_,g_]:=normalisedthetafunction[l,ml,g];(* finally normalised function *)
);

(* Case 2 : l is even, ml is zero *)
thetaconditionevenlzeroml:=
(
thetaBC=1.0;(* define boundary conditions *)
thetasolution[l_,ml_]:=NDSolve[{(1-g^2)theta''[g]-2g*theta'[g]+(l(l+1)-(ml^2)/(1-g^2))theta[g]==0,
theta[thetalowerB]==thetaBC,
theta[thetaupperB]==thetaBC},
theta,{g,thetalowerinterpolate,thetaupperinterpolate},PrecisionGoal->10];(* solve the 2nd order ODE *)
thetafunction[l_,ml_,g_]:=Evaluate[(theta[g])/.(thetasolution[l,ml])];(* obtain the Θ(x) function *)
thetaarea=NIntegrate[squarethetafunction[l,ml,g],{g,thetalowerinterpolate,thetaupperinterpolate}];(* find the area *)
THETA[l_,ml_,g_]:=normalisedthetafunction[l,ml,g];(* finally normalised function *)
);


(* Case 3 : ml is non-zero *)
thetaconditionnonzeroml:=
(
thetaBC=N[1*(10^(-decimalplace)),10];(* define boundary conditions *)
thetasolution[l_,ml_]:=NDSolve[{(1-g^2)theta''[g]-2g*theta'[g]+(l(l+1)-(ml^2)/(1-g^2))theta[g]==0,
theta[thetalowerB]==thetaBC,
theta[thetaupperB]==thetaBC},
theta,{g,thetalowerinterpolate,thetaupperinterpolate},PrecisionGoal->10];(* solve the 2nd order ODE *)
thetafunction[l_,ml_,g_]:=Evaluate[(theta[g])/.(thetasolution[l,ml])];(* obtain the Θ(x) function *)
thetaarea=NIntegrate[squarethetafunction[l,ml,g],{g,thetalowerinterpolate,thetaupperinterpolate}];(* find the area *)
zero:=0;
THETA[l_,ml_,g_]:=
Piecewise[
{{normalisedthetafunction[l,ml,g],thetalowerB≤g≤thetaupperB},
{zero,g<thetalowerB},
{zero,g>thetaupperB}}];(* finally normalised function *)
(*need to piecewise the function so that it is smooth at the interval of [-1,-value] and [value,1] *) 
);


(* a delayed function, so that different combination will get its corresponding boundary condition *)
thetacalculation:=
(
If[OddQ[l]&&ml==0,
thetaconditionoddlzeroml
];

If[EvenQ[l]&&ml==0,
thetaconditionevenlzeroml
];

If[(OddQ[l]||EvenQ[l])&&ml≠0,
thetaconditionnonzeroml
];

);


(*+++++++++++++++++++++++++++++++*)


(*****************)


(*****************)
(* RADIAL function *)

EigenEnergy[n_,l_]:=Abs[Energy[n,l]];(* delayed function, so that when energy is provided, it can be called into the radial function *)
shift=0.1;(* define an 'artificial' shift, △r *)
RupperB=100;(* approaches ∞ *)
approximatezero=0.0000001;
zero:=0;
RlowerB=approximatezero;
epsilon=N[10^-6];

Rupperinterpolate=RupperB+0.1;(* the interpolating range need to be slightly beyond the interval [RlowerB,RupperB] *)
Rlowerinterpolate=RlowerB/10;(* the interpolating range need to be slightly less than the interval [RlowerB,RupperB] *)

squareRfunction[n_,l_,r_]:=(Rfunction[n,l,r])^2;(* square the R function *)
normalisedRfunction[n_,l_,r_]:=(1/Sqrt[(Rarea)])*Rfunction[n,l,r];(* to divide R(r) with area *)

(* Case 1 : l≠0 *)
radialconditionnonzerol:=
(
Rsolution[n_,l_]:=NDSolve[{-(r^2)*(R''[r])-2r*R'[r]-2*r*R[r]+(l)*(l+1)*(R[r])==2*(r^2)*EigenEnergy[n,l]*R[r],
R[RupperB]==RupperB*Exp[-RupperB],
R[RlowerB]==approximatezero},
R,{r,zero,Rupperinterpolate},PrecisionGoal->10];
Rfunction[n_,l_,r_]:=Evaluate[((R[r]))/.(Rsolution[n,l])];(* to obtain the R(r) function *)
Rarea=NIntegrate[(r^2)*squareRfunction[n,l,r],{r,RlowerB,RupperB}];
RADIAL[n_,l_,r_]:=normalisedRfunction[n,l,r];
);

(* Case 2 : l=0 *)
(* need to transform into u(r)=r*R(r) *)
radialconditionzerol:=
(
Rsolution[n_,l_]:=NDSolve[{-(r^2)*(u''[r])-2*r*u[r]+(l)*(l+1)*(u[r])==2*(r^2)*EigenEnergy[n,l]*u[r],
u[RupperB]==RupperB*Exp[-RupperB],
u[RlowerB]==approximatezero},
u,{r,Rlowerinterpolate,Rupperinterpolate},PrecisionGoal->10];
Rfunction[n_,l_,r_]:=Evaluate[((u[r])/(r+shift))/.(Rsolution[n,l])];(* to obtain the R(r) function instead of u(x) *)
Rarea=NIntegrate[(r^2)*squareRfunction[n,l,r],{r,RlowerB,RupperB}];
RADIAL[n_,l_,r_]:=normalisedRfunction[n,l,r];
);


(* a delayed If-loop for solving R(r) *)
radialcalculation:=
(
If[l≠0,
radialconditionnonzerol;
];
If[l==0,
radialconditionzerol;
];
);


(*****************)


(*****************)
(* Final Wavefunction *)
(* a multiplication of all the Θ(x),Φ(ϕ) and R(r) *)
WAVEFUNCTION[n_,l_,ml_,ϕ_,g_,r_]:=(PHI[ml,ϕ])*(THETA[l,ml,g])*(RADIAL[n,l,r]);

(*****************)

PART III : OBTAINING THE PROBABILITY DENSITY (CARTESIAN)

In this part, hundreds or thousands of points will be transformed discretely into the Cartesian coordinate so that they can be visualised in the latter part.

Density Plot

(***********************)
(***********************)
(* NEED TO DEFINE *)
rlast=20.0;
(***********************)
(***********************)


(* ϕ plot range and ϕ[] function *)
numberofϕ=360;
ϕini=0.0;
ϕlast=2Pi;
ϕinterval=N[2Pi/numberofϕ];
ϕangle[a_]:=ϕini+(a-1)*ϕinterval;

(* r plot range and R[] function *)
rini=0.0;
rinterval=0.1;
numberofr=IntegerPart[(rlast-rini)/rinterval];
rlocation[a_]:=rini+(a-1)*rinterval;

(* θ plot range and Θ[] function *)(* but we need to convert x to θ first *)
twiceθupperB=2Pi;
θdegreeupperB=180;
numberofθ=θdegreeupperB+1;
gini=-1.0;(* lower boundary of substitution x *)
glast=1.0;(* upper boundary of substitution x *)
ginterval=N[(glast-gini)/numberofθ];(* to find the interval size of x *)
θini=ArcCos[glast];
θlast=ArcCos[gini];
θinterval=N[Abs[ArcCos[-1.0]-ArcCos[1.0]]/numberofθ];
θangle[a_]:=θini+(a-1)*θinterval;
gsubstitution[a_]:=Cos[θangle[a]];(* define function of x *)


(**********************************)
(* Cartesian coordinate *)
X[a_,b_,c_]:=rlocation[c]*Sin[θangle[b]]*Cos[ϕangle[a]];
Y[a_,b_,c_]:=rlocation[c]*Sin[θangle[b]]*Sin[ϕangle[a]];
Z[b_,c_]:=rlocation[c]*Cos[θangle[b]];
(**********************************)


(**********************************)
CartesianWavefunctionValue[n_,l_,ml_,a_,b_,c_]:=
WAVEFUNCTION[n,l,ml,ϕangle[a],gsubstitution[b],rlocation[c]];

squareCartesianWavefunctionValue[n_,l_,ml_,a_,b_,c_]:=
(CartesianWavefunctionValue[n,l,ml,a,b,c])^2;
(**********************************)


(* We will not plot the XY-plane cause it contains only the ϕ-information only *)
(*
XYCartesianProbabilityDensity[n_,l_,ml_,a_,b_,c_]:=
{X[a,b,c],Y[a,b,c],squareCartesianWavefunctionValue[n,l,ml,a,b,c][[1]]};
*)


(***************)
(* plot the probability density from 0 to Pi *)
rightXZCartesianProbabilityDensity[n_,l_,ml_,a_,b_,c_]:=
{X[a,b,c],Z[b,c],squareCartesianWavefunctionValue[n,l,ml,a,b,c][[1]]};
(***************)


(***************)
(* plot the probability density from Pi to 2Pi *)
mirrorX[a_,b_,c_]:=rlocation[c]*Sin[twiceθupperB-θangle[b]]*Cos[ϕangle[a]];
mirrorZ[b_,c_]:=rlocation[c]*Cos[twiceθupperB-θangle[b]];

leftXZCartesianProbabilityDensity[n_,l_,ml_,a_,b_,c_]:=
{mirrorX[a,b,c],mirrorZ[b,c],squareCartesianWavefunctionValue[n,l,ml,a,b,c][[1]]};
(***************)


(***************)
(* to combine the density plot from 0 to Pi and Pi to 2Pi *)
XZCartesianProbabilityDensity[n_,l_,ml_,a_,b_,c_]:=
{rightXZCartesianProbabilityDensity[n,l,ml,a,b,c],leftXZCartesianProbabilityDensity[n,l,ml,a,b,c]};
(***************)


(* information on XY-plane will not be used *)
(* ArcCos[g]==Pi/2 , X and Y only for RADIAL and PHI *) 
(*
gfixed=Cos[Pi/2];
gfixednumber=IntegerPart[Abs[(gfixed-gini)]/ginterval];
j=gfixednumber+1;
XYPlotList[n_,l_,ml_]:=
Flatten[Table[XYCartesianProbabilityDensity[n,l,ml,i,j,k],{i,1,numberofϕ},{k,1,numberofr+1}],1];
*)



(* ϕ==0 , X and Z only for RADIAL and THETA *) 
ϕfixed=0;
ϕfixednumber=IntegerPart[Abs[(ϕfixed-ϕini)/ϕinterval]];
i=ϕfixednumber+1;
XZPlotList[n_,l_,ml_]:=
Flatten[Table[XZCartesianProbabilityDensity[n,l,ml,i,j,k],{j,1,numberofθ},{k,1,numberofr+1}],2];(* to produce the list for density plot for 2Pi *)

PART IV : VISUALISING THE PROBABILITY DENSITY PLOT

The For-loop in this section is to import the data file saved in Part I and to execute all the delayed functions in Part II and Part III.

Compiling For-loop    (* n = 5 *)

(*********************)
(* important, must initialise *)
clast=5;(* define n-state *)
(*********************)


For[c=clast,c≥0,c--,

n=c;
filename="data"<>ToString[n]<>".dat";(* to import the energy into the kernel *)

For[d=c-1,d≥0,d--,

l=d;
Energy[n,l]=Import[filename,"Data"][[1]][[1]];(* to define the Energy *)

(*
Print[Energy[n,l]];
*)

(************************)
(* to the plot the R(r) only *)
(*
Print["  n is ",n];
Print[" ; l is ",l];
Print[" RADIAL function is"];
*)
radialcalculation;
(*
Print[Plot[RADIAL[n,l,r],{r,0,20}, PlotRange->All,AxesOrigin->{0,0}]];
*)
(************************)


For[e=d,e≥0,e--,

ml=e;

(************************)
(* to the plot the Θ(x) only *)
(*
Print["  n is ",n];
Print[" ; l is ",l];
Print["  ; ml is ",ml];
Print[" THETA function is"];
*)
thetacalculation;
(*
Print[Plot[THETA[l,ml,x],{x,-1.0,1.0}]];
*)
(************************)


(************************)
(* to the plot the Φ(ϕ) only *)
(*
Print["  n is ",n];
Print[" ; l is ",l];
Print["  ; ml is ",ml];
Print[" PHI function is"];
*)
phicalculation;
(*
Print[Plot[PHI[ml,ϕ],{ϕ,philowerB,phiupperB}, PlotRange->All]];
*)
(************************)


(*********)
(* calculate the discrete density plot points and export them *)
data[n,l,ml]=XZPlotList[n,l,ml];(* to generate the numerics for ListDensityPlot *)

(* Case A : to export the numerics for list density plot *)
(*
filename="data"<>ToString[n]<>ToString[l]<>ToString[ml]<>".dat";
Export[filename,data[n,l,ml],"Data"];(* export the numerics list *)
*)


(* Case B : to export directly the images of well-plotted list density plot *)
imagefilename="image"<>ToString[n]<>ToString[l]<>ToString[ml]<>".jpeg";
Export[imagefilename,ListDensityPlot[XZPlotList[n,l,ml],ColorFunction->"SunsetColors"]];
(*********)

];

];

];

If all the data file or image file has been saved properly in MyDocuments, one needs to execute this section only because its main purpose is to import all the data file or image file back into MATHEMATICA kernel.

List Density Plot

(************************)
(* NEED TO DEFINE *)
clast=5;
(************************)


For[c=clast,c≥0,c--,
n=c;
For[d=c-1,d≥0,d--,
l=d;
For[e=d,e≥0,e--,
ml=e;

(*********)
(* Case A : to import all the lists in MyDocuments and plot the using ListDensityPlot *)
(*
filename="data"<>ToString[n]<>ToString[l]<>ToString[ml]<>".dat";
data[n,l,ml]=Import[filename,"Data"];
Print["Density plot at (",n,",",l,",",ml,")"];
lplot[n,l,ml]=ListDensityPlot[data[n,l,ml],ColorFunction->"SunsetColors"];
*)
(*********)

(*********)
(* Case B : all the density plot has been saved as images, they will be imported back into the kernel *)
imagefilename="image"<>ToString[n]<>ToString[l]<>ToString[ml]<>".jpeg";
image[n,l,ml]=Import[imagefilename,"JPEG"];
(*********)
];

];

];

Visualisation

(* Case A : to import back all the lists and plot them using ListPlotDensity *)
(* 
Manipulate[lplot[n,l,ml],{n,1,clast,1},{l,0,n-1,1},{ml,0,l,1}] 
*)

(* Case B : to import back all the images and directly manipulate them *)
Manipulate[image[n,l,ml],{n,1,clast,1},{l,0,n-1,1},{ml,0,l,1}]