function sigemadB=SigmaCalc(f0,u,q)

%主函數
g=9.8;
c=3*10^8;
% u=24*0.5144;              %由用戶給出
% f0=10*10^6;               %由用戶給出
% q=5*pi/6;                 %由用戶給出
% w=4;                      %由用戶給出
f0=f0*10^6;
k0=2*pi*f0/c;
wB=sqrt(2*g*k0);
eta=-4:0.01:4;
N=length(eta);
sigema=zeros(1,N);
for i=1:N
    if eta(i)==-1
        sigema(i)=sigema1L(q,u,f0)/16;
    elseif eta(i)==1
        sigema(i)=sigema1R(q,u,f0)/16;
    else
        [sitaL,sitaR]=limitfunc(eta(i));
        sigema(i)=Rombergfunc(sitaL,sitaR,eta(i),u,q,f0)/32;
    end
end  
for i=1:N
    if sigema(i)<0.0000000001
       sigema(i)=0.0000000001;
    end
end

 x=ifft(sigema);
 ww=gausswin(801,1);
 y=x.*ww.';
 s2=fft(y);

 for i=1:N
    if abs(s2(i))<0.0000000001
       s2(i)=0.0000000001;
    end
 end
 
 
 nn = 20;
 n1 = 400-nn;
 n2 = 400+nn;
 NWild = 0;
 for i=n1:n2
     if i>n1
         ss(1:i-n1) = abs(s2(n1:i-1));
         ss(i-n1+1:n2-n1) = abs(s2(i+1:n2));
         mm = mean(ss);
         if abs(s2(i))>4*mm
             NWild = NWild+1;
             wildVal(NWild) = i;
         end
         if abs(s2(i))<mm/4
             NWild = NWild+1;
             wildVal(NWild) = i;
         end
     end
 end
 
 sumS = sum(abs(s2(n1:n2)));
 for i=1:NWild
     sumS = sumS-abs(s2(wildVal(i)));
 end
 mm = sumS/(n2-n1+1-NWild);
 for i=1:NWild
     s2(wildVal(i)) = mm;
 end
 
 nn = length(s2);
 for i=1:nn
     if i<nn-2
         s2(i) = mean(abs(s2(i:i+2)));
     else
         s2(i) = mean(abs(s2(i-2:i)));
     end
 end
 
% sigemadB= 10*log10(abs([sigema.' ])); 
% s2dB=10*log10(abs([ s2.']));
% %   plot(eta(101:700),sigemadB(101:700));
% %   plot(eta(101:700),s2dB(101:700));
%   figure;plot(eta(101:700),[sigemadB(101:700) s2dB(101:700)]);
% xlabel('歸一化多普勒頻率');ylabel('海浪散射系數(dB)');
% grid

sigemadB=10*log10(abs([ s2.']));
sigemadB=sigemadB';
% figure;
% plot(eta,sigemadB);
% grid;
% if (w<-4*wB)|(w>4*wB)
%     ans=0;
% else 
%     tem=floor(100*w/wB)+400;
%     ans=sigema(tem);
% end