function I = filon_cos(w,f,t,N)

wmin = w(1);
wmax = w(end);

N = N/2;

dw = (w(end)-w(1))/(2*N);
ww = wmin:dw:wmax;
ff = spline(w,f,ww);

%figure
%plot(ww,ff,w,f,'x');

I = zeros(length(t),1);
for i = 1:length(t)

C2n_1 = 0; S2n_1p = 0;
C2n = 0.5*ff(1)*cos(t(i)*ww(1))-0.5*ff(2*N+1)*cos(t(i)*ww(2*N+1)); 
for j = 1:N
  C2n_1 = C2n_1+ff(2*j)*cos(t(i)*ww(2*j));
  S2n_1p = S2n_1p+ff(2*j)*sin(t(i)*ww(2*j));  
  C2n = C2n+ff(2*j+1)*cos(t(i)*ww(2*j+1));
end
     
theta = t(i)*dw;
[alpha,beta,gamma] = filon_coeffs(theta);


I(i) = alpha*(-ff(1)*sin(t(i)*ww(1))+ff(2*N+1)*sin(t(i)*ww(2*N+1))) + beta*C2n + gamma*C2n_1 + 2*t(i)*dw^4*S2n_1p/45;
end

I = I*dw;




