球谐分析（如重力场）是将地球表面观测的某个物理量f(theta,lambda)展开成球谐函数的级数：
其中，theta为余纬，lambda：经度
一般地，Pnm为完全归一化的缔合勒让德多项式，其与无归一化的缔合勒让德多项式的Pnm0的关系为：
Pnm=(-1)^m*sqrt(k*(2N+1)(N-M)!/(N+M)!)Pnm0
其中k=1,当 m=0
k=2,当m>0
在Matlab中，有现成的缔合勒让德多项式：
-----------------------------------------------------------------------------------
LEGENDRE Associated Legendre function.
P = LEGENDRE(N,X) computes the associated Legendre functions
of degree N and order M = 0, 1, ..., N, evaluated for each element
of X. N must be a scalar integer and X must contain real values
between -1 <= X <= 1.
If X is a vector, P is an (N+1)-by-L matrix, where L = length(X).
The P(M+1,i) entry corresponds to the associated Legendre function
of degree N and order M evaluated at X(i).
There are three possible normalizations, LEGENDRE(N,X,normalize)
where normalize is 'unnorm','sch' or 'norm'.
The default, unnormalized associated Legendre functions are:
P(N,M;X) = (-1)^M * (1-X^2)^(M/2) * (d/dX)^M { P(N,X) },
where P(N,X) is the Legendre polynomial of degree N. Note that
the first row of P is the Legendre polynomial evaluated at X
(the M == 0 case).
SP = LEGENDRE(N,X,'sch') computes the Schmidt semi-normalized
associated Legendre functions SP(N,M;X). These functions are
related to the unnormalized associated Legendre functions
P(N,M;X) by:
SP(N,M;X) = P(N,X), M = 0
= (-1)^M * sqrt(2*(N-M)!/(N+M)!) * P(N,M;X), M > 0
因此，由sch正交化得到大地测量中的完全正交化，需要乘以renorm,其中renorm为
sqrt(2*N+1)
NP = LEGENDRE(N,X,'norm') computes the fully-normalized
associated Legendre functions NP(N,M;X). These functions are
associated Legendre functions NP(N,M;X). These functions are
normalized such that
/1
|
| [NP(N,M;X)]^2 dX = 1 ,
|
/-1
and are related to the unnormalized associated Legendre
functions P(N,M;X) by:
NP(N,M;X) = (-1)^M * sqrt((N+1/2)*(N-M)!/(N+M)!) * P(N,M;X)
因此，由norm正交化得到大地测量中的完全正交化，需要乘以renorm,其中
renorm=sqrt(2)，当m=0
renorm=2,当m>0
下面的函数使用matlab中的缔合勒让德函数legendre来由球谐系数计算球谐函数
-----------------------------------------------------------------------------------
例子1：确定CMB地形(数据来自Morelli,1987,nature)
下面%{ %}注释掉的分别为使用matlab legendre函数的sch正交化和norm正交化，然后进行modify的结果。
两者结果完全一致。
------------------------------------------------------------------------------
%clear
load('spherical_coeff.dat');
C=spherical_coeff(:,1);% coeff of cos
S=spherical_coeff(:,2);% coeff of sin
N=180;% grid resolution
Nmax=4;% max spherical degree of CMB topo
theta0=linspace(0,pi,N+1);% colatitude
phi0=linspace(0,2*pi,2*N+1);%longitude
Topo=zeros(length(theta0),length(phi0));
for i1=0:4
%{
% modified from 'norm ' option of matlab legendre function
P=legendre(i1,cos(theta0),'norm');%/sqrt(2*pi);
temp1=(i1+1)*i1/2+1;
for i2=1:i1+1
if (i2-1) == 0
renorm=sqrt(2);
else
renorm=2;
end
temp2=temp1+i2-1;
m=i2-1;
CS=C(temp2)*cos(m*phi0)+S(temp2)*sin(m*phi0);
Topo=Topo+renorm*P(i2,:).'*CS;
end
%}
%%{
% modified from 'sch ' option of matlab legendre function
P=legendre(i1,cos(theta0),'sch');%/sqrt(2*pi);
temp1=(i1+1)*i1/2+1;
for i2=1:i1+1
temp2=temp1+i2-1;
renorm=sqrt(2*i1+1);
m=i2-1;
CS=C(temp2)*cos(m*phi0)+S(temp2)*sin(m*phi0);
Topo=Topo+renorm*P(i2,:).'*CS;
end
%%}
end
Topo=Topo+2.5;%Notice that C00 is -2.5 relative PREM model
figure
axesm robinson
pcolorm((pi/2-theta0)/pi*180,phi0/pi*180,Topo)
%demcmap(caxis)
load coast
plotm(lat,long)
figure
v=-6:2:6;
contour(phi0/pi*180,theta0/pi*180,Topo,v),colorbar
set(gca,'ydir','reverse')
------------------------------------------------------------------------
其中的球谐系数文件spherical_coeff.dat如下：
-------------------------------------------------------------------
-2.5 0
0.34 0
0.07 -0.03
0.14 0
0.55 -1.02
0.29 0
-0.26 0
0.3 -0.04
1.15 -0.31
-0.11 -1.25
-0.11 0
0.22 0.36
-0.44 0.12
-0.66 -0.04
-0.23 -0.68
-------------------------------------------------------------------
所得结果与文章中图案一致，但是数值大小有差异，不知道其原因，有知道的，希望能告知
其结果见下面的伪彩图和等值线图
利用EGM球谐系数计算重力异常
下面是网上的利用EGM96球谐系数计算重力异常（相对于WGS-84椭球的重力场）
（http://gps.alaska.edu/jeff/Classes/GEOS602/downloads.html）
--------------------------------------------------------------------------------------------------------
% In order to refer heights to WGS-84, we must first subtract the even-numbered zonal
% coefficients(即次数为0的各偶数阶球谐系数，与经度无关，关于赤道对称) that go into the definition of the ellipsoid and gravity field.
%原因在于WGS-84中椭球项对应各偶数阶次数为0的球谐系数在EGM96中亦存在。因此，应减去之，才是差别
% See W. Torge (1980) Geodesy, section 5.2 for example, or page 20 of Lambeck (1988).
%For the WGS-84 ellipsoid, we remove the zonal coefficients J2 through J10
% based on the WGS-84 reference values.
%
% See http://cddis.nasa.gov/926/egm96/doc/S11.HTML for specific
% values of J2-J10. The values are just below equation 11-4.3. Note that there
% is an error in equation 11-4.3. The sqrt(2n+1) should be sqrt(4n+1).
%
%The equation for the gravity anomaly is
%
% delta_g = g0 * SUM(n=2,inf) [ SUM(m=0,n) [ (n-1) * Pnm(cos(theta) *
% { Cnm*cos(m*phi) + Snm*(sin(m*phi) } ] ]
%
% where the Cnm, Snm have had the zonal correction made, theta is
% colatitude and phi is longitude.
%注意：上面有一个(n-1)
%This program ignores the zero degree term, which amounts to a correction of
% -53 cm for EGM96 and WGS-84.
%
% This program makes two plots. The first is a rectangular plot of the
% potential (red=positive, blue=negative). The second plot shows
% the potential mapped onto a sphere.
%
%Author: Jeff Freymueller, Univ. of Alaska, Jan. 2004（http://gps.alaska.edu/jeff/jeff.html）
% Set the maximum degree and order here:
nmax = 360;
% Read EGM96. We read the nominal values of GM, a (earth radius), and the
% spherical harmonic coefficients Cnm and Snm. Each of these are
% dimensioned Cnm(360,361) for EGM96, which is complete to degree and order
% 360.
[GM, a, Cnm, Snm] = readEGM96('/home/hjs/yangting/elastic_thickness/EGM96/egm96_to360.ascii');
g0 = GM/a^2;
% Subtract the even zonals that are part of the ellipsoid definition for
% WGS-84
Cnm(2,1) = Cnm(2,1) + 0.108262982131*10^-2/sqrt(5); % J2
Cnm(4,1) = Cnm(4,1) - 0.237091120053*10^-5/sqrt(9); % J4
Cnm(6,1) = Cnm(6,1) + 0.608346498882*10^-8/sqrt(13); % J6
Cnm(8,1) = Cnm(8,1) - 0.142681087920*10^-10/sqrt(17); % J8
Cnm(10,1) = Cnm(10,1) + 0.121439275882*10^-13/sqrt(21); % J10
% Define a regular grid of theta (colatitude) and phi (longitude)
theta = linspace(0, pi, 2*180+1);
phi = linspace(0, 2*pi, 2*360+1);
% Calculate the Pnm. Note MATLAB's normalization is for the associated
% Legendre functions only: (Pnm, Pnm) = 1. When multiplied by the
% cos(m*phi) and sin(m*phi) terms, MATLAB's version does not result in a
% fully normalized spherical harmonic. Thus we add a renormalization to
% give the standard geodetic usage. It seems easiest to calculate the
% Schmidt normalization and renormalize that.
% Initialize vals with zeros. The summation of terms for gravity anomalies
% begins with n=2. The n=0 term (GM/r^2) is part of the
% definition of the ellipsoid, so it is subtracted implicitly along with
% the even zonals above.
delta_g = zeros(length(theta), length(phi));
% Do the summation over degree and order
% Note that if nmax<2, MATLAB will not execute this loop. Some languages
% execute a loop once no matter what, so beware if this is translated.
for in = 2:nmax,
Pnm = legendre(in, cos(theta), 'sch');
for m = 0:in,
% Renormalize the function in accord with Kaula standard used in geodesy
% This is the difference between Kaula's fully normalized and the Schmidt
% normalization
if ( m == 0 )
renorm = sqrt(2*(2*in+1));
%注意：源程序是错的，应该改为sqrt(2*in+1);
else
renorm = sqrt(in+1/2);
end
% Add the contribution from this n, m
delta_g = delta_g + renorm*(in-1)*g0*Pnm(m+1,:)'*( Cnm(in,m+1)*cos(m*phi) ...
+ Snm(in,m+1)*sin(m*phi) );
end
end
% Print out the range of geoid heights
disp('Minimum and maximum gravity anomaly')
disp( [max(max(delta_g)) min(min(delta_g)) ])
% Now make figures
% first figure: a rectangular plot
% Note the transformation of the axes.
% x-axis: longitude given in degrees
% y-axis: latitude in degrees
% recall that latitude = 90 - colatitude
figure
pcolor(phi*180/pi, 90-theta*180/pi, delta_g*1e5);
shading('flat');
% second figure: map the harmonic onto a sphere.
% We need to reverse the order of the rows, so that the northern
% hemisphere (theta = 0 to pi/2) ends up on the top of the sphere.
[nx ny] = size(delta_g);
for i=1:nx,
delta_g2(i,:) = delta_g(nx+1-i,:);
end
figure
sphere; h = findobj(gcf, 'Type', 'surface');
set(h, 'CData', delta_g2, 'FaceColor', 'texturemap')
axis equal;
ylabel('long = 0')
%%%%%%%%%%%%%%%自己修改的程序如下
nmax = 120;% Set the maximum degree and order here:
filename='/home/hjs/yangting/elastic_thickness/EGM96/egm96_to360.ascii';
GM = 3.986004415*10^14; % m^2/sec^3
a = 6378136.3; % m
% Read EGM96. We read the nominal values of GM, a (earth radius), and the
% spherical harmonic coefficients Cnm and Snm. Each of these are
% dimensioned Cnm(360,361) for EGM96, which is complete to degree and order
% 360.
[fid, err] = fopen(filename, 'r');
if ( fid == -1 )
disp(err);
end
N=sum(3:nmax+1);
temp1=dlmread(filename,'',[0 2 N-1 3]);
Cnm=temp1(:,1);
Snm=temp1(:,2);
g0 = GM/a^2;
% Subtract the even zonals that are part of the ellipsoid definition for
% WGS-84
Cnm(1) = Cnm(1) + 0.108262982131*10^-2/sqrt(5); % J2
Cnm(8) = Cnm(8) - 0.237091120053*10^-5/sqrt(9); % J4
Cnm(19) = Cnm(19) + 0.608346498882*10^-8/sqrt(13); % J6
Cnm(34) = Cnm(34) - 0.142681087920*10^-10/sqrt(17); % J8
Cnm(53) = Cnm(53) + 0.121439275882*10^-13/sqrt(21); % J10
% Define a regular grid of theta (colatitude) and phi (longitude)
theta = linspace(0, pi, nmax+1);
phi = linspace(0, 2*pi, 2*nmax+1);
% Calculate the Pnm. Note MATLAB's normalization is for the associated
% Legendre functions only: (Pnm, Pnm) = 1. When multiplied by the
% cos(m*phi) and sin(m*phi) terms, MATLAB's version does not result in a
% fully normalized spherical harmonic. Thus we add a renormalization to
% give the standard geodetic usage. It seems easiest to calculate the
% Schmidt normalization and renormalize that.
% Initialize vals with zeros. The summation of terms for gravity anomalies
% begins with n=2. The n=0 term (GM/r^2) is part of the
% definition of the ellipsoid, so it is subtracted implicitly along with
% the even zonals above.
delta_g = zeros(length(theta), length(phi));
% Do the summation over degree and order
% Note that if nmax<2, MATLAB will not execute this loop. Some languages
% execute a loop once no matter what, so be ware if this is translated.
k=1;
for in = 2:nmax,
Pnm = legendre(in, cos(theta), 'sch');
for m = 0:in,
% Renormalize the function in accord with Kaula standard used in geodesy
% This is the difference between Kaula's fully normalized and the Schmidt
% normalization
if ( m == 0 )
renorm = sqrt(2*(2*in+1));
else
renorm = sqrt(in+1/2);
end
% Add the contribution from this n, m
delta_g = delta_g + renorm*(in-1)*g0*Pnm(m+1,:)'*( Cnm(k)*cos(m*phi) ...
+ Snm(k)*sin(m*phi) );
k=k+1;
end
end
% Print out the range of geoid heights
disp('Minimum and maximum gravity anomaly')
disp( [max(max(delta_g)) min(min(delta_g)) ])
% Now make figures
% first figure: a rectangular plot
% Note the transformation of the axes.
% x-axis: longitude given in degrees
% y-axis: latitude in degrees
% recall that latitude = 90 - colatitude
figure
worldmap world
load coast
plotm(lat,long)
[long1,lat1]=meshgrid(phi*180/pi,90-theta*180/pi);
pcolorm(lat1,long1,delta_g*1e5);
figure
sphere;
h = findobj(gcf, 'Type', 'surface');
set(h, 'CData', delta_g2, 'FaceColor', 'texturemap')
axis equal;
axis off;
ylabel('long = 0')
for i1=0:360
view(i1,30);
drawnow
end