MATLAB Code | ece6390-tesseract

%plot satellite coverage

clear all;

close all;

clc;

%========= initialize constants===================

earthrad = 6371; % radius of earth in kilometers

fov = 28; % field of view of imaging equipment(degrees)

Timage = 13 ; % seconds in between images

xpixels = 10328;

ypixels = 7760;

numsats = 25;

gshpbw = 16;

latgs = 0;

longs = 0;

%========= Plot Satellite Track for Sat 1==========================

%---- load sat track----------------

load Sentinel5p.mat;

latsubpoint1 = esplat;

latsubpoint1 = interp1(1:length(latsubpoint1),latsubpoint1,1:.01:length(latsubpoint1));

lonsubpoint1 = esplon;

indexcorrections = find(abs(lonsubpoint1(1:length(lonsubpoint1)-1)-lonsubpoint1(2:length(lonsubpoint1))) >5);

for m = 1:length(indexcorrections)

lonsubpoint1(indexcorrections(m)+1:end) = lonsubpoint1(indexcorrections(m)+1:end) - 360;

end

lonsubpoint1(find(lonsubpoint1 < 0)) = lonsubpoint1(find(lonsubpoint1 < 0))+360;

lonsubpoint1 = interp1(1:length(lonsubpoint1),lonsubpoint1,1:.01:length(lonsubpoint1));

altsat1 = altitude; % altitude of satellite in km

clear esplat esplon altsat;

%------------------------------------

%----------- calculate longitude shift between satellites ------------

Torbit = 2*pi*sqrt(((altsat1+earthrad))^3/(398600.44)); % period of orbit in seconds

londeg = (Torbit/86164.09)*360;

lonoffset = londeg/numsats;

%---------------------------------------------------------------------

%========== calculate image window radius ================

windowradius = altsat1*tand(fov/2); % radius of fov of view on earth's surface in km

resolution = (2*windowradius*1000)/sqrt(xpixels^2 + ypixels^2); % resolution of imagery (meters/pixel)

%===================================================

datapoints = length(latsubpoint1); % total number of points in one sidereal day

secperpoint = 86164.09/datapoints; % calculate seconds between position measurements

Timagesindex = 1:floor(Timage/secperpoint):datapoints; % how many data points between appropriately spaced images

figure1 = figure();

worldmap world;

load coast;

plotm(lat,long);

hold on;

%-------------------- determine sat positions that are in the sun light

daytime = zeros(1,datapoints); % initialize logical indicator of satellite being of daylit area

db = 360; % initialize position of day break longitude

nf = 360 - .5*360; % initilize position of night fall longitude

increment = 360/datapoints;

for d = 1:datapoints

position = 180+rem((lonsubpoint1(d) + 180),360);

if db > nf && position > nf && position < db

daytime(d) = 1;

end

if db < nf && (position > nf || position < db)

daytime(d) = 1;

end

db = db - increment;

nf = nf - increment;

if nf < -180

nf = nf + 360;

end

if db < -180

db = db + 180;

end

%-------------------------------------------------------------------

for s = 1:numsats

plotm( latsubpoint1,lonsubpoint1+s*lonoffset,'m-'); % total track

%plotm( latsubpoint1(find(daytime)),lonsubpoint1(find(daytime))+s*lonoffset,'r.'); % daylit track

end

title('Satellite Track');

clear lat long ;

hold on;

%=--------- plot ground station coverage ----------------------

% gswindowradius = altitude*tand(gshpbw/2);

% angle = linspace(0,359,360); % for circular window

% xmap = gswindowradius*cosd(angle);

% ymap = gswindowradius*sind(angle);

% latgs = latgs + (180/pi) * (ymap./2)./earthrad;

% longs = longs + (180/pi) * (xmap./2)./ (earthrad * cosd(latgs));

% plotm(latgs,longs,'red-');

%---------------------------------------------------------------

%--------- Plot circles around Image locations-------------------

imagelats = latsubpoint1(Timagesindex);

imagelons = lonsubpoint1(Timagesindex);

for b = 1:numsats % all sats

for n = 2000:2025

% angle = linspace(0,359,360); % for circular window

% xmap = windowradius*cosd(angle);

% ymap = windowradius*sind(angle);

xmap = [ones(1,10)*windowradius*cosd(atand(1536/2048)),linspace(windowradius*cosd(atand(1536/2048)),-1*windowradius*cosd(atand(1536/2048)),10) ,-1*ones(1,10)*windowradius*cosd(atand(1536/2048)) ,linspace(windowradius*cosd(atand(1536/2048)),-1*windowradius*cosd(atand(1536/2048)),10)];

ymap = [linspace(windowradius*sind(atand(1536/2048)),-1*windowradius*sind(atand(1536/2048)),10) ,ones(1,10)*windowradius*sind(atand(1536/2048)),linspace(windowradius*sind(atand(1536/2048)),-1*windowradius*sind(atand(1536/2048)),10),-1*ones(1,10)*windowradius*sind(atand(1536/2048)) ];

latmap = imagelats(n) + (180/pi) * (ymap./2)./earthrad;

lonmap = imagelons(n) + b*lonoffset + (180/pi) * (xmap./2)./ (earthrad * cosd(imagelats(n)));

%figure(1)

%plotm(imagelats(n),imagelons(n),'r*');

%hold on;

plotm(latmap,lonmap,'black-');

end

%-----------------------------------------------------------------

%====================================================================

%=================== Clear Workspace =================================

clear all;

close all;

clc;

%=====================================================================

%=================== Initialization ==================================

%-------------- Constants ------------------------------------

earthrad = 6371E3; % earth's radius in m

G = 6.67384E-11; % gravitational constant (N*m^2/kg^2)

Mp = 5.9722E24; % Mass of earth in kg

c = 2.99E8;

%-------------------- Input Variables-----------------------------

days = 1; % set number of simulation days

delt = 20; % time step in seconds

steps = round(days*86164.09/delt); %total number of timesteps

altitude = 400; % altitude of satellite in km

dperigee = earthrad + altitude*1000; % Distance at apogee (m)

dapogee = earthrad + altitude*1000; % Distance at perigee (m)

esplat0 = 85; % initial earth sub point lat

esplon0 = 0; % initial earth sub point lon

a = (dperigee + dapogee)/2;

e = (a-dperigee)/(a);

r0 = dapogee; % initial radial distance from center of earth at apogee (m)

vr0 = 0; % Initial radial velocity at apogee is zero (m/sec)

absolutevel0 = sqrt((G*Mp*(1-e))/(a*(1+e))); % absolute initial velocity at apogee

theta0 = 0; % Initial angular position (degrees)

vtheta0 = absolutevel0/dapogee; % Initial angular velocity (radians/sec)

atllat0 = 33.7758; % latitude of Van Leer

atllon0 = -84.39738; % Longitude of Van Leer

%---------------Fixed Parameters and Initial Calculations---------------

r(1:steps) = 0;

r(1) = r0;

vr(1:steps) = 0;

vr(1) = vr0;

theta(1:steps) = 0;

theta(1) = theta0;

vtheta(1:steps) = 0;

vtheta(1) = vtheta0;

delr(1:steps) = 0;

delr(1) = vr(1)*delt;

deltheta(1:steps) = 0;

deltheta(1) = vtheta0*delt;

satx(1:steps) = 0;

saty(1:steps) = 0;

satz(1:steps) = 0;

atlx(1:steps) = 0;

atly(1:steps) = 0;

atlz(1:steps) = 0;

%----------- Initialize Earth -----------------------------

earthx = earthrad*cosd(1:361); % initialize equator

earthy = earthrad*sind(1:361); % initialize equator

earthz(1:361) = 0; % initialize equator

earthx = [earthx earthrad*cosd(1:361)]; % initialize prime meridian

earthy = [earthy zeros(1, 361)]; % initialize prime meridian

earthz = [earthz earthrad * sind(1:361)]; % initialize prime meridian

earthx = [earthx earthrad*cosd(1:90)]; % initialize equator

earthy = [earthy earthrad*sind(1:90)]; % initialize equator

earthz = [earthz zeros(1,90)]; % initialize equator

earthx = [earthx zeros(1,361)]; % initialize lon = +/- 90deg

earthy = [earthy (earthrad * cosd(1:361))]; % initialize lon = +/- 90deg

earthz = [earthz (earthrad * sind(1:361))]; % initialize lon = +/- 90deg

espx(1:steps) = 0;

espy(1:steps) = 0;

espz(1:steps) = 0;

earthzerolongitude(1:steps) = 0;

load coast;

coastx = (earthrad*cosd(long).*sind(lat+90))';

coasty = (earthrad*sind(long).*sind(lat + 90))';

coastz = (earthrad*cosd(90 - lat))';

atlx(1) = (earthrad*cosd(atllon0).*sind(atllat0+90))'; % x position of Atlanta time zero

atly(1) = (earthrad*sind(atllon0).*sind(atllat0 + 90))'; % y position of atlanta time zero

atlz(1) = (earthrad*cosd(90 - atllat0))'; % z position of altanta time zero

R = rotz(360*days/steps); % proper earth rotation for each time step (degrees)

%------- Rotate Earth to Proper Initial angle -------------------

earthzerolongitude(1) = -1*esplon0;

R0 = rotz( -1*esplon0); % initial rotation

XYZearth = [earthx;earthy;earthz];

XYZcoast = [coastx;coasty;coastz];

XYZatl = [atlx(1);atly(1);atlz(1)];

XYZearth = R0*XYZearth;

XYZcoast = R0*XYZcoast;

XYZatl = R0*XYZatl;

earthx = XYZearth(1,:);

earthy = XYZearth(2,:);

earthz = XYZearth(3,:);

coastx = XYZcoast(1,:);

coasty = XYZcoast(2,:);

coastz = XYZcoast(3,:);

atlx(1) = XYZatl(1,:);

atly(1) = XYZatl(2,:);

atlz(1) = XYZatl(3,:);

%-----------------------------------------------------------------------

%======================================================================

%==================== Main Iterative Loop =============================

mov = avifile('hw1.avi'); % Open Movie file

for n = 1:steps

if n ~= steps

r(n+1) = r(n) + delr(n); % update satellite position

%theta(n+1) = theta(n) - deltheta(n); %retrograde orbit

theta(n+1) = theta(n) + deltheta(n); %non-retrograde orbit

delr(n+1) = delr(n) + (r(n)+ 0.5*delr(n))*(deltheta(n))^2 -(G*Mp*delt^2)/(r(n))^2;

deltheta(n+1) = deltheta(n) - (2*delr(n)*deltheta(n))/(r(n) + 0.5*delr(n));

XYZatl = [atlx(n);atly(n);atlz(n)];

XYZatl = R*XYZatl; % rotates atlanta with earth revolution

atlx(n+1) = XYZatl(1);

atly(n+1) = XYZatl(2);

atlz(n+1) = XYZatl(3);

end

satx(n) = r(n)*cos(theta(n)); % calculates x,y,z coordinates of satellite before inclination

saty(n) = r(n)*sin(theta(n));

satz(n) = 0;

Rinc = roty(-1*esplat0); % rotate xyz position of satellite around y axis to initial latitude

XYZsat = [satx(n);saty(n);satz(n)];

XYZsat = Rinc*XYZsat;

satx(n) = XYZsat(1,1); % calculates x, y, z coordinates of satellite

saty(n) = XYZsat(2,1);

satz(n) = XYZsat(3,1);

espx(n) = earthrad * (satx(n)/sqrt(satx(n)^2 + saty(n)^2 + satz(n)^2)); % calculate x,y,z satellite sub point

espy(n) = earthrad * (saty(n)/sqrt(satx(n)^2 + saty(n)^2 + satz(n)^2));

espz(n) = earthrad * (satz(n)/sqrt(satx(n)^2 + saty(n)^2 + satz(n)^2));

if n~= 1

earthzerolongitude(n) = earthzerolongitude(n-1) + 360*days/steps ;

XYZearth = [earthx;earthy;earthz];

XYZcoast = [coastx;coasty;coastz];

XYZearth = R*XYZearth; % rotates equator and prime meridian with earth revolution

XYZcoast = R*XYZcoast; % rotates coast with earth revolution

earthx = XYZearth(1,:);

earthy = XYZearth(2,:);

earthz = XYZearth(3,:);

coastx = XYZcoast(1,:);

coasty = XYZcoast(2,:);

coastz = XYZcoast(3,:);

end

figure1 = figure();

plot3(earthx,earthy,earthz,'m:');

xlim([-1E7, 1E7]);

ylim([-1E7, 1E7]);

zlim([-1E7, 1E7]);

% xlim([-10000, 10000]);

% ylim([-10000, 10000]);

% zlim([-10000, 10000]);

hold on;

plot3(espx(n), espy(n),espz(n),'r*'); % plots red '*' on earth subpoint

plot3(coastx,coasty,coastz,'b-'); % plot coast

plot3(satx(n),saty(n),satz(n),'b*'); % plots blue '*' on satellite

plot3([0 satx(n)],[0 saty(n)],[0 satz(n)],'b:') % plots dotted line between center of earth and satellite

plot3(atlx(n), atly(n), atlz(n),'g*'); % plots red '*' at the location of atlanta

xlabel('x');

ylabel('y');

zlabel('z');

f = getframe(gcf);

mov = addframe(mov,f); % adds current updated frame to movie file

end

mov = close(mov); % close movie file

%======================================================================

%------------ Calculate Earth Sub Point in theta, phi and then lat lon

NH = find(espz >0);

esplon = atand(espy./espx) - earthzerolongitude ;

esplon(NH) = esplon(NH) + 180;

esplat = atand(espz./sqrt(espx.^2 + espy.^2));

%================ Calculate Azmith and Elevation from Atlanta =========

esplon = rem(esplon,360);

esplon(find(esplon<-180)) = esplon(find(esplon<-180)) + 360;

sigma = acosd((satx.*atlx + saty.*atly + satz.*atlz)./(sqrt(satx.^2 + saty.^2 + satz.^2).*sqrt(atlx.^2 + atly.^2 + atlz.^2)));

alpha = asind(sind(abs(atllon0-esplon)).*((cosd(esplat))./(sind(sigma))));

XYZesp = [espx;espy;espz];

for m = 1:steps

Rz = rotz(-1*rem(earthzerolongitude(m),360)); %align zerolongitude with x = 0 axis

XYZesp(:,m) = Rz*XYZesp(:,m);

end

Rz = rotz(-1*atllon0); %align atlanta longitude with x = 0

XYZesp = Rz*XYZesp;

Ry = roty(atllat0); %align atlanta latitude with x = 0

XYZesp = Ry*XYZesp;

x = XYZesp(1,:);

y = XYZesp(2,:);

z = XYZesp(3,:);

NE = find(y >= 0 & z >= 0);

NW = find(y < 0 & z > 0);

SE = find(y > 0 & z <0);

SW = find(y <= 0 & z <=0);

alpha(NE) = alpha(NE); % no change

alpha(SE) = 180 - alpha(SE);%

alpha(SW) = 180 + alpha(SW);

alpha(NW) = 360 - alpha(NW);

El = acosd((sind(sigma))./(sqrt(1 + (earthrad./r).^2 - 2*(earthrad./r).*cosd(sigma))));

visible = find(sigma <= acosd(earthrad./r));

figure()

%plot(alpha,El,'r.');

%hold on;

plot(alpha(visible),El(visible),'b.');

ylim([0,90]);

xlim([0,360]);

xlabel('Azimuth (deg)');

ylabel('Elevation (deg)');

%======================================================================

%================Plot satellite position over time==================

figure1 = figure();

worldmap world;

load coast;

plotm(lat,long);

hold on;

plotm( esplat,esplon,'r-');

%plotm(esplat(visible),esplon(visible),'b.');

title('Satellite Track');

clear lat long ;

figure2 = figure();

plot3(earthx,earthy,earthz,'m:');

% xlim([-1E7, 1E7]);

% ylim([-1E7, 1E7]);

% zlim([-1E7, 1E7]);

xlim([-2E7, 7E7]);

ylim([-6E7, 3E7]);

zlim([-8E7, 1E7]);

hold on;

plot3(coastx,coasty,coastz,'b-'); % plot coast

plot3(satx,saty,satz,'m-');

plot3(satx(500),saty(500),satz(500),'m*'); % plots blue '*' on satellite

Tesseract Inc.