MATLAB實作地球表面上兩點之間的仰角和方位角計算

背景：兩個天線分別坐落在兩個山頭上（或者其他什么地方均可，中間要保證無阻礙物），它們的波束角度都很小，大概為一度，且相距很遠（幾十公里），無法直接用望遠鏡觀察位置來手動調整，因此需要通過計算調整各自的方位和仰角，兩個點的經緯度和海拔均為已知條件，

首先，地球并不是完美的球體，它更接近一個兩極比赤道更扁的橢球，查閱資料可得，赤道半徑約為6378137米，而極半徑約為6356752米，中間大概相差了22000米，因此其造成的誤差已經無法忽略，對于這種情況，即使兩個點的海拔一樣高，在實際空間上也會有一定的高度差，會需要仰角的小幅度偏移，

我的思路：

第一步，將經緯度和海拔轉化成常用的笛卡爾坐標（xyz），這樣計算起來比較熟悉和方便，通過對橢球中任意一個子午面的分析可得出換算公式：

圖中的y對應笛卡爾坐標系內P點的z軸，x是坐標系內P點投影向量在（x，y）與原點之間連接的線段，原理類似于球體極坐標轉換成直角坐標系，根據勾股定理等可推出：

P(x) = ( N + 海拔 ) .* cos( 緯度) .* cos( 經度)
P(y) = ( N + 海拔 ) .* cos( 緯度) .* sin( 經度 )
P(z)= ( ( 1.0 - e2) .* N + 海拔 ) .* sin( 緯度 )

其中e2 = (a^2-b^2)/(a^2)
W ^2 =(1.0-e2*((sin(緯度))^2)) （這是一個常用的輔助函式，了解即可）
N = a/W;

根據這樣的思路即可給出代碼：

clear;
clc;
%%
%從橢球坐標（經緯度和海拔）轉化成笛卡爾坐標
 plh = [30.76,104.08,503.67;
        30.58,104.08,503.67]; %隨便給兩個點的坐標，分別對應緯度經度和海拔
 
 lat = plh(:,1)*pi/180; %latitude
 lon = plh(:,2)*pi/180; %longitude
 hgt = plh(:,3); %altitude
   a = 6378137.0;
   b = 6356752.3;
       
if abs(lat) > pi/2 | abs(lon) > 2*pi
   warning( 'In plh2xyz(), the input latitude and/or longitude may not be in units of radians.' );
end

e2 = (a^2-b^2)/(a^2);
W  = sqrt(1.0-e2*((sin(lat)).^2));
N  = a./W;

pos(:,1) = ( N + hgt ) .* cos( lat ) .* cos( lon );   %x
pos(:,2) = ( N + hgt ) .* cos( lat ) .* sin( lon );   %y
pos(:,3) = ( ( 1.0 - e2 ) .* N + hgt ) .* sin( lat ); %z

Raim=sqrt(pos(:,1).^2+pos(:,2).^2+pos(:,3).^2);
pos

這樣得到的pos就是一個2x3的矩陣，包含了兩個點的xyz坐標，

第二步，分別計算各自的仰角：

仰角很容易理解，如果是P0對P1點的仰角，可以看成是向量P0P1和P0點對應的橢球切面的夾角，那么只要計算出平面的法向量和向量P0P1的夾角，再用90°減去它即可，然后P1對P0的仰角同理，

根據這個思路，可以寫出如下代碼：

%%
%求任意一點P0(x0,y0,z0)對于點P1(x1,y1,z1)的仰角w1
%過P0點的切平面：x0*x/a^2+y0*y/a^2+z0*z/b^2-1=0
%過P0P1的直線：(x-x0)/(x1-x0)=(y-y0)/(y1-y0)=(z-z0)/(z1-z0)
x0=pos(1,1);y0=pos(1,2);z0=pos(1,3);
x1=pos(2,1);y1=pos(2,2);z1=pos(2,3);
upmn=(x1-x0)*x0/a^2+(y1-y0)*y0/a^2+(z1-z0)*z0/b^2;
downmn=sqrt(((x1-x0).^2+(y1-y0).^2+(z1-z0).^2)*((x0^2+y0^2)/a^4+z0^2/b^4));

if pos(1,3)>=pos(2,3)
    w1=90-acosd(upmn/downmn);
else
    w1=acosd(upmn/downmn)-90;
end
w1

需要注意的是，因為在P0和P1兩個點上，它們各自的切面斜率不同，因此兩個點的仰角絕對值有細微差別，并且因為方向一定是相反的（一個向上一個向下），所以得到的數值一個是正數一個是負數（要么就都是0），

第三步，計算各自的方向角（如：北偏東45°）：

思路可參考：https://blog.csdn.net/u014531567/article/details/70767732?utm_medium=distribute.pc_relevant.none-task-blog-blogcommendfrommachinelearnpai2-1.edu_weight&depth_1-utm_source=distribute.pc_relevant.none-task-blog-blogcommendfrommachinelearnpai2-1.edu_weight，只要求出這個二面角即可，

根據這個思路，可以寫出下列代碼：

%%
%調整方向角度w2
P1OC=acosd(abs(z1)/sqrt(x1^2+y1^2+z1^2));
EOF=plh(2,2)-plh(1,2);
P0OP1=acosd((x0*x1+y0*y1+z0*z1)/sqrt((x0^2+y0^2+z0^2)*(x1^2+y1^2+z1^2)));
C0P0P1=asind(sind(P1OC)*sind(EOF)/sind(P0OP1));
w2=abs(C0P0P1);
if plh(2,1)>=plh(1,1) & plh(2,2)>=plh(1,2)
    ww=char('The azimuth is from North to East,degree is', num2str(w2));
elseif plh(2,1)>=plh(1,1) & plh(2,2)<plh(1,2)
    ww=char('The azimuth is from North to West,degree is', num2str(w2));
elseif plh(2,1)<plh(1,1) & plh(2,2)<plh(1,2)
    ww=char('The azimuth is from South to West,degree is', num2str(w2));
else
    ww=char('The azimuth is from South to East,degree is', num2str(w2));
end
    w2 %方位角w2
    ww %列印指導：從南（北）向東（西）轉向

這樣之后，天線的水平轉向角即可確定，至此，天線的波束理論上即可相互連接，只需要在之后做一些較準即可，

附上運行結果：

以兩個點P0(30.76,104.08,1503.67)和P1(30.58,104.04,503.67)舉例（分別是緯度，經度和海拔）：

結果如下：

如圖所示，P0對P1的仰角為w1=-2.9086°，P0對P1的方位角為S to W，w2=10.8822°，驗證成功，