How to find minimum value for a loop iteration if we have two variable at some min value i need to save time and velocity ?

2 visualizzazioni (ultimi 30 giorni)
%file first
function [r_b,v_b,zenit,Elevation,Azi,close,minimum_time]=main(a,e,eta,omega_1,omega_2)
GM= 398600.44;
n = sqrt(GM/(a*a*a)) %mean angular velocity/mean motion
T= 2*pi/(n);
k=1;
for i=0:10:2*T
M= n *(i-0); %mean anomly
E0=M;
err=1;
while(err>10^-12)
E= M+e*sin(E0);
err=abs(E-E0); % I want to understand this step
E0=E; % eccentric anomly
end
v=2*atan(sqrt(1+e)*tan(E/2)/sqrt(1-e)); %atan2 is not working % true anomly
r = a*(1-e*cos(E)); %distance of the satellite
r_b= [r*cos(v);r*sin(v);0] ; %position vector in the orbit system
r_b1(k,:)=[a*(cos(E)-e);a*(sqrt(1-(e*e)))*sin(E);0]; %position vector in the orbit system
v_b(k,:)=[-a*n*sin(E)*a/r;a*n*sqrt(1-(e*e))*cos(E)*a/r]; %velocity vector in the orbit system
r_3_omega_1= [cos(-omega_1) sin(-omega_1) 0;-sin(-omega_1) cos(-omega_1) 0;0 0 1];
r_3_omega_2= [cos(-omega_2) sin(-omega_2) 0;-sin(-omega_2) cos(-omega_2) 0;0 0 1];
r_1=[1 0 0;0 cos(-eta) sin(-eta);0 -sin(-eta) cos(-eta)];
r_i= r_3_omega_1*r_1*r_3_omega_2*r_b;
r_i1(k,:)=transpose(r_i);
w_a=((2*pi)/86164);
phi=w_a*(i-0);
r_phi= [cos(phi) sin(phi) 0;-sin(phi) cos(phi) 0;0 0 1];
r_e=r_phi*r_i ;
r_e1(k,:)=r_phi*r_i ;
z_lat(k,:)= r_e1(k,3);
y_lat(k,:)=r_e1(k,2);
x_lat(k,:)=r_e1(k,1);
lat(k,:)= atan2((z_lat(k,:)),(sqrt((x_lat(k,:)*x_lat(k,:))+(y_lat(k,:)*y_lat(k,:)))));
long(k,:)= atan2(y_lat(k,:),x_lat(k,:)) ;
r_w=[4075.53022;931.78130;4801.61819];
lat_w(k,:)=atan2(r_w(3,1),(sqrt(r_w(1,1)*r_w(1,1)+r_w(2,1)*r_w(2,1))));
long_w(k,:)=atan2(r_w(2,1),r_w(1,1));
A = [-sin(lat_w(k,:))*cos(long_w(k,:)) -sin(long_w(k,:)) cos(lat_w(k,:))*cos(long_w(k,:));-sin(lat_w(k,:))*sin(long_w(k,:)) cos(long_w(k,:)) cos(lat_w(k,:))*sin(long_w(k,:));cos(lat_w(k,:)) 0 sin(lat_w(k,:))];
r_top = transpose(A)*(r_e-r_w);
r_top1(k,:)=transpose(A)*(r_e-r_w);
close=min(r_e-r_w,T)
S=sqrt(r_top(1,1).^2+r_top(2,1).^2+r_top(3,1).^2);
S_1(k,:)=sqrt(r_top1(k,1).^2+r_top1(k,2).^2+r_top1(k,3).^2);
zenit(k,:)=acos(r_top(3,1)/S);
Elevation(k,:)=(pi/2) -zenit(k,:);
Azi(k,:)= atan2(-r_top1(k,2),r_top1(k,1));
k=k+1;
end
% subplot(1,1,1);
figure()
plot3(r_e1(:,1),r_e1(:,2),r_e1(:,3));
save('fig(1)');
% subplot(1,1,2);
figure()
plot3(r_i1(:,1),r_i1(:,2),r_i1(:,3));
save('fig(2)');
%subplot(1,1,3)
figure()
pax=polaraxes;
pax.ThetaDir = 'clockwise';
polarplot(Azi,rad2deg(Elevation))
pax.RDir = 'reverse';
%fig4=polarplot(zenit,rad2deg(Elevation))
% subplot(1,1,4)
figure()
plot(rad2deg(long),rad2deg(lat))
hold on
load coastlines
geoshow(coastlat,coastlon)
hold off
save('fig(5)')
minimum_closest_distance=min(r_e-r_w) ;
minimum_time=(T),
end
%file second, run this file
a_gps=26500;
e_gps=0.01;
eta_gps=deg2rad(55);
omega_1_gps=deg2rad(0);
omega_2_gps=deg2rad(0);
[r_b,v_b,zenit,Elevation,Azi,close]= GPSex1original(a_gps,e_gps,eta_gps,omega_1_gps,omega_2_gps);
a_gnss=42164.140;
e_gnss=0;
eta_gnss=deg2rad(63);
omega_1_gnss=deg2rad(0);
omega_2_gnss=deg2rad(0);
[gnss_r_b,gnss_v_b,gnss_zenit,gnss_Elevation,gnss_Azi]= GPSex1original(a_gnss,e_gnss,eta_gnss,omega_1_gnss,omega_2_gnss);
a=29994;
e=0;
eta=deg2rad(56);
omega_1=deg2rad(0);
omega_2=deg2rad(0);
[galileo_r_b,galileo_v_b,galileo_zenit,galileo_Elevation,galileo_Azi]= GPSex1original(a,e,eta,omega_1,omega_2);
a=6836;
e=0.004;
eta=deg2rad(87);
omega_1=deg2rad(0);
omega_2=deg2rad(0);
[champ_r_b,champ_v_b,champ_zenit,champ_Elevation,champ_Azi]= GPSex1original(a,e,eta,omega_1,omega_2);
a=26554;
e=0.7;
eta=deg2rad(65);
omega_1=deg2rad(245);
omega_2=deg2rad(270);
[molniya_r_b,molniya_v_b,molniya_zenit,molniya_Elevation,molniya_Azi]= GPSex1original(a,e,eta,omega_1,omega_2);
% I have attached files in the attachment please look into it and help me.

Risposte (0)

Categorie

Scopri di più su CubeSat and Satellites in Help Center e File Exchange

Prodotti


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by