Welcome to our Mobius STEM page. This page contains information that will help you explore coding using Octave, an open source code. Follow the instructions to assess a make-believe constellation sensing an make-believe threat.

Octave Example Instructions

The open-source language “Octave” will run the code in this example. Download it from https://octave.org/download.

See if you can use Octave, code that Mobius wrote, and sample trajectories to assess the distance from satellites to missiles.

  1. Download and install Octave

  2. Create a New Folder on your C:\ Drive and label it “MobiusSpaceMissileSensor” without quotes

  3. Click on “Trajectory Data Download” below and save the file with name “trajectory_J2000PosVel” in the MobiusSpaceMissileSensor folder. Make sure these is no suffix on the file - it won’t load it has a .dat, .txt or any othe file suffix.

  4. Open Octave and Copy/Paste the “Mobius Octave Code” from below

  5. Copy the file Trajectory Data file location and paste it into the “Load Trajectory” section of the Octave code. Example: “C:\Users\[your user name]\Documents\MobiusSpaceMissileSensor\”

  6. Notice the notes and constants.  See if you can recognize the constants, equations, and functions.

  7.  Select “Run” in Octave to see the code run

Assumptions: spherical earth, circular satellite orbits, Keplerian gravity.

Hint: To manage computation time the analysis can be done at 10 second timesteps.

A given Walker constellation is defined by 50deg inclination, 60 satellites, 6 planes, and a Walker number spacing of 3. The constellation altitude is 700km Satellite constellation - Wikipedia.

Define the position of all the satellites in the constellation as a function of time in the J2000 Earth-Centered Inertial (ECI) frame Earth-centered inertial - Wikipedia.

Note: Depending on the phasing chosen the student may have slightly different results than the ones shown below.

The 3-D rotation matrices are provided Rotation Matrix -- from Wolfram MathWorld. In the J2000 ECI frame, satellite constellation inclination can be represented by rotation about the x-axis, and the right ascension of the ascending node for the angle between the planes Longitude of the ascending node - Wikipedia can be represented by rotation about the z-axis.

A missile trajectory is also provided in the J2000 ECI frame. [see data below]

Plot the satellite constellation orbits, and the missile trajectory data, on a sphere representing the earth.

Generate a figure showing the distance of each satellite to the missile trajectory at each time step. Set the y-axis range to 0 to 5000km to exclude the satellites that would be on the opposite side of the earth from the missile. Include a horizontal line for 1500km range to allow for visual inspection of how often the target is within the sensor range.

The distance equation is provided Distance - Wikipedia, where (x1, y1, z1) is the position of the missile in J2000 ECI coordinates and (x2, y2, z2) is the position of a satellite in J2000 ECI coordinates.

The Distance Equation

Rotation Matrix

Trajectory Plot

Distance Plot - from satellites to object

 

Mobius Code

% Alison Yu, Mobius, 2023

%% Initial Code

close all

clear

%%%%%%%%%%%%%%%% Inputs %%%%%%%%%%%%%%%%

% Constellation inputs

conType     = 'Walker';

incl        = 50*pi/180;    %rad

tot_sats    = 6*10; % the total number of satellites = orbits x satellites per orbit

num_orbs    = 6;    % number of orbits

f_spacing   = 3;

alt         = 700e3;        % meters

num_sats    = tot_sats/num_orbs; % number of satellites per orbit

% Additional model inputs

traj        = 'trajectory_J2000PosVel'

sensorDist  = 1500e3; %m

tstep       = 10;     %sec

% Constants

G    = 6.67e-11; %N m^2 / kg^2, gravitational constant

M    = 5.972e24; %kg, mass of earth

Re   = 6371e3;   %m, radius of earth

% Plot options

markers = {'o','+',   'd','x',   's','<',   'h','*',   'p','^',   'v','>'};

    colors  = {...

    [0    0   0   ]   % 1 black

    [0.5  0   0   ]   % 2

    [1    0   0   ]   % 3 red

    [1    0.5 0   ]   % 4 orange

    [1    1   0   ]   % 5

    [0.5  1   0   ]   % 6

    [0    1   0   ]   % 7 green

    [0    0.5 0   ]   % 8

    [0    0.5 1   ]   % 9

    [0    1   0.5 ]   % 10

    [0    1   1   ]   % 11 cyan

    [0    0   1   ]   % 12 blue

    [0    0   0.5 ]   % 13

    [0.5  0   1   ]   % 14

    [1    0   0.5 ]   % 15

    [1    0   1   ]   % 16 magenta

    };

%%%%%%%%%%%%%% Load Trajectory %%%%%%%%%%%%%%%%%%%%%%%

load(['C:\Users\[your user name]\Documents\MobiusSpaceMissileSensor\' traj]) %ttraj, xtraj, ytraj, ztraj

ttraj_orig = ttraj;

xtraj_orig = xtraj;

ytraj_orig = ytraj;

ztraj_orig = ztraj;

% interpolate

ttraj = 0:tstep:ttraj_orig(end)-ttraj_orig(1); %seconds

xtraj = interp1(ttraj_orig-ttraj_orig(1),xtraj_orig,ttraj);

ytraj = interp1(ttraj_orig-ttraj_orig(1),ytraj_orig,ttraj);

ztraj = interp1(ttraj_orig-ttraj_orig(1),ztraj_orig,ttraj);

t_idx = 1:length(ttraj);

%%%%%%%%%%%%%%%%%%%%% Initialize figures %%%%%%%%%%%%%%%%%%%%

%% Figure 1

figure(1, 'position',[100,100,1000,1000]); hold on; grid on

xlabel('x [m]')

ylabel('y [m]')

zlabel('z [m]')

title('J2000 ECI Position of Satellites and Trajectory')

set(gca,'FontSize',20)

lgd_str_fig1 = {};

% plot trajectory

figure(1)

plot3(xtraj,ytraj,ztraj,'color','r','LineStyle','none','Marker','o','MarkerSize',1)

%lgd_str_fig1{end+1} = 'trajectory';

lgd_str_fig1{end+1} = 'trajectory timesteps';

% plot a sphere to represent earth

figure(3)

[xsphere, ysphere, zsphere] = sphere(20);

close

figure(1)

hSurface = surf(Re*xsphere,Re*ysphere,Re*zsphere);

set(hSurface,'FaceColor',[0.9 0.9 1],'EdgeColor',[1 1 1])

lgd_str_fig1{end+1} = 'earth';

view([1 1 1])

%% Figure 2

figure(2, 'position',[800,200,1000,1000]); grid on; hold on

xlabel('Threat Trajectory Time After Lift Off [sec]')

ylabel('Distance from Satellite to Target [km]')

title('Distance from Satellite to Target vs Time')

set(gca,'FontSize',20)

lgd_str_fig2 = {};

axis([0 ttraj(end) 0 5e3])

plot([ttraj(1) ttraj(end)],[sensorDist/1e3 sensorDist/1e3],'Color','k','LineWidth',5)

lgd_str_fig2{end+1} = 'Sensor Required Distance';

%%%%%%%%%% Calculate Orbit Coordinates %%%%%%%%%%%%%%%

'Orbit Coordinates'

r        = Re + alt;

grav     = G*M/r^2; %m/s^2, function of altitude

v        = sqrt(grav*r);

thetaDot = v/r; %rad/sec

for i_orbit = 1:num_orbs

  i_orbit

  if strcmp(conType,'Walker')

    RAAN  = (i_orbit-1)*2*pi/num_orbs; %rad

  elseif strcmp(conType,'Streets')

    RAAN  = (i_orbit-1)*pi/num_orbs; %rad

  endif

  phase = i_orbit*f_spacing*2*pi/tot_sats; %rad

  for i_sat = 1:num_sats

    %i_sat

    theta_vec = phase + 2*pi*(i_sat-1)/num_sats + thetaDot*ttraj; %rad/sec

    x1 = r*cos(theta_vec);

    y1 = 0*theta_vec;

    z1 = r*sin(theta_vec);

    incl_calc = pi/2-incl;   %rad, angle used for calculations

    %% Rotations

    Rx = [1 0 0                 ;

          0  cos(incl_calc) sin(incl_calc);

          0 -sin(incl_calc) cos(incl_calc)];

    Rz = [cos(RAAN) sin(RAAN) 0;

          -sin(RAAN)  cos(RAAN) 0;

          0 0 1                 ];

    for t_var = t_idx

      % Inclination rotation

      x2(t_var) = Rx(1,1)*x1(t_var) + Rx(1,2)*y1(t_var) + Rx(1,3)*z1(t_var);

      y2(t_var) = Rx(2,1)*x1(t_var) + Rx(2,2)*y1(t_var) + Rx(2,3)*z1(t_var);

      z2(t_var) = Rx(3,1)*x1(t_var) + Rx(3,2)*y1(t_var) + Rx(3,3)*z1(t_var);

      % RAAN rotation

      x3(t_var) = Rz(1,1)*x2(t_var) + Rz(1,2)*y2(t_var) + Rz(1,3)*z2(t_var);

      y3(t_var) = Rz(2,1)*x2(t_var) + Rz(2,2)*y2(t_var) + Rz(2,3)*z2(t_var);

      z3(t_var) = Rz(3,1)*x2(t_var) + Rz(3,2)*y2(t_var) + Rz(3,3)*z2(t_var);

    endfor %for t_var

    %% Plot orbits

    figure(1)

    plot3(x3,y3,z3,'marker',markers{mod(i_orbit,12)+1},'MarkerSize',3,'color','k','linestyle','-','linewidth',1)

    %% Format output orbit data structures

    x_orbit{i_orbit,i_sat} = x3;

    y_orbit{i_orbit,i_sat} = y3;

    z_orbit{i_orbit,i_sat} = z3;

  endfor %for i_sat

endfor %for i_orbit

%%%%%%%%%%%%%%%%% Calculations on Orbits %%%%%%%%%%%%%%%%%%%

'Calclations on Orbits'

for i_orbit = 1:num_orbs

  i_orbit

  for i_sat = 1:num_sats

%%%%%%%%%%%%%% Conjunction Calculations %%%%%%%%%%%%%%

    x3 = x_orbit{i_orbit,i_sat};

    y3 = y_orbit{i_orbit,i_sat};

    z3 = z_orbit{i_orbit,i_sat};

    for i_t = t_idx

      %% Calculate distance between satellites and trajectory

      dist(i_orbit,i_sat,i_t) = sqrt( (xtraj(i_t)-x3(i_t))^2 + (ytraj(i_t)-y3(i_t))^2 + (ztraj(i_t)-z3(i_t))^2 ); %m

    endfor

    %check if distance is through the earth, dist<5000km

    figure(2)

    for i_vec = 1:length(dist(i_orbit,i_sat,:))

      dist_vec(i_vec) = dist(i_orbit,i_sat,i_vec); %m

    endfor

    plot(ttraj,dist_vec/1e3,'color',colors{mod(i_sat,16)+1},'Marker',markers{mod(i_orbit,12)+1},'MarkerSize',3);

    %legend

    #lgd_str_sats{i_lgd} = ['orbit ' num2str(i_orbit) ', sat ' num2str(i_sat)];

    #i_lgd = i_lgd+1;

  endfor %for i_sat

endfor %for i_orbit

STK Multi-Hop Communication Example

More challenging with less help

Visit the Mobius STEM website to see an STK video and example: [https:\\www.mobius-llc.com\STEM}

Download a free version of STK as a student

Task : Create a 24 hour scenario that calculate the duration of access link between a target, a 6x6 walker constellation of SAR satellites (Capella Whitney), a 6x11 Walker constellation of comms satellites( Iridium) , and a ground station

STK Scenario: the video below shows 30 seconds of the simulation of the scenario in STK

Legend

  • Yellow cone: a SAR satellite collecting data at the target location (Uplink comms)

  • Red cone: a SAR satellite communicating with a comms satellite ( Crosslink comms) and a comms satellite delivering data to the ground station (Downlink comms)

  • Red square block: Target location

  • Green square block: Ground Station location

Results

In 24 hours, the communication was successful twice:

  • 1st time: March 14th 2023 at 4 pm for 457.38 seconds

  • 2nd time: March 14th 2023 at 5:42 pm for 140 seconds