diff --git a/config/plot_tle_catalog.json.empty b/config/plot_tle_catalog.json.empty new file mode 100644 index 0000000..3d56e16 --- /dev/null +++ b/config/plot_tle_catalog.json.empty @@ -0,0 +1,22 @@ +// Copyright (c) 2014-2016 Kartik Kumar, Dinamica Srl (me@kartikkumar.com) +// Copyright (c) 2016 Enne Hekma, Delft University of Technology (ennehekma@gmail.com) +// Distributed under the MIT License. +// See accompanying file LICENSE.md or copy at http://opensource.org/licenses/MIT +{ + // Path to SQLite database containing scan data. Maximum of + "databases" : [ "", + "", + "", + ""], + // Directory where scan figure is stored. + "output_directory" : ".", + + // Set parameters for figures. + "figure_dpi" : 300, + + // # Set number of lines per entry in TLE catalog (2 or 3). + "tleEntryNumberOfLines" : 3, + + // # Set font size for axes labels. + "fontSize" : 22 +} diff --git a/python/plot_tle_catalog.py b/python/plot_tle_catalog.py new file mode 100644 index 0000000..0decc43 --- /dev/null +++ b/python/plot_tle_catalog.py @@ -0,0 +1,265 @@ +''' +Copyright (c) 2014-2016, K. Kumar (me@kartikkumar.com) +Copyright (c) 2016, E.J. Hekma (ennehekma@gmail.com) +All rights reserved. +''' + +################################################################################################### +# Set up input deck +################################################################################################### + +# Set path to TLE catalog file. + + +################################################################################################### + +''' + DO NOT EDIT PARAMETERS BEYOND THIS POINT!!! +''' + +################################################################################################### +# Set up modules and packages +################################################################################################### + +from sgp4.earth_gravity import wgs72 +from sgp4.io import twoline2rv +from sgp4.propagation import getgravconst + +from matplotlib import rcParams +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from twoBodyMethods import convertMeanMotionToSemiMajorAxis +import sys + +# I/O +import commentjson +import json +from pprint import pprint +import sqlite3 + + +# Parse JSON configuration file +# Raise exception if wrong number of inputs are provided to script +if len(sys.argv) != 2: + raise Exception("Only provide a JSON config file as input!") + +json_data = open(sys.argv[1]) +config = commentjson.load(json_data) +json_data.close() +pprint(config) + +fontSize = config['fontSize'] +figureDPI = config['figure_dpi'] +# exit() + +################################################################################################### + +################################################################################################### +# Read and store TLE catalog +################################################################################################### + +sma = [] +ecc5 = [] +incl5 = [] + +raan3 = [] +ecc3 = [] +aop3 = [] +inclinations3 = [] + +# order = ['all','SSO','GEO','HEO'] +for x in xrange(0,len(config['databases'])): + # tleCatalogFilePathNew = eval(str('tleCatalogFilePath' + str(x))) + tleCatalogFilePathNew = config['databases'][x] + # Read in catalog and store lines in list. + fileHandle = open(tleCatalogFilePathNew) + catalogLines = fileHandle.readlines() + fileHandle.close() + + # Strip newline and return carriage characters. + for i in xrange(len(catalogLines)): + catalogLines[i] = catalogLines[i].strip('\r\n') + + # Parse TLE entries and store debris objects. + debrisObjects = [] + for tleEntry in xrange(0,len(catalogLines),config['tleEntryNumberOfLines']): + debrisObjects.append(twoline2rv(catalogLines[tleEntry+1], catalogLines[tleEntry+2], wgs72)) + + # Sort list of debris objects based on inclination. + inclinationSortedObjects = sorted(debrisObjects, key=lambda x: x.inclo, reverse=False) + + inclinations = [] + raan = [] + ecc = [] + aop = [] + for i in xrange(len(inclinationSortedObjects)): + inclinations.append(inclinationSortedObjects[i].inclo) + raan.append(inclinationSortedObjects[i].nodeo) + ecc.append(inclinationSortedObjects[i].ecco) + aop.append(inclinationSortedObjects[i].argpo) + + smatemp = [] + smatemp = [convertMeanMotionToSemiMajorAxis(debrisObject.no/60.0, \ + getgravconst('wgs72')[1])-6373 \ + for debrisObject in debrisObjects] + smatemp2 = [] + smatemp2 = pd.DataFrame(smatemp, columns=[str(x)]) + sma.append(smatemp2) + + ecctemp = [] + ecctemp = [debrisObject.ecco for debrisObject in debrisObjects] + ecc4 = [] + ecc4 = pd.DataFrame(ecctemp, columns=[str(x)]) + ecc5.append(ecc4) + + incltemp = [] + incltemp = [np.rad2deg(debrisObject.inclo) for debrisObject in debrisObjects] + incl4 = [] + incl4 = pd.DataFrame(incltemp, columns=[str(x)]) + incl5.append(incl4) + + raan2 = [] + raan2 = pd.DataFrame(raan, columns=[str(x)]) + raan3.append(raan2) + + ecc2 = [] + ecc2 = pd.DataFrame(ecc, columns=[str(x)]) + ecc3.append(ecc2) + + aop2 = [] + aop2 = pd.DataFrame(aop, columns=[str(x)]) + aop3.append(aop2) + + inclinations2 = [] + inclinations2 = pd.DataFrame(inclinations, columns=[str(x)]) + inclinations3.append(inclinations2) + +sma = pd.concat(sma, axis=1) +ecc5 = pd.concat(ecc5, axis=1) +incl5 = pd.concat(incl5, axis=1) + +raan3 = pd.concat(raan3, axis=1) +ecc3 = pd.concat(ecc3, axis=1) +aop3 = pd.concat(aop3, axis=1) +inclinations3 = pd.concat(inclinations3, axis=1) + +################################################################################################### + +################################################################################################### +# Generate plots +################################################################################################### + +# Set font size for plot labels. +rcParams.update({'font.size': fontSize}) + +# Plot distribution of eccentricity [-] against semi-major axis [km]. +figure = plt.figure() +axis = figure.add_subplot(111) +plt.xlabel("Semi-major axis altitude [km]") +plt.ylabel("Eccentricity [-]") +plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0)) +plt.plot(sma['0'],ecc3['0'], \ + marker='.', markersize=1, color='k', linestyle='none') +if len(config['databases']) > 1: + plt.plot(sma['1'],ecc3['1'], \ + marker='s', markersize=10, color='c', linestyle='none') + if len(config['databases']) > 2: + plt.plot(sma['2'],ecc3['2'], \ + marker='^', markersize=10, color='g', linestyle='none') + if len(config['databases']) > 3: + plt.plot(sma['3'],ecc3['3'], \ + marker='D', markersize=6, color='r', linestyle='none') +# axis.set_xlim(xmax=0.4e5) +figure.set_tight_layout(True) +plt.savefig(config['output_directory'] + + "/figure1_debrisPopulation_eccentricityVsSemiMajorAxis.pdf", + dpi = figureDPI) +plt.close() + +# Plot components of eccentricity vector [-]. +figure = plt.figure() +axis = figure.add_subplot(111) +plt.xlabel("$e \cos{\omega}$ [-]") +plt.ylabel("$e \sin{\omega}$ [-]") +plt.plot(ecc3['0']*np.cos(aop3['0']),ecc3['0']*np.sin(aop3['0']), \ + marker='.', markersize=1, color='k', linestyle='none') +if len(config['databases']) > 1: + plt.plot(ecc3['1']*np.cos(aop3['1']),ecc3['1']*np.sin(aop3['1']), \ + marker='s', markersize=10, color='c', linestyle='none') + if len(config['databases']) > 2: + plt.plot(ecc3['2']*np.cos(aop3['2']),ecc3['2']*np.sin(aop3['2']), \ + marker='^', markersize=10, color='g', linestyle='none') + if len(config['databases']) > 3: + plt.plot(ecc3['3']*np.cos(aop3['3']),ecc3['3']*np.sin(aop3['3']), \ + marker='D', markersize=6, color='r', linestyle='none') +plt.axis('equal') +# axis.set_xlim(xmin=-.82, xmax=.82) +# axis.set_ylim(ymin=-.82, ymax=.82) +# axis.set(xticks=[-.8,-.4,0,.4,.8]) +# axis.set(yticks=[-.8,-.4,0,.4,.8]) +figure.set_tight_layout(True) +plt.savefig(config['output_directory'] + + "/figure2_debrisPopulation_eccentricityVector.pdf", + dpi = figureDPI) +plt.close() + + +# Plot distribution of inclination [deg] against semi-major axis [km]. +figure = plt.figure() +axis = figure.add_subplot(111) +plt.xlabel("Semi-major axis altitude [km]") +plt.ylabel("Inclination [deg]") +plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0)) +plt.plot(sma['0'],incl5['0'], \ + marker='.', markersize=1, color='k', linestyle='none') +if len(config['databases']) > 1: + plt.plot(sma['1'],incl5['1'], \ + marker='s', markersize=10, color='c', linestyle='none') + if len(config['databases']) > 2: + plt.plot(sma['2'],incl5['2'], \ + marker='^', markersize=10, color='g', linestyle='none') + if len(config['databases']) > 3: + plt.plot(sma['3'],incl5['3'], \ + marker='D', markersize=6, color='r', linestyle='none') +# axis.set_xlim(xmax=0.4e5) +figure.set_tight_layout(True) +plt.savefig(config['output_directory'] + + "/figure3_debrisPopulation_inclinationVsSemiMajorAxis.pdf", + dpi = figureDPI) +plt.close() + +# Plot components of inclination vector [deg]. +figure = plt.figure() +axis = figure.add_subplot(111) +plt.xlabel("$i \cos{\Omega}$ [deg]") +plt.ylabel("$i \sin{\Omega}$ [deg]") +plt.plot(np.rad2deg(inclinations3['0'])*np.cos(raan3['0']), + np.rad2deg(inclinations3['0'])*np.sin(raan3['0']), + marker='.', markersize=1, color='k', linestyle='none') +if len(config['databases']) > 1: + plt.plot(np.rad2deg(inclinations3['1'])*np.cos(raan3['1']), \ + np.rad2deg(inclinations3['1'])*np.sin(raan3['1']), \ + marker='s', markersize=10, color='c', linestyle='none') + if len(config['databases']) > 2: + plt.plot(np.rad2deg(inclinations3['2'])*np.cos(raan3['2']), \ + np.rad2deg(inclinations3['2'])*np.sin(raan3['2']), \ + marker='^', markersize=10, color='g', linestyle='none') + if len(config['databases']) > 3: + plt.plot(np.rad2deg(inclinations3['3'])*np.cos(raan3['3']), \ + np.rad2deg(inclinations3['3'])*np.sin(raan3['3']), \ + marker='D', markersize=6, color='r', linestyle='none') +plt.axis('equal') +# axis.set_xlim(xmin=-110.0, xmax=110.0) +# axis.set_ylim(ymin=-110.0, ymax=110.0) +figure.set_tight_layout(True) +plt.savefig(config['output_directory'] + + "/figure4_debrisPopulation_inclinationVector.pdf", + dpi = figureDPI) +plt.close() + + + + +################################################################################################### +################################################################################################### \ No newline at end of file diff --git a/python/twoBodyMethods.py b/python/twoBodyMethods.py new file mode 100644 index 0000000..748532e --- /dev/null +++ b/python/twoBodyMethods.py @@ -0,0 +1,33 @@ +''' +Copyright (c) 2014, K. Kumar (me@kartikkumar.com) +All rights reserved. +''' + +################################################################################################### +# Set up modules and packages +################################################################################################### + +import numpy as np + +################################################################################################### + +################################################################################################### + +################################################################################################### +# Convert mean-motion to semi-major axis. +################################################################################################### + +def convertMeanMotionToSemiMajorAxis(meanMotion, gravitionalParameter): + return ( gravitionalParameter / ( meanMotion ** 2.0 ) ) ** ( 1.0 / 3.0 ) + +################################################################################################### + +################################################################################################### +# Convert mean-motion to semi-major axis. +################################################################################################### + +# Convert semi-major axis to mean-motion. +def convertSemiMajorAxisToMeanMotion(semiMajorAxis, gravitionalParameter): + return np.sqrt( gravitionalParameter / ( semiMajorAxis ** 3.0 ) ) + +################################################################################################### \ No newline at end of file diff --git a/src/atomScanner.cpp b/src/atomScanner.cpp index 5984f15..b10156b 100644 --- a/src/atomScanner.cpp +++ b/src/atomScanner.cpp @@ -23,6 +23,8 @@ #include #include +#include + #include #include diff --git a/src/j2Analysis.cpp b/src/j2Analysis.cpp index 8ed768c..6288763 100644 --- a/src/j2Analysis.cpp +++ b/src/j2Analysis.cpp @@ -17,6 +17,8 @@ #include #include +#include + #include #include diff --git a/src/lambertFetch.cpp b/src/lambertFetch.cpp index d7e153c..187027a 100644 --- a/src/lambertFetch.cpp +++ b/src/lambertFetch.cpp @@ -10,6 +10,8 @@ #include +#include + #include #include diff --git a/src/lambertScanner.cpp b/src/lambertScanner.cpp index 073e8e0..13bcb76 100644 --- a/src/lambertScanner.cpp +++ b/src/lambertScanner.cpp @@ -21,6 +21,8 @@ #include #include +#include + #include #include @@ -191,7 +193,7 @@ void executeLambertScanner( const rapidjson::Document& config ) { DateTime departureEpoch = input.departureEpochInitial; departureEpoch = departureEpoch.AddSeconds( input.departureEpochStepSize * m ); - + const Eci tleDepartureState = sgp4Departure.FindPosition( departureEpoch ); const Vector6 departureState = getStateVector( tleDepartureState ); @@ -209,7 +211,7 @@ void executeLambertScanner( const rapidjson::Document& config ) = astro::convertCartesianToKeplerianElements( departureState, earthGravitationalParameter ); const int departureObjectId = static_cast< int >( departureObject.NoradNumber( ) ); - + // Loop over time-of-flight grid. for ( int k = 0; k < input.timeOfFlightSteps; k++ ) { @@ -354,7 +356,7 @@ void executeLambertScanner( const rapidjson::Document& config ) // Reset SQL insert query. query.reset( ); } - } + } } ++showProgress; diff --git a/src/sgp4Scanner.cpp b/src/sgp4Scanner.cpp index dc53987..1738731 100644 --- a/src/sgp4Scanner.cpp +++ b/src/sgp4Scanner.cpp @@ -27,6 +27,8 @@ #include #include +#include + #include #include