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_]:=-((p myNorm myU)/(k Abs[mom3He])) (((Sin[((k/p)-((Sqrt[2 mass (-Esep+myU)])/p)) range])/(2 ((k/p)-((Sqrt[2 mass (-Esep+myU)])/p))))-((Sin[((k/p)+((Sqrt[2 mass (-Esep+myU)])/p)) range])/(2 ((k/p)+((Sqrt[2 mass (-Esep+myU)])/p)))))

Plot[TMatrix0[k],{k,0,50},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_]:= (1/q)(Sqrt[(1/(p ))])myNorm1 (((Sin[((q/p)-((Sqrt[2 mass1 (-Esep1+myU1)])/p)) range1])/(2 ((q/p)-((Sqrt[2 mass1 (-Esep1+myU1)])/p))))-((Sin[((q/p)+((Sqrt[2 mass1 (-Esep1+myU1)])/p)) range1])/(2 ((q/p)+((Sqrt[2 mass1 (-Esep1+myU1)])/p))))+myCoeff1(((E^(-((Sqrt[2 mass1 Esep1])/p) range1)) ((q/p) Cos[(q/p) range1]+((Sqrt[2 mass1 Esep1])/p) Sin[(q/p) range1]))/(((((q/p)))^2)+(((((Sqrt[2 mass1 Esep1])/p)))^2))))

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

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

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

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

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

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

Ms = 6567.18;
Mb = 7482.538;
Mp = 938.272;
Mt = 1875.612;
Mtp = ((Mt Mp)/(Mt + Mp));
Mstp = ((Ms*(Mt+Mp))/(Ms+Mt+Mp));
(*K = 357.789;*)
Kin = 1785.15; (* 8He lab momentum *)
K = (Mt/(Mt+Mb))Kin; (* 8He CM momentum *)
(*Q = -19.320;*)
(*Ek = 42.682;*)
Ek = ((K^2)/(2 ))*(((Mb+Mt))/(Mb * Mt));
K1[Et_]:=Sqrt[Mstp*(Ek+Q-Et)]
q1[Et_,thetaCM_]:=Sqrt[(((K1[Et] Sin[thetaCM]))^2)+(((K1[Et]Cos[thetaCM]-(Ms/Mb)K))^2)];
q2[Et_,thetaCM_]:=Sqrt[((((Mt/(Mt+Mp))K1[Et]Sin[thetaCM]))^2)+((((Mt/(Mt+Mp))K1[Et]Cos[thetaCM]-K))^2)];
Q = Mt+Mb - Ms - 2808.391;

dSigma[Et_,thetaCM_]:= Mstp Mp ((K1[Et])/K)  ((Abs[PsiP[q1[Et,thetaCM]]])^2)((Abs[TMatrix0[q2[Et,thetaCM]]])^2)

ContourPlot[dSigma[Et,(thetaCM/180)Pi],{Et,0,100},{thetaCM,0,45},PlotPoints->50,Contours->10,PlotRange->All]
(*Plot[dSigma[Et,(0/180)Pi],{Et,0,25}]*)

ContourPlot[q1[Et,(thetaCM/180)Pi],{Et,0,50},{thetaCM,0,180 },PlotPoints->50,PlotRange->All,ContourLabels->True]

ContourPlot[q2[Et,(thetaCM/180)Pi],{Et,0,50},{thetaCM,0,180},PlotPoints->50,Contours->10,PlotRange->All,ContourLabels->True]

Plot[{q1[1,thetaCM],q1[2,thetaCM],q1[10,thetaCM],q1[15,thetaCM]},{thetaCM,0,Pi 45/180},PlotRange->{0,350}]

Plot[{q2[1,thetaCM],q2[2,thetaCM],q2[10,thetaCM],q2[15,thetaCM],q2[25,thetaCM]},{thetaCM,0,Pi 45/180},PlotRange->{0,350}]

Plot[K1[Et],{Et,-5,50}]