MATLAB球谐系数的可视化

球谐展开是众多地球科学应用中经常使用的一种方法。球谐系数根据其阶次可以分为三类：（1）带谐系数：m=0对应的系数。此时带谐函数与经度λ无关并将地球球面分成许多纬度带。（2）扇谐系数：L=m≠0对应的系数。扇谐函数将地球球面分成正和负的扇形，相应在经度0~360区间内谐函数变换符号2L次。（3）田谐系数：L≠m≠0对应的系数。田谐函数将地球球面分为若干格，各个格之间正负交替变换。

下面我分享一个利用matlab对球谐系数进行可视化的程序：

% Draw a stretch sphereifext e n d =1 ,
% or a shaded s p h e r e i f e x t e n d=0
extend = 0 ;
% C a l c u l a t e t h e f i r s t few r e a l s p h e r i c a l harmon ics Y lm
for l = 0: 3
for m = -l : l % l ow e r c a s e L , no t one
% Se t up a g r i d
phis = linspace ( 0 , 2*pi , 200) ;
thetas = linspace ( 0 , 1*pi , 200) ;
[ phi , theta ] = ndgrid ( phis , thetas ) ;
% F i r s t c a l c u l a t e t h e a s s o c i a t e d Legendre f u n c t i o n o f
% c o s ( t h e t a ) w i t h d e g r e e l and o r de r ab s (m)
P = legendre ( l , cos ( theta ) ) ;
if ( l > 0 )
% I f l == 0 , l e g e n d r e r e t u r n s an o r de r 3 t e n s o r
% c o n t a i n i n g r e s u l t s f o r m=0 , m=1 , . . .
Plm = reshape (P( abs (m) + 1 , : , : ) , size ( phi ) ) ;
else
% I f l == 0 , l e g e n d r e r e t u r n s j u s t
% t h e m a t r ix f o r m=0
Plm = P;
end
% C a l c u l a t e t h e n o rm a l i z a t i o n c o n s t a n t N
a = (2*l +1)*factorial(l-abs (m)) ;
b = 4*pi*factorial ( l+abs (m) ) ;
N = sqrt ( a/b ) ;
% F i n a l l y , c a l c u l a t e t h e r e a l s p h e r i c a l harmon ic
if (m == 0 )
H = N.*Plm ;
elseif (m < 0 )
H = sqrt( 2 )*N*Plm .*sin ( abs (m)*phi ) ;
elseif (m > 0 )
H = sqrt( 2 )*N *Plm .*cos (m*phi ) ;
end

if ( extend ) % Draw by s t r e t c h i n g t h e s p h e r e
C = sign(H ) ; % Color v a l u e
S = abs (2*H) ; % s t r e t c h f a c t o r
Xm = cos(phi ).*sin (theta ) .*S ;
Ym = sin(phi ).*sin (theta ) .*S ;
Zm = cos(theta ).*S ;
else % Draw as a shaded s p h e r e
C = H; % Color v a l u e
Xm = cos(phi).*sin(theta ) ;
Ym = sin(phi).*sin(theta ) ;
Zm = cos(theta ) ;
end
% T r a n sl a t e t o p o s i t i o n in l a r g e r g r i d
X = Xm;
Y = Ym + m*2.9 ;
Z = Zm + 3 -l*2.9 ;
% Draw i t
colormap default ;
oldcmap = colormap ;
colormap( flipud ( oldcmap ) ) ;
h = surf (X, Y, Z , C ) ;
h . AmbientStrength = 0.6 ;
h . DiffuseStrength = 0.6 ;
h . SpecularStrength = 0.8 ;
h . SpecularExponent = 25;
shading flat ;
hold on ;
end
end
axis ([ -10 , 10 , -10, 10 , -10, 10 ] )
view (90 , 7 ) ;
grid off;
axis off
axis equal
colormap jet
xlabel ('x') ;
ylabel ('y') ;
zlabel ('z') ;
light ( 'Position', [ 10 -10 20 ] , 'Style', 'local ') ;
colorbar
% Reverse t h e X a x i s (MATLAB u s u a l l y draws x g o i n g i n t o t h e page )
gca.XDir = 'reverse ';
hold off

运行结果：

对应一般论文中经常出现的结果图：

（Chen et al., 2019,JGR-SE）

参考文献：

reeves_lee_apm_502_project_2_math_ma_portfolio_spring_2020-publish.pdf