From octave-sources-request at bevo dot che dot wisc dot edu  Fri Oct 23 04:04:04 1998
Subject: Legendre functions
From: Xavier Calbet <xcalbet at inm dot es>
To: octave-sources at bevo dot che dot wisc dot edu
Cc: xcalbet at inm dot es
Date: Fri, 23 Oct 1998 09:59:24 +0100

This is a multi-part message in MIME format.
--------------1790046A7D8C6D452CBEA343
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Hello,

  I believe there are no Legendre polinomials right now in Octave.

  Mari Carmen Romero and myself, Xavier Calbet, have
translated the Legendre functions written in Yorick to Octave.

  I enclose the function if it is of any use for anybody.

  This function seems to work for the Legendre polinomials (m=0)
--------------1790046A7D8C6D452CBEA343
Content-Type: text/plain; charset=us-ascii; name="legen.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline; filename="legen.m"


%   Compute associated Legendre functions.
%
%   Translated from the corresponding Yorick function
%   to Octave by
%              Mari Carmern Romero
%                     and
%              Xavier Calbet
%


function retval = legndr(l,m, x)
% DOCUMENT legndr(l,m, x)
%     return the associated Legendre function Plm(x).  The X may
%     be an array (-1<=x<=1), but L and M (0<=M<=L) must be scalar
%     values.  For m=0, these are the Legendre polynomials Pl(x).
%     Relation of Plm(x) to Pl(x):
%       Plm(x) = (-1)^m (1-x^2)^(m/2) d^m/dx^m(Pl(x))
%     Relation of Plm(x) to spherical harmonics Ylm:
%       Ylm(theta,phi)= sqrt((2*l+1)(l-m)!/(4*pi*(l+m)!)) *
%                           Plm(cos(theta)) * exp(1i*m*phi)
%   SEE ALSO: ylm_coef
%

  if (nargin ~= 3)
    usage("legndr(integer,integer,vector)");
  endif

%  m+= 0;
%  l+= 0;

  if (m<0 || m>l) 
    usage( "l, m integers with 0<=m<=l for Plm");
  endif

  if (m > 0)
      % sqrt blows up if x out of range
      if ( fix(m/2)*2 == m )
	pmm = exp(sum(log(1:2:2*m))) * (sqrt((1.-x).*(1.+x)).^m);
      else
	pmm =-exp(sum(log(1:2:2*m))) * (sqrt((1.-x).*(1.+x)).^m);
      endif
    else 
    if (any(abs(x)>1.0)) 
      usage( "-1<=x<=1 for Plm");
    endif
    pmm= ones(1, length(x));  
  endif

  if (l==m)
    retval=pmm;
    return;
  endif

  pmmp1= (2*m+1)*x.*pmm;
  
  
  if (l==m+1) 
    retval= pmmp1;
    return;
  endif
  ll=m+2;
  while(1)
    
    c1= ll-m;
    c= (ll+m-1)/c1;
    c1= (2*ll-1)/c1;
    pll= c1*x.*pmmp1 - c*pmm;
    if (ll==l) 
      retval= pll;
      return;
    endif
    pmm= pmmp1;
    pmmp1= pll;
    ll++; 
  endwhile

endfunction
  





















--------------1790046A7D8C6D452CBEA343--