function [A]=Romberg(sitaL,sitaR,eta,u,q,f0)
%Romberg積分法
%sitaL,sitaR表示被積區間；
%A表示用Romberg加速算法所求得的積分值；
%eta,u,q,f0風別為歸一化多普勒頻率，風速，風向及雷達工作頻率。
 eps2=10^-7;
 h=sitaR-sitaL;
 R(1,1)=h*(integrandfunc(eta,sitaL,u,f0,q)+integrandfunc(eta,sitaR,u,f0,q))/2;
 M=1;
 J=0;
 err=1;
 while err>eps2
     J=J+1;
     h=h/2;
     S=0;
     for p=1:M
         sita=sitaL+h*(2*p-1);
         S=S+integrandfunc(eta,sita,u,f0,q);
     end
     R(J+1,1)=R(J,1)/2+h*S;
     M=M*2;
     for K=1:J
         R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);
     end
     err=abs(R(J+1,J)-R(J+1,J+1));
 end
 A=R(J+1,J+1);