(* 
Solution of the spherical 8-band PB Hamiltonian        
limited to |X>, |Y>, |Z> basis, spin-orbit splitting D=0
and to F=0.
Infinite barrier is set outside the dot region;
a combination of two linear independent solutions 
inside the dot (represented by spherical Bessel two-component
functions) are matched with two linearly independent Henkel
functions (decying asymptotically for r-> inf).
*)


a=-1.0;
b=-1.0;
p=1.23008;
g=0.1048; 
vc=0.0*g;
vv=1.0*g;


bb=b*(g+vc-e)+a*(e+vv)+p*p;
dd=Sqrt[Abs[bb*bb-4*a*b*(e+vv)*(g+vc-e)]];
k1=(-bb-dd)/(2*a*b);
k2=(-bb+dd)/(2*a*b);
kk1=Sqrt[k1];
kk2=Sqrt[k2];
akk1=Sqrt[Abs[k1]];
akk2=Sqrt[Abs[k2]];
ak1=Abs[k1];
ak2=Abs[k2];

bx=b*(g+v-e)+a*(e+v)+p*p;
dx=Sqrt[Abs[bx*bx-4*a*b*(e+v)*(g+v-e)]];
x1=(-bx-dx)/(2*a*b);
x2=(-bx+dx)/(2*a*b);
xx1=Sqrt[x1];
xx2=Sqrt[x2];
axx1=Sqrt[Abs[x1]];
axx2=Sqrt[Abs[x2]];
ax1=Abs[x1];
ax2=Abs[x2];

(*
Plot[k1, {e, g, 2.5*g}]
Plot[k2, {e, g, 2.5*g}]
Plot[y1, {e, 0, 2.5*g}]
Plot[y2, {e, 0, 2.5*g}]
*)

l=40;
ll=90;

s1=-p*kk1/(b*k1+e+vv);
s2=-p*kk2/(b*k2+e+vv);

z1=-p*xx1/(b*x1+e);
z2=-p*xx2/(b*x2+e);


(*
rr=BesselJ[0,I*Sqrt[Abs[k1]]*l]*s2*BesselJ[1,Sqrt[Abs[k2]]*l]-BesselJ[0,Sqrt[Abs[k2]]*l]*s1*I*BesselJ[1,I*Sqrt[Abs[k1]]*l]

Plot[rr, {e, g, 2.5*g}]

yy=BesselJ[0,I*Sqrt[Abs[k1]]*l]*s2*I*BesselJ[1,I*Sqrt[Abs[k2]]*l]-BesselJ[0,I*Sqrt[Abs[k2]]*l]*s1*I*BesselJ[1,I*Sqrt[Abs[k1]]*l]

Plot[yy, {e, 0.001, g-0.001}]

xx=BesselJ[0,I*Sqrt[Abs[x1]]*ll]*z2*I*BesselJ[1,I*Sqrt[Abs[x2]]*ll]-BesselJ[0,I*Sqrt[Abs[x2]]*ll]*z1*I*BesselJ[1,I*Sqrt[Abs[x1]]*ll]

Plot[xx, {e, 0.16, 0.25}]
*)


a11=BesselJ[0,kk1*l];
a12=BesselJ[0,kk2*l];
a21=s1*BesselJ[1,kk1*l];
a22=s2*BesselJ[1,kk2*l];



tt={{a11, a12}, {a21, a22}};
td=Re[Det[tt]];
Plot[td, {e, -0.1, 0.2}]
(*
FindRoot[td==0, {e, 0, g}]
*)
(*
c=-a11/a12
f=BesselJ[0, kk1*x]+c*BesselJ[0, kk2*x]
Plot[f, {x, 0, l}]
*)