輸入地理經緯度，怎樣用matlab計算地磁經緯度？

% 將大地經緯度轉化為地磁經緯度

% 磁極點取：經度72.2°W，緯度80.0°N（IGRF，2011）

% 輸入及輸出值為弧度制單位（rad）

% by chensd HHU 2016.01.02

% 磁極點經緯度

b=80.0*pi/180;

l=-72.2*pi/180;

% 以下公式基於劉長建，2011

sb=sin(Bg)*sin(b)+cos(Bg)*cos(b)*cos(Lg-l);

Bm=asin(sb);

sl=cos(Bg)*sin(Lg-l)/cos(Bm);

cl=(sin(b)*sb-sin(Bg))/(cos(b)*cos(Bm));

Lm=atan2(sl,cl);