dFdEgal.jl

function r(l::T, psi::T) where { T <: Number }
    return sqrt(RSOLAR^2 - 2*RSOLAR*l*cos(psi) + l^2)
end

function lmax(psi::T) where { T <: Number }
    return sqrt(RH^2 - RSOLAR^2*sin(psi)^2) + RSOLAR * cos(psi)
end

function rho_nfw(l::T, psi::T) where { T<: Number }
    return 0.4 * (RSOLAR/r(l,psi)) * ((1+RSOLAR/RS) / (1+r(l,psi)/RS))^2
end

function J(psi::T) where { T<: Number }
    return quadgk(l -> rho_nfw(l, psi), 0, lmax(psi), rtol=1e-10)[1]
end

dFdE_gal = nu_total[:,:dNdE] * FPBH/(MPBH*G2GEV) * 0.5 * quadgk(psi -> sin(psi) * J(psi), 0, pi,rtol=1e-10)[1] * KPC2CM