How can I measure the number of areas covering a point at one time?

3 次查看（过去 30 天）
Scott Hunter
Scott Hunter 2018-10-30
编辑： Scott Hunter 2018-10-30
This is quite an open question unfortunately because I have no clue where to begin. Essentially, I need to be able to tell how many satellites are providing coverage for a point on the surface of a planet at any given time. I have the satellites orbiting the planet plotted with time, along with the surface caps that each satellite can see. I have not yet found a way of producing the point on the surface of the planet (which rotates with time) that I need to be able to see the satellites from. The second thing I can't yet think of a way to do is measure how many of these surface caps cover my point at each timestep in my simulation.
The code I am using is shown below. The plotting section is at the end, the other code above is for anyone interested in running it.
I currently have the viewcones shown by edge colour, however would like to eventually have them invisible. Plotting all the satellites also isn't crucial so I have only plotted the lead satellite of each orbit and their viewcones at the moment.
Thanks lot for any help or ideas!
clear
close
%%%Things to consider:
%Satellite Constellation Idea 1
%MARSMAP
%Select Timestep
timestep = 600;
% Timescale
% Time will need calculating and conversion involving mean julian day
% for now using arbritrary values of t for example
%--------------------------------------------------------------------------
% General data
G = 6.67408*10^-11; %[m^3kg^-1s^-2]
% Mars data
Mars.m = 6.39*10^23; %[kg]
Mars.mu = Mars.m*G; %[m^3s^-2]
Mars.r = 3389.5*10^3; %[m]
Mars.v = 241.17; %[ms^-1]
Mars.T = (2*pi*Mars.r)/Mars.v; %[s]
Mars.o = Mars.v/Mars.r; %[rads^-1]
Mars.th = Mars.o*timestep; %[rad] for every second
t = (0:timestep:Mars.T); %Mars.T, Length of martian day in seconds
% Orbit data
% NOTE: All orbits are circular and equal to each other therefore a=r
% except RAAN, each orbit is spaced by 45 deg
% Orbits
i = deg2rad([90 90 90 90]); %[rad] Inclinations for each orbit
H = 9448.61*10^3; %[m] Altitude of each orbit
a = Mars.r + H; %[m] Semi-major axis of each orbit
o = deg2rad([0 0 0 0]); %[rad] Small omega of each orbit
e = [0 0 0 0]; % Eccentricity of each orbit
T = 2*pi*sqrt(a.^3/Mars.mu); %[s] Period of each orbit
n = (2*pi)./T; %[kms^-1] Mean motion of each orbit
h = sqrt(Mars.mu*a); %[kgm^2s^-1] Angular momentum on each orbit
v = sqrt(Mars.mu*(1./a)); %[ms^-1] Velocity of sats on each orbit
%Since it's circular M = n(t-tau)
% is actually just M = n(t)
M = n*t'; %[rad] Mean anomaly of each orbit
O = deg2rad([0 45 90 135]); %[rad] RAANs for each orbit
%--------------------------------------------------------------------------
%IN-PLANE ORBITAL PARAMETERS
%Set up satellite arrays (Equivalent to writing for Sat 1 through 32)
th0 = zeros(1,32); %Initial angle of each satellite for each timestep
%This step is essentially filling th0 array with values depending upon the
%orbit
for j = 1:8
th0(:,j) = deg2rad(45*(j-1)); %For orbit 1, sats start at 0, spaced 45 deg
end
for j = 9:16
th0(:,j) = deg2rad(22.5+(45*(j-9))); %For orbit 2, sats start at 22.5, spaced 45
end
for j = 17:24
th0(:,j) = deg2rad(45*(j-17)); %For orbit 3, sats start at 0, spaced 45 deg
end
for j = 25:32
th0(:,j) = deg2rad(22.5+(45*(j-25))); %For orbit 4, sats start at 22.5, spaced 45
end
th = th0 + 2*atan(tan(M/2)); %[rad] Angle of each sat at each timestep
xp = a*cos(th); %[m] x pos of each sat at each timestep
subx = Mars.r*cos(th); %[m] sub sat x pos
yp = a*sin(th); %[m] y pos of each sat at each timestep
suby = Mars.r*sin(th); %[m] sub sat y pos
zp = xp*0; %[m] z pos is 0 as it is in plane motion
subz = subx*0;
xp_dot = -sqrt(Mars.mu/a)*sin(th); %[ms^-1] x vel of each sat at each timestep
yp_dot = sqrt(Mars.mu/a)*cos(th); %[ms^-1] y vel of each sat at each timestep
zp_dot = xp_dot*0; % z vel is 0 as it is in plane motion
pos = [xp yp zp]'; % pos matrix (takes form of [X Y Z]') of each sat at each timestep
sub = [subx suby subz]'; % sub sat pos matrix of each sat at each timestep
vel = [xp_dot yp_dot zp_dot]'; % vel matrix of each sat at each timestep
%--------------------------------------------------------------------------
%TRANSFORMATION FROM KEP TO CART
%Rotation Matrix for orbits
R1 = [cos(o).*cos(O)-sin(o).*cos(i).*sin(O) % Column 1 of rotation matrix for the 4 orbits
cos(o).*sin(O)+sin(o).*cos(i).*cos(O)
sin(o).*sin(i)];
R2 = [-sin(o).*cos(O)-cos(o).*cos(i).*sin(O) % Column 1 of rotation matrix for the 4 orbits
-sin(o).*sin(O)+cos(o).*cos(i).*cos(O)
cos(o).*sin(i)];
R3 = [sin(i).*sin(O) % Column 1 of rotation matrix for the 4 orbits
-sin(i).*cos(O)
cos(i)];
R = [R1 R2 R3];
%Select individual rotation matrices for each orbit, relates to the columns
%of R
Ro1 = [R(:,1) R(:,5) R(:,9)];
Ro2 = [R(:,2) R(:,6) R(:,10)];
Ro3 = [R(:,3) R(:,7) R(:,11)];
Ro4 = [R(:,4) R(:,8) R(:,12)];
%Covert from Kep position matrices to Cart position matrices
%This can be done for sub sat point and velocity matrices in similar ways
%if required
POS = zeros(96,148);
SUB = zeros(96,148);
%Fill array with values depending on orbit
for j = 1:8
POS([j,j+32,j+64],:)= Ro1*[pos(j,:); pos(j+32,:); pos(j+64,:);];
SUB([j,j+32,j+64],:)= Ro1*[sub(j,:); sub(j+32,:); sub(j+64,:);];
end
for j = 9:16
POS([j,j+32,j+64],:)= Ro2*[pos(j,:); pos(j+32,:); pos(j+64,:);];
SUB([j,j+32,j+64],:)= Ro2*[sub(j,:); sub(j+32,:); sub(j+64,:);];
end
for j = 17:24
POS([j,j+32,j+64],:)= Ro3*[pos(j,:); pos(j+32,:); pos(j+64,:);];
SUB([j,j+32,j+64],:)= Ro3*[sub(j,:); sub(j+32,:); sub(j+64,:);];
end
for j = 25:32
POS([j,j+32,j+64],:)= Ro4*[pos(j,:); pos(j+32,:); pos(j+64,:);];
SUB([j,j+32,j+64],:)= Ro4*[sub(j,:); sub(j+32,:); sub(j+64,:);];
end
%--------------------------------------------------------------------------
%PLOT
%Create Mars as a surface
marsterrain = imread('mars.jpg');
[X, Y, Z] = sphere(20);
mars = mesh(X*Mars.r, Y*Mars.r, Z*Mars.r);
set(mars,'FaceColor','texturemap','cdata',marsterrain,'edgecolor','none');
hold on
%Plot viewcones as spherical caps
%Each cone is 'cut' at a distance of 1507.09Km from the centre of mars
%(used hand calcs to figure out), call d
cap = zeros(1,32);
d = 1507.09*10^3;
[x, y, z] = sphere(20);
xcap = Mars.r*x;
ycap = Mars.r*y;
zcap = Mars.r*z;
zcap(zcap < d) = NaN;
%Each cap has an initial position rotated by the initial angle of it's
%satellite. Since the first satellite starts at 0 on the y and z axis i.e
%on the equator at 0 deg latitude, and the cap is created on the top of the sphere,
%the first one is rotated by 90 deg to begin with. Call it thcap0. thcap
%will be the constant rotation of the cap per second, i.e angular speed of
%the satellites. The axis about which each cap rotates is
%rotated by the RAAN initially. The vector direction which represents the
%RAAN(square brackets in rotate function) simply needs to create the angle
%given by the RAAN Anticlockwise.
thcap0 = zeros(1,32);
thcap = th(1,1) - th(2,1);
for j = 1
thcap0(:,j) = (pi/4)+th0(:,j);
cap(:,j) = surf(xcap, ycap, zcap, 'FaceColor', 'none', 'edgecolor', 'r');
rotate(cap(:,j), [0 1 0], rad2deg(thcap0(:,j)));
end
for j = 9
thcap0(:,j) = (pi/4)+th0(:,j);
cap(:,j) = surf(xcap, ycap, zcap, 'FaceColor', 'none', 'edgecolor', 'g');
rotate(cap(:,j), [-0.5 0.5 0], rad2deg(thcap0(:,j)));
end
for j = 17
thcap0(:,j) = (pi/4)+th0(:,j);
cap(:,j) = surf(xcap, ycap, zcap, 'FaceColor', 'none', 'edgecolor', 'b');
rotate(cap(:,j), [-1 0 0], rad2deg(thcap0(:,j)));
end
for j = 25
thcap0(:,j) = (pi/4)+th0(:,j);
cap(:,j) = surf(xcap, ycap, zcap, 'FaceColor', 'none', 'edgecolor', 'y');
rotate(cap(:,j), [-0.5 -0.5 0], rad2deg(thcap0(:,j)));
end
hold on
%Draw Orbits
Orbit1 = animatedline(POS(1,1), POS(33,1), POS(65,1), 'color', 'r');
Orbit2 = animatedline(POS(9,1), POS(41,1), POS(73,1), 'color', 'g');
Orbit3 = animatedline(POS(17,1), POS(49,1), POS(81,1), 'color', 'b');
Orbit4 = animatedline(POS(25,1), POS(57,1), POS(89,1), 'color', 'y');
%Set axis limits
axis equal
xlim([-a, a]);
ylim([-a, a]);
zlim([-a, a]);
xlabel('x');
ylabel('y');
zlabel('z');
set(gcf, 'color', [0 0 0])
axis off
hold on
%Plot over time
for j = 1:length(t)
%Draw orbits from leading satellite
addpoints(Orbit1, POS(1,j), POS(33,j), POS(65,j));
addpoints(Orbit2, POS(9,j), POS(41,j), POS(73,j));
addpoints(Orbit3, POS(17,j), POS(49,j), POS(81,j));
addpoints(Orbit4, POS(25,j), POS(57,j), POS(89,j));
%Plot satellites (currently only plotting lead sats due to lag)
S = zeros(1,32);
for k = 1
S(:,k) = scatter3(POS(k,j), POS(k+32,j), POS(k+64,j), 'Filled', 'MarkerFaceColor', 'r');
end
for k = 9
S(:,k) = scatter3(POS(k,j), POS(k+32,j), POS(k+64,j), 'Filled', 'MarkerFaceColor', 'g');
end
for k = 17
S(:,k) = scatter3(POS(k,j), POS(k+32,j), POS(k+64,j), 'Filled', 'MarkerFaceColor', 'b');
end
for k = 25
S(:,k) = scatter3(POS(k,j), POS(k+32,j), POS(k+64,j), 'Filled', 'MarkerFaceColor', 'y');
end
drawnow
%Rotate Mars
rotate(mars, [0 0 1], (rad2deg(Mars.th)));
%Rotate Caps
for k = 1
rotate(cap(:,k), [0 1 0], (rad2deg(thcap)));
end
for k = 9
rotate(cap(:,k), [-0.5 0.5 0], (rad2deg(thcap)));
end
for k = 17
rotate(cap(:,k), [-1 0 0], (rad2deg(thcap)));
end
for k = 25
rotate(cap(:,k), [-0.5 -0.5 0], (rad2deg(thcap)));
end
%Delete satellite from current frame
delete(S)
end

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Satellite Mission Analysis 的更多信息

标签

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by