ВНИМАНИЕ! На форуме началось голосование в конкурсе - астрофотография месяца СЕНТЯБРЬ!
0 Пользователей и 1 Гость просматривают эту тему.
%начальные условия tb=11700;ti=12300;S0=7500;JD0=2456177.625;xtb=7003.008789; ytb=-12206.626953;ztb=21280.765625;Vxtb=0.7835417; Vytb=-2.8042530;Vztb=1.3525150;%расчет начальных улсовияwz=7.2921151467*10^(-5);Stb=S0+wz*(tb-10800);x0tb=xtb*cos(Stb)-ytb*sin(Stb);y0tb=xtb*sin(Stb)+ytb*cos(Stb);z0tb=ztb;xt0tb=Vxtb*cos(Stb)-Vytb*sin(Stb)-wz*y0tb;yt0tb=Vxtb*sin(Stb)+Vytb*cos(Stb)+wz*x0tb;zt0tb=Vztb;%направляющие косинусы и радиус-вектор на лунуT=(JD0+(tb-10800)/86400-2451545)/36525;e=0.4090926006-0.0002270711*T;G=1.4547885346 + 71.0176852437*T-0.0001801481*(T^2);sigmal=2.1824391966-33.7570459536*T + 0.0000362262*(T^2);ql=2.3555557435 + 8328.6914257190*T + 0.0001545547*(T^2);al=3.84385243e+05;% большая полуось орбиты Луныel=0.054900489;% эксцентриситет лунной орбитыil=0.0898041080;%среднее наклонение орбиты Луны к плоскости эклиптикиSezv=1-(cos(sigmal))^2 *(1-cos(il));Etazv=sin(sigmal)*sin(il);Jzv=cos(sigmal)*sin(il);Se11=sin(sigmal)*cos(sigmal)*(1-cos(il));Se12=1-(sin(sigmal))^2 * (1-cos(il));Eta11=Sezv*cos(e)-Jzv*sin(e);Eta12=Se11*cos(e)+Etazv*sin(e);J11=Sezv*sin(e)+Jzv*cos(e);J12=Se11*sin(e)-Etazv*cos(e);%Решение уравнения кеплера решается методом итераций, пока Е Е (пред) л л ? не будет меньше 10–8,eps = 10^(-8); % точностьx0 = 0.00; % начальное приближение N = 100; % количество итераций, чтобы не было зацикливанийg = inline('ql+el*sin(x0)','x0','ql','el'); x1 = g(x0,ql,el); % первое значение for i = 1 : N % делаем максимум 100 итераций if abs(x1 - x0) <= eps break end x0 = x1; x1 = g(x0,ql,el);endEl=x1;sinVl=(sqrt(1-el^2)*sin(El))/(1-el*cos(El));cosVl=(cos(El)-el)/(1-el*cos(El));Sel=(sinVl*cos(G)+cosVl*sin(G))*Se11+(cosVl*cos(G)-sinVl*sin(G))*Se12;Etal=(sinVl*cos(G)+cosVl*sin(G))*Eta11+(cosVl*cos(G)-sinVl*sin(G))*Eta12;Jl=(sinVl*cos(G)+cosVl*sin(G))*J11+(cosVl*cos(G)-sinVl*sin(G))*J12;rl=al*(1-el*cos(El));%направляющие косинусы и радиус-вектор на солнцеws = -7.6281824375 + 0.0300101976*T + 0.0000079741*(T^2);qs=6.2400601269+628.3019551714*T- 0.0000026820*(T^2);es=0.016719;as=1.49598*10^8;%Решение уравнения кеплера решается методом итераций, пока Е Е (пред) л л ? не будет меньше 10–8,eps = 10^(-8); % точностьy0 = 0.00; % начальное приближение N = 100; % количество итераций, чтобы не было зацикливанийg = inline('qs+es*sin(y0)','y0','qs','es'); y1 = g(y0,qs,es); % первое значение for i = 1 : N % делаем максимум 100 итераций if abs(y1 - y0) <= eps break end y0 = y1; y1 = g(y0,qs,es);endEs=y1;sinVs=(sqrt(1-es^2)*sin(Es))/(1-es*cos(Es));cosVs=(cos(Es)-es)/(1-es*cos(Es));Ses=cosVs*cos(ws)-sinVs*sin(ws);Etas=(sinVs*cos(ws)+cosVs*sin(ws))*cos(e);Js=(sinVs*cos(ws)+cosVs*sin(ws))*sin(e);;rs=as*(1-es*cos(Es));%Ускорения от лунных и солнечных гравитационных возмущенийGl=4902.799*10^(9);Gs=13271244*10^(13);xkl=x0tb/rl;ykl=y0tb/rl;zkl=z0tb/rl;Gkl=Gl/(rl^2);deltal=((Sel-xkl)^2+(Etal-ykl)^2+(Jl-zkl)^2)^(3/2);jx0l=Gkl*(((Sel-xkl)/deltal)-Sel);jy0l=Gkl*(((Etal-ykl)/deltal)-Etal);jz0l=Gkl*(((Jl-zkl)/deltal)-Jl);xks=x0tb/rs;yks=y0tb/rs;zks=z0tb/rs;Gks=Gs/(rs^2);deltas=((Ses-xks)^2+(Etas-yks)^2+(Js-zks)^2)^(3/2);jx0s=Gks*(((Ses-xks)/deltas)-Ses);jy0s=Gks*(((Etas-yks)/deltas)-Etas);jz0s=Gks*(((Js-zks)/deltas)-Js);J02=1082625.75*10^(-9);GM=398600441.8*10^6;ae=6378136;r0=sqrt(x0tb^2+y0tb^2+z0tb^2);GMk=GM/(r0)^2;xk0=x0tb/r0;yk0=y0tb/r0;zk0=z0tb/r0;ro=ae/r0;fun1=-GMk*xk0 - (3/2)*J02*GMk*xk0*ro^2 *(1-5*zk0^2)+jx0s+jx0l;fun2=-GMk*yk0 - (3/2)*J02*GMk*yk0*ro^2 *(1-5*zk0^2)+jy0s+jy0l;fun3=-GMk*zk0 - (3/2)*J02*GMk*zk0*ro^2 *(3-5*zk0^2)+jz0s+jz0l;