[GSoC2019|PlaneSpotting|Shoumik|3] Decoding of secondary flight parameters and recording overlap checker.

In the previous blog:

I discussed about how I am decoding the position of the aircraft from 2 ADS-B frames and also from 1 ADS-B frame and a reference position and how I found both handsome and erractic readings.

Let’s look into how other airborne parameters, such as ‚aircraft identity‘ and ‚velocity‘ are being decoded and calculated.

There is a special type of frame which is known as ‚Enhanced Mode-S‘ frames. They also contain similar information at times like the non-Enhanced Mode-S frames. I will elaborate on both in this blog.

Frames with downlink format 20 & 21 are known as Enhanced Mode-S frames.

The following methods were implemented in the code for decoding the information from the ADS-B frames.

Aircraft Identification

Frames with Downlink Format 17-18 and Type Code 1-4 are known as aircraft identification messages.

The structure of the 56-bit DATA field of any such aircraft identification frame will be as follows:

TCECC1C2C3C4C5C6C7C8
5366666666
TC: Type Code
EC: Emitter Code
Second row signifies the bit-length

The EC value with regard to the TC value determines the type of aircraft (Such as : Super, Heavy, light etc). When the EC is set to ‚zero‘ signifies such information is not being transmitted by the aircraft.

An Airbus A380 is known as 'Super' and a Boeing 777-300ER is known as 'Heavy'. Small, private or training aircrafts such as Cessna 172s are known as 'Light'.

For determining the callsign a lookup table is used:

#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######

‚#‘ signifies a null value and ‚_‘ signifies a space between two words.

Let us take an Aircraft identifier frame recorded by one of our ground stations:

a0001910200490f1df2820700716

Thus the DATA part of this frame is:

200490f1df2820
BIN00100000000001001001000011110001110111110010100000100000
INT41934955503232
CHAR[TYPE CODE]AIC172__

The integer number corresponds to the corresponding element in the lookup table array.

Thus, the decoded callsign is „AIC172“ which belongs to Air India on the route 172

BDS 2,0 frames are also known as Aircraft Identification Messages. In these frames the first 8 bit of the DATA segment contains the BDS number and the rest is same as the previous one.

Airborne Velocity

ADS-B frames with Downlink Format 17-18 and TC 19 contain the velocity information of the aircraft.

Subtype 1 (Ground Speed)

The typical structure of the frame:

Position in frame Position in data block Length Abbreviation Data
33 – 37 1 – 5 5 TC Type code
38 – 40 6 – 8 3 ST Subtype
41 9 1 IC Intent flag
42 10 1 RESV_A Reserved-A
43 – 45 11 – 13 3 NAC Velocity uncertainty (NAC)
46 14 1 S_ew East-West velocity sign
47 – 56 15 – 24 change 10 V_ew East-West velocity
57 25 1 S_ns North-South velocity sign
58 – 67 26 – 35 10 V_ns North-South velocity
68 36 1 VrSrc Vertical rate source
69 37 1 S_vr Vertical rate sign
70 – 78 38 – 46 9 Vr Vertical rate
79 – 80 47 – 48 2 RESV_B Reserved-B
81 49 1 S_Dif Diff from baro alt, sign
82 – 88 50 – 56 7 Dif Diff from baro alt

These two bits determine the direction at which the aircraft was flying:

S_ns:
     1 -> flying North to South
     0 -> flying South to North
 S_ew:
     1 -> flying East to West
     0 -> flying West to East

The speed and the heading can thus be computed by:





If the calculated heading is negative, then we add 360 to it.


Vertical Speed

The S_Vr bit (69) determines the direction(ascent/descent).

0- Up
1- Down

Therefore the actual speed is calculated by:

V = (int(Vr) – 1) * 64

Therefore, V is the speed and S_Vr states whether the aircraft is climbing or descending.

The VrSrc field (Vertical Rate Source) determines the measured altitude type:

0 - Barometric altitude rate
1 - Geometric altitude rate

Subtype 3 (AirSpeed)

Here let us see how the structure of Subtype 3 differs from Subtype1:

Position in frame Position in data block Length Abbreviation Data
33 – 37 1 – 5 5 TC Type code
38 – 40 6 – 8 3 ST Subtype
41 9 1 IC Intent change flag
42 10 1 RESV_A Reserved-A
43 – 45 11 – 13 3 NAC Velocity uncertainty (NAC)
46 14 1 S_hdg Heading status
47 – 56 15 – 24 10 Hdg Heading (proportion)
57 25 1 AS-t Airspeed Type
58 – 67 26 – 35 10 AS Airspeed
68 36 1 VrSrc Vertical rate source
69 37 1 S_vr Vertical rate sign
70 – 78 38 – 46 9 Vr Vertical rate
79 – 80 47 – 48 2 RESV_B Reserved-B
81 49 1 S_Dif Difference from baro alt, sign
82 – 88 50 – 66 7 Dif Difference from baro alt

Heading

The S_hdg bit shows whether heading data is available or not:

0 - Heading data not available
1 - Heading data available

If heading data is available then it is calculated by:

Heading = Hdg(10 bit data) * 360/1024

Velocity (Airspeed)

Simply converting the AS bits gives the airspeed in knots.

Vertical Rate is calculated in the same manner as in Subtype 1

Enhanced Mode-S:

Aircraft Intention (BDS4,0)

The structure:

FieldPositionLength
Status11
MCP selected altitude
Precision = 16ft
2-1312
Status141
FMS selected altitude
Precision = 16ft
15-2612
Status271
Barometric setting (in mb and not hPa)
Value is actual minus 800
Precision = 0.1mb
28-3912
Reserved40-478
Status481
VNAV hold state491
ALT hold state501
APP hold state511
Reserved52-532
Status541
Altitude source:
00: Unknown
01:Aircraft altitude
10: MCP set altitude
11: FMS set altitude9
55-562
Note: MCP stands for Mode Control Panel where the pilot feeds in the data to the autopilot
FMS stands for Flight Management Computer where the data is fed to calculate various other flight parameters.

Track and Turn (BDS 5,0):

Many more information is given such as roll angle, true track angle etc:

FieldPositionLength
Status11
Sign, 111
Roll angle
range = [-90, 90] degrees
Precision = 45/256 degree
3-119
Status121
Sign131
True Track angle
Precision: 90/512 degree-
14-2310
Status241
Ground Speed
Precision: 2 knots
25-3410
Status241
Sign361
Track angle rate
Precision: 8/256 degree/second
37-458
Status461
True Airspeed
Precision: 2 knots
47-5610
Note: The sign bit determines whether to use 2's complement to determine the value present after it. It it is set to '1' then we have to use 2's complement.

Airspeed and Heading (BDS 6,0)

A lot more information is conveyed in this kind of frames:

FieldPositionLength
Status11
Sign11
Magnetic Heading
Precision = 90/512 degree
3-1210
Status131
Indicated Airspeed
Precision: 1 knots
14-2310
Status241
Mach (speed)
Precision: 2.048/512 Mach
25-3410
Status351
Sign361
Barometric altitude rate
Precision: 32 feet/ minute
37-459
Status461
Sign471
Inertial altitude rate
Precision: 32 ft/minute
48-569
Unfortunately the percentage of decodable Enhanced Mode-S message received was very low in both short recordings and long recordings which was worth hours.

File overlap checker:

The ground stations were programed to span each recorded file for 2 minutes(approx.). Now this means that for a long duration, there will be an array of files from several ground stations.

There has to be a logic, where the code finds all the files from several ground stations recorded at the same or has overlapping recording times and club them together for the next step.

I identified few such cases of overlaps:

The overlapping files are clubbed together into a list for further processing.

Next step is to find the same frames present in different station recordings, so we need to know the files that overlap and process it accordingly.

Feel free to get in touch with me:

[GSoC2019|DirectDemod|Vladyslav M] Final Blog Report

In this blog, I will summarize all the work, that has been done during my three months of work at AerospaceResearch.net.

The resulting product is a fully functional web application and a set of supporting logic, that allows uploading, georeferencing, combining and displaying NOAA weather images. Below, more specific features are listed.

Work done

  • implemented functionality to extract and store satellite metadata
  • designed and implemented an accurate image georeferencing algorithm
  • designed and implemented an algorithm to merge several satellite images, taking care of overlapping areas
  • created CLI interfaces for each part of package API
  • implemented generation of web map and virtual globe
  • built custom web application, which includes web server backend with Flask and frontend
  • deployed the working application to the web (demo version at http://206.81.21.147:5000/map)
  • wrote tests, documentation, installation instructions and usage tutorials

Links to work

Getting started

For installation please see the installation guide in the repository description. To start using specific features please refer to the tutorials page, it contains usage examples with corresponding sample files for each use case. If you want to know more about how the software is working, please read my blogs. They describe in sufficient detail how the software is working.

Further work

  • Set up automatic data uploading from recording stations.
  • Set up a DNS domain.

Support

If you have any questions or issues concerning the usage of the software you can ask for help / open issues on GitHub or contact me directly.

[GSoC2019|DirectDemod|Vladyslav M] Web application.

This blog post will describe in general terms the architecture of web application,  main steps of deployment process, processing server and several additional features.

General App Architecture

The application is an end-to-end data processing pipeline starting from raw satellite signal recordings and ending with displaying of combined images on the map or virtual globe. The image below describes this in a more structured way. 

Figure 1. App architecture.

As showed on the image, the app consists of 4 main parts: 

  • Recording stations
  • Processing server
  • Web server
  • Frontend

Data is being passed along this processing chain. Functionality of each part will be described in following parts of the article.

Processing Server

The task of processing server is to gather the raw satellite signal recordings and transform it to georeferenced images. This transformation happens using directdemod API:

  1. DirectDemod decoder extracts the image from recording.
  2. Image is preprocessed, two channels are extracted.
  3. Each part is georeferenced and saved.
  4. Processed data is being sent to web server.

Server uses ssh connection to securely transmit data, namely scp command is used. SCP stands for Secure Copy, it allows secure transmission of arbitrary files from local computer to remote machine.

Web Application

Web application is implemented using Flask. Web server has the following structure:

Figure 2. Web server.

It contains two main directories images/ – stores uploaded georeferenced image and tms/ – stores Tile Map Service representation of processed images. Images could be uploaded in two ways, either via ssh or throught /upload page (could be removed in production for safety reasons). After upload is done and the images are properly stored, then server will automatically merge this images using directdemod.merge and create TMS, which will be saved in tms/ directory. After processing is done, metadata will be updated and  images will be displayed on website. Website itself contains 3 pages:

  1. map.html
  2. globe.html
  3. upload.html (may be removed)

Map page displays images on Leaflet webmap, globe page uses WebGLEarth library to display images on Virtual Globe.

Should be noted that final implementation is somewhat different then showed on figure 2. Ability of viewing different channels was added, therefore tms stored in 2 directories corresponding to each channel.

Figure 3. Image of second channel.

Another major feature is called „Scroll in time“. Images are stored in chronological order and could be viewed sequentially using special slider. For each entry both channels could be viewed using a toggle button.

Deployment

Application has been deployed to the web http://206.81.21.147:5000/map. Appropriate domain name will be set up shortly.

Conclusion

This is the final part of the GSoC project. All main functionality and features have been implemented, web application is up and running.

[GSoC2019|PlaneSpotting|Shoumik|2] Extraction, decoding of data & findings

Previously:

In the previous blog, I penned down a brief introduction to ADS-B(Automatic Dependant Survelliance-Broadcast) and how this project has planned to have all the ground stations self-synchronised.

We now have the main data in our hand :- The ADS-B frames itself! Lets dive deep and see how I am extracting the data and what are the, both normal and non-normal, findings that I’ve made.

Findings

Starting with the first version of the my decoder and data extractor code, I tested it with a rather small, approximately 2 minutes long recording from one of the ground stations.

Before I discuss about what I found in the recording, here is an overview of what I had expected to be present in the recording and what are the capabilities of the current decoder file:

  • The structure of a typical ADS-B message with Downlink Format 17 when converted to its binary equivalent.:
Length Position Abbreviation Name
5 1 – 5 DF Downlink Format
3 6 – 8 CA Capability
24 9 – 32 ICAO ICAO aircraft address
56 33 – 88 DATA Data
[33 – 37] [TC] Type code
24 89 – 112 PI Parity/Interrogator ID
  • Length of the messages according to Downlink Format(DF):
    • 112 bit (long): DF 16, 17, 18, 19, 20, 21, 22, 24
    • 56 bit (short): DF 0, 4, 5, 11
  • Not all data is contained in a single message, here’s how to know what data is contained in a particular message based on the Type Code:
Type Code(TC) Data
1 – 4 Aircraft identification
5 – 8 Surface position
9 – 18 Airborne position (Baro Altitude)
19 Airborne velocities
20 – 22 Airborne position (GNSS Altitude)
23 – 27 Reserved
28 Aircraft status
29 Target state and status information
31 Aircraft operation status

Now, the percentage of each type of message that I found in the recording:

  • Aircraft Identifiers: 0.77 %
  • Messages with surface position: 0.00 %
  • Airborne Position (Barometric Altitude): 8.21 %
  • Airborne Velocity messages: 7.60 %
  • Airborne Position (GNSS altitude): 0.00 %
  • Operation Status messages: 0.82 %
  • Enhanced Mode-S ADS-B data: 34.46 %
  • Squawk/Ident data: 3.09 %
  • Short Air to Air ACAS messages: 14.47 %
  • Encoded Altitude only: 29.22 %
  • Long ACAS messages: 1.34 %

A total of 13172 messages were detected in this recording.

After a brief analysis what I observed is, the velocity of each aircraft is sent almost after every 4-5 position messages. The reason for the absense of any surface position messages can be that the ground station was not close enough to any airfield. And, by looking at the statistics, usage of barometric altitude seems to be a standard than GNSS altitude.

Decoding the Data

Our prime point of focus are the Airborne Position messages (DF 17-18, TC 9-18) as the title of the project suggests. Let’s discuss about how my code goes about decoding the location from each frame.

Pre-Step : Creating ‚identifier‘ functions which determines the message type. One such identifier function would be an ‚Airborne Position‘ message identifier.

Typical structure of an ‚Airborne Position‘ Message is given as:

Airborne position message bits explained
Position Position in Data block Bit Length Content
33 – 37 1 – 5 5 Type code
38 – 39 6 – 7 2 Surveillance status
40 8 1 NIC supplement-B
41 – 52 9 – 20 12 Altitude
53 21 1 Time
54 22 1 CPR odd/even frame flag
55 – 71 23 – 39 17 Latitude in CPR format
72 – 88 40 – 56 17 Longitude in CPR format

There are two methods of determining the location:

  • Globally Unambiguous method: This method requires two frames for the calculation. The frames should be a pair of an odd frame and an even frame, received in succession. Two even frames or two odd frames will not work in this method.

Bit-54 in the binary equivalent of the message is the flag bit that tells us whether the frame is even or odd.

0 – Even Frame

1 – Odd Frame

Bits 55-71 contain latitude in CPR format, and bits 72-88 contain longitude in CPR format. Thus the CPR is calculated by

LAT_CPR = decimal(bits 55-71) / 131072

LON_CPR = decimal(bits 72-88)/ 131072

The decimal equivalent of the CPR is encoded in 17 bits. So the maximum value that it can have is 2^17 or 131072, thus the CPR is mapped to a floating value below 1 by dividing it by 131072.

Calculation of Latitude:

Latitude Index(j):

We need two constants to compute the relative latitudes


Where NZ is the abbreviation of number of geographic latitude zones present between the equator and the poles. NZ = 15 is a standard for Mode-S position encoding.

Thus, the calculation of the relative latitude coordinate:



For locations in southern hemisphere, the latitude lies between 270 – 360 degrees, but we need to restrict the range to [-90, 90]



The final latitude is chosen based on the newest frame received. If the newest received frame has an even CPR flag, then the even latitude (LatEven) will be the final latitude and vice-versa.

We now have one coordinate in our hand !!

Before we place foot in computing the longitude, we should check if both the latitudes in the two frames are consistent, i.e. they lie in the same latitude zone

  • Check for Latitude zone consistency:

Let’s define a function which returns the number of latitude zones. ‚NL(latitude)‘ is how it is named.

NL():


 The returned value lies in the closed interval [1, 59]

Both the calculated latitude of the even and the odd frames are passed to this function, if the returned value is equal for both the cases, then both of them lies in the same zone, and the code has successfully removed the ambiguity in them.

If the returned values are unequal, then skip the calculation and wait for a new message.

Calculation of Longitude:

If the newer message among the pair of frames is even:




When the newer message to be received is odd:





For both cases, if Lon is greater than 180 degrees

Now we have the coordinates of the aircraft (lat, lon) calculated using the ‚Globally unambiguous method‘

  • Locally unambiguous method:

This method requires a reference location for calculation. This location must be within 180NM from the location to be calculated.

Calculation of dLon


Latitude index (j):


Latitude is thus given by:


Calculation of longitude:

Here NL function is the same as it has been defined for the 'Globally Unambiguous method'.

Longitude index (m):


Thus, longitude is given by:


  • Calculation of Altitude:

Bits 41-52 contain the altitude in a DF 17 TC 9-18 ADS-B frame

Bit 48 is known as the Q-Bit which determines the encoding multiple of the altitude

If Q-Bit = 1; 25 feet

else if Q-Bit = 0; 100 feet

The binary string spans from bit 41-52 excluding Q-bit(48)

Thus the Altitude is given by:


Where N is the decimal equivelent of bit 41-52 excluding the Q-bit

It was now time to execute this implementation !!

And here are the first set of results:

Here is a comparison of the results obtained from my code and data acquired from trusted 3rd party sources
The blue plots are the locations calculated by the current code. 
Whereas, the yellow plots are the ones obtained from other sources.
The 3rd party sources do not give us with raw calculated data, that's why the plots do not match even though the flight path overlaps.

Here is a plot where you can identify a cruising aircraft:

And also an aircraft which is beginning is decent towards Leipzig Airport:

Unusual findings:

The first stage of the code where I had implemented the decoding of position only by ‚Globally Unambiguous Method‘, I found a few stray points in the plots. A good term for this would be removal of ‚Vomit Comets‘. For example:

As you can see, it does not define a proper flight path.

Here, I now introduce a checking system, where the position of the newest frame is calculated using the ‚Locally unambiguous method‘ and it then compares both the positions.

If there is a difference, then the position calculated by the latter method is taken as ultimate.

Also, many singular frames were skipped from the intial position calculation logic as there wasn’t a compatible succesive frame. For them, the location calculated by the latter method was final.

Up Next:

I will dig further into how I am decoding and calculating other flight parameters such as, airspeed, ground speed, heading, aircraft identify codes (such as SQUAWK and Callsign).

Also, the initial phases after the completion of the decoding layer: Searching of overlapping files from different ground stations.

Feel free to get in touch with me:

[GSoC2019|Ksat|Themistoklis] 2nd Coding Period

Hi there!

This is the second in a series of blog posts, where I am documenting the progress of my project and my experience working with aerospaceresearch.net for Google Summer of Code 2019.

Last time I went into details about the often synonymous concepts, in the engineering world, of designing and optimizing a system. The provided context served as a prerequisite for explaining the purpose and the motivation behind both the KSat project, as well as my contribution in the project, which is the development of an automated visualization system to explore different design points, gain insights about the optimization process and the genetic algorithm’s performance, aid in debugging and more. In case you missed it and you would like to have more information, you can find the first blog post here: https://aerospaceresearch.net/?p=1542 . In fact, I highly recommend that you have read the first blog post before moving on with the second one.

From this point on, this blog post will focus on the progress during the second coding period. Before I get into technical topics, since I am documenting the whole experience, I would like to write a couple of words about the collaboration and the communication within aerospaceresearch.net. Like mentioned in the first blog post, everyone inside the organization has been extremely welcoming and helpful even from the pre-application stage. As it is to be expected during the second coding period, I was even more familiar with my project as well as with the exact fit of my work in the bigger picture. The communication with my main mentor, Manfred, continues to take place in an almost daily base for anything related to my work and I would like to thank him one more time for his valuable input and our great collaboration.

Regarding the technical aspects of my project, attention was first given in making the visualizations more presentable and readable. For this purpose, the structure of the xml that serves as an input to the visualization system, was extended to accommodate the following functionality:

-) Application of custom axes limits or limits that are calculated from the optimization data.

-) Definition of viewing angle for the 3d plot.

-) Definition of plot title.

-) Definition of axes labels that correspond to degrees-of-freedom (DoFs) that are assigned in x-, y- and z- axis.

-) Definition of legends that correspond to DoFs that are assigned in line style and line color as well as marker and marker color.

Some additional options like defining the font size for the title, the labels and the legends are available. The above functionality is more or less self-explanatory, but in order to make things more tangible the following visualizations from the first blog post are presented. The difference here is, that the visualizations have been properly annotated by activating the corresponding options through the xml file.

The first visualization from the first blog post, would take the following form:

In the figure above, I have defined “Propulsion Unit Optimization” as the title to applied on the plot as well as “Effective Exhaust Velocity / m/s”, “Thrust / N”, and “Total Mass Fraction / -“ as the labels of the x-, y- and z- axis correspondingly. In addition, I have specified that the axes limits should be calculated from the optimization data and that the viewing angle should be [azimuth, elevation] = [30, 60]. I have also enabled the legend for the remaining of the activated DoFs, which in this case are propulsion type and propellant. As it can be seen on the legend, propulsion type (arc jet or grid ion thruster) has been assigned to line style (solid or dashed), while propellant (He, Xe or NH3) has been assigned to marker (crosshair, circle or star). The chosen font sizes in this example are 14 for the title, 12 for the labels and 9 for the legends. At this point it should be noted that some of the above options, concerning the appearance of the visualization, were also manually applied at the corresponding visualization of the first blog post. The difference here is, that their automatic application has now been implemented as a series of activatable options in the xml file.

Now that we have an idea about the way in which the majority of a visualization’s elements are annotated, lets also see how potential additional continuous DoFs that are assigned to line and marker color are presented. For this purpose, we are going to use the last visualization from the first blog post, which would take the following form:

The difference in this visualization compared to the previous one is, that two additional hypothetical continuous DoFs have been assigned to line and marker colour. To indicate this, the corresponding DoFs “Hypothetical Continuous DoF 1” and “Hypothetical Continuous DoF 2” are included in the legend with a colored line and mark. Notice that are all other elements in the legend are black. The colored elements indicate that the value of these continuous dofs corresponds to the line and marker color in the visualization.

Now that we have covered the annotation of the visualizations, lets move to the remaining implemented functionality. A second type of visualization, a 2d visualization, was added in the visualization system. The motivation for the 2d plot is similar to the motivation for the 3d plot. Different system DoFs can be assigned to y-axis, line style and color as well as marker style and color. Notice that the x- and z-axis are no longer available to the user. At this point you may ask what the purpose of the 2d plot is then. In the 2d plot the x-axis corresponds to the generations of the system’s evolution, thus the 2d plot can be used to acquire a clearer picture regarding the chronological order in which different lineages evolve through different design points.

As usual I am going to give a practical example. I assign total mass fraction to y-axis, propulsion type to line style, propellant to marker, effective exhaust velocity to line color and thrust to marker color. I deactivate the visualization of failed mutations, I sort the lineages by total mass fraction and I keep only the 5 best lineages. I also activate the use of the implemented annotations. The plotting system returns the following figure:

As we can see from the figure above, all lineages converge to almost identical final design points. One lineage achieves this in just 3 generations, while another one need as much as 10 generations. In between of improved design points, we can also observe the number of generations where failed mutations occurred. For the lineage that converges on the 10th generation we can see that the successful mutations take place at the 3rd, 6th and 10th generation. This means that the mutations of the generations 2, 4-5 and 7-9 were failed mutations.

Another implemented feature is the ability to save the generated visualizations in a desired file format. This is required in order to review the progress of the optimization process after the execution of an optimization cycle. For this purpose, the following functionality was implemented as an activatable option through the xml file.

-) Saving of the visualization in one or more file formats. The available file formats are: *.ps, *.eps, *.jpg, *.png, *.emf, *.pdf

-) Generating filenames according to input case number, input/mission parameters and plot case number.

In order to avoid loss of any generated visualization it is important that the generated filenames always remain different regardless of the number of visualizations or the activated filename options. For this purpose, a hard-coded fail-safe check was also included. A sample generated filename could be “input_case_1_totalimpulse_112670_deltav_686_plot_case_1”. In this case, the input case number, the mission parameters (total impulse and velocity increment) along with their corresponding values and the plot case number, were all included in the title. For distinction purposes between different visualizations the plot case number is always included in the title, while in the absence of mission parameters the input case number is also automatically included.

Finally, another way to reap more potential insights from the visualizations is to include the dimension of time in them. This was achieved to a certain extend with the use of the 2d plot where the x-axis corresponds to the generations. There is another way to do so, where the graphics objects that form the visualizations are animated through a series of frames that follow the evolution of the system that is being optimized. The output is a series of frames that form a gif animation, rather than a single snapshot. The graphics objects can be animated in two different ways.

According to the first approach, the graphics objects are animated separately for each lineage. The sequence of the design points that correspond to the best lineage are animated first, then those that correspond to the second best and so on. This animated version of the last visualization of the first post, is the following:

In the animation above, only the 10 best lineages have been included. The frame rate of the animation is defined at 2 fps through the xml file.

The other way in which the graphics objects of a visualization can be animated, is for all lineages together. In this case, the sequence of the design points of all selected lineages is animated according to the progression of the generations.  Design points that appear on the same frame correspond to mutations of different lineages that took place on the same generation during the optimization process. This animated version is the following:

It should be noted that the examples for which the animations were applied, are simplified, thus a certain visual overlap occurs between different lineages. This causes the animations to appear stagnant at times, while what is really happening is that some lineages follow an identical sequence of design points as others. These design points have already appeared in the visualization due to already visualized lineages or due to lineages that have advanced through them at an earlier generation, giving the illusion that the animation remains stagnant at times, but this is not the case.

One final thing to mention is, that not only a 3d plot but also a 2d plot can be animated in a similar fashion.

Summing up this second blog post, we got a brief overview as well as a quick demonstration of the newly added features. We saw how these features build on top of the existing ones and open the possibilities for further exploration of the evolution data. There are some cases where minor bug fixes are required, and potential improvements and extensions have already been identified. Documentation is also going to be developed, but all that is work for the following weeks …

Congrats for making it through the second blog post and many thanks for reading!

[GSOC2019|VisMa|Mayank] Week 5 to 7 (Phase II)

It’s been almost two months I started to work on VisMa as my GSoC-19 Project (although I have been working with the community since December 2018) and it has been a quiet learning experience. This post will focus on technical details of what VisMa team has improved in the project roughly during the span of Phase II of GSoC.

Week 5: During this time, my goal was to complete the work regarding the Discrete Mathematics module. As I have mentioned in the previous Blog, we had added some basic Combinatorics modules (Factorial, Permutation & Combinations), this week our plan was to take this expedition forward and add some more to this Module. Also, we intended to implement the combinatorics module in CLI/GUI, which has not been done yet. Firstly I added the comments and animations in the above-mentioned modules. Adding comments and animations always seems like a cakewalk, but as per me doing this is actually hard and time-consuming (reason being, you have to keep a track of all the equations which occur during any operation). However, once done it always gives a sweet sense of completion, as it did this case. Our next target during Week 5, was to add more to the Discrete Mathematics module. We decided on adding a Statistics Module, a Probability Module and a (bitwise) boolean algebra module. The Statistics module as of now contains basic functions to calculate mean, median and mode like measures. Statistics is a topic of prime importance, thus having a statistics module is useful for the project. The other reason behind adding this module is that VisMa already has a graph plotter, this in later versions when combined with Statistics module can be used for the analysis of user-entered data. The other major part of Discrete Mathematics module was (bitwise) boolean algebra modules. These modules are designed keeping the teaching purpose of the project in mind. The comments and animations in these modules are in such a way that student will be able to observe how each bit of a number is being operated with the subsequent bit of the other number to get the final result. This part has been solely implemented keeping teaching perspective in mind. Lastly, we added a simple Probability Module to the project. As of now, we have Combinatorics, Probability, Statistics and (bitwise) boolean Algebra module added in the project.

Week 6: This week our task was to improve the integration and differentiation module of the project. The earlier logic and code for these two modules was beautifully implemented. My task was to add integration and derivation function to all Function Classes of the Project. Function Classes, in very simple terms, is a name given to a super-class of tokens whose subclasses are Constant, Variable etc. I had to add differentiation and integration for all these subclasses. Many of these were implemented but some were missing. Among missing ones, were Trigonometric, Exponential and Logarithmic Classes. I wrote the respective functions for all these Function Classes, also adding comments and animations side by side. Also, I refactored the existing code in differentiation & integration modules to follow an object-oriented style. Also added some test cases for same.

Week 7: The task of this week was focussed on improving the tokenizing module, adding some corner cases in Expression Simplification (involving exponents) and fixing some potential bugs. The tokenizer module treated Variable raised to some power, as a single token of Variable type (with the value of pow parameter set to power), but it didn’t recognise power operator in any case, my task was to fix it for recognising power operator. The potential bugs with Expression Simplification could only be resolved after it was done. The Expression Simplification follows a recursive logic, thus adding even a small improvement in that module sometimes become much confusing. But finally it was done in a nice manner, and VisMa is now able to deal with almost all the cases involving Expressions and Exponents. I also added test cases to reflect the new behaviour of the project.

Current Goal: As of now, I am working on implementing Matrix Module in CLI/GUI, the matrix operations have been implemented and now the next goal is to enable users to enter Matrices interactively in CLI/GUI

Below images illustrate the GUI/CLI representation of „factorial“, „combination“ and „permutation“ features in action.

Lastly, the project development is going at a good rate :). Some times it becomes buggy and confusing too, but lastly, it is a learning process and each bug does teach something new. I will soon be updating all the logic and working of the project in the wiki so that future developers can be helped.

Chears to space, Open Source, Maths & Code!

An open-source program’s Call for Help: Coordinated team-up observation for StarLink and NOAA satellites in July. Contribute to the satellite observing community 2019-07-19 until 2019-07-21!

Dear SatObs/SeeSat-L community,

We are students, taking part in Google Summer of Code 2019, and we seek your help!

We are working on an open-source project called OrbitDeterminator [w] that aims to determine the orbital parameters of satellites based on your positional data in various formats. Currently, we prepared the software so far and our next step is to test it to see how accurately it can determine satellite orbits So we would like to ask for your help in a few points to test our code under real-life conditions and your use-cases.

Our Aim, tackling StarLink and NOAAs together

We would like to tackle data of special interest for us and you, namely the StarLink satellite train and NOAA satellites. With the much talked-about StarLink observation data we can hopefully help the community to get an additional set of improved TLEs. And NOAAs are interesting for us because the TLEs are already available and we will be able to compare our results with existing TLEs.
So it would be great if you could provide us with observational positional data produced by positional observers on SatObs in IOD, UK and RDE formats. And in case there are none available yet, we would be more than happy to organize a small “global” observation campaign within the next month.

With this call-for-help, we suggest the 19th to 21st of July for us all to hunt the StarLink and NOAA satellites (more info in the information box below)

Overall Aim

With your help, we will provide you with a tool that is freely available, is easy to use and produces accurate results and benefits the satellite observing community. This information and data that you will provide will give more purpose to our code, where as we will try our best to give meaning to your data.

It would be our pleasure to end our Google Summer of Code project with this “field testing”, bringing our OrbitDeterminator package to good use within the SatObs community, and reporting our test speculations as a technical paper for the next International Aeronautical Congress 2019 call for papers in Washington DC.

Below you will find our initial line of thoughts and we are open to all discussions that will promote improvements in the same.

So who would be kind enough to answer our “call for help”? Feel free to answer on the mailing list or contact [m] us directly. We are looking forward to your positive reply!

Best regards,

Krishna, Rakshit & Vidhan (GSoC students)
Nilesh, Arya & Aakash (GSoC mentors)

m: gsoc2019@aerospaceresearch.net
c: https://aerospaceresearch.zulipchat.com/#narrow/stream/147024-OrbitDeterminator
w: https://github.com/aerospaceresearch/orbitdeterminator

Here are the areas where we need you and you people can help us:

Coordinating the Field Test:

  • Proposed Time-Span until third week of July 2019
  • Coordinated measurements during 48hours
  • We suggest to team up during 2019-07-19 18:00 UTC until 2019-07-21 18:00 UTC
  • Team-up together or single people. We organize via this SatObs Mailinglist!
  • As many single measurements as possible send via the SatObs mailing list for EVERYONE, not just us.
  • Deliverables by us: new code to OrbitDeterminator and also the results (under public domain licence)

Questions:

  1. Since a lot of data that is reported in the mailing list does not contain detailed location information about the observer, we would like the location of all observers ordered by site/station number that contains: [site/station number, observer code, observer name, latitude, longitude, elevation, active/inactive, preferred reporting format].
  2. Since all data is reported with an intrinsic measurement uncertainty, there will be some degree of error in the determined orbits as well. Can you please specify the what a typical accepted margin of error is for you?
  3. We currently produce the 7 keplerian elements and 3D orbit plots as output, in a specified format using a command line user interface. Do you people have any other preferences regarding the user interface, input options, file extensions and the output format?
  4. Can you provide the historical data for some satellites in all 3 (UK, IOD, RDE) formats? It will help us test our code and remove possible errors and exceptions.
  5. We are very much interested in the Starlink satellite train. Would you please provide us with current observations for the same?
  6. We would like to compare NOAA satellite orbits as well. But we could not find a lot of reports on the mailing list regarding the same. Could you provide some?
  7. Lastly, can you provide some reports with the results already known? This will help us compare our results with the expected results and improve our algorithms accordingly.

If we will get more questions during our discussion on the mailinglist, we will send it to the ML. And if we need to share some data, we can use this google drive folder.
https://drive.google.com/drive/folders/17K1yleBHyZw_c_vOYpbwx5IJRmoQfz2i?usp=sharing

PS: I am forwarding this email on behalf of the GSOC students. They subsribed a few weeks ago and just received the „Your subscription request has been forwarded to the list administrator at seesat-l-owner@satobs.org for review.“ notice and no further confirmation. Would you please be so kind and subscribe them to the list? It will make conversion much easier for them.

[GSoC2019|Ksat|Themistoklis] 1st Coding Period

Hello everyone!

With these series of blogs posts I will be documenting the progress of my project and my experience working with aerospaceresearch.net for Google Summer of Code 2019.

First, a few words about my mentors from aerospaceresearch.net and the communication between us. The main mentor of my project is Manfred Ehresmann, with whom we discuss almost daily for everything regarding my project. Since there is also a lot of development going on from Manfred’s side, effective communication for defining my project’s requirements, coming to common ground regarding technical matters that affect the work of both and solving any issues that arise has been key to the progress of my project and I am very happy that we have effortlessly achieved this. My second mentor is Andreas Hornig and since he is the main mentor of other projects, we communicate on a less regular basis for everything regarding my work with aerospaceresearch.net and my participation in GSoC. Andreas has also been incredibly helpful throughout this time, from guiding us through the application process to taking care of our onboarding in the organization and making sure that we have everything set up to be able to work undistracted on our projects. The platform that we use for our communications is Zulip.

Now it’s time for an introduction to the GSoC project and my contribution. Manfred is developing a software for designing space craft systems, currently limited to electric propulsion units for satellites. Since in the engineering world, the word “design” is most of the times synonyms to the word “optimize”, the later is actually the goal of the software. That is, to use data from existing systems and additional resources along with system modelling and scaling rules, in order to apply a genetic algorithm that can design an optimized space craft subsystem, which is an electric propulsion system currently, for different mission scenarios. This optimization is not straight forward, especially when there are a lot of sub-systems with non-linear interactions between them to be included in future development. Additionally, the algorithm itself has a lot of hyperparameters, whose tuning is a non-trivial task. For this reason, it is useful to have a way to evaluate the progress of the optimization process, in order to gain insights that can either lead in improvements of the algorithm or be directly utilized my human designers as well as aid in debugging of newly added features. This is where my contribution to the code comes. I am developing a visualization system which can be used to automatically generate useful visualizations from the evolution data of the genetic algorithm according to user defined parameters.

Let’s say a few words about the work that has been completed during the 1st coding period. During the previous weeks I have been laying out the architecture of the visualization system, defining the structure of the xml file that is used by the user to communicate with the system and implementing the first type of plot, the 3d plot. By this time, you may wonder what this plot is, in which way it is using the evolution data and how is it even useful in the first place! To get a better understanding, let’s talk about optimization for a second.

What does it mean to optimize an electric propulsion system? Well, first let’s define what it means to design a system. To design a system means to select the appropriate values for some variables of the system. These are here the effective exhaust velocity, thrust, propulsion type, propellant and more. For initial system design these are enough to fully define the rest of the systems parameters by using known system modelling equations and scaling laws equations as well as any system or mission requirements, for example total impulse, velocity increment, propulsion power. We can name the first set of variables, independent degrees-of-freedom (dofs) of the system, and the second set of variables that are calculated from the first, dependent dofs. The goal of the design process is to come up with a system that achieves a certain performance. If we define a specific performance criterion (ex. mass fraction of the electric propulsion unit), then by systematically searching for the best (optimal) values of the independent dofs we can come up with a system that best meets this performance criterion (ex. has the lowest mass fraction). This process is called optimization. In our case, a genetic algorithm is used for searching the optimal values of the independent dofs, which can be either continuous variables (ex. effective exhaust velocity, thrust) or discrete variables (ex. propulsion type, propellant). The electric propulsion system has also additional dependent variables/parameters (ex. mass fractions and efficiencies of sub-systems) which are calculated from the independent variables.

By visualizing the evolution data, we want to observe how the system moves along the design space during the optimization process. We want to see whether the mutations of the system are successful or not and we want to identify the design points where the system converges and possibly observe any other piece of information or pattern that can help improve the optimization and better understand the optimal designs. For this purpose, we want to plot different system variables (ex. effective exhaust velocity, thrust, propulsion type, propellant etc.) against the objective criterion (ex. system mass fraction). Using a 3d plot we can assign the objective criterion to the z-axis, the effective exhaust velocity to the x-axis and the thrust to the y-axis. In order to include additional dofs into the visualization, we can encode them into other aspects of the visualization, specifically line style, line color, marker and marker color. Continuous or discrete dofs can be assigned to line color and marker color by mapping the values of the dofs to the colors of a corresponding color map. Additionally, discrete dofs can be assigned to line style and marker by mapping the values of the dofs to a corresponding set of line styles and markers.

Using the xml file, we can define which system dofs (ex. system mass fraction, effective exhaust velocity etc.) are assigned to each plot dof (x-axis, y-axis, line style, line color etc.). The key capabilities of the visualization system can be summed up in the following points:

-) Flexible assignment of system dofs to plot dofs.

-) Selection of color maps, line styles, markers etc. to be used in the visualization.

-) Definition of default visualization style.

-) Activating or deactivating the visualization of failed mutations.

-) Defining custom visualization style for failed mutations.

-) Defining custom visualization size for seed points.

-) Defining multiple plot cases.

-) Specifying active input cases.

-) Sorting lineages by a specified dof.

-) Selective plotting of lineages.

In order to get a better understanding of the capabilities of the visualization system, let’s go over a few potential use cases. Let’s say that an optimization campaign has been completed and we would like to start exploring the evolution data that has been created by the genetic algorithm.

First, I assign the system dofs, for which I am interested, to corresponding plot dofs. In these examples I am going to assign effective exhaust velocity to x-axis, thrust to y-axis, total mass fraction to z-axis, propulsion type to line style and propellant to marker. I will also define that the possible propulsion systems (arcjet, grid ion thruster) are going to map to line styles (solid, dashed) and the possible propellants (He, Xe, NH3) are going to map to markers (crosshair, circle, star). I am also going to activate the visualization of failed mutations, but I am going to leave deactivated the custom visualization style for failed mutations, since with this first visualization I just want to observe how the design space has been explored by the algorithm. The visualization system returns the following figure:

In the figure above, we can see that only a certain portion of the design space has been explored. The seed points are represented by the large blue markers. Since line color and marker color has not been assigned to a system dof, they are black, which is the default visualization style that I have specified.

Now we are going to use the capabilities of the visualization system to gradually get a deeper understanding of our data. It would be interesting to see which portion of the mutations are successful. For this purpose, I am going to activate custom visualization style for the failed mutations, and I will specify red color for line color and marker color. The visualization system returns the following figure:

From the figure above it is clear that the majority of the mutation are failed mutation. For this reason, I would like to visualize only the successful mutation. To do this, I deactivate the visualization of failed mutations. The visualization system returns the following figure:

The trend of the optimal design points to converge to high effective exhaust velocity, low thrust, grid ion thruster propulsion (dashed line) and Xe propellant (circle marker) is easily recognizable. Note that the seeds points have now different colors, which is currently a bug rather than an implemented feature. Next, I want to visualize only the three best designs. For this purpose, I am going to define the total mass fraction as the sorting criterion of lineages and select “increasing” as the sorting direction, since the best designs are the ones with the lowest total mass fractions. The visualization system returns the following figure:

We can see that all 3 lineages converge to the same design point. It would be interesting to also explore the behavior of the 3 worst lineages. The visualization system returns the following figure:

By surprise, the three worst lineages also converge to the same optimal design point. This means that our algorithm is very robust to the selection of initial seed points, at least for this simpler example. To confirm that all lineages converge to the same design point, I am also going to visualize only the worst lineage. The visualization system returns the following figure:

Indeed, all lineages converge to the same design point because the worst one does so. At this point we could come up with more visualizations to explore, but for the purpose of briefly demonstrating the flexibility of the visualization system these examples are considered enough.

One last thing to demonstrate is what the visualization would look like when line color and marker color are also assigned to a system dof. That would be the case if the model of the electric propulsion system had additional independent dofs. For this example, we are going to assume that two hypothetical continuous dofs are assigned to line color and marker color respectively. I am going to define “jet” as the color map for both line color and marker color. I am also going to activate the visualization of failed mutations and include all lineages in the visualization. The visualization system returns the following figure:

As we can see from the figure above, the two additional dofs have been encoded in the visualization.

Summing up this first blog post, there are still a few things that need to be implemented in the visualization system regarding the presentation of the 3d plot (labels, legends, viewing angle, limits etc.) which will be added soon. An additional type of plot as well as animation capabilities are going to be added during the following coding periods. But for now, …

Many thanks for reading and cheers to the space and GSoC community!

[GSOC2019|VisMa|Mayank] Week 1 to 3

This blog post is intended to pen down my experience of working with aerospaceresearch.net and on VisMa – Visual Math as a part of Google Summer of Code 2019. A part of this blog post will be about how it feels to be working here & other parts will focus on technicalities of the project (about how VisMa has been improved over the past three weeks).

Firstly, working with aerospaceresearch.net has been a great experience. The mentors of the project Shantanu, Manfred & Siddharth are very helpful. We almost regularly have chat on Zulip (& sometimes on WhatsApp) wherein we discuss the plans for the project and about implementation. Owing to these healthy discussions the project development is going at a fast rate (which is indeed a good thing!). The code base of VisMa is very well written and the complex task of simplifying mathematical equations has been beautifully coded. About the workload, I work regularly for 5 – 6 hours, but sometimes when a crucial bug finds its way in, the work time has to be increased. Overall it’s a Fun & Learn experience.

Now let’s talk about what we (the VisMa Team) have achieved so far. During the community bonding period, I spend my time to get equipped with basics of PyQt (on which the GUI of VisMa is based) and spend the remaining time studying the tokenizing and simplify modules. Honestly, these modules were hard to understand. But the effort turned fruitful during the coding period. Also, I completed the task of simultaneous equation solver during this time. The logic was easy to frame, but the inclusion of comments and animations was tricky and involved the entire redesigning of the LaTex generator module (the part which created the final latex to be displayed in GUI).

Over to the first week of the coding period, the task which was decided was to re-design all the simplify module in an object-oriented fashion, which before coding period I believed would be most difficult to write (after all more than 2000 lines of complex Pythonic logic, to be honest, was scary). But, efforts to understand module during the Community Bonding Period turned out to be fruitful, and I was happy to achieve this task in the desired time. Like any other programming project, there are some minor bugs out there, which will be fixed in the upcoming time!

Coming to week second, we decided to work on Expression Simplification. Now, what expression is, briefly (& vaguely) any part of inputted equation enclosed in brackets can be termed as an Expression. So like when we solve an equation by hand, we solve the bracket part first, our task was to make visma do same! Earlier it didn’t support Expression simplification. For that, we decided to go for a complex recursive logic (solving the inner expression first then the outer expressions & so on). Recursive logic was complex & time consuming to implement, taking care of all comments & animations during recursion was a challenge. It took somewhat 9 days to come up with something fruitful. But yes it was worth it, VisMa could easily solve Expressions now. A happy moment indeed! 🙂

Next phase was to include a new module to VisMa. It was decided it has to be a Discrete Math module (‚cause why the heck not!). We decided to go for Combinatorics stuff first, thus implementing Combinations & Permutation stuff. These modules require factorial of a number, I got this idea why not to implement factorial as a separate module and then use in Permutations & Combinations stuff. This added another factorial feature to VisMa.

For now, I am working of adding triple degree equation solvers to project. We as of now have decided to do so using Cardano Algorithm. If this turns out well it will also be done for four-degree equations. After this will be done, VisMa will be capable of solving & showing detailed steps for up to 4-degree equations.

Below is a short GIF showing some of new features in action:

Hope you enjoyed reading Chapter One of my Developer Diary. From now on, I will write these blogs bi-weekly.

Chears to Code, Space & Maths!

[GSoC2019|PlaneSpotting|Shoumik] Introduction: Self – Synchronisation

Hi there!

My name is Shoumik Dey. I am currently a second year student at Manipal Institute of Technology, India. I preserve strong interest in aviation, aeromodelling and maybe this is the reason that I have chosen to work on this project, this summer, for Google Summer of Code 2019.

You can also find me here:

This is the first blog post for this project. This post will describe the work done so far, the current outcomes, hurdles faced and also they were/can be be solved.

Introduction

‚ADS-B‘ stands for Automatic Dependent Surveillance–Broadcast. The aircraft automatically brodcasts each frame on the 1090MHz frequency periodically which contains navigational data, such as, location, altitude, airspeed etc.

All ADS-B frames do not contain location information in them, therefore the ultimate goal of this project can be described in two parts:

  • Self-syncronisation of the ground stations in focus without an external agent.
  • Determining the location of the aircraft when location data is not received using multi-lateral positioning

Self-synchronisation

The current and widely used method of synchronising the clock of any ground based station is by using GNSS(Global Navigation Satellite System)-receiver interrogation. The local clocks gets aligned with the atomic clocks in the satellite.

But in this case, the synchronisation takes place by calculating the actual point of time when an aircraft broadcasts a message. This serves as the reference time at all ground stations, and since this time value is same all over, hence all the stations become self-synchronised.

Data required for self synchronisation:

  • Time of arrival of the message at the ground station
  • Time difference of arrival at that station.

Calculating the time difference of arrival(TDOA):

  1. The receiver provides the time at which the frame(with location) is received.
  2. The time of travel of the frame is calculated by using simple math and that is subtracted from the time of arrival.
sqrt((x2-x1)^2 + (y2 - y1)^2 + (z2 - z1)^2)

where x1, y1, z1: Latitude, longitude and altitude of ground station; x2, y2, z2: Latitude, longitude and altitude reported by aircraft.

Data that we have for self-synchronisation:

We have decided on using dump1090 for generating the ADS-B messages from the IQ stream as recorded by the SDR radio. Dump1090 provides the mlat(multilateration) data in avr format. The first 6 bytes of the message provides the sample position of the last bit of that message.

Example:

Raw: @0000000929E28e3ff6e6990c4684000011548194;
Sample Position: 0000000929E2
Message: 8e3ff6e6990c4684000011548194

Next, we shift the sample position record from the last bit of the message to the first bit.

Typical recording of a 112 bit ADS-B message
Typical recording of a 56 bit ADS-B message

As we can see, the length of the messages is exactly half of each other.

The sample position reported by dump1090 is exactly 240 samples ahead of the preamble (You can see the preamble at the beginning, by 4 spikes in the signal).

So, the sample position is taken back by 240 samples for both the cases to reach the beginning of the preamble.

Now, data in hand:

  • Start time of the stream
  • Sample position of the beginning of the message in the stream.

Therefore the time of arrival(TOA) can thus be calcualted.

Work done so far

  • A plugin to import the mlat output of dump1090 (ADS-B raw) has been created, which imports the message and the sample position into the main program
  • The raw data is reorganised in our suitable format and terms.
  • With each frame present;( for 112 bit frames as of now)
    • Determination of the type of frame and it’s content
    • Decoding those data present in that frame
  • Output dump file which contains the frame and the data contained in it. (Only extraction, no calcualation)

Next up….

The next blog post will contain the implementation of the calculation layer of the project. In this layer, the data in the frames will be processed and information such as position, velocity, altitude will be covered.

The post will also contain the unusual findings and irregularities that were found in the stream and how we have decided to deal with it.