% gauleg.m     P145 Numerical Recipes in Fortran               
% compute x(i) and w(i)  i=1,n  Legendre ordinates and weights 
% on interval  -1.0 to 1.0 (length is 2.0)                     
% use ordinates and weights for Gauss Legendre integration     
%  needs  gaulegf.m  available

function gauleg
  clear
  format compact
  diary 'gauleg_m.out'
  disp('test gauleg.m on interval -1.0 to 1.0  ordinates, weights')
  for i1=1:15
    [x1 w1 ] = gaulegf(-1.0, 1.0, i1);
    sum = 0.0;
    for j1=1:i1
      disp(sprintf('x1(%d)=%21.13E, w1(%d)=%21.13E', j1, x1(j1), j, w1(j1)));
      sum = sum + w1(j1);
    end
    disp(sprintf('              integral(1.0, -1.0..1.0)=%21.13E', sum));
  end

  a2 = 0.5;
  b2 = 1.0;
  for i2=2:10
    [x2 w2] = gaulegf(a2, b2, i2);
    sum = 0.0;
    for  j2=1:i2
      sum = sum + w2(j2)*sin(x2(j2));
    end
    disp(sprintf('integral (0.5,1.0) sin(x) dx = %21.13E', sum))
  end
  disp(sprintf('-cos(1.0)+cos(0.5) = %21.13E', (-cos(1.0)+cos(0.5))))
  disp(sprintf('Maple says 3.372802560E-001'))
  diary off
  
end % gauleg.m  test program