calculate spherical harmonics G [usage] [G ,Itbl, Dim, Rmax] = vb_spherical_harmo_G( 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] G : 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 [G ,Itbl, Dim, Rmax] = vb_spherical_harmo_G(X,N,Rmax) 0002 % calculate spherical harmonics G 0003 % [usage] 0004 % [G ,Itbl, Dim, Rmax] = vb_spherical_harmo_G( 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 % G : 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=1~N 0036 G = 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 : ¶Ä³Ñ (elevation angle) : [-pi/2 pi/2 ] (0:¿åÊ¿ pi/2:¿¿¾å) 0043 % phi : Êý°Ì³Ñ(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 nn = 0; 0070 0071 for n = 0:N 0072 0073 % G = r .^ (-(n+1)) .* Y 0074 Rmulti = -(n+1); 0075 Rn = r .^ Rmulti; 0076 0077 % Cos-part 0078 for m = 0:n 0079 % Y(X,{n,m,cos}) 0080 nn = nn + 1; 0081 k = ix_tbl(n+1,m+1); 0082 G(:,nn) = Rn .* cos_Phi(:,m+1) .* P(:,k); 0083 Dim(:,nn) = [n; m; 0]; 0084 Itbl(n+1,m+1,1) = nn; 0085 end 0086 0087 % Sin-part 0088 for m = 1:n 0089 % Y(X,{n,m,sin}) 0090 nn = nn + 1; 0091 k = ix_tbl(n+1,m+1); 0092 G(:,nn) = Rn .* sin_Phi(:,m+1) .* P(:,k); 0093 Dim(:,nn) = [n; m; 1]; 0094 Itbl(n+1,m+1,2) = nn; 0095 end 0096 0097 % G = r .^ (-(n+1)) .* Y 0098 % Rn = Rn .* r; 0099 end 0100 0101 if nn ~= D, fprintf('Number of function = %d\n',nn); end;