Mp=938.272; (*mass of proton, MeV*)
Mn=939.565;(*mass of neutron, MeV*)
Ebind=3.925 ;(*binding energy ((of )^8)He, MeV*)
Esep=24.81432; (*1p separation energy ((for )^8)He, MeV*)
p=197.327;
mass=821;
range=3.735;

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,50}]

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}]

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.2}]

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 для гелия8 == 1.688 фм *)

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

(*PsiP[q_]:=Simplify[(Sqrt[(2/π)])myNorm Integrate[ Sin[((Sqrt[2 mass (-Esep+myU)])/p) r]Sin[((q r)/p)],{r,0,range},Assumptions->q>0 && q∈Reals] + (Sqrt[(2/π)])myNorm myCoeff Integrate[ Exp[-((Sqrt[2 mass Esep])/p) r]Sin[((q r)/p)],{r,range,Infinity},Assumptions->q>0 && q∈Reals]]
Print[PsiP[q]]*)

Integrate[Sin[a x] Exp[-b x],{x,r,Infinity},Assumptions->a>0 && b>0]

Integrate[Sin[a x] Sin[b x],{x,0,r},Assumptions->a>0 && b>0]

PsiP[q_]:= (Sqrt[(2/(π p))])myNorm (((((Sqrt[2 mass (-Esep+myU)])/p) Cos[((Sqrt[2 mass (-Esep+myU)])/p) range] Sin[(q/p) range]-(q/p) Cos[(q/p) range] Sin[((Sqrt[2 mass (-Esep+myU)])/p) range])/(((((q/p)))^2)-(((((Sqrt[2 mass (-Esep+myU)])/p)))^2)))+myCoeff(((E^(-((Sqrt[2 mass Esep])/p) range)) ((q/p) Cos[(q/p) range]+((Sqrt[2 mass Esep])/p) Sin[(q/p) range]))/(((((q/p)))^2)+(((((Sqrt[2 mass Esep])/p)))^2))))

(* Проверим нормировку откуда взялся 1/p??*)
NIntegrate[(((PsiP[q]))^2),{q,0,Infinity},AccuracyGoal->4]