function [I,F,S] =  f3_hc_ff(q,scaling,F_scal,F_rad,F_sigma,F_height,S_phi)

global glob

glob.H=F_height;

glob.R=F_rad;

glob.sigma=F_sigma;

h = waitbar(0,'Please wait...');

for z=1:max(size(q))

   waitbar(z/max(size(q)))

   glob.q=q(z);    

   accuracy=1e-12;

   F(z)=0;

   numsteps=300;

   for alpha=pi/100:pi/numsteps:pi/2

       glob.alpha=alpha;

       FF= scaling.*F_scal.*S_phi.*quad(@integrandHCM2c,max([F_rad-3*F_sigma*F_rad 0]),F_rad+3*F_sigma*F_rad,accuracy);

       F(z)=F(z)+FF;

   end

   F(z)=F(z)/size([pi/100:pi/numsteps:pi/2],2);

end

close(h);

%-------------------------------------------------------------------

function P=integrandHCM2c(R)

   global glob

   q=glob.q;

   H=glob.H;

   sigma=glob.sigma;

   RM=glob.R;

   alpha=glob.alpha;

   Vc =pi*R.^2.0.*2.0.*H;

   P=weight(R,sigma,RM) .* ...

     ((2.0.*Vc.*sin(q.*H.*cos(alpha))./(q.*H.*cos(alpha)).*besselj(1,q.*R.*sin(alpha))./(q.*R.*sin(alpha))).^2.0.*sin(alpha)./Vc);

%--------------------------------------------------------------------