Matlab code

function [I] = f3_ella_2y(q,scaling,F_scaling,F_radius1,F_radius2,S_phi)

nue=F_radius1/F_radius2;

FF1=(quadv(@(t) (3.*(sin((q.*F_radius2.*sqrt(1+(t).^2.*(nue.^2-1))))-(q.*F_radius2.*sqrt(1+(t).^2.*(nue.^2-1))).*cos((q.*F_radius2.*sqrt(1+(t).^2.*(nue.^2-1)))))./((q.*F_radius2.*sqrt(1+(t).^2*(nue.^2-1))).^3)).^2, 0, 1));

vol=((4/3)*pi*F_radius1*F_radius2^2) * (10^-6)^3;

I=F_scaling.*scaling.*S_phi.*vol.*FF1;