FFTLog — See MouseOver

Spherical Bessel/Hankel Transforms (∫_0^∞)(Subscript[j, l])(k y)f(y)((dy)/y) over logarithmically spaced array of k

DoTest=True

(* trial function that has analytic Bessel/Hankel transform *)
Clear[z,l,k,trialfunc,trialfuncSBT];
If[DoTest,
trialfunc[y_]=(y^z)(E^(-(y^2)/2));
(*trialfuncSBT[l_,k_,z_]=Integrate[SphericalBesselJ[l,k y]((trialfunc[y])/y),{y,0,∞},Assumptions->k>0&&l+z>0];
Simplify[trialfuncSBT[l,k,z]]*)
trialfuncSBT[l_,k_,z_]=(2^((1/2) (-4-l+z))) (k^l) (Sqrt[π]) Gamma[((l+z)/2)] Hypergeometric1F1Regularized[((l+z)/2),(3/2)+l,-((k^2)/2)]
]

(* define ln[(Subscript[U, l])(z)], (Subscript[U, l])(z) ≡ (∫_0^∞)(Subscript[j, l])(y)(y^z)((\[DifferentialD]y)/y) *)
Clear[lnUl];
lnUl[l_,z_]=(1/2)Log[π]+(z-2) Log[2]+(Piecewise[{{0, l+z==0}, {LogGamma[(1/2)(l+z)], True}}])-LogGamma[(1/2)(l-z+3)];
Clear[lnUlplain];
lnUlplain[l_,z_]=(1/2)Log[π]+(z-2) Log[2]+LogGamma[(1/2)(l+z)]-LogGamma[(1/2)(l-z+3)];

(* for real expr, Mathematica's Fourier[expr] has second half of array complex conjugate of first half; clip to the first half *)
Clear[FourierPars,ClipFourier];
FourierPars={-1,-1};
ClipFourier[expr__]:=Fourier[expr,FourierParameters->FourierPars][[1;;1+IntegerPart[((Length[expr])/2)]]]

(* ((∫_0)^∞)f(y)(Subscript[j, l])(k y)(dy/y) over logarithmically spaced table of k of width lnymax-lnymin centered at ln((Subscript[k, c])) *)
Clear[SphericalBesselTransform];
Clear[DoPrint,lnyc,Δlny,lnyTab,funcTab,ClipFourierfuncTab,lnkcyc,lnUlTab,ClipFourierfuncUlTab,FourierfuncUlTab,funcSBTTab];
SphericalBesselTransform[func_,ny_Integer,lnymin_?NumericQ,lnymax_?NumericQ,q_?NumericQ,l_?NumericQ]:=(*SphericalBesselTransform[func,ny,lnymin,lnymax,q,l]=*)(
DoPrint=False;
lnyc=N[(1/2)(lnymin+lnymax)];
Δlny=N[lnymax-lnymin];
lnyTab=Table[lnyc+(iy/ny)Δlny,{iy,-IntegerPart[(ny/2)],IntegerPart[((ny-1)/2)]}];
If[DoPrint,Print["lnyTab = ",Short[lnyTab,3]]];
funcTab=Table[(E^(-q lny))func[(E^lny)],{lny,lnyTab}];
If[DoPrint,Print["funcTab = ",Short[funcTab,3]]];
ClipFourierfuncTab=ClipFourier[funcTab];
If[DoPrint,Print["ClipFourierfuncTab = ",Short[ClipFourierfuncTab,3]]];
(* ln((Subscript[k, c])(Subscript[y, c])) where (Subscript[k, c]) and (Subscript[y, c]) are centres of logarithmic intervals *)
lnkcyc=N[((Δlny)/ny)Mod[(1/π)Im[lnUl[l,q+((π I ny)/(Δlny))]],1,-(1/2)]];
If[DoPrint,Print["ln((Subscript[k, c])(Subscript[y, c])) = ",lnkcyc]];
lnUlTab=Table[-((2π I iy)/(Δlny))lnkcyc+lnUl[l,q+((2π I iy)/(Δlny))],{iy,0,(ny/2)}];
If[DoPrint,Print["lnUlTab = ",Short[lnUlTab,3]]];
ClipFourierfuncUlTab=ClipFourierfuncTab Exp[lnUlTab];
If[DoPrint,Print["ClipFourierfuncUlTab = ",Short[ClipFourierfuncUlTab,3]]];
FourierfuncUlTab=Flatten[(Piecewise[{{{ClipFourierfuncUlTab[[;;-2]],Re[ClipFourierfuncUlTab[[-1]]],Conjugate[Reverse[ClipFourierfuncUlTab[[2;;-2]]]]}, EvenQ[ny]}, {{ClipFourierfuncUlTab,Conjugate[Reverse[ClipFourierfuncUlTab[[2;;]]]]}, True}}])];
funcSBTTab=Exp[-q(lnkcyc-lnyTab)]Re[InverseFourier[FourierfuncUlTab,FourierParameters->FourierPars]];
If[DoPrint,Print["funcSBTTab = ",Short[funcSBTTab,3]]];
Table[{lnkcyc-lnyTab[[i]],funcSBTTab[[i]]},{i,1,ny}]
)

(* plot Spherical Bessel Transform of trial function *)
Clear[fSBT];
If[DoTest,
Block[{l=5,z=l+3},
fSBT=SphericalBesselTransform[trialfunc,1024,-10,8,0,l];
ListLogPlot[{Table[{lnk,trialfuncSBT[l,(E^lnk),z]},{lnk,lnyTab}],fSBT},Joined->True,PlotRange->{(10^(-20)),All}]]
]