(*Задание Хюльтеновского анзаца, энергий связи и масс частиц*) 
m1=938.2723;
m2=939.5656;
Eb1=2.224;
Eb2=1.296;
Eb3=7.77;
m=((m1 m2)/(m1+m2));
j=197.327;

Psi[r_, q_, k_]:=n*((E^(-k r))-(E^(-q r)))

(*Поиск нормировочного коэффициента, подстановка в ВФ*)

A=Integrate[((Psi[r, q, k])^2),{r, 0, Infinity},Assumptions->Re[k]>0&&Re[q]>0];

K=Simplify[Solve[A==1, n]];

K[[2]][[1]]

Psi1[r_,q_,k_]=Psi[r,q,k]/.K[[2]][[1]];

Psi1[r,q,k]

(*Задание численного значения параметра q для Eb=2.224 МэВ. Поиск точки, в окрестности которой будем искать решение уравнения относительно q*)
Plot[Integrate[((Psi1[r,(q/j),((Sqrt[2 m Eb1])/j)])^2)*(r^2),{r,0,Infinity}]-16,{q,238,238.4}]

Q1=FindRoot[Integrate[((Psi1[r,(q/j),((Sqrt[2 m Eb1])/j)])^2)*(r^2),{r,0,Infinity},Assumptions->Re[q]>0]==16,{q,238.2}]

(*Задание численного значения параметра q для Eb=1.296 МэВ. Поиск точки, в окрестности которой будем искать решение уравнения относительно q*)

Plot[Integrate[((Psi1[r,(q/j),((Sqrt[2 m Eb2])/j)])^2)*(r^2),{r,0,Infinity}]-16,{q,4613250, 4613300}]

Q2=FindRoot[Integrate[((Psi1[r,(q/j),((Sqrt[2 m Eb2])/j)])^2)*(r^2),{r,0,Infinity},Assumptions->Re[q]>0]==16,{q,4613280}]

(*Задание численного значения параметра q для Eb=77.7 МэВ. Поиск точки, в окрестности которой будем искать решение уравнения относительно q*)

Plot[Integrate[((Psi1[r,(q/j),((Sqrt[2 m Eb3])/j)])^2)*(r^2),{r,0,Infinity}]-16,{q,80, 100}]

Q3=FindRoot[Integrate[((Psi1[r,(q/j),((Sqrt[2 m Eb3])/j)])^2)*(r^2),{r,0,Infinity},Assumptions->Re[q]>0]==16,{q,85}]

(*Графики полученных ВФ для энергий связи Eb1=2.224 МэВ, Eb2=1.296 МэВ, Eb3=7.77 МэВ*)

V=Psi1[r,(q/j),((Sqrt[2 m Eb1])/j)]/.Q1[[1]];

L=Psi1[r,(q/j),((Sqrt[2 m Eb2])/j)]/.Q2[[1]];

S=Psi1[r,(q/j),((Sqrt[2 m Eb3])/j)]/.Q3[[1]];

Plot[{L,V,S},{r,0,30},PlotLegends->Automatic]
LogPlot[{L,V,S},{r,0,100},PlotLegends->Automatic]

(Legended[, (LineLegend[1,296МэВ, 2,224  МэВ, 7,7 МэВ])])

(Legended[, (LineLegend[1,296МэВ, 2,224  МэВ, 7,7 МэВ])])

Psi1[r,q,k]

(*Фурье-преобразование Хюльтеновского анзаца*)
Psi2[k_,q_,p_]:=Simplify[((Sqrt[4 π])/p)Integrate[Psi1[r,q,k] Sin[p r],{r,0,Infinity}],Assumptions->Re[k]>Im[p]&&Re[q]>Im[p]]
Print[Psi2[k,q,p]]

(*Проверка нормировки получившейся ВФ*)

Simplify[Integrate[4 π (p^2) ((Psi2[k,q,p])^2) (1/(((2 π))^3)),{p,0,Infinity}],Assumptions->Re[k]>0&&Re[q]>0]

(*Плотность распределения в зависимости от модуля импульса для энергий связи Eb1=2.224 МэВ, Eb2=1.296 МэВ, Eb3=7.77 МэВ*)

M[k_,q_,p_]=((((1/(Sqrt[4 π])) Psi2[k,q,p] p))^2);

V1=M[(Sqrt[2 m Eb1]),q,p]/.Q1[[1]];

L1=M[(Sqrt[2 m Eb2]),q,p]/.Q2[[1]];

S1=M[(Sqrt[2 m Eb3]),q,p]/.Q3[[1]];

Plot[{L1,V1,S1},{p,0,500},PlotLegends->Automatic]

(Legended[, (LineLegend[1,296 МэВ, 2,224  МэВ, 7,7 МэВ])])

(*Задание распределений по перпендикулярным и параллельным импульсам, проверка нормировки функции, где х - параллельный импульс, а у - перпендикулярный импульс*)

Integrate[Integrate[((Psi2[k,q,(Sqrt[(x^2)+(y^2)])])^2) (1/(((2 π))^3)) 2 π y,{y,0, Infinity}, Assumptions->(Sqrt[-(k^2)-(x^2)])≠0&&(Sqrt[-(q^2)-(x^2)])≠0&&((Re[(Sqrt[-(q^2)-(x^2)])]≤0&&Re[(q^2)+(x^2)]>0)||(Sqrt[-(q^2)-(x^2)])∉(Reals))&&((Re[(k^2)+(x^2)]>0&&Re[(Sqrt[-(k^2)-(x^2)])]≤0)||(Sqrt[-(k^2)-(x^2)])∉(Reals))],{x,-Infinity, Infinity}, Assumptions->(Re[(k^2)]≥0||(k^2)∉(Reals))&&(Re[(q^2)]≥0||(q^2)∉(Reals))&&Re[k]>0&&Re[q]>0]

P=(1/π)((((1/(Sqrt[4 π])) Psi2[k,q,(Sqrt[(x^2)+(y^2)])]))^2) y;

(*Построение распределения по параллельным импульсам*)

P1[k_,q_,x_]=Integrate[P,{y,0,Infinity},Assumptions->(Sqrt[-(k^2)-(x^2)])≠0&&(Sqrt[-(q^2)-(x^2)])≠0&&((Re[(Sqrt[-(q^2)-(x^2)])]≤0&&Re[(q^2)+(x^2)]>0)||(Sqrt[-(q^2)-(x^2)])∉(Reals))&&((Re[(k^2)+(x^2)]>0&&Re[(Sqrt[-(k^2)-(x^2)])]≤0)||(Sqrt[-(k^2)-(x^2)])∉(Reals))];

Print[P1[k,q,x]]

V2=P1[(Sqrt[2 m Eb1]),q,x]/.Q1[[1]];

L2=P1[(Sqrt[2 m Eb2]),q,x]/.Q2[[1]];

S2=P1[(Sqrt[2 m Eb3]),q,x]/.Q3[[1]];

Plot[{L2,V2,S2},{x,0,50},PlotLegends->Automatic]

(Legended[, (LineLegend[1,296 МэВ, 2,224 МэВ, 7,7 МэВ])])

(*Построение распределения по перпендикулярным импульсам*)

P2[k_,q_,y_]=Integrate[P,{x,0,Infinity},Assumptions->Im[(Sqrt[-(k^2)-(y^2)])]>0&&Im[(Sqrt[-(q^2)-(y^2)])]>0];

V3=P2[(Sqrt[2 m Eb1]),q,y]/.Q1[[1]];

L3=P2[(Sqrt[2 m Eb2]),q,y]/.Q2[[1]];

S3=P2[(Sqrt[2 m Eb3]),q,y]/.Q3[[1]];

Plot[{L3,V3,S3},{y,0,300},PlotLegends->Automatic]

(Legended[, (LineLegend[1,296  МэВ, 2,224  МэВ, 7,7 МэВ])])