issatellitevisiable.m

根据卫星星历及测者位置

function isSatVis = isSatelliteVisiable(UserPosition,Satellite,t)
%　验证在某一时刻，某一位置，某一颗卫星的可见性，卫星可见则返回1,否则返回0
%　UserPosition = [latitude longitude];
%    其中-180 < longitude < 180 东经为正，西经为负；-90< latitude <90, 北纬为正，南纬为负
%  Satellite = [ i Omiga a e w M ]; 　　　　其中6个分量为初始时刻t0时的初始值
%　　6个分量分别是轨道面倾角，升交点赤经，长半轴，轨道偏心率，近地点角中，平近地角
%　求出t时刻卫星位置，将t时刻用户地理坐标转换成空间坐标，验证卫星的可见性

DtoR = 2*pi/360;
CommandAngle = 90 * 2*pi/360; % 地心－用户－卫星　间夹角大于该角度时，卫星可见
                              %　如果要求10度仰角，则最夹角最低为100度
%%%%%%%%%%%%%%%%%%%%% 求解卫星位置
%  GASTw : 星期六午夜至星期天子夜的交换时刻格林尼治视恒星时
% SatPos = [Xk Xk Xk]';
GASTw = pi/3;
SatPos = SatellitePos(Satellite,GASTw,t);

%%%%%%%%%%%%%%%%%%%%%　用户坐标转换
Re = 6366707;   % 地球圆体半径
long = UserPosition(2) * DtoR;
lat = UserPosition(1) * DtoR;
% 用户的地心直角坐标
Xu = Re * cos(long);
Yu = Re * sin(long);
Zu = Re * sin(lat);
UserPos = [Xu Yu Zu]';

Delta_Pos = SatPos - UserPos;
RR = [-sin(lat)*cos(long)  -sin(lat)*sin(long)  cos(lat);
      -sin(long)           cos(long)            0;
      cos(lat)*cos(long)   cos(lat)*sin(long)   sin(lat);];  % 转换距阵
PosUser = RR *  Delta_Pos;  %　卫星在观测站的站心直角坐标系坐标
Xps = PosUser(1);
Yps = PosUser(2);
Zps = PosUser(3);
%%%%%%%%%%%%%%%%%%%% 　验证卫星是否可见
% 用站心极坐标表示卫星位置
% Rp - 测者到卫星间的距离；Ap - 卫星方位角；Ep - 卫星高度角
Rp = sqrt(Xps^2+Yps^2+Zps^2);
Ap = atan(Yps/Xps);
Ep = atan(Zps/sqrt(Xps^2+Yps^2));
if Ep>0
    isSatVis = 1;
else
    isSatVis=0;
end