Matlab code

function [I] =  f3_s_hs(q,scaling,F_rad,S_phi)

sigma=2*F_rad;

phi=S_phi;

A0=-((1+2.*phi).^2)./((1-phi).^4);

A1=(6.*phi.*(1+phi./2).^2)./((1-phi).^4);

A2=-(0.5.*phi.*(1+2.*phi).^2)./((1-phi).^4);

S1=24.*phi.*((-2.*A1./(q.*sigma).^4)+(24.*A2./(q.*sigma).^6)-(cos(q.*sigma)./(q.*sigma).^2).*(A0+A1+A2)+ ...

   (sin(q.*sigma)./(q.*sigma).^3).*(A0+2.*A1+4.*A2)+(cos(q.*sigma)./(q.*sigma).^4).*(2.*A1+12.*A2)-24.*A2.* ...

   (sin(q.*sigma)./(q.*sigma).^5)-24.*A2.*(cos(q.*sigma)./(q.*sigma).^6));

S=1./(1-S1);