%
%  code05.m -- oscillatory quadrature (djm - 23 sept 01)
%
%    - needs code05a.m & code05b.m for integrand
%

choice = 2;

switch (choice)
  case {1}
    R = 1.0;  integrand = 'code05a'; tol = [1e-3 1e-5];
    T = 6;  dt = T/6000;  tt = 0:dt:T;
  case {2}
    R = 1.0;  integrand = 'code05b'; tol = [1e-3 1e-5];
    T = 18;  dt = T/6000;  tt = 0:dt:T;
end

%  plot the integrand

clf;  subplot(2,1,1)

reR = feval(integrand,tt,R,0);
imR = feval(integrand,tt,R,1);

plot(tt,reR,'r');  hold on;
plot(tt,imR,'b');

title('\bf oscillatory integrand')
  
xlabel('\bf t-axis')
ylabel('\bf real (red) & imag (blue)')

%  single quadrature to x = 6

[reInt work1] = quadl(integrand,0,T,tol,0,R,0);
[imInt work2] = quadl(integrand,0,T,tol,0,R,1);

Int = reInt + i*imInt;

disp(' ')
disp(['quadrature = ' num2str(Int) '     (' num2str([work1 work2]) ')'])
disp(' ')

%  cumulative quadrature

ti = 0:T/60:T;  Qd = zeros(size(ti));

reInt = 0;  imInt = 0;
for j = 2:size(ti,2)
  reInt = reInt + quadl(integrand,ti(j-1),ti(j),tol,0,R,0);
  imInt = imInt + quadl(integrand,ti(j-1),ti(j),tol,0,R,1);

  Qd(j) = reInt + i*imInt;
end

subplot(2,1,2)
plot(ti,real(Qd),'r');  hold on
plot(ti,imag(Qd),'b');
plot(ti,real(Qd),'k.');
plot(ti,imag(Qd),'k.');

title('\bf cumulative quadratures')
xlabel('\bf T-axis')
ylabel('\bf integral 0->T')