(*
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 at radius equal l;
Infinite barrier is set for ll>l in the external region.
a combination of two linear independent solutions
inside the dot (represented by spherical Bessel two-component
functions) are matched with four linearly independent Henkel
functions (decying asymptotically for r-> inf).

loop over dot radius; fixed width of the external barrier

[Note:
matching with two outside-Bessel functions vanishing at ll
is impossible]
*)

a=-0.5;
b=-2.0;
p=1.23008;
g=0.05 ;
vtt=2*g;
vmi=0.5*g;
dv=0.05*g;

vc=3.68*g;
vv=1.32*g;

vic=-0.25*g;
viv=-0.25*g;
lm=1.0;
(*
For[ii=1, ii<18, ii++,
*)
ii=14
(*
l=65-(ii-1)*3;
*)
l=55;
ll=l+22;     
(*
vc=vtt-vmi+2*dv;
vv=vmi-2*dv;
*)
tdx=0;
ex=0;
(*
For[i=1, i<51, i++,
*) 
(*
e=0.0 + g + (i-1)*0.02*vc;
*)
e=0.161635;

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+vc-e)+a*(e+vv)+p*p;
dx=Sqrt[Abs[bx*bx-4*a*b*(e+vv)*(g+vc-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];


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

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


(*
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);

w01=Exp[axx1*l]/(2*axx1*l);
w02=Exp[axx2*l]/(2*axx2*l);
w11=(Exp[axx1*l])*(1/axx1/l-1)/(2*I)/(axx1*l);
w12=(Exp[axx2*l])*(1/axx2/l-1)/(2*I)/(axx2*l);

(* w(ll) *)
c01=-Exp[-axx1*ll]/(2*axx1*ll);
c02=-Exp[-axx2*ll]/(2*axx2*ll);
c11=-(Exp[-axx1*ll])*(1+1/axx1/ll)/(2*I)/(axx1*ll);
c12=-(Exp[-axx2*ll])*(1+1/axx2/ll)/(2*I)/(axx2*ll);

u01=Exp[axx1*ll]/(2*axx1*ll);
u02=Exp[axx2*ll]/(2*axx2*ll);
u11=(Exp[axx1*ll])*(1/axx1/ll-1)/(2*I)/(axx1*ll);
u12=(Exp[axx2*ll])*(1/axx2/ll-1)/(2*I)/(axx2*ll);


a11=BesselJ[0,kk1*l];
a12=BesselJ[0,kk2*l];
a13=-h01;
a14=-h02;
a15=-w01;
a16=-w02;
a21=s1*BesselJ[1,kk1*l];
a22=s2*BesselJ[1,kk2*l];
a23=-z1*h11;
a24=-z2*h12;
a25=-z1*w11;
a26=-z2*w12;
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;
a35=-(-a*xx1-p*z1)*w11;
a36=-(-a*xx2-p*z2)*w12;
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);
a45=-(-2*b*z1*w11/l + (b*z1*xx1+p)*w01);
a46=-(-2*b*z2*w12/l + (b*z2*xx2+p)*w02);
a51=0;
a52=0;
a53=c01;
a54=c02;
a55=u01;
a56=u02;
a61=0;
a62=0;
a63=z1*c11;
a64=z2*c12;
a65=z1*u11;
a66=z2*u12;

tt={{Re[a11], Re[a12], Re[a13], Re[a14], Re[a15], Re[a16]}, {Re[a21], Re[a22], Re[a23], Re[a24], Re[a25], Re[a26]}, {Re[a31], Re[a32], Re[a33], Re[a34], Re[a35], Re[a36]}, {Re[a41], Re[a42], Re[a43], Re[a44], Re[a45], Re[a46]}, {Re[a51], Re[a52], Re[a53], Re[a54], Re[a55], Re[a56]}, {Re[a61], Re[a62], Re[a63], Re[a64], Re[a65], Re[a66]}};

td=N[Det[tt]/10^30];

(*
If[ (i>2) && (td*tdx<0), Print[ll,"  ",ex,"   ",e,"   ",((e+ex)/2)*13.6] ];
*)

tdx=td;
ex=e;

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

(*
]
*)
(*
]
*)

Print[e]
m21=Re[a21];
m22=Re[a22];
m23=Re[a23];
m24=Re[a24];
m25=Re[a25];
m26=Re[a26];
m31=Re[a31];
m32=Re[a32];
m33=Re[a33];
m34=Re[a34];
m35=Re[a35];
m36=Re[a36];
m41=Re[a41];
m42=Re[a42];
m43=Re[a43];
m44=Re[a44];
m45=Re[a45];
m46=Re[a46];
m51=Re[a51];
m52=Re[a52];
m53=Re[a53];
m54=Re[a54];
m55=Re[a55];
m56=Re[a56];
m61=Re[a61];
m62=Re[a62];
m63=Re[a63];
m64=Re[a64];
m65=Re[a65];
m66=Re[a66];
mm={{m22,m23,m24,m25,m26},{m32,m33,m34,m35,m36},{m42,m43,m44,m45,m46},{m52,m53,m54,m55,m56},{m62,m63,m64,m65,m66}};
bb={-m21, -m31, -m41, -m51, -m61};
LinearSolve[mm, bb]

{1.34349649471457*10^31, 8.98499261235014*10^56, 2.263283582908738*10^34, -0.0003701506357174367, 9.76562040915672*10^28}

Print[e]
ff=Re[BesselJ[0, kk1*x]]+(1.34349649471457*10^31)*Re[BesselJ[0, kk2*x]];
Plot[ff*10^(-30), {x, 0, l}]

Print[e]
ffo=(8.98499261235014*10^56)*Re[-Exp[-axx1*x]/(2*axx1*x)] + (2.263283582908738*10^34)*Re[-Exp[-axx2*x]/(2*axx2*x)] + (-0.0003701506357174367)*Re[Exp[axx1*x]/(2*axx1*x)] + (9.76562040915672*10^28)*Re[Exp[axx2*x]/(2*axx2*x)]
Plot[ffo*10^(-30), {x, l, ll}]
Print[l,"   ",ll]

Print[e]
{2.447106502258125*10^13, 1.02833091633915*10^26, -(3.34047019212494*10^15), -101.9115640250068, 
 
  4.327846021088476*10^12}