(* Input *)
(*kin = 20.;*)
mass= 939.;
range = 3;
t = 197.327;

Integrate[((Sin[x a])^2),x]
Integrate[Sin[x a]Sin[x b],x]
D[((Sin[a x])/(a x)),x]
D[((Cos[a x])/(a x)),x]

sps={{Dashing[None],RGBColor[0,0,0]}, {Dashing[Large],RGBColor[1,0,0]},{Dashing[Medium],RGBColor[0,0,1]},{Dashing[None],RGBColor[0,1,0]},{Dashing[None],RGBColor[0,1,1]}};

kOut[e_] := ((Sqrt[2*mass*e])/197.327)

kIn[e_,U0_] := ((Sqrt[2 mass (e- U0)])/197.327)

ang = 0;
Plot[{e*SphericalBesselJ[ang,e],Sin[e-ang*(Pi/2)],Sqrt[((Pi e)/2)]BesselJ[ang+0.5,e]},{e,0,50},PlotStyle->sps]

Plot[{-e*SphericalBesselY[ang,e],Cos[e-ang*(Pi/2)],-Sqrt[((Pi e)/2)]BesselY[ang+0.5,e]},{e,0,50},PlotStyle->sps]

(*δ[e_,ang_,r_,u0_]:= ArcCot[(SphericalBesselY[ang,kOut[e] r]D[SphericalBesselJ[ang,kIn[e,u0]r],r] - SphericalBesselJ[ang,kIn[e,u0]r]D[SphericalBesselY[ang,kOut[e] r],r])/(SphericalBesselJ[ang,kOut[e]r]D[SphericalBesselJ[ang,kIn[e,u0]r],r] - SphericalBesselJ[ang,kIn[e,u0]r]D[SphericalBesselJ[ang,kOut[e]r],r]) ]*)
(*δ[pIn_,pOut_,ang_,r_,u0_]:= ArcCot[(SphericalBesselY[ang,pOut r]D[SphericalBesselJ[ang,pIn r],r] - SphericalBesselJ[ang,pIn r]D[SphericalBesselY[ang,pOut r],r])/(SphericalBesselJ[ang,pOut r]D[SphericalBesselJ[ang,pIn r],r] - SphericalBesselJ[ang,pIn r]D[SphericalBesselJ[ang,pOut r],r]) ]*)
δ[q_,ang_,r_,u0_]:= ArcCot[(SphericalBesselY[ang,q r]D[SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r],r] - SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r]D[SphericalBesselY[ang,q r],r])/(SphericalBesselJ[ang,q r]D[SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r],r] - SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r]D[SphericalBesselJ[ang,q r],r]) ]

(*myDelta[e_,ang_,u0_]:= Simplify[δ[e,ang,r,u0]/.r->range]*)
myDelta[q_,ang_,u0_]:= Simplify[δ[q,ang,r,u0]/.r->range]

myU = -20;
myAng = 0;
(*Plot[{myDelta[kIn[e,-2],kOut[e],myAng,-2],myDelta[kIn[e,-5],kOut[e],myAng,-5],myDelta[kIn[e,-10],kOut[e],myAng,-10]},{e,0,40},PlotRange->All,PlotStyle->sps]*)
Plot[{myDelta[kOut[e],myAng,-2],myDelta[kOut[e],myAng,-5],myDelta[kOut[e],myAng,-10]},{e,0,40},PlotRange->All,PlotStyle->sps]

(*CS[pIn_,pOut_,ang_,u0_]:=(1/100)((4 π )/(pOut^2)) (2ang +1) (((Sin[myDelta[pIn,pOut,ang,u0]]))^2)*)
CS[q_,ang_,u0_]:=(1/100)((4 π )/(q^2)) (2ang +1) (((Sin[myDelta[q,ang,u0]]))^2)

csAng = 1;
(*Plot[{CS[kIn[ee,-5],kOut[ee],csAng,-5],CS[kIn[ee,-10],kOut[ee],csAng,-10],CS[kIn[ee,-15],kOut[ee],csAng,-15],CS[kIn[ee,-20],kOut[ee],csAng,-20]},{ee,0,40},PlotRange->All,PlotStyle->sps]*)
Plot[{CS[kOut[ee],csAng,-5],CS[kOut[ee],csAng,-10],CS[kOut[ee],csAng,-15],CS[kOut[ee],csAng,-20]},{ee,0,40},PlotRange->All,PlotStyle->sps]

(*A[e_,ang_,u0_]:=((kOut[e] SphericalBesselJ[ang,kOut[e]range]Cos[myDelta[e,ang,u0]] -kOut[e] SphericalBesselY[ang,kOut[e] range]Sin[myDelta[e,ang,u0]])/(kIn[e,u0]SphericalBesselJ[ang,kIn[e,u0]range]))*)
(*A[pIn_,pOut_,ang_,u0_]:=(1/(pIn SphericalBesselJ[ang,pIn range]))(pOut SphericalBesselJ[ang,pOut range]Cos[myDelta[pIn,pOut,ang,u0]] -pOut SphericalBesselY[ang,pOut range]Sin[myDelta[pIn,pOut,ang,u0]])*)
A[q_,ang_,u0_]:=(1/(Sqrt[(q^2)-((2mass u0)/(t^2))] SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] range]))(q SphericalBesselJ[ang,q range]Cos[myDelta[q,ang,u0]] -q SphericalBesselY[ang,q range]Sin[myDelta[q,ang,u0]])

(*fIn[e_,r_,ang_,u0_]:=A[e,ang,u0]kIn[e,u0]r SphericalBesselJ[ang,kIn[e,u0]r]*)
(*fIn[pIn_,pOut_,r_,ang_,u0_]:=A[pIn,pOut,ang,u0]pIn r SphericalBesselJ[ang,pIn r]*)
fIn[q_,r_,ang_,u0_]:=A[q,ang,u0]Sqrt[(q^2)-((2mass u0)/(t^2))] r SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r]
(*fIn[q_,r_,ang_,u0_]:=A[q,ang,u0]Sqrt[((q^2)-2 mass u0)] r SphericalBesselJ[ang,Sqrt[((q^2)-2 mass u0)] r]*)

(*fOut[e_,r_,ang_,u0_]:=kOut[e]r SphericalBesselJ[ang,kOut[e]r]Cos[myDelta[e,ang,u0]] - kOut[e]r SphericalBesselY[ang,kOut[e]r]Sin[myDelta[e,ang,u0]]*)
fOut[q_,r_,ang_,u0_]:=q r SphericalBesselJ[ang,q r]Cos[myDelta[q,ang,u0]] - q r SphericalBesselY[ang,q r]Sin[myDelta[q,ang,u0]]
(*fOut[q_,r_,ang_,u0_]:=q r SphericalBesselJ[ang,q r]Cos[myDelta[Sqrt[((q^2)-2 mass u0)],q,ang,u0]] - q r SphericalBesselY[ang,q r]Sin[myDelta[Sqrt[((q^2)-2 mass u0)],q,ang,u0]]*)

myE = 30;
myU = -20;
myAng = 1;
Plot[{fIn[kOut[myE],rr,myAng,myU],fOut[kOut[myE],rr,myAng,myU]},{rr,0.1,15},PlotRange->All,PlotStyle->sps]
(*Plot[{fIn[kOut[myE],rr,myAng,myU],fOut[kOut[myE],rr,myAng,myU]},{rr,0.1,15},PlotRange->All,PlotStyle->sps]*)

myR = 15;
myU = -15;
myAng = 1;
Plot[{fIn[kOut[ee],myR,myAng,myU],fOut[kOut[ee],myR,myAng,myU]},{ee,0.1,20},PlotRange->All,PlotStyle->sps]

(*fOut1[e_,r_,ang_]:=kOut[e] r  SphericalBesselJ[ang,kOut[e]r]*)

(* l == 0*)

δZero[q_,r_,u0_]:=ArcTan[ ((((Sin[q r])/(q r))D[((Sin[Sqrt[(q^2)-((2mass u0)/(t^2))] r])/(Sqrt[(q^2)-((2mass u0)/(t^2))] r)),r]-((Sin[Sqrt[(q^2)-((2mass u0)/(t^2))] r])/(Sqrt[(q^2)-((2mass u0)/(t^2))] r))D[Sin[q r],r])/(((Sin[Sqrt[(q^2)-((2mass u0)/(t^2))] r])/(Sqrt[(q^2)-((2mass u0)/(t^2))] r))D[((Cos[q r])/(q r)),r]- D[((Sin[Sqrt[(q^2)-((2mass u0)/(t^2))] r])/(Sqrt[(q^2)-((2mass u0)/(t^2))] r))]((Cos[q r])/(q r))))]
δZero1[q_,r_,u0_]:= ArcCot[(SphericalBesselY[ang,q r]D[SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r],r] - SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r]D[SphericalBesselY[ang,q r],r])/(SphericalBesselJ[ang,q r]D[SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r],r] - SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r]D[SphericalBesselJ[ang,q r],r]) ]
(*δZero1[q_,r_,u0_]:= ArcCot[(Cos D[SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r],r] - SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r]D[SphericalBesselY[ang,q r],r])/(SphericalBesselJ[ang,q r]D[SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r],r] - SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r]D[SphericalBesselJ[ang,q r],r]) ]*)

myDeltaZero[q_,u0_]:= Simplify[δZero[q,r,u0]/.r->range]

Azero[q_,u0_]:= ((Sin[q range]Cos[myDeltaZero[q,u0]]+Cos[q range] Sin[myDeltaZero[q,u0]])/(Sin[Sqrt[(q^2)-((2mass u0)/(t^2))]range]))

TMatriXZero[k_,q_,u0_,rr_]:= ((u0 Azero[q,u0])/197.327) Integrate[Sin[k rr]  Sin[q rr],rr]

myTMatriXZero[k_,q_,u0_]:= Simplify[TMatriXZero[k,q,u0,rr]/.rr->range]

TElasticZero[q_,u0_]:=myTMatriXZero[q,q,u0]

Plot[((mass^2)/(4 (Pi^2)(197.327^2)))((TElasticZero[(q/197.327),-5])^2),{q,0.01,150}]

TMatriX[k_,q_,r_,ang_,u0_]:= u0 Integrate[Sin[k r]A[q,ang,u0] Sin[q r],r]

Print[TMatriX[k,q,range,0,u]]

intT0[pIn_,pOut_,rr_,u0_]:=Integrate[Sin[pIn rr]Sin[pOut rr],rr]
Print[intT0[pIn,pOut,rr,u0]]

TElastic[q_,u_,r_]:= ((u 197.327)/(q^2))((r/2)-197.327((Sin[((2 q r)/197.327)])/(4 q)))

Plot[((mass^2)/(4 (Pi^2)(197.327^2)))((TElastic[q,-20,3])^2),{q,0.01,250}]

Clear[TMatriX0]
(* L=0 *)
TMatriX0[pIn_,pOut_,rr_,u_]:=((u A[pIn,pOut,0,u](197.327^2))/(2 pIn pOut))(((Sin[(((pIn-pOut) rr)/197.327)])/(pIn-pOut))-((Sin[(((pIn+pOut) rr)/197.327)])/(pIn+pOut)))
(*Print[Simplify[TMatriX0[pIn,pOut,rr,u0]]]*)

ContourPlot[((TMatriX0[p1,p2,3,20])^2),{p1,0.1,15},{p2,0.1,15},PlotPoints->20]

(* L=1 *)

TMatriX1[e_,aa_,bb_]:=Integrate[(((Sin[aa rr])/(aa rr)) - Cos[aa rr])(((Sin[bb rr])/(bb rr)) - Cos[bb rr]),rr]

Print[Simplify[TMatriX1[e,aa,bb]]]

tmpT1[e_,u_,aa_,bb_,rr_]:= ((u A[e,0,u])/(aa bb  ))(((aa^2) bb rr Cos[bb rr] Sin[aa rr]-(aa (bb^2) rr Cos[aa rr]+((aa^2)-(bb^2)) Sin[aa rr]) Sin[bb rr])/(bb ((aa^3)-aa (bb^2)) rr))

u0 = -20;
T1[e_]:= tmpT1[e,u0,kIn[e,u0],kOut[e],3]-tmpT1[e,u0,kIn[e,u0],kOut[e],0.0001]

Plot[((T1[e])^2),{e,0,50}]