legendre 球諧函數的正交性與系數

物理大地測量的各項球諧函數滿足正交性， 因此可以用多階球諧函數的和擬合地球重力位， 將各階球諧函數視作是總合函數的基底。 由於其正交性質，使用內積，（對函數也就是積分。） 就能由總合中分離出各項系數。 本次作業練習由 legendre 函數定義出各階球諧函數， 再以任意系數組合出一組重力位， 最後再嘗試從合成的重力位裡分離出各項系數。

legendre 函數

legendre 函數是一系列的遞迴函數， 總合稱為 legendre 級數。 用遞迴的定義方法， 可以由第一二階定義出任意階數， 而第 n 階的第 m 次，則需要以微分計算， 這裡使用數值微分。

global p0
global p1
p0 = @(t) ones(1,columns(t));
p1 = @(t) t;
function pn = p_nth(n)
  global p0
  global p1
  if (n == 0)
    pn = p0;
  elseif (n == 1)
    pn = p1;
  else
    pn_m1 = p_nth(n-1);
    pn_m2 = p_nth(n-2);
    pn = @(t) -(n-1)/n .* pn_m2(t) + (2*n-1) / n .* t .* pn_m1(t);
  end
end

function dydx = derive_at (f, x, dx = 0.005)
  dydx = ( f(x.+dx) .- f(x.-dx) ) / (dx*2);
end

function dfdx = derive (f, dx = 0.005)
  dfdx = @(x) derive_at(f, x, dx);
end

function pnm = p_nmth(n,m)
  dpnm = p_nth(n);
  i = 0;
  while (i < m)
    dpnm = derive(dpnm);
    i++;
  end
  pnm = @(t) (1-t.^2).^(m/2) .* dpnm(t);
end

呼叫 p_nth 可以取得第 n 階的 legendre 函數， 呼叫 p_nmth 則可以取得第 n 階第 m 次的 legendre 函數， 其返回值都是函數，而且各次呼叫間獨立； 如果要計算高階項效率會很差， 但本次作業最高階只用到 6，所以還能接受。

球諧函數

球諧函數就是將 legendre 函數 輸入定義在 cos(ψ) ，再乘上正弦或餘弦函數。 這裡也用高階函數的型式，重新包裝前面定義的 p_nmth ， 使其接受 ψ 作為輸入，再分別定義 rnm snm 為原函數乘上正弦或餘弦的結果。

p_nmth_cos = @(n,m) @(r) p_nmth(n,m)(cos(r));

global rnm
rnm = @(n,m) @(th,l) p_nmth_cos(n,m)(th) * cos(m*l);
global snm
snm = @(n,m) @(th,l) p_nmth_cos(n,m)(th) * sin(m*l);

各階各次球諧函數圖型

在繪圖時，發現如果用 180 * 360 的球去畫， 我的程式會 buffer overflow， 所以就簡化成 45 * 45 的球， 每 8° 經度 4° 緯度取樣一次。 網格會變比較粗糙，但色階大致正確。

q1_rnm = [2 0
          3 1
          4 0
          4 4
          5 0
          6 0];
q1_snm = [2 1
          3 3
          4 2
          4 2
          5 1
          6 1];

global x_degree_range = 0:359;
global x_radian_range = x_degree_range /180 *pi;
global y_degree_range = 0:179;
global y_radian_range = y_degree_range /180 *pi;

function hdl = plot_earth_grid_180x360(grid)
  % grid (latitude, longitude) = (1:180, 1:360)
  grid_square = grid(1:4:end,1:8:end);
  [x y z] = sphere(44);
  hdl = surf(x,y,z, 'cdata', grid_square);
  colorbar;
  shading interp;
end

function plot_grid_180x360_tofile(grid, name)
  fig = figure('visible', 'off');
  plot_earth_grid_180x360(grid);
  print('-dpng', [name '-sphere.png']);
  close(fig);
end

for i=1:rows(q1_rnm)
  n = q1_rnm(i,1);
  m = q1_rnm(i,2);
  frnm = rnm(n, m);
  gridi = frnm(y_radian_range', x_radian_range);
  plot_grid_180x360_tofile(gridi, ['r' num2str(n) num2str(m)]);

  n = q1_snm(i,1);
  m = q1_snm(i,2);
  fsnm = snm(n, m);
  gridi = fsnm(y_radian_range', x_radian_range);
  plot_grid_180x360_tofile(gridi, ['s' num2str(n) num2str(m)]);
end

帶諧

扇諧

田諧

合成球諧函數

先隨機產生 n=0:6 m=0:n 的系數， 再將這些系數乘上各球諧函數， 加總得到合成的球諧函數網格。 系數是隨機產生的，在 0 1 之間。

global coefficient_data
coefficient_data = {};

function xnm = coefficient_xnm(type, n, m, set_value = [])
  global coefficient_data
  if (type == 'r')
    type_id = 1;
  else
    type_id = 2;
  end

  if (isempty(set_value))
    xnm = coefficient_data{type_id, n+1, m+1};
  else
    xnm = coefficient_data{type_id, n+1, m+1} = set_value;
  end
end

coefficient_sum_grid = zeros(columns(y_radian_range), columns(x_radian_range));
for n = 0:6
  for m = 0:n
    % cell start from 1
    rnmc = rand(1);
    coefficient_xnm('r', n, m, rnmc);
    coefficient_sum_grid += rnmc * rnm(n, m)(y_radian_range', x_radian_range);

    snmc = rand(1);
    coefficient_xnm('s', n, m, snmc);
    coefficient_sum_grid += snmc * snm(n, m)(y_radian_range', x_radian_range);
  end
end
plot_grid_180x360_tofile(coefficient_sum_grid, 'coefficient-sum');

估計各項球諧系數的振幅

由之後分離系數所使用的積分，其積分結果產生的系數 2*pi/(2*n+1)*factorial(n+m)/factorial(n-m) 可以得知， 雖然各項球諧的系數都介於 0-1 之間， 但 m 的次數越高，其結果的振幅也會越大， 產生的系數也會越大，二者約是階乘的關係。

我選取了 S42 項作實驗，其圖形如上述圖九。 在維持其它系數不變，獨立調整 S42 係數， 繪出 B42 為 500 與 999 的情況如下。 在 500 時就能在原本合成的彩度中看出 S42 的形狀， 到 999 時則是二者混合，再增加就和 S42 的獨立圖形無太大區別了， 也就不畫了。

以正交性質分離系數

因為球諧函數各項間為正交， 故在函數上積分後如同內積， 可以消去其它所有正交的不同次項， 只留下原本的項。 可以運用此一性質從合成的函數結果中， 分離出各項的系數。

function coefficient = coefficient_devide(grid_all, sh_type, n, m)
  global x_radian_range
  global y_radian_range
  global rnm
  global snm
  if (sh_type == 'r')
    fnm = rnm(n,m);
  else % 's'
    fnm = snm(n,m);
  end

  if (m == 0)
    coefficient_integral = 4*pi/(2*n+1);
  else
    coefficient_integral = 2*pi/(2*n+1)*factorial(n+m)/factorial(n-m);
  end

  grid_single = fnm(y_radian_range', x_radian_range);
  grid_sin = sin(y_radian_range' * ones(size(x_radian_range)));
  sum2d = @(grid) sum(sum(grid));
  grid_integral = sum2d(grid_single .* grid_all .* grid_sin) * (
                    pi*2*pi /180/360
                  );
  coefficient = grid_integral / coefficient_integral;
end


for n=0:6
  for m = 0:n
    printf('R%d%d %d\n',
           n, m, coefficient_devide(coefficient_sum_grid, 'r', n,m));
    printf('S%d%d %d\n',
           n, m, coefficient_devide(coefficient_sum_grid, 's', n,m));
  end
end

其中 Snm 中 m=0 的項是 0， 所以該項的系數其實是沒有意義的； 但我在產生系數時沒有特別濾掉， 所以這裡也才會有數值，但分離出來還是 0。 分離出的系數都相當準確，可以準確到小數下第二位。