(* 
Solution of the spherical 8-band PB Hamiltonian        
limited to |X>, |Y>, |Z> basis, spin-orbit splitting D=0
and to F=0.
Finit barrier V 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=-0.8;
b=-1.7;
p=1.23008;
g=0.1048; 
v=1.5*g;

For[i=1, i<50, i++,
e=0.0+i*0.02*(g+v);
bb=b*(g-e)+a*e+p*p;
dd=Sqrt[Abs[bb*bb-4*a*b*e*(g-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);
s2=-p*kk2/(b*k2+e);

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


(*
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}]
*)


h01=-Exp[-axx1*l]/(2*axx1*l);
h02=-Exp[-axx2*l]/(2*axx2*l);
h11=-(Exp[-axx1*l])*(1+1/axx1/l)/(2*I)/(axx1*l);
h12=-(Exp[-axx2*l])*(1+1/axx2/l)/(2*I)/(axx2*l);


a11=BesselJ[0,kk1*l];
a12=BesselJ[0,kk2*l];
a13=-h01;
a14=-h02;
a21=s1*BesselJ[1,kk1*l];
a22=s2*BesselJ[1,kk2*l];
a23=-z1*h11;
a24=-z2*h12;
a31=(-a*kk1-p*s1)*BesselJ[1,kk1*l];
a32=(-a*kk2-p*s2)*BesselJ[1,kk2*l];
a33=-(-a*xx1-p*z1)*h11;
a34=-(-a*xx2-p*z2)*h12;
a41=-2*b*s1*BesselJ[1,kk1*l]/l + (b*s1*kk1+p)*BesselJ[0,kk1*l];
a42=-2*b*s2*BesselJ[1,kk2*l]/l + (b*s2*kk2+p)*BesselJ[0,kk2*l];
a43=-(-2*b*z1*h11/l + (b*z1*xx1+p)*h01);
a44=-(-2*b*z2*h12/l + (b*z2*xx2+p)*h02);


tt={{Re[a11], Re[a12], Re[a13], Re[a14]}, {Re[a21], Re[a22], Re[a23], Re[a24]}, {Re[a31], Re[a32], Re[a33], Re[a34]}, {Re[a41], Re[a42], Re[a43], Re[a44]}};


td=N[Det[tt]*10^5];

Print[e, "   ",td, "   ", 0] ]