Matlab code

function [I] = f3_ulv_ff(q,scaling,F_scal,F_rad,F_epsilon1,F_cont1,F_sigma1, ...

                        F_epsilon2,F_cont2,F_sigma2,F_epsilon3,F_cont3,F_sigma3,S_phi)

eps = [F_epsilon1, F_epsilon2, F_epsilon3];

rho = [F_cont1, F_cont2, F_cont3];

sig = [F_sigma1, F_sigma2, F_sigma3];

I = 0;

for i=1:3,

   for j=1:3,

       I = I + ((R0 + eps(i)) * (R0 + eps(j)) * rho(i)*rho(j)*sig(i)*sig(j)).* ...

       (exp( (-q.^2).*(sig(i).^2 + sig(j).^2)/2 )).* ...

       cos( q.*(eps(i) - eps(j)) );

   end

end

I = F_scal.*scaling.*S_phi.*I.*q.^(-2);