function [f]=integrand(eta,sita,u,f0,q)
  g=9.8;
  c=3*10^8;
  k0=2*pi*f0/c;
  [m1,m2]=mfunc(eta);
  y=calculateyfunc(eta,sita,m1,m2);
  f1=P1func(y,sita,m1,m2);
  f2=S1func(y,sita,m1,k0,u,q);
  f3=S2func(y,sita,m2,k0,u,q);
  f4=TaoHfunc(y,sita,eta,m1,m2);
  f5=TaoEMfunc(y,sita);
  r=f4+f5;
  r2=(abs(r))^2;
  f=16*pi*f1*f2*f3*y^3*r2;