Wolfram Computation Meets Knowledge

Visualize a Satellite Path with Wolfram SystemModeler and Mathematica

Explore the contents of this article with a free Wolfram SystemModeler trial.Today we rely heavily on satellites orbiting Earth for a variety of purposes. Mapping satellites are used to collect satellite images used in maps. Communication satellites are used for both telecommunication and internet access or for navigation services like GPS and GLONASS. Other usage areas are weather study, scientific observation, and reconnaissance.

The following model, created in Wolfram SystemModeler, is of a geocentric, inclined circular Low Earth Orbit (LEO) satellite. Geocentric means that it orbits around the Earth. An inclined circular orbit means that the orbit follows a circle, but is not aligned with the equator of the Earth. LEO is the name given to the altitude range below 2,000 kilometers (1,200 miles).

Suppose you are considering using this geocentric LEO satellite to collect image data. To achieve this, you would want to know where it is at the moment, how high it is, and how fast it’s going. If you want images of cities, you want to know over which cities it currently is. A SystemModeler model combined with data and computational resources in Mathematica can answer all of these questions.

Diagram of a geocentric, inclined circular Low Earth Orbit (LEO) satellite

Creating such a model is straightforward in SystemModeler. Using drag-and-drop, create three subsystems. Model the Earth using a mass with constant rotation, the satellite using a mass with propulsion forces, and the control logic using two proportional derivative (PD) controllers.

This blog post focuses on illustrating the orbit and flight of the satellite in the above model.

Simulation results from a modeled system can often be hard to grasp and hard to understand without further processing and reasoning. The following plot shows the satellite position in space:

Plot showing the satellite position in space

The three curves show the position in the directions x, y, and z. This plot alone will not tell you a lot about the satellite. It cannot, for example, answer questions like “Over which part of the globe is the satellite?“, “What is the velocity of the satellite?“, or “What large city is the satellite closest to?” It could answer the question “At what altitude is the satellite?“, but that would require some manual computation.

The goal of this post is to illustrate how Mathematica and SystemModeler can work together to answer the questions above in a visual and clear way.

First, connect Mathematica to SystemModeler by loading the WSMLink package:

Needs["WSMLink`"]

Then load and simulate the satellite model for 100,000 seconds:

Import[FileNameJoin[{NotebookDirectory[], "Satellite.mo"}], "ModelicaModel"];

sim = WSMSimulate["Satellite.ExampleModel", {0, 100000}];

Then pick out some data needed for the visualizations from the simulation result object:

position = sim[{"satellite_position_over_earth.r[1]", "satellite_position_over_earth.r[2]", "satellite_position_over_earth.r[3]"}, t]; velocity = sim[{"satellite_velocity.v[1]", "satellite_velocity.v[2]", "satellite_velocity.v[3]"}, t]; {startTime, stopTime} = sim["SimulationInterval"];

As we work in Mathematica, we have access to continuously updated data sources for a variety of things. We are going to use that data for two purposes: retrieving a list of large cities and their coordinates and getting the radius of the Earth. Lets start with the cities, where we use CountryData and CityData:

(* A list of all countries *) countries = CountryData[]; (* Get a list of all large cities *)  largeCities = Join @@ Map[Function[co, CityData[{Large, co}]], countries]; (* restrict the list to cities with over 3 million inhabitants *)  hugeCities = Select[largeCities, CityData[#, "Population"] > 3000000 &]; (* Get coordinates for the cities *) cityCoords = Map[{#, CityData[#, "Coordinates"]} &, hugeCities];

The Earth radius is available in AstronomicalData:

earthRadius = AstronomicalData["Earth", "Radius"]

6.3675*10^6

In the interest of conciseness, we use a ready-made 3D plot of the Earth. An example of how to create such a plot can be found on the Texture reference page.

3D plot of the Earth

Now we are ready to start creating the pieces of our visualization. Start by creating a 3D plot with a point at the satellite’s current position. Define a function that returns that plot:

satellite3D[time_] := ListPointPlot3D[{position} /. t → time, PlotStyle → Directive[Orange, PointSize[0.03]]];

Show the position together with the Earth. The ViewVector option defines from which direction the plot is viewed. Also define a variable for the time at which we test all our visualizations:

showTime = 21000;

options3D[   time_] := {ViewVector → {2 position /. t → time, {0, 0, 0}}, Boxed → False, PlotRange → 2 earthRadius, ViewAngle → 70°}

Show[earthplot, satellite3D[showTime], options3D[showTime]]

3D plot of the Earth with the location of the satellite

The path of the satellite would be clearer if we also plot that path in the same 3D plot. Let’s do that using a 3D parametric plot, again defining a function to calculate it. We also add some color to the path, using the ColorFunction option, and increase the number of PlotPoints to get a smooth look:

satelliteTrail3D[time_] := ParametricPlot3D[position, {t, startTime, time}, ColorFunction → (ColorData["Rainbow"][#3] &), PlotPoints → 500]

Adding this to the previous plot is simple; we just add it to the list of things to show:

Show[earthplot, satellite3D[showTime], satelliteTrail3D[showTime], options3D[showTime]]

3D Plot of the satellite with paths

That is a nice 3D view of the satellite, and we can see where the satellite is above the Earth at any time. Such a 3D plot makes for a great interactive view, allowing you to rotate and manipulate the display.

Another way to display the path of the satellite is on a 2D map of the world. This is how Wolfram|Alpha displays the current position of the International Space Station.

The function CountryData that we used before also has graphic representations of countries and regions. We use this as our source for a world map. Choose an equirectangular projection:

projection = "Equirectangular"; worldMap = Show[CountryData["World", {"Shape", projection}], ImageSize → 400]

World map

Place the satellite and its path on this map using the graphic primitives Disk and Line. Because the equirectangular projection maps directly to longitude and latitude, we can use LatitudeLongitude to pick out the graphics coordinates from a GeoPositionXYZ object:

satellite2D[time_] := satellite2D[time] = Graphics[{Orange, Disk[Reverse[LatitudeLongitude[GeoPositionXYZ[position /. t → time]]], 3.5]}]

We split the position list into multiple lines, starting a new line each time the satellite crosses the edge of the graphic.

satelliteTrail2D[time_] := satelliteTrail2D[time] = Graphics[{Orange, Line /@ Split[Table[Reverse[LatitudeLongitude[GeoPositionXYZ[position /. t → x]]], {x, startTime, time, 100}], Abs[#1[[1]] - #2[[1]]] < 100 &]}]

World map showing the satellite paths

Another interesting thing to know is which cities are closest to the satellite. Use the list of cities we retrieved earlier in the blog post. Create a table of all the cities and their distances to the satellite. The distance is calculated from the point directly under the satellite using GeoDistance.

Then sort the table by distance and take the first three cities:

nearestCities[time_] := nearestCities[time] = Module[{cityDistances}, cityDistances = Table[{city[[1]], GeoDistance[LatitudeLongitude[GeoPositionXYZ[position /. t → time]], city[[2]]]/1000}, {city, cityCoords}]; Take[Sort[cityDistances, #1[[2]] < #2[[2]] &], 3]    ]

Get the nearest cities and show them in a grid:

nearestCities[showTime] // Grid

{Riyadh, Riyadh, SaudiArabia} 390.016 {Jiddah, Mecca, SaudiArabia} 1177.92 {Baghdad, Baghdad, Iraq} 1231.97

Also create functions for the velocity and altitude of the satellite:

satelliteAltitude[time_] :=   satelliteAltitude[time] = GeoPosition[GeoPositionXYZ[position /. t → time]][[1, 3]]/1000

satelliteAltitude[showTime]

1124.7

satelliteSpeed[time_] := satelliteSpeed[time] = Norm[velocity /. t → time]

satelliteSpeed[showTime]

7287.79

Now we have all the visualizations needed to create a nice report:

showReport[time_, imSize_] := Grid[{{"Altitude", satelliteAltitude[time]}, {"Speed", satelliteSpeed[time]}, {"Closest cities", Grid[nearestCities[time]]}, {Show[earthplot, satellite3D[time], satelliteTrail3D[time], options3D[time], ImageSize → imSize], SpanFromLeft}, {Show[worldMap, satellite2D[time], satelliteTrail2D[time], ImageSize → imSize], SpanFromLeft}}, Frame → All, FrameStyle → LightGray]

showReport[showTime, 450]

Complete report showing nearest cities, 3D plot of Earth with satellite paths, and world map with satellite paths

A way to further visualize the simulation results and the satellite trajectory is to make an animation over time. This is done in the following video by rearranging into a horizontal layout and exporting as video:

With these visualization functions and the model, which are downloadable below, you can do further analysis, answering and studying a range of different questions. For example, what happens if different altitude controllers are used? How much energy is needed to launch the satellite or to keep it at a stable altitude? What about if a sudden avoidance maneuver is required to avoid space debris or another satellite? How many satellites like this would be required to give “coverage” to the whole Earth or to a certain percentage of the population or area?

Download the full satellite model, including documentation and examples, from the industry examples page and start exploring!

Download this post as a Computable Document Format (CDF) file.

Comments

Join the discussion

!Please enter your comment (at least 5 characters).

!Please enter your name.

!Please enter a valid email address.