clc
clear
clf
% Physical constants
mu = 398600.4418; % km^3/s^2
R_E = 6378.137; % km
omega_E = 7.2921158553e-5; % rad/s
nu = 0;
fpa = 0;
% Satellite parameters
e = 0; % Eccentricity
i = 30; % Inclination [deg]
RAAN = 20; % RAAN [deg]
n = sqrt(mu/R_E^3); % Mean motion [rad/s]
h = R_E*sqrt(1-e^2); % Specific angular momentum
alt= 35786; %km
r0m = R_E + alt %semi major axis
v0m = 3.11;
En=v0m^2/2-mu/r0m;
a=-mu/(2*En)
r0=a*[cosd(nu) sind(nu) 0];
v0=v0m*[cosd(nu-90+fpa),sind(nu-90+fpa), 0 ];
y0=[r0 v0];
Period=2*pi*sqrt(a^3/mu);
fprintf(1,'Period = %6.0fsecs, %6.1fmins, %6.2fhrs\n',Period,Period/60,Period/3600)
t0=0;
nrevs = 100;
tmax = nrevs*Period;
% Sensor parameters
nadir_angle = 4; % [deg]
nadir_angle_rad = deg2rad(nadir_angle);
distance_to_earth = R_E + (a*sin(nadir_angle_rad)); % Distance to the point where the sensor line of sight intersects the earth
% NMW Inertial coordinate systemclc
x_NMW = [cosd(RAAN), sind(RAAN), 0]; % x-axis
z_NMW = [0, 0, 1]; % z-axis
y_NMW = cross(z_NMW, x_NMW); % y-axis
rref=a; %Set tolerances relative to r0
vref=v0m; % and v0
reltol=(1e-6); %Sets tolerance
abstol=reltol*[rref rref rref vref vref vref]; %this make each tol relate to r0m, v0m
settings=odeset('RelTol',reltol,'AbsTol',abstol); %sets relative and absolute tolerances
[~,Y]=ode45(@twoBody,[0 tmax],y0,settings);
Yend=Y(end,:);
tolset=1e-6;
for k=1:length(tolset)
fprintf(1,'Tolerance= %4.0e',tolset(k))
[~,Y]=ode45(@twoBody,[0 tmax],y0,settings);
r_NMWt=Y(:,1:3)' ;
v_NMWt=Y(:,4:6)';
end
r0_NMW=[a, 0, 0];
v0_NMW=[0, sqrt(mu)/a, 0];
figure
hold off
[ex,ey,ez]=sphere(24);surf(R_E*ex,R_E*ey,R_E*ez)
hold on; axis equal;
rplot=8000*1.1; rlabel=rplot*1.1; %Dimension for placing axes labels
sizec=8000*1.1; %Sizes length of axes relative to semimajor axis
axis([-sizec sizec -sizec sizec -sizec sizec]); %Sizes the plot space -x +x -y +y -z +z
plot3([0 rplot],[0 0],[0 0],'k','LineWidth',2);text(rlabel,0,0,'X NMW km','FontWeight','bold')%x axis
plot3([0,0],[0,rplot],[0,0],'k','LineWidth',2);text(0,rlabel,0,'Y NMW km','FontWeight','bold')%y axis
plot3([0,0],[0,0],[0,rplot],'k','LineWidth',2);text(0,0,rlabel,'Z NMW km','FontWeight','bold')%z axis
hold on
plot3(r_NMWt(1,:),r_NMWt(2,:),r_NMWt(3,:),'- g','LineWidth',2,'MarkerSize',1);
axis equal
hold off
grid on
% Line-of-sight point of intersection with the Earth as the satellite makes one orbit revolution
theta = linspace(0, 2*pi, 1000);
x = distance_to_earth*cos(theta);
y = (acosd(i) + distance_to_earth*sin(nadir_angle_rad))*sin(theta);
z = (asind(i))*sin(theta);
% Earth Centered Inertial (ECI) frame
t = linspace(0, 2*pi/n, 1000);
r = [x', y', z'];
n_vec = repmat([0, 0, n], size(r, 1), 1);
v = cross(n_vec, r);
eci_states = [r, v];
figure
hold off
[ex,ey,ez]=sphere(24);surf(R_E*ex,R_E*ey,R_E*ez)
hold on; axis equal;
rplot=8000*1.1; rlabel=rplot*1.1; %Dimension for placing axes labels
sizec=8000*1.1; %Sizes length of axes relative to semimajor axis
axis([-sizec sizec -sizec sizec -sizec sizec]); %Sizes the plot space -x +x -y +y -z +z
plot3([0 rplot],[0 0],[0 0],'k','LineWidth',2);text(rlabel,0,0,'X ECI km','FontWeight','bold')%x axis
plot3([0,0],[0,rplot],[0,0],'k','LineWidth',2);text(0,rlabel,0,'Y ECI km','FontWeight','bold')%y axis
plot3([0,0],[0,0],[0,rplot],'k','LineWidth',2);text(0,0,rlabel,'Z ECI km','FontWeight','bold')%z axis
hold on
%ECI_path = eci_states(1:3,:);
ECI_path = eci_states(:,1:3);
plot3(ECI_path(:,1), ECI_path(:,2), ECI_path(:,3), '-o','MarkerSize',1,'linewidth',1)
xlabel('x (km)')
ylabel('y (km)')
zlabel('z (km)')
title('Path of the satellite in ECI frame')
grid on
axis equal
hold off
% Argument of Latitude
lat_30 = 30;
u = atan2d(tand(lat_30), cosd(RAAN));
fprintf('Argument of latitude: %.2f degrees\n', u);
% Earth Centered Earth Fixed (ECEF) frame
GMST = @(t) mod(omega_E*t, 2*pi);
ecef_states = zeros(size(eci_states));
for i = 1:size(eci_states, 1)
T = [cos(GMST(t(i))), sin(GMST(t(i))), 0;
-sin(GMST(t(i))), cos(GMST(t(i))), 0;
0, 0, 1];
ecef_states(i, :) = eci_states(i, :)';
end
ECEF_path = ecef_states(:, 1:3);
figure
hold off
[ex,ey,ez]=sphere(24);surf(R_E*ex,R_E*ey,R_E*ez)
hold on; axis equal;
rplot=8000*1.1; rlabel=rplot*1.1;
sizec=8000*1.1;
axis([-sizec sizec -sizec sizec -sizec sizec]);
plot3([0 rplot],[0 0],[0 0],'k','LineWidth',2);text(rlabel,0,0,'X ECEF km','FontWeight','bold')%x axis
plot3([0,0],[0,rplot],[0,0],'k','LineWidth',2);text(0,rlabel,0,'Y ECEF km','FontWeight','bold')%y axis
plot3([0,0],[0,0],[0,rplot],'k','LineWidth',2);text(0,0,rlabel,'Z ECEF km','FontWeight','bold')%z axi
plot3(ECEF_path(1,:), ECEF_path(2,:), ECEF_path(3,:))
xlabel('Longitude (deg)')
ylabel('Latitude (deg)')
title('Ground track of the satellite in ECEF frame')
grid on
axis equal
hold on
plot(u, lat_30, 'r')
axis equal
hold off
legend('Ground track', 'Spot beam')
grid on
% Part of the world where the spot beam is hitting
fprintf('The spot beam is hitting the part of the world near %.2f degrees N, %.2f degrees E.\n', lat_30, u);
