Matlab code

function [I] =  f3_s_sw(q,scaling,F_rad,S_phi,S_u,S_epsil)

phi=S_phi;

sigma=2*F_rad;

u=S_u;

epsilon=S_epsil;

x=q.*sigma;

delta=(1./(12.*epsilon)).*exp(-u)+phi./(1-phi);

lambda=6.*(delta-(delta.^2-phi).^0.5)./phi;

mu=lambda.*phi.*(1-phi);

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

beta=-(3.*phi.*(2+phi).^2-2.*mu.*(1+7.*phi+phi.^2)+(2+phi).*mu.^2)./(2.*(1-phi).^4);

C=-24.*phi.*(alph.*(sin(x)-x.*cos(x))./x.^3+beta.*(2.*x.*sin(x)-(x.^2-2).*cos(x) ...

   -2)./x.^4+(phi./2).*alph.*((4.*x.^3-24.*x).*sin(x)-(x.^4-12.*x.^2+24).*cos(x) ...

   +24)./x.^6)-(4.*(phi.*lambda.*epsilon).^2.*(sin(epsilon.*x)- ...

   epsilon.*x.*cos(epsilon*x))./(epsilon.*x).^3-(phi./2).*(2.*epsilon.*x.*sin(epsilon.*x) ...

   -((epsilon.*x).^2-2).*cos(epsilon.*x)-2)./(epsilon.*x).^4)-2.*phi.^2.*lambda.^2.*((1- ...

   cos(x))./x.^2-epsilon.^2.*(1-cos(epsilon.*x))./(epsilon.*x).^2)+(2.*phi.*lambda./epsilon).*((1- ...

   cos(x))./x.^2-(1-epsilon).^2.*(1-cos((1-epsilon).*x))./((1-epsilon).*x).^2)+ ...

   24.*phi.*((sin(x)-x.*cos(x))./x.^3-(1-epsilon).^3.*(sin((1-epsilon).*x)- ...

   (1-epsilon).*x.*cos((1-epsilon).*x))./((1-epsilon).*x).^3);

S=1./(1-C);