Диссертация (1091538), страница 16

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 16)

2.6function fp_n0i% Система (n,i)global p time;p=0.375;time=6e+14;fp_n0i.mfigure,hold on;% Именуем осиxlabel(’i’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);ylabel(’n_0’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);% Строим на графике стационарные точкиxc=GetRoot(1,0,0,-1,p); %Находим корни уз уравненияplot([-0.01 0.01], [xc(1) xc(1)],’-k’);plot([-0.01 0.01], [xc(2) xc(2)],’-k’);% Рисуем фазовые линии137˓→w=[[1e-3 3.14]; [1e-3 2.8]; [1e-3 2.3]; [6.5e-4 1.8]; [7e-4 1.7]; [5.8e-4 1.8];[3e-4 1.8]; [1e-3 1.6]; [1e-3 1];];for j=w’y=plotFP(j);end%--- График точек, в коротых производная равна 0 (зеленый)otv1=[];z=0:0.005:3;z1=0;for i=1:length(z);%находим корниotv1=[otv1; GetRoot(1,0,0,-cos(z(i)),p)];%Не учитываем комплексные корниif(imag(otv1(end))~=0 && z1==0)z1=i;endend%Если все корни действительныеif(z1==0)z1=length(z);endplot(z(1:z1), otv1(1:z1,:),’-g’);endfunction y=plotFP(y0)[t,y]=ode23s(@systNI, [0 1], y0);ind=find(y(:,1)>0.7);if(~isempty(ind))y(ind,:)=[];endplot(y(:,2), y(:,1), ’-b’);endfunction dx=systNI(t,x)%Система двух переменных n0, iglobal p time;n0=x(1);I=x(2);%Уравнения системыdn0 = -3*(n0^(16/3))*(cos(I) - p/(n0^(1/3)) - n0);di = -(n0^(13/3))*sin(I)/2;138dx=[dn0;di];dx=dx.*time; %Изменяем масштаб времениendФункции для построения графиков эволюции на рис.

1.6evol_EM.mfunction evol_EM% Построение эволюционных графиков для системы Земля-Лунаglobal gamma p time G0A delta f0 r;f=6.67*10^(-11); %Гравитационная постояннаяm = 5.9736*10^24;r = 6.378*10^6;Tp = 23.93419;mu = 734.9*10^20;a = 384.4*10^6;e = 0.0549;I = 5.15;Ts = 27.321661;disp(’Земля-Луна’);[p,n0,I]=getParam(m,mu,r,a,e,I,Tp,Ts)%p=0.17h = helpers;˓→%˓→%% Вычисляем параметр гамма для неизвестных параметровkoef = cos(I)-(p/(n0^(1/3)))*sqrt(1-e^2)dota = 1.204e-9;ur = -(mu*(1+mu/m)/((a^7)*(1-e^2)^(15/2)))*(h.F_3(e) (1/n0)*koef*((1-e^2)^1.5)*h.F_2(e));ur2=-(2*f0^(1/3)*(n0^(11/3))/((1-e^2)^(15/2)*G0A^(2/3)))*(n0*h.F_3(e) koef*((1-e^2)^1.5)*h.F_2(e))d2 = dota/ur2gamma =dota/uraa=gamma/r^7delta = gamma*mu*(G0A^(16/3))/(2*m*f*(f0^(5/3)))time=(1e+12)*31556926%1e+18;yearstime = (1e+12)%time/31556926year = 31556926˓→de = -delta*(n0^(13/3)/((1-e^2)^(13/2)))*(n0*h.F_4(e) koef*((1-e^2)^1.5)*h.F_5(e))139y0=[n0,e,I,a];t = linspace(0, 1, 100);y = lsode(@(y,t) systNEIA(t,y), y0,t);a2 = f0^(1/3)./(G0A*y(:,1)).^(2/3) ;hold on%plot(t, (a2./r), ’-g’);%plot(t, (y(:,4)./r), ’-r’);plot_a(t,y);plot_n0eia(t,y);% Дополнительные фазовые портреты%plot_n0e(t,y);%plot_n0i(t,y);%plot_ea(t,y);showEvolutionParams(t,y,yearstime);endfunction showEvolutionParams(t,y,yearstime)global r G0A;a_five = y(end,4)/1000;ar = y(:,4);k=0;for i=1:length(ar)-1;if(ar(i+1)<ar(i))k=i;break;endend% Значения параметров при максмальном удалении Луны от землиn0_max=y(k,1)e_max=y(k,2)I_max=y(k,3)*180/pia_max=ar(k)/1000a_max_r = ar(k)/rt_max=t(k)*yearstime %309.7505622712333e+9Ts_max=2*pi/(n0_max*G0A*86400)er=y(:,2);j=0;140for i=1:length(er);if(er(i)<=0)j=i-1break;endend% Значения параметров при достижении эксцентриситетом Луны нуляif(j~=0)n0_min=y(j,1)e_min=y(j,2)I_min=y(j,3)a_min=y(j,4)ar = y(:,4);a_min=ar(j)/1000a_min_r = ar(j)/rt_min=t(j)*yearstimeTs_min=2*pi/(n0_min*G0A*86400)endendfunction plot_a(t,y)% График эволюции большой полуосиglobal r;figure;plot(t, (y(:,4)./r), ’-b’);xlabel(’t’);ylabel(’a/R_{\oplus}’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);set(gca,’YLim’,[60 100]);box on;grid on;endfunction plot_n0eia(t,y)% График эволюции орбитальных параметровglobal r;figure;subplot(2, 2, 1);plot(t, y(:,1), ’-b’);xlabel(’t’);ylabel(’n_0’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;subplot(2, 2, 2);plot(t, y(:,2), ’-b’);xlabel(’t’);ylabel(’e’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);141grid on;subplot(2, 2, 3);plot(t, y(:,3).*180/pi, ’-b’);xlabel(’t’);ylabel(’i^{\circ}’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;subplot(2, 2, 4);plot(t, (y(:,4)./r), ’-b’);xlabel(’t’);ylabel(’a/R_{\oplus}’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;endfunction plot_ea(t,y)% Фазовый портрет системы (e,a)global r;figure;subplot(3, 1, 1);plot(t, (y(:,4)./r), ’-b’);xlabel(’t’);ylabel(’a/R_{\oplus}’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;subplot(3, 1, 2);plot(t, y(:,2), ’-b’);xlabel(’t’);ylabel(’e’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;subplot(3, 1, 3);plot((y(:,4)./r), y(:,2), ’-b’);ylabel(’e’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);xlabel(’a/R_{\oplus}’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;endfunction plot_n0e(t,y)% Фазовый портрет системы (n0,e)figure;subplot(3, 1, 1);plot(t, y(:,1), ’-b’);xlabel(’t’);142ylabel(’n_0grid on;’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);subplot(3, 1, 2);plot(t, y(:,2), ’-b’);xlabel(’t’);ylabel(’e’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;subplot(3, 1, 3);plot(y(:,2), y(:,1), ’-b’);xlabel(’e’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0)ylabel(’n_0’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0)grid on;endfunction plot_n0i(t,y)% Фазовый портрет системы (n0,i)figure;subplot(3, 1, 1);plot(t, y(:,1), ’-b’);xlabel(’t’);ylabel(’n_0’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;subplot(3, 1, 2);plot(t, y(:,3), ’-b’);xlabel(’t’);ylabel(’i’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;subplot(3, 1, 3);plot(y(:,3), y(:,1), ’-b’);xlabel(’i’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);ylabel(’n_0’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;endfunction dx=systNEIA(t,x)%Система переменных n0, e, i, aglobal time gamma G0A f0 delta p;n0=x(1);e=x(2);i=x(3);a=x(4);h=helpers;143koef = cos(i)-(p/(n0^(1/3)))*sqrt(1-e^2);dn0 = 3*(n0^(16/3)/((1-e^2)^(15/2)))*(n0*h.F_3(e) - koef*((1-e^2)^1.5)*h.F_2(e));de = -(e*n0^(13/3)/((1-e^2)^(13/2)))*(n0*h.F_4(e) - koef*((1-e^2)^1.5)*h.F_5(e));di = -(n0^(13/3)*sin(i)/(2*(1-e^2)^(5)))*h.F_1(e);da = -G0A*(2/3)*dn0*(a^(5/2))/(f0^(1/2));dx=[dn0;de;di;da];dx=dx.*delta*time; %Изменяем масштаб времениendfunction [p,n0,I]=getParam(m,mu,r,a,e,I,Tp,Ts)% Вычисление параметров p,n0,I для систем планета-спутник%m - масса планеты,%mu - масса спутника,%r - радиус планеты,%e - эксцентриситет,%I - наклонение,%Tp - период вращения планеты,%Ts- период обращения спутникаglobal G0A f0;I=I*pi/180; %Наклонение в радианахw=2*pi/(Tp*3600) %угловая скорость вращения планетыn=2*pi/(Ts*86400) %среднее движение спутникаf=6.67*10^(-11); %Гравитационная постояннаяf0=f*(m+mu)mr=m*mu/(m+mu)A=(2/5)*m*r^2G0=A*w + mr*f0^(2/3)*sqrt(1-e^2)*cos(I)/n^(1/3)p=((A^(1/3))*(f0^(2/3))*mr)/(G0^(4/3));n0=n*A/G0;G0A = G0/AendФункции для построения графиков эволюции Земля-Луна на рис.

1.6evol_n0ei.mfunction evol_n0ei%Построение графиков для системы уравнений (n0,e,i) 21.03.14% Для каждой системы строим графикиPlanetSys={{’EarthMoon’,’EM’},{’UpiterKalisto’,’UK’},{’UpiterIo’,’UI’},{’UpiterEurope’,’UE’},{’Ufor k=PlanetSysplotCurves(k);endend144function plotCurves(planetSys)global p time;p=1;time=1;% Получение начальных параметров[n0,e,i,p] = getP(planetSys{1}{1});y0=[n0,e,i];% Установке времениsetTime(planetSys{1}{1});% Интегрированиеt = linspace(0,1,100);%t = ode45(@systNEI, y0[0 1]);y = lsode("systNEI", y0,t);% Построение графиковhf = figure(’Units’, ’normalized’, ’OuterPosition’, [0.3 0.03 0.29 0.8]);grid on;subplot(3, 1, 1);plot(t, y(:,1), ’-b’);xlabel([’t_{’ planetSys{1}{2} ’} ’],’interpreter’,’tex’,’fontsize’,14,’rotation’,0);ylabel(’n_0’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);title([planetSys{1}{1}]);grid on;subplot(3, 1, 2);plot(t, y(:,2), ’-b’);xlabel([’t_{’ planetSys{1}{2} ’} ’],’interpreter’,’tex’,’fontsize’,14,’rotation’,0);ylabel(’e’,’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;subplot(3, 1, 3);plot(t, y(:,3), ’-b’);xlabel([’t_{’ planetSys{1}{2} ’} ’],’interpreter’,’tex’,’fontsize’,14,’rotation’,0);ylabel({’i,’,’rad’},’interpreter’,’tex’,’fontsize’,14,’rotation’,0);grid on;% Сохраняем график в векторним формате eps%s = [’/media/windows/Aspirant/Research/disser/FinalDoc_v3/images/’ planetSys{1}{1}˓→’.eps’]%saveas(hf,s, ’psc2’);%s = [’График сохранен в файл ’ s]145endfunction setTime(planetSys)% Устанавливаем время интегрирования% planetSys - название системыglobal time;switch planetSyscase ’EarthMoon’time=2e+12;case ’UpiterIo’time=1e+8;case ’UpiterEurope’time=1e+8;case ’UpiterGanimed’time=1e+11;case ’UpiterKalisto’time=1e+11;case ’MarsFobos’time=1.625e-4;case ’MarsDeimos’time=1e+8;case ’SaturnTitan’time=1e+10;endendfunction dx=systNEI(x,t)%Система трёх переменных n0, e, iglobal p time;n0=x(1);e=x(2);i=x(3);h=helpers;koef = cos(i)-(p/(n0^(1/3)))*sqrt(1-e^2);dn0 = 3*(n0^(16/3)/((1-e^2)^(15/2)))*(n0*h.F_3(e) - koef*((1-e^2)^1.5)*h.F_2(e));de = -(e*n0^(13/3)/((1-e^2)^(13/2)))*(n0*h.F_4(e) - koef*((1-e^2)^1.5)*h.F_5(e));di = -(n0^(13/3)*sin(i)/(2*(1-e^2)^(5)))*h.F_1(e);dx=[dn0;de;di];dx=dx.*time; %Изменяем масштаб времениend146function [n0,e,I,p] = getP(planetSys)global time;format long;switch planetSyscase ’EarthMoon’%а) Земля-Лунаm = 5.9736*10^24;r = 6.378*10^6;Tp = 23.93419;mu = 734.9*10^20;a = 384.4*10^6;e = 0.0549;I = 5.15;Ts = 27.321661;disp(’Земля-Луна’);[p,n0,I]=getParam(m,mu,r,a,e,I,Tp,Ts);case ’UpiterIo’%б) Юпитер-Иоm=1898.6*10^24;r = 71.398*10^6;Tp = 9.92425;mu = 893.3*10^20;a=421.6*10^6;e = 0.0041;I=0.04;Ts = 1.769138;disp(’Юпитер-Ио’);[p,n0,I]=getParam(m,mu,r,a,e,I,Tp,Ts);case ’UpiterEurope’%в) Юпитер-Европаm=1898.6*10^24;r = 71.398*10^6;Tp = 9.92425;mu = 479.7*10^20;a=670.9*10^6;e = 0.0101;I= 0.47;Ts = 3.55181;disp(’Юпитер-Европа’);[p,n0,I]=getParam(m,mu,r,a,e,I,Tp,Ts);case ’UpiterGanimed’%г) Юпитер-Ганимедm=1898.6*10^24;r = 71.398*10^6;Tp = 9.92425;147mu = 1482*10^20;a=1070*10^6;e = 0.0015;I=0.195;Ts = 7.154553;disp(’Юпитер-Ганимед’);[p,n0,I]=getParam(m,mu,r,a,e,I,Tp,Ts);case ’UpiterKalisto’%д) Юпитер-Каллистоm=1898.6*10^24;r = 71.398*10^6;Tp = 9.92425;mu = 1076*10^20;a=1883*10^6;e = 0.007;I=0.281;Ts = 16.689018;disp(’Юпитер-Каллисто’);[p,n0,I]=getParam(m,mu,r,a,e,I,Tp,Ts);case ’MarsFobos’%е) Марс-Фобосm=0.64185*10^24;r = 3.394*10^6;Tp = 24.622962;mu = 10.8*10^15;a=9337.2*10^3;e = 0.0151;I=1.082;Ts = 0.318910;disp(’Марс-Фобос’);[p,n0,I]=getParam(m,mu,r,a,e,I,Tp,Ts);case ’MarsDeimos’%ж) Марс-Деймосm=0.64185*10^24;r = 3.394*10^6;Tp = 24.622962;mu = 1.8*10^15;a=23463.2*10^3;e = 0.00033;I=1.791;Ts = 1.262441;disp(’Марс-Деймос’);[p,n0,I]=getParam(m,mu,r,a,e,I,Tp,Ts);case ’SaturnTitan’148%з) Сатурн-Титанm=568.46*10^24;r = 60.330*10^6;Tp = 10.65622;mu = 1345.5*10^20;a=1221.85*10^6;e = 0.0292;I=0.33;Ts = 15.945421;disp(’Сатурн-Титан’);[p,n0,I]=getParam(m,mu,r,a,e,I,Tp,Ts);endendfunction [p,n0,I]=getParam(m,mu,r,a,e,I,Tp,Ts)global time;% Вычисление параметров p,n0,I для систем планета-спутник%m - масса планеты,%mu - масса спутника,%r - радиус планеты,%e - эксцентриситет,%I - наклонение,%Tp - период вращения планеты,%Ts- период обращения спутникаI=I*pi/180; %Наклонение в радианахw=2*pi/(Tp*3600); %угловая скорость вращения планетыn=2*pi/(Ts*86400); %среднее движение спутникаf=6.67*10^(-11); %Гравитационная постояннаяf0=f*(m+mu);mr=m*mu/(m+mu);A=(2/5)*m*r^2;G0=A*w + mr*f0^(2/3)*sqrt(1-e^2)*cos(I)/n^(1/3);p=A^(1/3)*f0^(2/3)*mr/G0^(4/3);n0=n*A/G0;endВспомогательные функцииfunction h=helpersh.F_1 = @F_1;h.F_2 = @F_2;h.F_3 = @F_3;h.F_4 = @F_4;h.F_5 = @F_5;endhelpers.m149function otv=F_1(e)otv = 1 + 3*e^2 + 3*(e^4)/8;endfunction otv=F_2(e)otv = 1 + (15/2)*e^2 + (45/8)*e^4 + (5/16)*e^6;endfunction otv=F_3(e)otv = 1 + (31/2)*e^2 + (255/8)*e^4 + (185/16)*e^6 + (25/64)*e^8;endfunction otv=F_4(e)otv = 9 + (135/4)*e^2 + (135/8)*e^4 + (45/64)*e^6;endfunction otv=F_5(e)otv = (11/2) + (33/4)*e^2 + (11/16)*e^4;endGetRoot.mfunction x=GetRoot(A,B,C,D,E)% Нахождение корней многочлена 4-ой степени методом Феррари% http://ru.wikipedia.org/wiki/Метод_Феррариa = - (3/8)*(B/A)^2 + (C/A);b = (1/8)*(B/A)^3 - (B*C/A^2) + (D/A);g = -(3/256)*(B/A)^4 + (B^2*C/(16*A^3)) - (B*D)/(4*A^2) + (E/A);P=-a^2/12 - g;Q=- (a^3)/108 + a*g/3 - (b^2)/8;R= -Q/2 + sqrt(Q^2/4+P^3/27);U=R^(1/3);y=-5*a/6+U;if(U~=0)y=y-P/(3*U);elsey=y-(Q^(1/3));endW=sqrt(a+2*y);x1= -B/(4*A) + (W+sqrt(-(3*a+2*y+2*b/W)))/2;x2= -B/(4*A) + (W-sqrt(-(3*a+2*y+2*b/W)))/2;x=[x1 x2];x=x.^3;end§Б.2.

Характеристики

Список файлов диссертации