(* Input *) 
kProton = (Sqrt[2*2.225*((939.565*938.272)/(938.272+939.565))]);
radCh = 2.1421;

(* Найдём радиус *) 
Solve[(radii^2)/4 + 0.64 -0.125==(radCh^2),radii,Assumptions->{radii>0}]

(* Используем такое представление волновой функции *)
fl1[k_,q_,r_] := Exp[-k*r] - Exp[-q*r]

(* Нормируем волновую функцию *)
Norm1[k_,q_] := Integrate[fl1[k,q,r]^2,{r,0,Infinity},Assumptions->{k>0,q>0}]

Simplify[Norm1[k,q]]

PsiD[k_,q_,r_]:= ((fl1[k,q,r])/(Sqrt[Norm1[k,q]]))

PsiD[k,q,r]

(* Проверим нормировку, везде где у нас импульсы домножаем на коэффициент 1/hc==1/197.327 для размерности *)
Integrate[((((((Sqrt[2]) ((E^(-k r))-(E^(-q r))))/(Sqrt[(1/k)+(1/q)-(4/(k+q))]))))^2),{r,0,Infinity},Assumptions->{k>0,q>0}]

(* Поиск q как среднее (r^2) *)
(*Plot[Integrate[((PsiD[(kProton/197.327),(q/197.327),rad])^2)*(rad^2),{rad,0,Infinity}]-(4.04^2),{q,225,230}]*)
Plot[Integrate[((((((Sqrt[2]) ((E^(-k r))-(E^(-((q r)/197.327)))))/(Sqrt[(1/k)+(197.327/q)-(4/(k+(q/197.327)))]))))^2)*(r^2)/.k-> (kProton/197.327) ,{r,0,Infinity}]-(4.04^2),{q,225,230}]

FindRoot[Integrate[((PsiD[(kProton/197.327),(q/197.327),r])^2)*(r^2),{r,0,Infinity},Assumptions->Re[q]>0]==(4.04^2),{q,228,229}]

(*Plot[PsiD[(kProton/197.327),(228.465/197.327),rad],{rad,0,20}]*)

Plot[(((Sqrt[2]) ((E^(-k r))-(E^(-q r))))/(Sqrt[(1/k)+(1/q)-(4/(k+q))]))/.k-> (kProton/197.327) /.q->(228.465/197.327) ,{r,0,30}]

(*Фурье-преобразование*)
(*PsiP[k_,q_,p_]:=Simplify[((Sqrt[4 π])/p)Integrate[(((Sqrt[2]) ((E^(-k r))-(E^(-q r))))/(Sqrt[(1/k)+(1/q)-(4/(k+q))]))*Sin[p r],{r,0,Infinity}],Assumptions->k>0 && q>0]*)
PsiP[k_,q_,p_]:=Simplify[(Sqrt[(2/π)])Integrate[(((Sqrt[2]) ((E^(-k r))-(E^(-q r))))/(Sqrt[(1/k)+(1/q)-(4/(k+q))]))*Sin[p r],{r,0,Infinity}],Assumptions->k>0 && q>0]
Print[PsiP[k,q,p]]

(* Проверим нормировку*)

(*Integrate[4 π (p^2) ((PsiP[k,q,p])^2) (1/(((2 π))^3)),{p,0,Infinity},Assumptions->k>0 && q>0]*)
Integrate[(((((2 p ((1/((k^2)+(p^2)))-(1/((p^2)+(q^2)))))/((Sqrt[π]) (Sqrt[(1/k)+(1/q)-(4/(k+q))])))))^2),{p,0,Infinity},Assumptions->k>0 && q>0]

(* Строим Фурье образ. Коэффициент 197.327 не нужен, потому что мы не переводим в Ферми, а работаем в импульсном представлении  *)
Plot[(((((2 p ((1/((k^2)+(p^2)))-(1/((p^2)+(q^2)))))/((Sqrt[π]) (Sqrt[(1/k)+(1/q)-(4/(k+q))])))))^2)/.k-> kProton /.q->228.465,{p,0,300}]

Plot[{(((2  ((1/((k^2)+(p^2)))-(1/((p^2)+(q^2)))))/((Sqrt[π]) (Sqrt[(1/k)+(1/q)-(4/(k+q))]))))/.k-> kProton /.q->228.465,(0.0047/((((p/44))^2)+1))},{p,0,200},PlotRange->All]

Clear[x]

Clear[y]

PsiP1[k_,q_,px_,py_]:=((2 ((1/((k^2)+(((px+py))^2)))-(1/((((px+py))^2)+(q^2)))))/((Sqrt[π]) (Sqrt[(1/k)+(1/q)-(4/(k+q))])))

Clear[k]

PsiP2[k_,q_,px_,py_]:=Simplify[(((((2 ((1/((k^2)+(((px+py))^2)))-(1/((((px+py))^2)+(q^2)))))/((Sqrt[π]) (Sqrt[(1/k)+(1/q)-(4/(k+q))])))))^2)]

Print[PsiP2[k,q,px,py]]

(* Параллельный импульс - интегрирование по (PsiP^2)*py*d(py) *)

PsiPar[k_,q_,px_]=Simplify[Integrate[((4 k q (((k+q))^3))/(π ((((k^2)+(((px+py))^2)))^2) ((((px^2)+2 px py+(py^2)+(q^2)))^2)))*py,{py,0,Infinity},Assumptions->k>0 && q>0]]

paral[px_]:=Simplify[PsiPar[kProton,228.465,px]]

Print[paral[px]]

Plot[PsiPar[k,q,px_]/.k-> kProton /.q->228.465,{px,0,50}]

Plot[paral[px],{px,0,100}]

(* Поперечный импульс - интегрирование по (PsiP^2)*d(px) *)

PsiPerp[k_,q_,py_]=Simplify[Integrate[((4 k q (((k+q))^3))/(π ((((k^2)+(((px+py))^2)))^2) ((((px^2)+2 px py+(py^2)+(q^2)))^2))),{px,0,Infinity},Assumptions->k>0 && q>0]]

perp[py_]:=Simplify[PsiPerp[kProton,228.465,py]]

Print[perp[py]]

Plot[perp[py]py,{py,0,300}]