% test_aquad.m  1/x  eps = 0.001
function test_aquad
  fid = fopen('test_aquad_m.out', 'w');
  eps=0.001;
  fprintf(fid, 'test_aquad.m running  eps=%g\n', eps);
  sprintf('generating test_aquad_m.out \n')
  xmin = 0.1;
  xmax = 2.0;
  q = quad(@f,xmin,xmax,eps);
  e = log(xmax)-log(xmin);
  fprintf(fid,'xmin=%g, mmax=%g, area=%g, exact=%g, err=%g \n', xmin, xmax, e, q, e-q);
 
  xmin = 0.01;
  xmax = 2.0;
  q = quad(@f,xmin,xmax,eps);
  e = log(xmax)-log(xmin);
  fprintf(fid,'xmin=%g, mmax=%g, area=%g, exact=%g, err=%g \n', xmin, xmax, e, q, e-q); 

  xmin = 0.001;
  xmax = 2.0;
  q = quad(@f,xmin,xmax,eps);
  e = log(xmax)-log(xmin);
  fprintf(fid,'xmin=%g, mmax=%g, area=%g, exact=%g, err=%g \n', xmin, xmax, e, q, e-q);

  xmin = 0.0001;
  xmax = 2.0;
  q = quad(@f,xmin,xmax,eps);
  e = log(xmax)-log(xmin);
  fprintf(fid,'xmin=%g, mmax=%g, area=%g, exact=%g, err=%g \n', xmin, xmax, e, q, e-q);
 
  xmin = 0.000001;
  xmax = 2.0;
  q = quad(@f,xmin,xmax,eps);
  e = log(xmax)-log(xmin);
  fprintf(fid,'xmin=%g, mmax=%g, area=%g, exact=%g, err=%g \n', xmin, xmax, e, q, e-q);

  fprintf(fid,'1.0/((x-0.3)*(x-0.3)+0.01) + 1.0/((x-0.9).*(x-0.9)+0.04) -6\n');
  xmin = 0.0;
  xmax = 1.0;
  q = quad(@f1,xmin,xmax,eps);
  e = 29.8583;
  fprintf(fid,'xmin=%g, mmax=%g, area=%g, exact=%g, err=%g \n', xmin, xmax, e, q, e-q);
  fprintf(fid,'test_aquad.m finished\n'); 
  return;

  function y=f1(x)
    y=1.0 ./ ((x-0.3) .*(x-0.3)+0.01) + 1.0 ./ ((x-0.9).*(x-0.9)+0.04) -6.0; 
    return
  end

  function y=f(x)
    y = 1.0 ./ x;
    return
  end
end