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 views (last 30 days)
Show older comments
%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.
0 Comments
Answers (0)
See Also
Categories
Find more on CubeSat and Satellites in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!