通過兩個點的經緯度坐标，計算偏北角算法

以下是matlab代碼

clc

clear all

reflat = [121.5579868206571 29.901759261705514 0];%基準點經緯度

lla = [121.5575941023510 29.901770103264110 0];%待測點經緯度

% reflat = [121.5575660045000 29.901662429833333 0];

% lla = [121.5580470217427 29.901685205115914 0];

pointNum = size(lla,1);%待測點個數

enu = zeros(pointNum,3);%待測點東北天坐标

northAngle = zeros(pointNum,1);%待測點偏北角

for i = 1:pointNum

enu(i,:) = Fun.wgslla2enu(lla(i,2), lla(i,1), lla(i,3), reflat(2), reflat(1), reflat(3));%待測點東北天坐标求取

if abs(enu(i,1)) < 0.0001

if enu(i,2) > 0

northAngle(i) = 0;

else

northAngle(i) = 180;

end

else

if enu(i,1) > 0

northAngle(i) = 90 - atand(enu(i,2)/enu(i,1)); %目标偏北角

else

northAngle(i) = 270 - atand(enu(i,2)/enu(i,1)); %目标偏北角

end

end

end