Home > atmlab > math > gauss_laguerre.m

gauss_laguerre

PURPOSE ^

gauss_laguerre returns abscissas and weights for Laguerre Gauss quadrature

SYNOPSIS ^

function [x_i,w_i]=gauss_laguerre(alpha,N,varargin)

DESCRIPTION ^

 gauss_laguerre returns abscissas and weights for Laguerre Gauss quadrature
    
     Returns a vector with abscissas and weights for performing
     Laguerre Gauss quadrature to represent scattering properties
     of particle polydispersions. 

     To apply the integration weights use *gauss_laguerre_apply*.

     An example on usage (to generate a PSD for ice particles):
        xnorm = 50e-6;
        [x_i,w_i] = gauss_laguerre(0,10,xnorm);
        y = ice_psd_gamma( ... );
        y = gauss_laguerre_apply(y,x_i,w_i,xnorm);

 FORMAT   [x_i,w_i] = gauss_laguerre(alpha,N[,xnorm])     

 OUT      x_i is a vector with abscissas
          w_i is a vector with weights
 
 IN       alpha  exponent of weighting function x^alpha*exp(-x)
          N      number of points
 OPT      xnorm  Normalisation value for length scale. Original point of 1
                 will be mapped to this value. Default is 1.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

gauss_laguerre.m

SOURCE CODE ^

0001 % gauss_laguerre returns abscissas and weights for Laguerre Gauss quadrature
0002 %
0003 %     Returns a vector with abscissas and weights for performing
0004 %     Laguerre Gauss quadrature to represent scattering properties
0005 %     of particle polydispersions.
0006 %
0007 %     To apply the integration weights use *gauss_laguerre_apply*.
0008 %
0009 %     An example on usage (to generate a PSD for ice particles):
0010 %        xnorm = 50e-6;
0011 %        [x_i,w_i] = gauss_laguerre(0,10,xnorm);
0012 %        y = ice_psd_gamma( ... );
0013 %        y = gauss_laguerre_apply(y,x_i,w_i,xnorm);
0014 %
0015 % FORMAT   [x_i,w_i] = gauss_laguerre(alpha,N[,xnorm])
0016 %
0017 % OUT      x_i is a vector with abscissas
0018 %          w_i is a vector with weights
0019 %
0020 % IN       alpha  exponent of weighting function x^alpha*exp(-x)
0021 %          N      number of points
0022 % OPT      xnorm  Normalisation value for length scale. Original point of 1
0023 %                 will be mapped to this value. Default is 1.
0024 
0025 % History: 2005-05-13  Created by Bengt Rydberg
0026 
0027 function [x_i,w_i]=gauss_laguerre(alpha,N,varargin)
0028 %
0029 [xnorm] = optargs( varargin, { 1 } );
0030 
0031 
0032 L_size=100000;
0033 x=linspace(0,100,L_size)';
0034 % get Laguerre polynomials
0035 for i=0:N+1
0036     if i==0
0037        L(1:L_size,1)=1;
0038     elseif i==1
0039        L(1:L_size,2)=-x+alpha+1;
0040     else
0041       L(1:L_size,i+1)=((-x+2*(i-1)+alpha+1.0).*L(1:L_size,i)-(i-1.0+alpha)*L(1:L_size,i-1))/i;
0042     end
0043 end
0044 % find abscissas (roots of laguerre polynom of order N)
0045 s=1;
0046 for i=1:length(x)-1
0047     if L(i,N+1)*L(i+1,N+1)<=0
0048        x_i(s)=x(i);
0049        L_i(s)=L(i,N+1);
0050        L_i_a(s)=L(i,N+2);
0051        s=s+1;
0052     end
0053 end
0054 n_fak=1;
0055 for i=1:N
0056     n_fak=n_fak*i;
0057 end
0058 % get weights
0059 w_i=gamma(N+alpha+1)*x_i./n_fak/(N+1)^2./L_i_a.^2;
0060 
0061 
0062 
0063 x_i = x_i * xnorm;
0064

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005