function [I] =  f3_css_ff(q,scaling,F_scal,F_cont_core,F_cont_shell,F_rad_core,F_sigma,F_thick_shell,S_phi)

 

Io=F_scal;

phi=S_phi;

z=1/F_sigma^2-1;

Rc=F_rad_core;

tsh=F_thick_shell;

rho1=F_cont_core;

rho2=F_cont_shell;

rho3=0;

x=q.*Rc;

y=q.*tsh;

gamma=(rho3-rho2)./(rho1-rho2);

B=(z+1).^2./((z+1).^2+4*x.^2);

D=atan(2.*x./(z+1));

c1=0.5-gamma.*(cos(y)+y.*sin(y))+0.5.*gamma.^2.*(1+y.^2);

c2=gamma.*y.*(gamma-cos(y));

c3=((gamma.^2+1)./2)-gamma.*cos(y);

c4=gamma.^2.*(y.*cos(y)-sin(y)).^2-c1;

c5=2.*gamma.*sin(y).*(1-gamma.*(y.*sin(y)+cos(y)))+c2;

c6=c3-(gamma.*sin(y)).^2;

c7=gamma.*sin(y)-(gamma.^2/2).*(1+y.^2).*sin(2.*y)-c5;

c8=c4-0.5+gamma.*cos(y)-(gamma.^2./2).*(1+y.^2).*cos(2.*y);

c9=gamma.*sin(y).*(1-gamma.*cos(y));

F=S_phi.*scaling.*Io.*((16.*3.1416.^2)./q.^6).*((rho2-rho1).^2).*(c1+c2.*x+c3.*x.^2.*((z+2)./(z+1))+B.^((z+1)/2).*(c4.*cos((z+1).*D) ...

   +c7.*sin((z+1).*D))+x.*B.^((z+2)./2).*(c5.*cos((z+2).*D)+c8.*sin((z+2).*D))+((z+2)./(z+1)).*x.^2.*B.^((z+3)./2).* ...

   (c6.*cos((z+3).*D)+c9.*sin((z+3).*D)));