1.地球靜止軌道衛星，傾角分別為0，30，90度。

clc; clear;

t = 0:1:6;

we = 360/24;

u = we*t;

i = 30;

fai = asind( sind(i)*sind(u) );

deltalmd = atand( cosd(i)*tand(u) );

if(i==90)

deltalmd(end) = 90;

end

lmd = deltalmd - we*t;

% use symetry to generate the other data

for j = 1:6

lmd(j+7) = -lmd(7-j);

fai(j+7) = fai(7-j);

end

for j = 1:12

lmd(j+13) = lmd(13-j);

fai(j+13) = -fai(13-j);

end

h = geoshow('landareas.shp', 'FaceColor', [1 1 1]);

grid on

hold on

plot(lmd, fai); title(['GEO¹ìµÀ£¬Çã½Çi=', num2str(i)])

2.迴歸軌道衛星，迴歸週期1天，傾角分別為60度，週期為

clc; clear;

t = [0 1/3 1/2 2/3 4/5 1];

we = 360/24;

w = 180/2;

u = w*t;

i = 60;

fai = asind( sind(i)*sind(u) );

deltalmd = atand( cosd(i)*tand(u) );

lmd = deltalmd - we*t;

% use symetry to generate the other data

for j = 1:5

lmd(j+6) = lmd(6) + ( lmd(6) - lmd(6-j) );

fai(j+6) = fai(6-j);

end

for j = 1:10

if (lmd(11) + ( lmd(11) - lmd(11-j) )) > 180

lmd(j+11) = -180 + rem(lmd(11) + ( lmd(11) - lmd(11-j) ), 180);

else

lmd(j+11) = lmd(11) + ( lmd(11) - lmd(11-j) );

end

fai(j+11) = -fai(11-j);

end

cnt = 1;

for m = 1:5

for j = 1:21

if (lmd(j+21*(m-1)) + 60) > 180

lmd(j+21*m) = -180 + rem(lmd(j+21*(m-1)) + 60, 180);

record(m,cnt) = j; % record when tranverse from east to west

cnt = cnt + 1;

else

lmd(j+21*m) = lmd(j+21*(m-1)) + 60;

end

fai(j+21*m) = fai(j+21*(m-1));

end

cnt = 1;

end

load still

h = geoshow('landareas.shp', 'FaceColor', [1 1 1]);

grid on

hold on

plot(lmd1(2:20), fai1(2:20), 'b--'); % earth still

plot(lmd(1:6), fai(1:6), 'bo');

plot(lmd(21*6), fai(21*6), 'bo');

plot(lmd(1:13), fai(1:13)); plot(lmd(14:21), fai(14:21));

for m = 1:5

plot(lmd(21*m+1:record(m,1)+21*m-1), fai(21*m+1:record(m,1)+21*m-1)); plot(lmd(record(m,1)+21*m:21*(m+1)), fai(record(m,1)+21*m:21*(m+1)));

plot(lmd(21*m), fai(21*m), 'bo');

end

title(['ÐÇϵã¹ì¼££ºT=4h¹ìµÀ£¬Çã½Çi=', num2str(i)])

地球不轉時的星下點

clc; clear;
t = [0 1/3 1/2 2/3 4/5 1];
we = 360/24;
w = 180/2;
u = w*t;
i = 60;
fai = asind( sind(i)*sind(u) );
 
deltalmd = atand( cosd(i)*tand(u) );
lmd = deltalmd;    % earth still

% use symetry to generate the other data
for j = 1:5
    lmd(j+6) = lmd(6) + ( lmd(6) - lmd(6-j) );
    fai(j+6) = fai(6-j);
end
for j = 1:10
    if (lmd(11) + ( lmd(11) - lmd(11-j) )) > 180
        lmd(j+11) = -180 + rem(lmd(11) + ( lmd(11) - lmd(11-j) ), 180);
    else
        lmd(j+11) = lmd(11) + ( lmd(11) - lmd(11-j) );
    end
    fai(j+11) = -fai(11-j);
end
 
 
for j = 1:21
    if (lmd(j) + 180) > 180
        lmd(j) = -180 + rem(lmd(j) + 180, 180);
    else
        lmd(j) = lmd(j) + 180;
    end
    fai(j) = fai(j);
end
lmd(11) = 0;
 
lmd1 = lmd;
fai1 = fai;
save still lmd1 fai1
 
h = geoshow('landareas.shp', 'FaceColor', [1 1 1]);
grid on
hold on
plot(lmd(2:20), fai(2:20));
% plot(lmd(1:13), fai(1:13)); plot(lmd(14:21), fai(14:21));
% for m = 1:5
%     plot(lmd(21*m+1:record(m,1)+21*m-1), fai(21*m+1:record(m,1)+21*m-1)); plot(lmd(record(m,1)+21*m:21*(m+1)), fai(record(m,1)+21*m:21*(m+1)));
% end
title(['T=4h¹ìµÀ£¬Çã½Çi=', num2str(i)])