%
%   LPOLY_DEMO1.M    Sample script file to plot several low-order Legendre
%                  polynomials and to demonstrate their orthogonality property
%
%   This sample file generates plots of the P0(x) - P5(x) Legendre polynomials.  It
%   also uses Matlab's QUADL routine to evaluate integrals of Pm(x)Pn(x) over the 
%   interval [-1,1].  This should show the orthogonality property of the Legendre 
%   polynomials.
%  
%   The real purpose here is simply to demonstrate the use of Matlab's LEGENDRE 
%   function.
%
%   File written by J. R. White, UMass-Lowell (Aug. 2003)
%

%
%   getting started
      clear all, close all,  nfig = 0;
%
%   set color and marker code for creating plots
      Ncm = 6;
      scm = ['r-';      % red solid
             'g:';      % green dotted
             'b-';      % blue solid
             'm:';      % magenta dotted
             'c-';      % cyan solid
             'y:'];     % yellow dotted
%
%   set up independent variable
      Nx = 51;   x = linspace(-1,1,Nx);   
%
%   evaluate Pn(x) for n = 0:5  
%   (note that m = 0 in the associated Legendre polynomials gives the desired 
%    functions, and this is the first row of the variable returned from LEGENDRE)
      P = zeros(Nx,6);     % initialize space for storage of Pn(x)
      for n = 1:6
        AP = legendre(n-1,x);   P(:,n) = AP(1,:)';
      end
%
%   now let's plot all six curves
      nfig = nfig+1;  figure(nfig)
      for n = 1:6
        plot(x,P(:,n),scm(n,:),'LineWidth',2), grid on, hold on
        txt(n) = {['P',num2str(n-1),'(x)']};  
      end 
      title('LPoly\_Demo1:  Several Low-Order Legendre Polynomials')                               
      xlabel('x value'),ylabel('P_n(x)')
      legend(txt)
%
%   evaluate the orthogonality condition (use Matlab's QUADL routine)
      PmPn = zeros(6,6);   % initialize space for storage of PmPn integrals
      for n = 1:6
        for m = 1:6
          PmPn(m,n) = quadl('lpoly_demo1a',-1,1,[],[],m,n);
        end
      end
%
%   print out table of PmPn integrals
      fprintf('   n   m ->   0             1              2              3              4              5 \n')
      for n = 1:6
        fprintf(' %3i   %12.3e   %12.3e   %12.3e   %12.3e   %12.3e   %12.3e \n', ...
                 n-1,PmPn(:,n));
      end
%
%   end of demo
