Mp=938.272; (*mass of proton, MeV*)
Mn=939.565;(*mass of neutron, MeV*)
Ebind=7.718 ;(*binding energy ((of )^3)He, MeV*)
Esep=5.49351; (*1p separation energy ((for )^3)He, MeV*)
p=197.327;
mass=625.411;
range=1.6;

fIn[q_,r_,ang_]:=q r SphericalBesselJ[ang,q r]
fOut[k_,r_,ang_]:=((((2/Pi)k r))^(1/2)) BesselK[ang+(1/2),k r]
dfInR[q_,r_,ang_]:=D[fIn[q,r,ang],r]
dfIn[q_,ang_]:=Simplify[dfInR[q,r,ang]/.r->range]
dfOutR[k_,r_,ang_]:=D[((((2/Pi)k r))^(1/2)) BesselK[ang+(1/2),k r],r]
dfOut[q_,ang_]:=Simplify[dfOutR[q,r,ang]/.r->range]

Plot[((fIn[((Sqrt[2 mass (-Esep+U)])/p),range,0])/(dfIn[((Sqrt[2 mass (-Esep+U)])/p),0]))-((fOut[((Sqrt[2 mass Esep])/p),range,0])/(dfOut[((Sqrt[2 mass Esep])/p),0])),{U,0,100}]

myU=U/.FindRoot[((fIn[((Sqrt[2 mass (-Esep+U)])/p),range,0])/(dfIn[((Sqrt[2 mass (-Esep+U)])/p),0]))-((fOut[((Sqrt[2 mass Esep])/p),range,0])/(dfOut[((Sqrt[2 mass Esep])/p),0]))==0,{U,60}]

(*0.6401531831830206 2.5426173021971623 fOut[((Sqrt[2 mass Esep])/p),r,0]*)

myCoeff=coeff/.FindRoot[fIn[((Sqrt[2 mass (-Esep+myU)])/p),range,0]-coeff *fOut[((Sqrt[2 mass Esep])/p),range,0]==0,{coeff,1}]

myNorm=A/.FindRoot[(A^2) Integrate[((fIn[((Sqrt[2 mass (-Esep+myU)])/p),r,0])^2), {r,0,range}]+(((myCoeff*A))^2) Integrate[((fOut[((Sqrt[2 mass Esep])/p),r,0])^2), {r,range,Infinity}]==1, {A,0.5}]

Plot[Piecewise[{{myNorm fIn[((Sqrt[2 mass (-Esep+myU)])/p),r,0],r<range},{myNorm myCoeff fOut[((Sqrt[2 mass Esep])/p),r,0],r>range}}],{r,0,10}]

(*R для гелия3 == 1.561 фм *)

(2/3)Sqrt[Integrate[(myNorm^2)((fIn[((Sqrt[2 mass (-Esep+myU)])/p),r,0])^2)(r^2),{r,0,range}]+Integrate[(myNorm^2) (myCoeff^2) ((fOut[((Sqrt[2 mass Esep])/p),r,0])^2)(r^2),{r,range,Infinity}]]

PsiR[r_]:=Piecewise[{{myNorm fIn[((Sqrt[2 mass (-Esep+myU)])/p),r,0],r<range},{myNorm myCoeff fOut[((Sqrt[2 mass Esep])/p),r,0],r>range}}]

Plot[PsiR[r],{r,0,10}]

(* Импульс протона в 3He определяется энергией связи; то есть мы его знаем при постановке задачи *)

mom3He = Sqrt[2 mass (-Esep)]

(*М матрица для протона в прямоугольной яме 3He*)

Integrate[Sin[a x]Exp[- b x],x] (* интеграл в общем виде *)

Integrate[Sin[a x] Sin[b x],x]

TMatrix0[k_]:=-((myNorm myU)/(k mom3He)) (((Sin[(k-mom3He) range])/(2 (k-mom3He)))-((Sin[(k+mom3He) range])/(2 (k+mom3He))))

Plot[TMatrix0[k],{k,0,500},PlotRange->All]

(* Форм фактор протона в 8He *)

myNorm1 = 0.6514573663189586`;
myCoeff1 = 25.351749791847226`;
myU1=35.89034438287419`;
mass1 = 821;
range1=3.735;
Esep1=24.81432; (*1p separation energy ((for )^8)He, MeV*)
PsiP[q_]:= (Sqrt[(2/π)])myNorm1 (((Sin[(q-((Sqrt[2 mass1 (-Esep1+myU1)])/p)) range1])/(2 (q-((Sqrt[2 mass1 (-Esep1+myU1)])/p))))-((Sin[(q+((Sqrt[2 mass1 (-Esep1+myU1)])/p)) range1])/(2 (q+((Sqrt[2 mass1 (-Esep1+myU1)])/p))))+myCoeff1(((E^(-((Sqrt[2 mass1 Esep1])/p) range1)) (q Cos[q range1]+((Sqrt[2 mass1 Esep1])/p) Sin[q range1]))/((q^2)+(((((Sqrt[2 mass1 Esep1])/p)))^2))))

NIntegrate[(((PsiP[q]))^2),{q,0,Infinity}]

Plot[PsiP[q],{q,0,10},PlotRange->All]

(* Моя Т матрица *)

TMat0[k_,q_]:= TMatrix0[k]PsiP[q]

(* Строим что-то *)

ContourPlot[((Abs[TMat0[k,q]])^2),{k,0,100},{q,0,100}]