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}}]
PsiR1[r_]:=myNorm myCoeff fOut[((Sqrt[2 mass Esep])/p),r,0]

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

(* строю фуункцию ферми *)
funcFermi[x_,par_]:= (1/(Exp[((par-x)/0.03)]+1));
Plot[funcFermi[x,2],{x,0,3},PlotRange->All]

(* режу ВФ в радиальном представлении с помощью функции ферми *)
PsiRCut[r_,cuval_]:= PsiR[r] funcFermi[r,cuval]
PsiR1Cut[r_,cuval_]:= PsiR1[r] funcFermi[r,cuval]
Plot[{PsiR[r],PsiR1[r],PsiRCut[r,7],PsiR1Cut[r,7]},{r,0,10},PlotRange->All]

(* Форм фактор протона в 3He. Должно быть ((Фм)^1.5).  ((Фм)^0.5)/(MeV^0.5) *)
(*PsiP3He[q_]:= (1/q)(Sqrt[(1/(p ))])myNorm(((Sin[((q/p)-((Sqrt[2 mass (-Esep+myU)])/p)) range])/(2 ((q/p)-((Sqrt[2 mass (-Esep+myU)])/p))))-((Sin[((q/p)+((Sqrt[2 mass (-Esep+myU)])/p)) range])/(2 ((q/p)+((Sqrt[2 mass (-Esep+myU)])/p))))+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))))*)
PsiP3He[q_]:= (1/q) myNorm(((((Sqrt[2 mass (-Esep+myU)])/p) Cos[((Sqrt[2 mass (-Esep+myU)])/p) range] Sin[q range]-q Cos[q range] Sin[((Sqrt[2 mass (-Esep+myU)])/p) range])/((((q))^2)-(((((Sqrt[2 mass (-Esep+myU)])/p)))^2)))+myCoeff(((E^(-((Sqrt[2 mass Esep])/p) range)) (q Cos[q range]+((Sqrt[2 mass Esep])/p) Sin[q range]))/((((q))^2)+(((((Sqrt[2 mass Esep])/p)))^2))))

(* Форм фактор из обрезанной ВФ *)
PsiP3HeCut[q_,cutVal_]:= (p/q) NIntegrate[PsiR1Cut[r,cutVal]Sin[((q r)/p)],{r,0,500}]
PsiP3HeCut[500,7]

NIntegrate[(((q PsiP3He[q]))^2),{q,0,Infinity},AccuracyGoal->5]

Plot[PsiP3HeCut[q,7],{q,0,2000},PlotRange->All]

data3He10 =Table[{q,PsiP3HeCut[q,10]},{q,1,1800}];
data3He12 =Table[{q,PsiP3HeCut[q,12]},{q,1,1800}];
data3He14 =Table[{q,PsiP3HeCut[q,14]},{q,1,1800}];
data3He18 =Table[{q,PsiP3HeCut[q,18]},{q,1,1800}];
data3He20 =Table[{q,PsiP3HeCut[q,20]},{q,1,1800}];
data3He22 =Table[{q,PsiP3HeCut[q,22]},{q,1,1800}];

ListPlot[{data3He10},PlotRange->All]

(* Data interpolation *)
factor3He10=Interpolation[data3He10,InterpolationOrder->2]
factor3He12=Interpolation[data3He12,InterpolationOrder->2]
factor3He14=Interpolation[data3He14,InterpolationOrder->2]
factor3He18=Interpolation[data3He18,InterpolationOrder->2]
factor3He22=Interpolation[data3He22,InterpolationOrder->2]
Plot[{factor3He10[q],PsiP3He[(q/p)]},{q,1,500},PlotRange->All]

/

NIntegrate[(((q PsiP3He[q]))^2),{q,0,Infinity},AccuracyGoal->5] (* нормировка на (Pi/2)*)

Plot[PsiP3He[(q/p)],{q,0,500},PlotRange->All]

(* численная оценка в нуле *)
myNorm(((((Sqrt[2 mass (-Esep+myU)])/p) Cos[((Sqrt[2 mass (-Esep+myU)])/p) range]range-Sin[((Sqrt[2 mass (-Esep+myU)])/p) range])/(-(((((Sqrt[2 mass (-Esep+myU)])/p)))^2)))+myCoeff(((E^(-((Sqrt[2 mass Esep])/p) range)) (1+((Sqrt[2 mass Esep])/p) range))/(((((Sqrt[2 mass Esep])/p)))^2)))

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

mom3He = Sqrt[2 mass (-Esep)]

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

Integrate[Sin[a x]Exp[- b x],{x,myRan,Infinity},Assumptions->a>0 && b>0] (* интеграл в общем виде *)

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

(*TMatrix3He[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))))) (* (Fm^2) , надо в (Fm^(3 ))MeV*)*)
TMatrix3He[k_]:=-((p myNorm myU)/(k Abs[mom3He])) (((((Sqrt[2 mass (-Esep+myU)])/p) Cos[((Sqrt[2 mass (-Esep+myU)])/p) range] Sin[(k/p) range]-(k/p) Cos[(k/p) range] Sin[((Sqrt[2 mass (-Esep+myU)])/p) range])/(((((k/p)))^2)-(((((Sqrt[2 mass (-Esep+myU)])/p)))^2)))) (* (Fm^2) , надо в (Fm^(3 ))MeV*)

(* численная оценка в нуле *)
-((myNorm myU)/(Abs[mom3He]))(((((Sqrt[2 mass (-Esep+myU)])/p)Cos[((Sqrt[2 mass (-Esep+myU)])/p)range]range-Sin[((Sqrt[2 mass (-Esep+myU)])/p)range])/(-(((((Sqrt[2 mass (-Esep+myU)])/p)))^2))))

Plot[TMatrix3He[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*)

PsiR8He[r_]:=Piecewise[{{myNorm1 fIn[((Sqrt[2 mass1 (-Esep1+myU1)])/p),r,0],r<range1},{myNorm1 myCoeff1 fOut[((Sqrt[2 mass1 Esep1])/p),r,0],r>range1}}]
PsiR8He1[r_]:=myNorm1 myCoeff1 fOut[((Sqrt[2 mass1 Esep1])/p),r,0]

Plot[{PsiR8He[r],PsiR8He1[r]},{r,0,10},PlotRange->All]

(* режу ВФ в радиальном представлении с помощью функции ферми *)
PsiR8HeCut[r_,cutVal_]:= PsiR8He[r] funcFermi[r,cutVal]
PsiR8HeCut1[r_,cutVal_]:= PsiR8He1[r] funcFermi[r,cutVal]
Plot[{PsiR8He[r],PsiR8He1[r],PsiR8HeCut[r,5],PsiR8HeCut1[r,5]},{r,0,15},PlotRange->All]

(* Форм фактор из обрезанной ВФ *)
PsiP8HeCut[q_,cutVal_]:= (p/q)NIntegrate[PsiR8HeCut1[r,cutVal]Sin[((q r)/p)],{r,0,500}]
PsiP8HeCut[500,5]

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

data8He10 =Table[{q,PsiP8HeCut[q,10]},{q,1,1000}];
data8He12 =Table[{q,PsiP8HeCut[q,12]},{q,1,1000}];
data8He14 =Table[{q,PsiP8HeCut[q,14]},{q,1,1000}];
data8He18 =Table[{q,PsiP8HeCut[q,18]},{q,1,1000}];
data8He22 =Table[{q,PsiP8HeCut[q,22]},{q,1,1000}];

PsiP8He[q_]:= (1/q)myNorm1(((((Sqrt[2 mass1 (-Esep1+myU1)])/p) Cos[((Sqrt[2 mass1 (-Esep1+myU1)])/p) range1] Sin[q range1]-q Cos[q range1] Sin[((Sqrt[2 mass1 (-Esep1+myU1)])/p) range1])/((((q))^2)-(((((Sqrt[2 mass1 (-Esep1+myU1)])/p)))^2)))+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))))

ListPlot[{data8He14},PlotRange->All]

(* Data interpolation *)
factor8He10=Interpolation[data8He10,InterpolationOrder->2]
factor8He12=Interpolation[data8He12,InterpolationOrder->2]
factor8He14=Interpolation[data8He14,InterpolationOrder->2]
factor8He18=Interpolation[data8He18,InterpolationOrder->2]
factor8He22=Interpolation[data8He22,InterpolationOrder->2]
Plot[{factor8He22[q]},{q,1,500},PlotRange->All]

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

Plot[PsiP8He[q],{q,0,5},PlotRange->All]

mom8He = Sqrt[2 mass1 (-Esep1)]

TMatrix8He[k_]:=-((p myNorm1 myU1)/(k Abs[mom8He])) (((Sin[((k/p)-((Sqrt[2 mass1 (-Esep1+myU1)])/p)) range1])/(2 ((k/p)-((Sqrt[2 mass1 (-Esep1+myU1)])/p))))-((Sin[((k/p)+((Sqrt[2 mass1 (-Esep1+myU1)])/p)) range1])/(2 ((k/p)+((Sqrt[2 mass1 (-Esep1+myU1)])/p)))))

Plot[TMatrix8He[k],{k,0,50},PlotRange->All]

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

TMat0[k_,q_]:= TMatrix3He[k]PsiP8He[(q/p)]

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

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));
Mbt = ((Mb Mt)/(Mb+Mt));
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[2 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_]:= (1/100)((Mstp Mbt)/(p^2)) ((2 K1[Et])/(3 K))  ((Abs[PsiP8He[((q1[Et,thetaCM])/p)]TMatrix3He[q2[Et,thetaCM]]+PsiP3He[((q2[Et,thetaCM])/p)]TMatrix8He[q1[Et,thetaCM]]])^2)

(*  обрезание 3He  *)
dSigmaHe3Cut10[Et_,thetaCM_]:= (1/100)((Mstp Mbt)/(p^2)) ((2 K1[Et])/(3 K))  ((Abs[factor3He10[((q2[Et,thetaCM])/p)]TMatrix8He[q1[Et,thetaCM]]])^2) 
dSigmaHe3Cut12[Et_,thetaCM_]:= (1/100)((Mstp Mbt)/(p^2)) ((2 K1[Et])/(3 K))  ((Abs[factor3He12[((q2[Et,thetaCM])/p)]TMatrix8He[q1[Et,thetaCM]]])^2) 
dSigmaHe3Cut14[Et_,thetaCM_]:= (1/100)((Mstp Mbt)/(p^2)) ((2 K1[Et])/(3 K))  ((Abs[factor3He14[((q2[Et,thetaCM])/p)]TMatrix8He[q1[Et,thetaCM]]])^2) 
dSigmaHe3Cut18[Et_,thetaCM_]:= (1/100)((Mstp Mbt)/(p^2)) ((2 K1[Et])/(3 K))  ((Abs[factor3He18[((q2[Et,thetaCM])/p)]TMatrix8He[q1[Et,thetaCM]]])^2) 
dSigmaHe3Cut22[Et_,thetaCM_]:= (1/100)((Mstp Mbt)/(p^2)) ((2 K1[Et])/(3 K))  ((Abs[factor3He22[((q2[Et,thetaCM])/p)]TMatrix8He[q1[Et,thetaCM]]])^2)

(*  обрезание 8He  *)
dSigmaHe8Cut10[Et_,thetaCM_]:= (1/100)((Mstp Mbt)/(p^2)) ((2 K1[Et])/(3 K))  ((Abs[factor8He10[((q1[Et,thetaCM])/p)]TMatrix3He[q2[Et,thetaCM]]])^2) 
dSigmaHe8Cut12[Et_,thetaCM_]:= (1/100)((Mstp Mbt)/(p^2)) ((2 K1[Et])/(3 K))  ((Abs[factor8He12[((q1[Et,thetaCM])/p)]TMatrix3He[q2[Et,thetaCM]]])^2) 
dSigmaHe8Cut14[Et_,thetaCM_]:= (1/100)((Mstp Mbt)/(p^2)) ((2 K1[Et])/(3 K))  ((Abs[factor8He14[((q1[Et,thetaCM])/p)]TMatrix3He[q2[Et,thetaCM]]])^2) 
dSigmaHe8Cut18[Et_,thetaCM_]:= (1/100)((Mstp Mbt)/(p^2)) ((2 K1[Et])/(3 K))  ((Abs[factor8He18[((q1[Et,thetaCM])/p)]TMatrix3He[q2[Et,thetaCM]]])^2) 
dSigmaHe8Cut22[Et_,thetaCM_]:= (1/100)((Mstp Mbt)/(p^2)) ((2 K1[Et])/(3 K))  ((Abs[factor8He22[((q1[Et,thetaCM])/p)]TMatrix3He[q2[Et,thetaCM]]])^2)

ContourPlot[dSigmaHe8Cut22[Et,(thetaCM/180)Pi],{Et,0,18},{thetaCM,0,50},Contours->20,PlotRange->All]
Plot[{dSigmaHe8Cut22[2.2,(thetaCM/180)Pi]},{thetaCM,0,35},PlotRange->All]

|

(* 1-Dit Cross section* )

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

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

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