solar_pos.txt

利用matlab实现太阳在地心惯性坐标系的位置计算

clear;
time=2500000;
nnn=0;
for time=2500000:0.001:2500002.4;
nnn=nnn+1;
mjd=time-2415020;
secday=1;
AU=1.5*10^11;              %日地平均距离
year=1900+mjd/365.25;
DeltaET_year=26.465+0.747622*(year-1950)+1.886913*sin(2.0*pi*(year-1975.0)/33.0);
T=(mjd+DeltaET_year/secday)/36525.0;
M=(mod(358.47583+mod(35999.04975*T,360.0)-(0.000150...
    +0.0000033*T)*sqrt(T),360.0))*pi/180;
L=(mod(279.69668+mod(36000.76892*T,360.0)+0.0003025*sqrt(T),360.0))*pi/180;
e=0.01675104-(0.0000418+0.000000126*T)*T;
C=((1.919460-90.004789+0.000014*T)*T)*sin(M)...
    +(0.020094-0.000100*T)*sin(2*M)+0.000293*sin(3*M)*pi/180;
O=mod(259.18-1934.142*T,360.0)*pi/180;
Lsa=mod(L+C-(0.00569-0.00479*sin(O))*pi/180,pi*2.0);
nu=mod(M+C,pi*2);
R=1.0000002*(1-sqrt(e))/(1+e*cos(nu));
eps=(23.452294-(0.0130125+(0.00000164-0.000000503*T)*T)*T+0.00256*cos(O))*pi/180;
R=AU*R;
solar_vector_x(nnn)=R*cos(Lsa);
solar_vector_y(nnn)=R*sin(Lsa)*cos(eps);
solar_vector_z(nnn)=R*sin(Lsa)*sin(eps);
solar_vector_w=R;

end
[x0,y0,z0]=sphere(20);
x=637500000*x0;y=y0*637500000;z=z0*637500000;
plot3(solar_vector_x,solar_vector_y,solar_vector_z);
hold on;
surf(x,y,z);