根据开普勒元素快速计算卫星 LLA matlab代码.rar

function [elon nlat height]=eph2lla(filename) %{ README: text file format (delimited by space) -------------------- INTELSAT ZCZC 1002COE Epoch (UTCG): 24 Feb 2016 11:18:19.49 Semi-major Axis (km): 42167.061724 Eccentricity: 0.00015155449517 Inclination (deg): 0.039235178024 RAAN (deg): 136.08472194 AOP (deg): 198.01134961 Mean Anomaly (deg): 348.27154761 --------------------- author: Jonathon Y. Cheah, PhD. email: jcheah@ieee.org %} %------------------ % common constants: EarthR=6335.43; rad=pi/180; %------------------ %get ephermeris data from file %s{1} picks up the file by words delimited by space fid = fopen(filename); satID=fgets(fid); s = textscan(fid,'%s','delimiter',' '); fclose(fid); %------------------ %process the file data S=s{1}; %------------------ %get Kepler elements a=str2num(S{10}); e=str2num(S{12}); i=str2num(S{15})*rad; o=str2num(S{18})*rad; w=str2num(S{21})*rad; M=str2num(S{25})*rad; %------------------- %get ephermeris time stamp tstart=strcat(S{3},'-',S{4},'-',S{5},{' '},S{6}); [yr,mm,dy,h,m,s] = datevec(tstart); %------------------- %covert to Julian date if (mm<3) yr = yr-1; mm = mm+12; end A = fix(yr/100); B = 2 - A + fix(A/4); jd = floor(365.25*(yr+4716))+floor(30.6001*(mm+1))+ ... dy+B-1524.5+h/24+ m/1440+s/86400; %------------------- %The angle between vernal equinox and the Greenwich Meridian T=(jd-2451545)/36525; az=280.46061837+13185000.77005374*T+0.000387933*T^2-T^3/38710000; %keep az within 360 degrees. azG=(az-floor(az/360)*360)*rad; %------------------- %calculate true anomaly: %Eccentric anomaly using Mean anomaly by Kepler's equation err = 1; E0 = M; itt = 0; while(err > 1e-6) E = M + e*sin(E0); err=abs(E-E0); E0 = E; end v=2*atan2(sqrt(1+e)*sin(E/2),sqrt(1-e)*cos(E/2)); %-------------------- %The radius of satellite from true anomaly r=a*(1-e^2)/(1-e*cos(v)); %-------------------- %perifocal coordinates rpqw=[r*cos(v);r*sin(v);0]; %-------------------- % make ECI rotations % rotation to w rw=[cos(w) -sin(w) 0 ;sin(w) cos(w) 0; 0 0 1]; % rotation to i ri=[1 0 0; 0 cos(i) -sin(i); 0 sin(i) cos(i)]; %rotation to o ro=[cos(o) -sin(o) 0 ;sin(o) cos(o) 0; 0 0 1]; %-------------------- %Earth fixed frame convert LLA Ef=ro*ri*rw*rpqw; Az=(atan2(Ef(2),Ef(1))-azG)/rad; %Take care of >360 deg longitude rotation if abs(Az)>360 Az=Az-floor(Az/360)*360; end if abs(Az)>180 Az=Az-sign(Az)*360; end long=Az; lat=asin(Ef(3)/sqrt(Ef(1)^2+Ef(2)^2))/rad; ht=sqrt(Ef(1)^2+Ef(2)^2+Ef(3)^2)-EarthR; %--------------------- disp ('------------------------------------') disp ([' longitude = ',num2str(long),'�']) disp ([' latidue = ',num2str(lat),'�']) disp ([' height= ',num2str(ht), ' Km']) disp ('-------------------------------------') end