Home > vbmeg > functions > leadfield > spherical_harmo > vb_spherical_harmo_F.m

vb_spherical_harmo_F

PURPOSE ^

calculate spherical harmonics

SYNOPSIS ^

function [Y ,Itbl, Dim, Rmax] = vb_spherical_harmo_F(X,N,Rmax)

DESCRIPTION ^

 calculate spherical harmonics
 [usage]
   [Y ,Itbl, Dim, Rmax] = vb_spherical_harmo_F( X, N, Rmax )
 [input]
   X : (NP*3) three-D coordinates of NP points
   N : maximum order of spherical harmonics (generally 10~20)
   Rmax : reference radius for normalization [m]
 [output]
   Y   : Spherical harmonics (NP x D),   D = N*(N+2)
  Itbl : Index table to map (n,m,0/1) to 'k' of Y(:,k)
  Dim  : Degree of Legendre function
  Rmax : reference radius for normalization [m]

 --- Example to get Spherical harmonics for (X ,n, m, c/s)
 k = Itbl(n+1,m+1,i+1)    for ( n = 0:N ; m = 0:N ; i=0,1:(c,s) )
 
 Y(:,k) = Sharmonics(X ,n, m, i) for ( n = 0:N ; m = 0:n ; i=0,1:(c,s) )
  if m > n, k = Itbl(n+1,m+1,i+1) = 0, which means Y(:,k) = 0

 Y(:,1) = Legendre(z ,n=0 ,m=0 ,i=0 ) = 1 
 Y(:,2) = Legendre(z ,n=1 ,m=0 ,i=0 )
 Y(:,3) = Legendre(z ,n=1 ,m=1 ,i=0 )

 Dim(1,k) = n : degree of polynomial r^n
 Dim(2,k) = m : degree of azimuthal angle rotation
 Dim(3,k) = 0 (:cos) or 1 (:sin)

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Y ,Itbl, Dim, Rmax] = vb_spherical_harmo_F(X,N,Rmax)
0002 % calculate spherical harmonics
0003 % [usage]
0004 %   [Y ,Itbl, Dim, Rmax] = vb_spherical_harmo_F( X, N, Rmax )
0005 % [input]
0006 %   X : (NP*3) three-D coordinates of NP points
0007 %   N : maximum order of spherical harmonics (generally 10~20)
0008 %   Rmax : reference radius for normalization [m]
0009 % [output]
0010 %   Y   : Spherical harmonics (NP x D),   D = N*(N+2)
0011 %  Itbl : Index table to map (n,m,0/1) to 'k' of Y(:,k)
0012 %  Dim  : Degree of Legendre function
0013 %  Rmax : reference radius for normalization [m]
0014 %
0015 % --- Example to get Spherical harmonics for (X ,n, m, c/s)
0016 % k = Itbl(n+1,m+1,i+1)    for ( n = 0:N ; m = 0:N ; i=0,1:(c,s) )
0017 %
0018 % Y(:,k) = Sharmonics(X ,n, m, i) for ( n = 0:N ; m = 0:n ; i=0,1:(c,s) )
0019 %  if m > n, k = Itbl(n+1,m+1,i+1) = 0, which means Y(:,k) = 0
0020 %
0021 % Y(:,1) = Legendre(z ,n=0 ,m=0 ,i=0 ) = 1
0022 % Y(:,2) = Legendre(z ,n=1 ,m=0 ,i=0 )
0023 % Y(:,3) = Legendre(z ,n=1 ,m=1 ,i=0 )
0024 %
0025 % Dim(1,k) = n : degree of polynomial r^n
0026 % Dim(2,k) = m : degree of azimuthal angle rotation
0027 % Dim(3,k) = 0 (:cos) or 1 (:sin)
0028 %
0029 % Copyright (C) 2011, ATR All Rights Reserved.
0030 % License : New BSD License(see VBMEG_LICENSE.txt)
0031 
0032 if ~exist('N','var'), N=10; end;
0033 
0034 NP = size(X,1);
0035 D  = (N + 1)^2 ; % !)(2n+1) : n=0:N
0036 Y  = zeros(NP,D);
0037 
0038 Dim  = zeros(3,D);
0039 Itbl = zeros(N+1,N+1,2);
0040 
0041 % convert coordinates from cartisian to spherical
0042 % theta  : $B6D3Q(B  (elevation angle)  : [-pi/2  pi/2 ] (0:$B?eJ?(B pi/2:$B??>e(B)
0043 % phi    : $BJ}0L3Q(B(azimuthal angle)  : [   0   2*pi ]
0044 [phi, theta, r] = cart2sph(X(:,1),X(:,2),X(:,3));
0045 
0046 % --- Legendre basis function
0047 % k = ix_tbl(n+1,m+1)        for ( n = 0:N ; m = 0:N )
0048 % P(:,k) = Legendre(z ,n, m)
0049 %    If m > n,  P(:,k) = 0
0050 
0051 z  = X(:,3)./r; % = cos(pi/2 - theta) (NP * 1)
0052 
0053 [P ,ix_tbl] = vb_legendre_basis(z,N);
0054 
0055 % --- Sin/Cos function
0056 % cos( m * phi ) = cos_Phi( m+1 )  for ( m = 0:N )
0057 % sin( m * phi ) = sin_Phi( m+1 )  for ( m = 0:N )
0058 
0059 Phi = phi * (0:N);  % (NP * (N+1))
0060 cos_Phi = cos(Phi); % (NP * (N+1))
0061 sin_Phi = sin(Phi); % (NP * (N+1))
0062 
0063 
0064 if ~exist('Rmax','var'), Rmax=sum(r)/NP; end;
0065 
0066 r  = r / Rmax; % (NP * 1)
0067 Rn = r;;
0068 
0069 % Sharmonics(X ,n=0 ,m=0 , i=0:(cos) ) = 1
0070 Itbl(1,1,1) = 1;
0071 Dim(:,1)    = [0; 0; 0];
0072 
0073 nn = 1;
0074 Y(:,1) = 1;
0075 
0076 for n = 1:N
0077 
0078     % Cos-part
0079     for m = 0:n
0080         % Y(X,{n,m,cos})
0081         nn = nn + 1;
0082         k  = ix_tbl(n+1,m+1);
0083         Y(:,nn) = Rn .* cos_Phi(:,m+1) .* P(:,k); 
0084         Dim(:,nn) = [n; m; 0];
0085         Itbl(n+1,m+1,1) = nn;
0086     end
0087     
0088     % Sin-part
0089     for m = 1:n
0090         % Y(X,{n,m,sin})
0091         nn = nn + 1;
0092         k  = ix_tbl(n+1,m+1);
0093         Y(:,nn) = Rn .* sin_Phi(:,m+1) .* P(:,k); 
0094         Dim(:,nn) = [n; m; 1];
0095         Itbl(n+1,m+1,2) = nn;
0096     end
0097     
0098     Rn = Rn .* r;
0099 end
0100 
0101 if nn ~= D, fprintf('Number of function = %d\n',nn); end;

Generated on Mon 22-May-2023 06:53:56 by m2html © 2005