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

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 = 1;
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]

δ[q_,ang_,r_,u0_]:= ArcTan[(D[SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r],r]SphericalBesselJ[ang,q r] - D[SphericalBesselJ[ang,q r],r]SphericalBesselJ[ang,Sqrt[(q^2)-((2mass u0)/(t^2))] r])/(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]) ]

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

myAng = 1;
Plot[{myDelta[(q/t),myAng,-16],myDelta[(q/t),myAng,-18],myDelta[(q/t),myAng,-22.737]},{q,0,200},PlotRange->All,PlotStyle->sps]

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

csAng = 1;
Plot[{CS[(q/t),csAng,-16],CS[(q/t),csAng,-18],CS[(q/t),csAng,-20]},{q,0,200},PlotRange->All,PlotStyle->sps]

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

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

myQ = 10;
(*myU = -22.736;*)
myU = -22.;
myAng = 1;
Plot[{fIn[(myQ/t),rr,myAng,myU],fOut[(myQ/t),rr,myAng,myU]},{rr,0.1,10},PlotStyle->sps]

myR = 3;
myU = -16;
myAng = 1;
Plot[{fIn[(q/t),myR,myAng,myU],fOut[(q/t),myR,myAng,myU]},{q,0.1,1000},PlotRange->All,PlotStyle->sps]

(*l=0*)
δ0[q_,r_,u0_]:= ArcTan[(D[((Sin[Sqrt[(q^2)-((2mass u0)/(t^2))]r])/(Sqrt[(q^2)-((2mass u0)/(t^2))]r)),r]((Sin[q r])/(q r)) - D[((Sin[q r])/(q r)),r]((Sin[Sqrt[(q^2)-((2mass u0)/(t^2))]r])/(Sqrt[(q^2)-((2mass u0)/(t^2))]r)))/(-((Cos[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[((Cos[q r])/(q r)),r]) ]

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

Plot[{myDelta0[((Sqrt[2*mass*e])/t),-2],myDelta0[((Sqrt[2*mass*e])/t),-4],myDelta0[((Sqrt[2*mass*e])/t),-5.6842]},{e,0,50},PlotRange->All,PlotStyle->sps]

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

Plot[{((A0[q,-5.6842])/q),((A[q,0,-5.6842])/q)},{q,0,0.000001},PlotStyle->sps]

TMatrix0[q_,u0_]:= ((u0 A0[q,u0])/(197.327 q q))(((Sin[(q-Sqrt[(q^2)-((2mass u0)/(t^2))]) range])/(2 (q-Sqrt[(q^2)-((2mass u0)/(t^2))])))-((Sin[(Sqrt[(q^2)-((2mass u0)/(t^2))]+q) range])/(2 (Sqrt[(q^2)-((2mass u0)/(t^2))]+q))))

CS0[q_,u0_]:=(1/100)((mass^2)/((197.327^2)  ))(4^2) Pi ((Abs[TMatrix0[q,u0]])^2)

Plot[{CS0[(q/t),-16.],CS[(q/t),0,-16.]},{q,0,5},PlotRange->All,PlotStyle->sps]

TMat0[q_,k_,u0_]:= ((u0 A0[q,u0])/(197.327 q q))(((Sin[(q-k) range])/(2 (q-k)))-((Sin[(k+q) range])/(2 (k+q))))

ContourPlot[TMat0[q,k,-5.6842],{q,0.1,20},{k,0.1,20}]

(* L=1 *)

δ1[q_,r_,u0_]:= -ArcTan[(D[((Sin[Sqrt[(q^2)-((2mass u0)/(t^2))]r])/(((q^2)-((2mass u0)/(t^2)))(r^2)))-((Cos[Sqrt[(q^2)-((2mass u0)/(t^2))]r])/(Sqrt[(q^2)-((2mass u0)/(t^2))]r)) ,r](((Sin[q r])/((q^2)(r^2)))-((Cos[q r])/(q r)) )-D[((Sin[q r])/((q^2)(r^2)))-((Cos[q r])/(q r)) ,r](((Sin[Sqrt[(q^2)-((2mass u0)/(t^2))]r])/(((q^2)-((2mass u0)/(t^2)))(r^2)))-((Cos[Sqrt[(q^2)-((2mass u0)/(t^2))]r])/(Sqrt[(q^2)-((2mass u0)/(t^2))]r))))/(D[((Sin[Sqrt[(q^2)-((2mass u0)/(t^2))]r])/(((q^2)-((2mass u0)/(t^2)))(r^2)))-((Cos[Sqrt[(q^2)-((2mass u0)/(t^2))]r])/(Sqrt[(q^2)-((2mass u0)/(t^2))]r)) ,r](((Cos[q r])/((q^2)(r^2)))+((Sin[q r])/(q r)))-D[((Cos[q r])/((q^2)(r^2)))+((Sin[q r])/(q r)),r](((Sin[Sqrt[(q^2)-((2mass u0)/(t^2))]r])/(((q^2)-((2mass u0)/(t^2)))(r^2)))-((Cos[Sqrt[(q^2)-((2mass u0)/(t^2))]r])/(Sqrt[(q^2)-((2mass u0)/(t^2))]r))))]

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

Plot[{myDelta1[(q/t),-16],myDelta1[(q/t),-18],myDelta1[(q/t),-22.737]},{q,0,200},PlotRange->All,PlotStyle->sps]

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

Plot[{A1[q,-20],A[q,1,-20]},{q,0,10},PlotStyle->sps,PlotRange->All]

Integrate[(((Sin[a r])/(a r))-Cos[a r])(((Sin[b r])/(b r))-Cos[b r]),r]

TMatrix1[q_,u0_]:= ((u0 A1[q,u0])/(197.327 q q))((1/(q (((Sqrt[(q^2)-((2mass u0)/(t^2))])^3)-Sqrt[(q^2)-((2mass u0)/(t^2))] (q^2)) range))(((q^2)-((2mass u0)/(t^2))) q range Cos[q range] Sin[Sqrt[(q^2)-((2mass u0)/(t^2))] range]-(Sqrt[(q^2)-((2mass u0)/(t^2))] (q^2) range Cos[Sqrt[(q^2)-((2mass u0)/(t^2))] range]+(((q^2)-((2mass u0)/(t^2)))-(q^2)) Sin[Sqrt[(q^2)-((2mass u0)/(t^2))] range]) Sin[q range]))

CS1[q_,u0_]:=(1/100)((mass^2)/((197.327^2)  ))(4^2) Pi ((Abs[TMatrix1[q,u0]])^2)3

Plot[{CS1[(q/t),-16.],CS[(q/t),1,-16.]},{q,0,200},PlotRange->All,PlotStyle->sps]

TMat1[q_,k_,u0_]:=((u0 A1[q,u0])/(197.327 q q))((1/(q ((k^3)-k (q^2)) range))((k^2) q range Cos[q range] Sin[k range]-(k (q^2) range Cos[k range]+((k^2)-(q^2)) Sin[k range]) Sin[q range]))

ContourPlot[TMat1[q,k,-20],{q,0.1,20},{k,0.1,20}]