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)
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;