(* Created in 2003 by author(s) Tran The Trung and Nguyen Hong Quang, Institute of Physics, NCST, 46 Nguyen Van Ngoc, Hanoi, Vietnam *) (* Version 1.1, 10/2003 *) (* Distributable free of charge for educational purposes provided the following acknowledgement is printed in a standard place: "SWSED.m is distributed with permission by Nguyen Hong Quang and Tran The Trung." ; otherwise contact author(s) *) (* Compatible to Mathematica version 3 or later *) Off[General::spell]; Off[General::spell1]; BeginPackage["SolidState`SWSED`","NumericalMath`BesselZeros`"]; Unprotect[EnergyLevel,EnergyStructure,StructurePlot, RadialWaveFunction, DensityFunction, PlotPDF, LevelLabel, StructureOrder,EnergyStructureQ]; Remove[EnergyLevel,EnergyStructure,StructurePlot, RadialWaveFunction, DensityFunction, PlotPDF, LevelLabel, StructureOrder,EnergyStructureQ]; EnergyLevel::usage = "EnergyLevel[{R,V},{l,n}] gives energy in eV of state of the single electron indicated by l and n within the square well quantum dot specified by {R,V} with default material as Si inside and SiO2 outside. Other material can be specified by using options EMass and Permeability. Type ?EMass and ?Permeability for more help on these options."; EnergyStructureQ::usage="EnergyStructure[data] gives True if data is output of EnergyStructure function, and False ortherwise. EnergyStructureQ has no option."; EnergyStructure::usage = "EnergyStructure[R,V] gives energies (in eV) of all bound states of the single electron within the square well quantum dot specified by {R,V} with default material as Si inside and SiO2 outside. Other material can be specified by using options EMass and Permeability. Type ?EMass and ?Permeability for more help on these options."; StructurePlot::usage = "StructurePlot[output of EnergyStructure] plots all bound states calculated by EnergyStructure of the single electron within the square well quantum dot input to EnergyStructure. To check if a data is output of EnergyStructure, use the function EnergyStructureQ. StructurePlot has the following options: . To get more information about each option, type ? infront of the option name."; RadialWaveFunction::usage = "RadialWaveFunction[{R,V},{l,n}] gives wave function in standard Mathematica function format of the single electron inside square well quantum dot {R,V} at the state {l,n}. The dot has default material as Si inside and SiO2 outside. Other material can be specified by using options EMass and Permeability. Type ?EMass and ?Permeability for more help on these options."; DensityFunction::usage = "DensityFunction[{R,V},{l,n}] gives probability density function in standard Mathematica function format of the random distribution of the single electron inside the space of square well quantum dot {R,V} at the state {l,n}. The dot has default material as Si inside and SiO2 outside. Other material can be specified by using options EMass and Permeability. Type ?EMass and ?Permeability for more help on these options."; EnergyLevel::badin = "Please check input. Correct form should be EnergyLevel[{R,V},{l,n}] with l = 0,1, ..., n = 1,2,.., R is positive radius of barrier in Angstrom, V is positive barrier potential in eV and Min, Mout are effective electron mass inside and out side barrier."; EnergyLevel::unbound = "EnergyLevel encounters level demanded unbound or potential well too shallow."; EnergyLevel::noroot = "EnergyLevel finds no bound energy level, output Null."; EnergyLevel::badfun = "EnergyLevel encounters numerically unstable matching equation to solve."; EnergyLevel::bessel = "Matching equation badly behaves near root, EnergyLevel assumes root close to Bessel root."; PlotPDF::error = "Option PlotPDF for StructurePlot must be either True or False (default value is False)."; StructurePlot::nobound = "StructurePlot receives no bound state to plot."; Options[StructurePlot] = {PlotPDF -> False, LevelLabel->Automatic,Frame->True, PlotStyle->{GrayLevel[0]},FrameLabel->{"Quantum number l","Energy in eV"},PlotLabel->None,AxesLabel->None,Axes->False,AspectRatio->1/GoldenRatio,DisplayFunction->$DisplayFunction,ImageSize->Automatic,LevelShape->Automatic}; PlotPDF::usage = "PlotPDF is an option for StructurePlot. Setting PlotPDF -> True cause the StructurePlot to plot probability density function graph above each energy level. Default value is False which causes StructurePlot to plot only the energy levels."; LevelLabel::usage = "LevelLabel is an option for StructurePlot. It specifies the text to be printed next to each energy level. Default values is Automatic, which causes StructurePlot to print standard atomic label for each level. You can set LevelLabel to any expression or list of expressions, StructurePlot will try to print them next to each energy level in order from bottom to top and from left to right."; LevelShape::usage = "LevelShape is an option for StructurePlot. It is effective only when option PlotPDF of StructurePlot is set to False. The setting of LevelShape -> None causes StructurePlot to draw no shape at the energy levels. Default value is Automatic, which makes StructurePlot to plot a horizontal line at each energy level. When PlotPDF -> True, each energy level will have a graph of probability density function associated with the level."; StructureOrder::usage = "StrutureOrder[output of EnergyStructure] gives {{R,V},{Order List}} where OrderList is the list of string represent the energy level's name, under ascending order of level, for all bound states calculated by EnergyStructure of the single electron within the square well quantum dot {R,V}. To check if a data is output of EnergyStructure use the function EnergyStructureQ. StructureOrder has no option."; StructureOrder::badin = "StructureOrder receives no energy structure level."; LevelShape::error = "LevelShape option for StructurePlot can take value of only either Automatic or None. For unrecognized input LevelShape, Automatic value is then assumed."; LevelShape::usage = "Option LevelShape specify the shape of level. LevelShape option for StructurePlot can take value of only either Automatic or None."; Options[EnergyLevel] = {EMass->{0.27,0.5}, Permeability->{11.7,11.7}}; Options[EnergyStructure] = {EMass->{0.27,0.5}, Permeability->{11.7,11.7}}; Options[RadialWaveFunction] = {EMass->{0.27,0.5}, Permeability->{11.7,11.7}}; Options[DensityFunction] = {EMass->{0.27,0.5}, Permeability->{11.7,11.7}}; EMass::error = "Options EMass must be a list of two real positive numbers."; Permeability::error = "Options Permeability must be a list of two real positive numbers."; EMass::usage = "EMass is an option for EnergyLevel, EnergyStructure, RadialWaveFunction and DensityFunction functions, it can be specified by a list of two positive numbers, representing electron effective mass inside and outside of SWSED, in unit of free electron mass. Default value is {0.27, 0.5}, which is correct for Si-SiO2 dot."; Permeability::usage = "Permeabilitys is an option for EnergyLevel, EnergyStructure, RadialWaveFunction and DensityFunction functions, it can be specified by a list of two positive numbers, representing electric permeability inside and outside of SWSED. Default value is {11.7, 11.7}, which is correct for Si-SiO2 dot." StructureOrder::noop = "StructureOrder has no option. Type ?StructureOrder for the usage of this function. Also type ? infront of the option name for the usage of the option you intend to use."; EnergyStructureQ::error="Input energy structure is not valid."; Begin["`Private`"]; (* Private definitions *) (* Length conversion to special unit: input length: x in angstrom output length: x Minside / aBohr Energy conversion to special unit: input energy: e in eV output energy: e / (Minside EHart) kappa: kappa = alpha x Minside / aBohr kappa = beta x Minside / aBohr alpha = Sqrt[e Minside EHart] beta = Sqrt[Moutside/Minside(V Minside EHart - e Minside EHart)] Note that however, we take effective approximation: aBohrEffInside = aBohr EpsilonInside EHartEffInside = EHart / (2 EpsilonInside^2) aBohrEffOutside = aBohr EpsilonOutside EHartEffOutside = EHart / (2 EpsilonOutside^2) *) (* Initiations *) aBohr = 0.529177; (* Angstrom*) EHart = 13.6; (* halfed in eV *) LName = Join[{"s", "p", "d", "f", "g", "h", "i", "k", "l", "m", "n", "o","q","r","t","u","v","w","x","y","z"},Table["?",{100}]]; Remove[lhs,rhs,FunctionEquat,Bisect,Secand,BiSe,RootSearch,RootAll1,RootAll2]; Clear[l,kappa,R,alpha,Vo,gamma,f,x1,x2,tolerance,rm,RootLevel]; FunctionEquat[l_,R_,Vo_,gamma_,kappa_]:=Module[{ t1 = kappa R ,beta ,t2 ,t3=l-1/2 ,t4=l+1/2 }, beta = Sqrt[gamma (Vo-kappa^2)]; t2 = beta R; gamma kappa BesselJ[t3,t1]/BesselJ[t4,t1] + beta BesselK[t3,t2]/BesselK[t4,t2] + (l+1) (1 - gamma)/R ]; FunctionEquat[0,R_,Vo_,gamma_,kappa_]:=(1-gamma)/R+Sqrt[gamma (Vo-kappa^2)]+gamma kappa Cot[R kappa]; FunctionEquat[1, R_, Vo_, gamma_, kappa_] := Module[{t1, t2,t3,t4,t5}, t1 = Vo-kappa^2; t2 = R Sqrt[gamma t1]; t3 = R kappa; t4 = Cos[t3]; t5 = Sin[t3]; -(R kappa (-2 (1 + t2) + gamma (2 + 2 t2 - R^2 t1)) t4 + (2 + 2 t2 + gamma (-2 + R^2 Vo - 2 t2 + t3^2 t2 )) t5) / (1+t2)(t3 t4 - t5) ]; (* eliminated factor 1/R from equation *) FunctionEquat[2, R_, Vo_, gamma_, kappa_] := Module[{t2,t3,t4,t5,t6,t7}, t2 = gamma (Vo-kappa^2); t3 = Sqrt[t2]; t4 = kappa R; t5 = Cos[t4]; t6 = Sin[t4]; t7 = gamma-1; (t4*(-27*t7*R*t2 + 3 *gamma*R^3*Vo*t2 - 27*t7*t3 + gamma*t4^2*R^2*t3^3 - 3*gamma*R^2*t3*(-3*t7*kappa^2 + (-4 + 3*gamma)*Vo))*t5 + (27*t7*R*t2 + t4^2*R^3*t2^2 + 27*t7*t3 - 4*t7*t4^2*R^2*t3^3 - 3*R^3*t2*(3*t7*kappa^2 + gamma*Vo) + 3*R^2*t3*(-3*(-1 + gamma^2)*kappa^2 + gamma*(-4 + 3*gamma)*Vo))*t6)/ ((3*R*t2 + 3*t3 + R^2*t3^3)*(3*t4*t5 + (-3 + t4^2)*t6)) ]; (* eliminated factor 1/R from equation *) Bisect[f_,x1_,x2_,tolerance_,f1_,f2_]:=Module[{i = 0, imax = 5000, x, fx, f1t=f1, f2t=f2, x1t = x1, x2t = x2}, If[f1t f2t > 0, Print["Wrong starting x1 and x2 for Bisect"], While[i < imax, x = 0.5(x2t + x1t); fx = Re[f[x]]; If[Abs[(x2t - x1t)/x] < tolerance, i = imax, If[f2t fx < 0, x1t = x; f1t = fx; , x2t = x; f2t = fx; ] ]; i = i + 1; If[i == imax, Print["Reach Max iteration ", imax, " before solution found"]; Print["Intermediate result x = ", x]; Print[" x1 = ", x1t]; Print[" x2 = ", x2t]; Print[" f1 = ", f1t]; Print[" f2 = ", f2t]; Print[" f = ", fx]; ]; ]; {x1t, x2t, f1t, f2t, x}] ]; Secand[f_, x1_, x2_, f1_, f2_] := Module[{i = 0, imax = 5000, x, fx, f1t = f1, f2t = f2, x1t = x1, x2t = x2, tolerance}, If[f1t f2t > 0, Print["Wrong starting x1 and x2 for Secand"], tolerance = 10^(-10)(Abs[f1t]+Abs[f2t]); While[i < imax, x = x1t - (x2t - x1t)f1t/(f2t - f1t); fx = Re[f[x]]; If[Abs[fx] < tolerance, i = imax, If[f2t fx < 0, x1t = x; f1t = fx; , x2t = x; f2t = fx; ] ]; i = i + 1; If[i == imax, Print["Reach Max iteration ", imax, " before solution found"]; Print["Intermediate result x = ", x]; Print[" x1 = ", x1t]; Print[" x2 = ", x2t]; Print[" f1 = ", f1t]; Print[" f2 = ", f2t]; Print[" f = ", fx]; ]; ]; {x1t, x2t, f1t, f2t, x}] ]; BiSe[f_, x1_, x2_, f1_, f2_] := Module[{x, x1t, x2t, f1t, f2t}, {x1t, x2t, f1t, f2t, x} = Bisect[f, x1, x2,10^(-3), f1, f2]; Secand[f,x1t,x2t,f1t,f2t][[5]] ]; RootSearch[R_?Positive, Vo_?Positive, rm_?Positive, l_Integer, bin_Integer] := Module[{BesselRoot = Join[{0}, BesselJZeros[l +1/2, bin + 1]/R], Step, xmin, xmax, binmax = 0, i, f1, f2, hihi=Sqrt[Vo]}, Do[If[BesselRoot[[i]] < hihi, binmax = binmax + 1], {i, 2, bin + 2}]; (* If[FunctionEquat[l, R, Vo, rm, 0.2 BesselRoot[[2]]] < 0, bin = bin + 1]; *) Step = 0.2(BesselRoot[[bin + 1]] - BesselRoot[[bin]]); If[bin > binmax, If[bin == binmax + 1, xmax = hihi - 0.01Step; f2 = FunctionEquat[l, R, Vo, rm, xmax]; If[f2< 0, xmin = BesselRoot[[bin]] + Step; xmin=If[xmin>hihi,hihi,xmin]; f1 = FunctionEquat[l, R, Vo, rm, xmin]; If[f1 <0, xmin = xmin + Step; xmin=If[xmin>hihi,hihi,xmin]; f1 = FunctionEquat[l, R, Vo, rm, xmin]; If[f1 <0,Message[EnergyLevel::badfun]; Null, BiSe[FunctionEquat[l, R, Vo, rm, #1] &, xmin, xmax, f1, f2]] , BiSe[FunctionEquat[l, R, Vo, rm, #1] &, xmin, xmax, f1, f2]] , Message[EnergyLevel::unbound]; ] , Message[EnergyLevel::unbound]; ] , xmin = BesselRoot[[bin]] + Step; xmax = BesselRoot[[bin + 1]] - 0.01Step; f2 = FunctionEquat[l, R, Vo, rm, xmax]; If[f2 > 0, xmax = xmax + 0.005 Step; (* xmax=If[xmax>hihi,hihi,xmax]; *) f1 = FunctionEquat[l, R, Vo, rm, xmax]; If[f1 > 0, xmax = xmax - f1/(f1 - f2) 0.01Step; (* xmax=If[xmax>hihi,hihi,xmax]; *) f2 = FunctionEquat[l, R, Vo, rm, xmax]; , f2 = f1]; ]; If[f2 < 0 , f1 = FunctionEquat[l, R, Vo, rm, xmin]; If[f1<0, xmin = xmin + Step; (* xmin=If[xmin>hihi,hihi,xmin]; *) f1 = FunctionEquat[l, R, Vo, rm, xmin]; If[f1 < 0, xmin = xmin - 1.9 Step; f1 = FunctionEquat[l, R, Vo, rm, xmin]; If[f1 < 0, Message[EnergyLevel::badfun]; Plot[FunctionEquat[l, R, Vo, rm, xmin], {xmin, 0, BesselRoot[[bin + 1]]},PlotRange -> {-f1, f1},Frame->True,FrameLabel->{"x","Matching Equation"},PlotLabel->"Graph of matching equation near root"]; ,BiSe[FunctionEquat[l, R, Vo, rm, #1] &, xmin, xmax,f1,f2] ] , BiSe[FunctionEquat[l, R, Vo, rm, #1] &, xmin, xmax,f1,f2] ] , BiSe[FunctionEquat[l, R, Vo, rm, #1] &, xmin, xmax,f1,f2] ] , Message[EnergyLevel::bessel]; BesselRoot[[bin + 1]] ] ] ]; RootAll1[R_?NumberQ, Vo_?NumberQ, rm_?NumberQ, l_Integer] := Module[{Lr = {}, r, i = 1, continue = True}, While[continue, r = RootSearch[R, Vo, rm, l, i]; If[NumberQ[r], Lr = Append[Lr, r]; i=i+1;, continue = False]; ]; Lr ]; RootAll2[R_?NumberQ,Vo_?NumberQ, rm_?NumberQ]:= Module[{Lr={},tLr,continue=True,l=0}, Off[EnergyLevel::unbound]; While[continue, tLr = RootAll1[R,Vo,rm,l]; If[tLr == {},continue = False, Lr = Append[Lr,tLr];l=l+1;]; ]; On[EnergyLevel::unbound]; Lr ]; (* Public functions *) Clear[R,V,Minside,Moutside,ES]; EnergyLevel[{R_?Positive,V_?Positive},{l_Integer,n_Integer},opts___Rule] := Module[{XC, EC, r,min,mout,epin,epout}, {min,mout}=EMass/.{opts}/.Options[EnergyLevel]; {epin,epout}=Permeability/.{opts}/.Options[EnergyLevel]; If[NumericQ[min]&&Positive[min]&&NumericQ[mout]&&Positive[mout], If[NumericQ[epin]&&Positive[epin]&&NumericQ[epout]&&Positive[epout], XC = min / aBohr / epin; EC = min EHart / epin^2; If[l<0||n<1,Message[EnergyLevel::badin], If[V==Infinity,(BesselJZeros[l +1/2,n][[n]]/(R XC))^2 EC, r=RootSearch[R XC, V / EC, mout/min, l,n]; If[NumericQ[r],r^2 EC,Message[EnergyLevel::noroot];]]] ,Message[Permeability::error]] ,Message[EMass::error]] ]; EnergyStructure[R_?Positive,V_?Positive,opts___Rule]:=Module[ {XC, EC,min,mout,epin,epout}, {min,mout}=EMass/.{opts}/.Options[EnergyLevel]; {epin,epout}=Permeability/.{opts}/.Options[EnergyLevel]; If[NumericQ[min]&&Positive[min]&&NumericQ[mout]&&Positive[mout], If[NumericQ[epin]&&Positive[epin]&&NumericQ[epout]&&Positive[epout], XC = min / aBohr / epin; EC = min EHart / epin^2; {{R,V,min,mout,epin,epout},RootAll2[R XC, V / EC, mout/min]^2 EC} ,Message[Permeability::error]] ,Message[EMass::error]] ]; StructurePlot[{{R_?Positive,V_?Positive,min_?Positive,mout_?Positive,epin_?Positive,epout_?Positive},ES_List},opts___Rule]:=Module[{H,L1,pdf=PlotPDF/.{opts}/.Options[StructurePlot], t, graph, frm=Frame/.{opts}/.Options[StructurePlot], plst=PlotStyle/.{opts}/.Options[StructurePlot], frmlb=FrameLabel/.{opts}/.Options[StructurePlot], pllb=PlotLabel/.{opts}/.Options[StructurePlot], axlb=AxesLabel/.{opts}/.Options[StructurePlot], ax=Axes/.{opts}/.Options[StructurePlot], asra=AspectRatio/.{opts}/.Options[StructurePlot], dpfn=DisplayFunction/.{opts}/.Options[StructurePlot], imsi=ImageSize/.{opts}/.Options[StructurePlot], lvlb=LevelLabel/.{opts}/.Options[StructurePlot], lvsh=LevelShape/.{opts}/.Options[StructurePlot], dH,i,m,RS,lblist,LES,cnt,XC, EC }, If[EnergyStructureQ[{{R,V,min,mout,epin,epout},ES}], LES = Length[Flatten[ES]]; L1=Length[ES]; If[LES == 0,Message[StructurePlot::nobound] , Switch[lvlb, Automatic, lblist={}; Do[ Do[lblist=Append[lblist,ToString[m+i-1] <> LName[[i]]]; , {m, Length[ES[[i]]]}] ,{i,L1}]; , _, If[Depth[lvlb]==1, lblist=Table[lvlb,{LES}]; , lblist=Table["",{LES}]; Do[lblist[[i]]=lvlb[[i]],{i,Length[lvlb]}]; ] ]; H = 1.1 V; cnt=1; Switch[pdf, False, graph=Flatten[{ plst , Switch[lvsh ,Automatic ,Table[Map[Line[{{i - 0.45, #1}, {i + 0.2, #1}}] &, ES[[i]]], {i, L1}] ,None ,plst ,_ ,Message[LevelShape::error]; Table[Map[Line[{{i - 0.45, #1}, {i + 0.2, #1}}] &, ES[[i]]], {i, L1}] ] , Switch[lvlb, None, Dashing[{0.018}],_, { Table[ Table[Text[lblist[[cnt++]], {i - 0.65, ES[[i, m]]}] , {m, Length[ES[[i]]]}] ,{i,L1}] ,Dashing[{0.018}] } ], Line[{{0,V},{L1+0.3,V}}] }]; Show[Graphics[graph],Frame->frm, FrameLabel->frmlb,PlotLabel->pllb,AxesLabel->axlb,Axes->ax,AspectRatio->asra,DisplayFunction->dpfn,ImageSize->imsi, FrameTicks -> {None, Automatic}, Ticks -> {None, Automatic}] ,True, XC = min / aBohr / epin; EC = min EHart / epin^2; t= XC / 0.85; dH=If[Length[ES[[1]]]>1,0.7(ES[[1,2]]-ES[[1,1]]),0.7(V-ES[[1,1]])]; RS = Sqrt[ES / EC]; graph=Flatten[{ Graphics[Flatten[ { plst, Table[ Table[ Line[{{i - 0.85, ES[[i,m]]},{i - 0.85, ES[[i,m]]+dH},{i, ES[[i,m]]+dH}, {i, ES[[i,m]]},{i + 0.1, ES[[i,m]]},{i - 0.95, ES[[i,m]]}}] , {m, Length[ES[[i]]]}] , {i, L1}] , Switch[lvlb,None, Dashing[{0.018}] ,_, { Table[ Table[ Text[lblist[[cnt++]], {i - 0.4, ES[[i, m]]-dH/8}] , {m, Length[ES[[i]]]}] ,{i,L1}] ,Dashing[{0.018}] }], Line[{{0,V},{L1+0.1,V}}] } ]] , Table[ Table[ Plot[dH Pi BesselJ[i-1/2, RS[[i,m]] t R (x-i+0.85)]^2 /(2 RS[[i,m]] t R (x-i+0.85))+ ES[[i,m]], {x,i-0.85,i}, DisplayFunction->Identity, PlotStyle->plst] , {m, Length[ES[[i]]]}] ,{i,1,L1}] }]; Show[graph,Frame->frm, FrameLabel->frmlb,PlotLabel->pllb,AxesLabel->axlb,Axes->ax,AspectRatio->asra,DisplayFunction->dpfn,ImageSize->imsi, FrameTicks -> {None, Automatic}, Ticks -> {None, Automatic}] ,_,Message[PlotPDF::error] ]] ,Message[EnergyStructureQ::error]] ]; RadialWaveFunction[{R_?Positive,V_?Positive},{l_Integer, n_Integer},opts___Rule]:=Module[ {alpha, beta, A, B, r, v, g, XC, EC, T1, T2, temp1, temp2, temp3, temp4,min,mout,epin,epout,temp5}, If[l<0||n<1,Message[EnergyLevel::badin], {min,mout}=EMass/.{opts}/.Options[EnergyLevel]; {epin,epout}=Permeability/.{opts}/.Options[EnergyLevel]; If[NumericQ[min]&&Positive[min]&&NumericQ[mout]&&Positive[mout], If[NumericQ[epin]&&Positive[epin]&&NumericQ[epout]&&Positive[epout], XC = min / aBohr / epin; EC = min EHart / epin^2; r = R XC; v = V / EC; g = mout/min; alpha=If[V==Infinity,BesselJZeros[l +1/2,n][[n]]/r,RootSearch[r, v, g, l,n]]; If[NumericQ[alpha], beta = Sqrt[g(v-alpha^2)]; temp1 = BesselJ[l+1/2, alpha r]; temp2 = BesselK[l+1/2, beta r]; temp3 = temp1^2; temp4 = temp2^2; T1 = temp3 - BesselJ[l-1/2,alpha r]BesselJ[l+3/2,alpha r]; T2 = temp3 (BesselK[l+3/2,beta r]BesselK[l-1/2, beta r]/temp4 -1); B = temp1/temp2; A = Sqrt[(T1+ T2)/2] R; temp5 = alpha XC; If[l==0, ToExpression[ "If[#==0," <>ToString[InputForm[Sqrt[temp5 2/Pi]/A]] <>",If[# <" <>ToString[InputForm[R]] <>",BesselJ[" <>ToString[InputForm[2 l + 1]] <>"/2," <>ToString[InputForm[temp5]] <>" #]/Sqrt[#]," <>ToString[InputForm[B]] <>" BesselK[" <>ToString[InputForm[2 l + 1]] <>"/2," <>ToString[InputForm[beta XC]] <>" #]/Sqrt[#]]/" <>ToString[InputForm[A]] <>"]&" ] , ToExpression[ "If[#==0,0" <>",If[# <" <>ToString[InputForm[R]] <>",BesselJ[" <>ToString[InputForm[2 l + 1]] <>"/2," <>ToString[InputForm[temp5]] <>" #]/Sqrt[#]," <>ToString[InputForm[B]] <>" BesselK[" <>ToString[InputForm[2 l + 1]] <>"/2," <>ToString[InputForm[beta XC]] <>" #]/Sqrt[#]]/" <>ToString[InputForm[A]] <>"]&" ] ] ,Message[EnergyLevel::noroot]] ,Message[Permeability::error]] ,Message[EMass::error]] ]]; DensityFunction[{R_?Positive,V_?Positive},{l_Integer, n_Integer},opts___Rule]:=Module[ {alpha, beta, A2, B2, r, v, g, XC, EC , T1, T2, temp1, temp2,min,mout,epin,epout,temp5}, If[l<0||n<1,Message[EnergyLevel::badin], {min,mout}=EMass/.{opts}/.Options[EnergyLevel]; {epin,epout}=Permeability/.{opts}/.Options[EnergyLevel]; If[NumericQ[min]&&Positive[min]&&NumericQ[mout]&&Positive[mout], If[NumericQ[epin]&&Positive[epin]&&NumericQ[epout]&&Positive[epout], XC = min / aBohr / epin; EC = min EHart / epin^2; r = R XC; v = V / EC; g = mout/min; alpha=If[V==Infinity,BesselJZeros[l +1/2,n][[n]]/r,RootSearch[r, v, g, l,n]]; If[NumericQ[alpha], beta = Sqrt[g(v-alpha^2)]; temp1 = BesselJ[l+1/2, alpha r]^2; temp2 = BesselK[l+1/2, beta r]^2; T1 = temp1 - BesselJ[l-1/2,alpha r]BesselJ[l+3/2,alpha r]; T2 = temp1 (BesselK[l+3/2,beta r]BesselK[l-1/2, beta r]/temp2 -1); B2 = temp1 / temp2; A2 = (T1+ T2) 2 Pi R^2; temp5 = alpha XC; If[l==0, ToExpression[ "If[#==0," <>ToString[InputForm[temp5 2/Pi/A2]] <>",If[# < " <>ToString[InputForm[R]] <>", BesselJ[" <>ToString[InputForm[2 l+1]] <>"/2," <>ToString[InputForm[temp5]] <>" #" <>"]^2/#, " <>ToString[InputForm[B2]] <>" BesselK[" <>ToString[InputForm[2 l+1]] <>"/2," <>ToString[InputForm[beta XC]] <>" #" <>"]^2/#]/" <>ToString[InputForm[A2]] <>"]&" ] , ToExpression[ "If[#==0,0" <>",If[# < " <>ToString[InputForm[R]] <>", BesselJ[" <>ToString[InputForm[2 l+1]] <>"/2," <>ToString[InputForm[temp5]] <>" #" <>"]^2/#, " <>ToString[InputForm[B2]] <>" BesselK[" <>ToString[InputForm[2 l+1]] <>"/2," <>ToString[InputForm[beta XC]] <>" #" <>"]^2/#]/" <>ToString[InputForm[A2]] <>"]&" ] ] ,Message[EnergyLevel::noroot]] ,Message[Permeability::error]] ,Message[EMass::error]] ]]; StructureOrder[{{R_?Positive,V_?Positive,min_?Positive,mout_?Positive,epin_?Positive,epout_?Positive},ES_List},opts___Rule]:=Module[{temp,l1=Length[ES],i,m,temp2,temp3}, If[Length[{opts}]>0, Message[StructureOrder::noop]; , If[l1 == 0, Message[StructureOrder::badin]; , If[EnergyStructureQ[{{R,V,min,mout,epin,epout},ES}], temp=ES; Do[ Do[ temp[[i,m]] = {temp[[i,m]], ToString[m+i-1] <> LName[[i]]}; ,{m,Length[ES[[i]]]}] ,{i,l1}]; {{R,V,min,mout,epin,epout},Transpose[Sort[Flatten[temp, 1], (#1[[1]] < #2[[1]])&]][[2]]} ,Message[EnergyStructureQ::error]] ] ] ]; Remove[data]; EnergyStructureQ[data_]:=Module[{ok=False}, If[ToString[Head[data]]=="List", If[Length[data]==2, If[VectorQ[data[[1]]]&&Length[data[[1]]]==6, If[MatchQ[data[[1]],{_?Positive,_?Positive,_?Positive,_?Positive,_?Positive,_?Positive}], If[And@@Map[NumericQ,Flatten[data[[2]]]], ok=True; ] ]] ]]; ok ]; End[]; Protect[EnergyLevel,EnergyStructure,StructurePlot, RadialWaveFunction, DensityFunction, PlotPDF, LevelLabel, StructureOrder, EnergyStructureQ]; EndPackage[]; On[General::spell]; On[General::spell1];