Matlab code

function [I] =  f3_hs_ff(q,scaling,F_scal,F_rad,F_sigma,S_phi)

Z=(1./F_sigma).^2-1;

alpha=(Z+1)./(q.*F_rad);

G1=alpha.^6-(1+4./alpha.^2).^(-(Z+1)./2).*alpha.^6.*cos((Z+1).*atan(2./alpha))+ ...

       (Z+2).*(Z+1).*(alpha.^(4)+(1+4./alpha.^2).^(-(Z+3)./2).*alpha.^4.*cos((Z+3).*atan(2./alpha)))- ...

       2.*(Z+1).*(1+4./alpha.^2).^(-(Z+2)./2).*alpha.^5.*sin((Z+2).*atan(2./alpha));

G2=sin((Z+1).*atan(1./alpha))-(Z+1).*(1+alpha.^2).^(-.5).*cos((Z+2).*atan(1./alpha));  

F=S_phi*scaling*F_scal.*(F_rad.^6.*(Z+1).^(-6).*G1);