(* 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 = 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]

δ[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 = 0;
Plot[{myDelta[((Sqrt[2*mass*e])/t),myAng,-2],myDelta[((Sqrt[2*mass*e])/t),myAng,-4],myDelta[((Sqrt[2*mass*e])/t),myAng,-6]},{e,0,20},PlotRange->All,PlotStyle->sps]

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

csAng = 0;
Plot[{CS[(q/t),csAng,-5],CS[(q/t),csAng,-10],CS[(q/t),csAng,-15],CS[(q/t),csAng,-20]},{q,0,50},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 = 200;
myU = -4;
myAng = 1;
Plot[{fIn[(myQ/t),rr,myAng,myU],fOut[(myQ/t),rr,myAng,myU]},{rr,0.1,15},PlotRange->All,PlotStyle->sps]

myR = 4;
myU = -4;
myAng = 0;
Plot[{fIn[(q/t),myR,myAng,myU],fOut[(q/t),myR,myAng,myU]},{q,0.1,500},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),-6]},{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,-4],A[q,0,-4]},{q,0,4},PlotStyle->sps]

TMatrix0[q_,u0_]:= ((u0 A0[q,u0])/(197.327 q))((range/(2 ))-((Sin[2q range])/(2 q)))

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

crossSection[q_,thetaCos_]:=0.01* (Pi/(q^2))(((Abs[Sum[((2 l+1)/(Sqrt[2]))*(Exp[2*I*myDelta0[q,u0]] -1)*LegendreP[l,thetaCos],{l,0,lmax}]]))^2)

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

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)/(2(q^2))) (1-197.327 ((Sin[((2q r)/197.327)])/(2q )))

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

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

ContourPlot[((TMatriX0[p1,p2])^2),{p1,0.1,10},{p2,0.1,10},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}]