// /**************************************************************************************** // * * // * OpenSpace * // * * // * Copyright (c) 2014-2018 * // * * // * Permission is hereby granted, free of charge, to any person obtaining a copy of this * // * software and associated documentation files (the "Software"), to deal in the Software * // * without restriction, including without limitation the rights to use, copy, modify, * // * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to * // * permit persons to whom the Software is furnished to do so, subject to the following * // * conditions: * // * * // * The above copyright notice and this permission notice shall be included in all copies * // * or substantial portions of the Software. * // * * // * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, * // * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A * // * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT * // * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF * // * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE * // * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. * // ****************************************************************************************/ // #include // #include // #include // #include // #include // #include // #include // #include // #include // #include // #include // #include // #include // #include // #include // #include // #include // #include // #include // // Todo: // // Parse epoch correctly // // read distances using correct unit // // ... // namespace { // constexpr const char* ProgramName = "KeplerTrails"; // constexpr const char* KeyFile = "File"; // constexpr const char* KeyLineNum = "LineNumber"; // static const openspace::properties::Property::PropertyInfo PathInfo = { // "Path", // "Path", // "The file path to the CSV file to read" // }; // static const openspace::properties::Property::PropertyInfo SegmentsInfo = { // "Segments", // "Segments", // "The number of segments to use for each orbit ellipse" // }; // static const openspace::properties::Property::PropertyInfo EccentricityColumnInfo = { // "EccentricityColumn", // "EccentricityColumn", // "The header of the column where the eccentricity is stored" // }; // static const openspace::properties::Property::PropertyInfo SemiMajorAxisColumnInfo = { // "SemiMajorAxisColumn", // "SemiMajorAxisColumn", // "The header of the column where the semi-major axis is stored" // }; // static const openspace::properties::Property::PropertyInfo SemiMajorAxisUnitInfo = { // "SemiMajorAxisUnit", // "SemiMajorAxisUnit", // "The unit of the semi major axis. For example: If specified in km, set this to 1000." // }; // static const openspace::properties::Property::PropertyInfo InclinationColumnInfo = { // "InclinationColumn", // "InclinationColumn", // "The header of the column where the inclination is stored" // }; // static const openspace::properties::Property::PropertyInfo AscendingNodeColumnInfo = { // "AscendingNodeColumn", // "AscendingNodeColumn", // "The header of the column where the ascending node is stored" // }; // static const openspace::properties::Property::PropertyInfo ArgumentOfPeriapsisColumnInfo = { // "ArgumentOfPeriapsisColumn", // "ArgumentOfPeriapsisColumn", // "The header of the column where the argument of periapsis is stored" // }; // static const openspace::properties::Property::PropertyInfo MeanAnomalyAtEpochColumnInfo = { // "MeanAnomalyAtEpochColumn", // "MeanAnomalyAtEpochColumn", // "The header of the column where the mean anomaly at epoch is stored" // }; // static const openspace::properties::Property::PropertyInfo EpochColumnInfo = { // "EpochColumn", // "EpochColumn", // "The header of the column where the epoch is stored" // }; // } // namespace openspace { // documentation::Documentation RenderableSatellites::Documentation() { // using namespace documentation; // return { // "Renderable Kepler Orbits", // "space_renderable_kepler_orbits", // { // { // SegmentsInfo.identifier, // new DoubleVerifier, // Optional::No, // SegmentsInfo.description // }, // { // PathInfo.identifier, // new StringVerifier, // Optional::No, // PathInfo.description // }, // { // EccentricityColumnInfo.identifier, // new StringVerifier, // Optional::No, // EccentricityColumnInfo.description // }, // { // SemiMajorAxisColumnInfo.identifier, // new StringVerifier, // Optional::No, // SemiMajorAxisColumnInfo.description // }, // { // SemiMajorAxisUnitInfo.identifier, // new DoubleVerifier, // Optional::No, // SemiMajorAxisUnitInfo.description // }, // { // InclinationColumnInfo.identifier, // new StringVerifier, // Optional::No, // InclinationColumnInfo.description // }, // { // AscendingNodeColumnInfo.identifier, // new StringVerifier, // Optional::No, // AscendingNodeColumnInfo.description // }, // { // ArgumentOfPeriapsisColumnInfo.identifier, // new StringVerifier, // Optional::No, // ArgumentOfPeriapsisColumnInfo.description // }, // { // MeanAnomalyAtEpochColumnInfo.identifier, // new StringVerifier, // Optional::No, // MeanAnomalyAtEpochColumnInfo.description // }, // { // EpochColumnInfo.identifier, // new StringVerifier, // Optional::No, // EpochColumnInfo.description // } // } // }; // } // RenderableSatellites::RenderableSatellites(const ghoul::Dictionary& dictionary) // : Renderable(dictionary) // , _path(PathInfo) // , _nSegments(SegmentsInfo) // , _eccentricityColumnName(EccentricityColumnInfo) // , _semiMajorAxisColumnName(SemiMajorAxisColumnInfo) // , _semiMajorAxisUnit(SemiMajorAxisUnitInfo) // , _inclinationColumnName(InclinationColumnInfo) // , _ascendingNodeColumnName(AscendingNodeColumnInfo) // , _argumentOfPeriapsisColumnName(ArgumentOfPeriapsisColumnInfo) // , _meanAnomalyAtEpochColumnName(MeanAnomalyAtEpochColumnInfo) // , _epochColumnName(EpochColumnInfo) // { // documentation::testSpecificationAndThrow( // Documentation(), // dictionary, // "RenderableSatellites" // ); // _nSegments = // static_cast(dictionary.value(SegmentsInfo.identifier)); // _path = // dictionary.value(PathInfo.identifier); // _eccentricityColumnName = // dictionary.value(EccentricityColumnInfo.identifier); // _semiMajorAxisColumnName = // dictionary.value(SemiMajorAxisColumnInfo.identifier); // _inclinationColumnName = // dictionary.value(InclinationColumnInfo.identifier); // _ascendingNodeColumnName = // dictionary.value(AscendingNodeColumnInfo.identifier); // _argumentOfPeriapsisColumnName = // dictionary.value(ArgumentOfPeriapsisColumnInfo.identifier); // _meanAnomalyAtEpochColumnName = // dictionary.value(MeanAnomalyAtEpochColumnInfo.identifier); // _epochColumnName = // dictionary.value(EpochColumnInfo.identifier); // _semiMajorAxisUnit = // dictionary.value(SemiMajorAxisUnitInfo.identifier); // addPropertySubOwner(_appearance); // addProperty(_path); // addProperty(_nSegments); // addProperty(_semiMajorAxisUnit); // /* // * test // */ // const std::string& file = dictionary.value(KeyFile); // int lineNum = 1; // if (dictionary.hasKeyAndValue(KeyLineNum)) { // lineNum = static_cast(dictionary.value(KeyLineNum)); // readTLEFile(file, lineNum); // } // } // // The list of leap years only goes until 2056 as we need to touch this file then // // again anyway ;) // const std::vector LeapYears = { // 1956, 1960, 1964, 1968, 1972, 1976, 1980, 1984, 1988, 1992, 1996, // 2000, 2004, 2008, 2012, 2016, 2020, 2024, 2028, 2032, 2036, 2040, // 2044, 2048, 2052, 2056 // }; // // Count the number of full days since the beginning of 2000 to the beginning of // // the parameter 'year' // int countDays(int year) { // // Find the position of the current year in the vector, the difference // // between its position and the position of 2000 (for J2000) gives the // // number of leap years // constexpr const int Epoch = 2000; // constexpr const int DaysRegularYear = 365; // constexpr const int DaysLeapYear = 366; // if (year == Epoch) { // return 0; // } // // Get the position of the most recent leap year // const auto lb = std::lower_bound(LeapYears.begin(), LeapYears.end(), year); // // Get the position of the epoch // const auto y2000 = std::find(LeapYears.begin(), LeapYears.end(), Epoch); // // The distance between the two iterators gives us the number of leap years // const int nLeapYears = static_cast(std::abs(std::distance(y2000, lb))); // const int nYears = std::abs(year - Epoch); // const int nRegularYears = nYears - nLeapYears; // // Get the total number of days as the sum of leap years + non leap years // const int result = nRegularYears * DaysRegularYear + nLeapYears * DaysLeapYear; // return result; // } // // Returns the number of leap seconds that lie between the {year, dayOfYear} // // time point and { 2000, 1 } // int countLeapSeconds(int year, int dayOfYear) { // // Find the position of the current year in the vector; its position in // // the vector gives the number of leap seconds // struct LeapSecond { // int year; // int dayOfYear; // bool operator<(const LeapSecond& rhs) const { // return std::tie(year, dayOfYear) < std::tie(rhs.year, rhs.dayOfYear); // } // }; // const LeapSecond Epoch = { 2000, 1 }; // // List taken from: https://www.ietf.org/timezones/data/leap-seconds.list // static const std::vector LeapSeconds = { // { 1972, 1 }, // { 1972, 183 }, // { 1973, 1 }, // { 1974, 1 }, // { 1975, 1 }, // { 1976, 1 }, // { 1977, 1 }, // { 1978, 1 }, // { 1979, 1 }, // { 1980, 1 }, // { 1981, 182 }, // { 1982, 182 }, // { 1983, 182 }, // { 1985, 182 }, // { 1988, 1 }, // { 1990, 1 }, // { 1991, 1 }, // { 1992, 183 }, // { 1993, 182 }, // { 1994, 182 }, // { 1996, 1 }, // { 1997, 182 }, // { 1999, 1 }, // { 2006, 1 }, // { 2009, 1 }, // { 2012, 183 }, // { 2015, 182 }, // { 2017, 1 } // }; // // Get the position of the last leap second before the desired date // LeapSecond date { year, dayOfYear }; // const auto it = std::lower_bound(LeapSeconds.begin(), LeapSeconds.end(), date); // // Get the position of the Epoch // const auto y2000 = std::lower_bound( // LeapSeconds.begin(), // LeapSeconds.end(), // Epoch // ); // // The distance between the two iterators gives us the number of leap years // const int nLeapSeconds = static_cast(std::abs(std::distance(y2000, it))); // return nLeapSeconds; // } // double calculateSemiMajorAxis(double meanMotion) { // constexpr const double GravitationalConstant = 6.6740831e-11; // constexpr const double MassEarth = 5.9721986e24; // constexpr const double muEarth = GravitationalConstant * MassEarth; // // Use Kepler's 3rd law to calculate semimajor axis // // a^3 / P^2 = mu / (2pi)^2 // // <=> a = ((mu * P^2) / (2pi^2))^(1/3) // // with a = semimajor axis // // P = period in seconds // // mu = G*M_earth // double period = std::chrono::seconds(std::chrono::hours(24)).count() / meanMotion; // const double pisq = glm::pi() * glm::pi(); // double semiMajorAxis = pow((muEarth * period*period) / (4 * pisq), 1.0 / 3.0); // // We need the semi major axis in km instead of m // return semiMajorAxis / 1000.0; // } // double epochFromSubstring(const std::string& epochString) { // // The epochString is in the form: // // YYDDD.DDDDDDDD // // With YY being the last two years of the launch epoch, the first DDD the day // // of the year and the remaning a fractional part of the day // // The main overview of this function: // // 1. Reconstruct the full year from the YY part // // 2. Calculate the number of seconds since the beginning of the year // // 2.a Get the number of full days since the beginning of the year // // 2.b If the year is a leap year, modify the number of days // // 3. Convert the number of days to a number of seconds // // 4. Get the number of leap seconds since January 1st, 2000 and remove them // // 5. Adjust for the fact the epoch starts on 1st Januaray at 12:00:00, not // // midnight // // According to https://celestrak.com/columns/v04n03/ // // Apparently, US Space Command sees no need to change the two-line element // // set format yet since no artificial earth satellites existed prior to 1957. // // By their reasoning, two-digit years from 57-99 correspond to 1957-1999 and // // those from 00-56 correspond to 2000-2056. We'll see each other again in 2057! // // 1. Get the full year // std::string yearPrefix = [y = epochString.substr(0, 2)](){ // int year = std::atoi(y.c_str()); // return year >= 57 ? "19" : "20"; // }(); // const int year = std::atoi((yearPrefix + epochString.substr(0, 2)).c_str()); // const int daysSince2000 = countDays(year); // // 2. // // 2.a // double daysInYear = std::atof(epochString.substr(2).c_str()); // // 2.b // const bool isInLeapYear = std::find( // LeapYears.begin(), // LeapYears.end(), // year // ) != LeapYears.end(); // if (isInLeapYear && daysInYear >= 60) { // // We are in a leap year, so we have an effective day more if we are // // beyond the end of february (= 31+29 days) // --daysInYear; // } // // 3 // using namespace std::chrono; // const int SecondsPerDay = static_cast(seconds(hours(24)).count()); // //Need to subtract 1 from daysInYear since it is not a zero-based count // const double nSecondsSince2000 = (daysSince2000 + daysInYear - 1) * SecondsPerDay; // // 4 // // We need to remove additionbal leap seconds past 2000 and add them prior to // // 2000 to sync up the time zones // const double nLeapSecondsOffset = -countLeapSeconds( // year, // static_cast(std::floor(daysInYear)) // ); // // 5 // const double nSecondsEpochOffset = static_cast( // seconds(hours(12)).count() // ); // // Combine all of the values // const double epoch = nSecondsSince2000 + nLeapSecondsOffset - nSecondsEpochOffset; // return epoch; // } // void RenderableSatellites::readTLEFile(const std::string& filename, int lineNum){ // ghoul_assert(FileSys.fileExists(filename), "The filename must exist"); // std::ifstream file; // file.exceptions(std::ofstream::failbit | std::ofstream::badbit); // file.open(filename); // // All of the Kepler element information // struct { // double inclination = 0.0; // double semiMajorAxis = 0.0; // double ascendingNode = 0.0; // double eccentricity = 0.0; // double argumentOfPeriapsis = 0.0; // double meanAnomaly = 0.0; // double meanMotion = 0.0; // double epoch = 0.0; // } keplerElements; // std::string line; // // Loop through and throw out lines until getting to the linNum of interest // for (int i = 1; i < lineNum; ++i) { // std::getline(file, line); // } // std::getline(file, line); // Throw out the TLE title line (1st) // std::getline(file, line); // Get line 1 of TLE format // if (line[0] == '1') { // // First line // // Field Columns Content // // 1 01-01 Line number // // 2 03-07 Satellite number // // 3 08-08 Classification (U = Unclassified) // // 4 10-11 International Designator (Last two digits of launch year) // // 5 12-14 International Designator (Launch number of the year) // // 6 15-17 International Designator(piece of the launch) A // // 7 19-20 Epoch Year(last two digits of year) // // 8 21-32 Epoch(day of the year and fractional portion of the day) // // 9 34-43 First Time Derivative of the Mean Motion divided by two // // 10 45-52 Second Time Derivative of Mean Motion divided by six // // 11 54-61 BSTAR drag term(decimal point assumed)[10] - 11606 - 4 // // 12 63-63 The "Ephemeris type" // // 13 65-68 Element set number.Incremented when a new TLE is generated // // 14 69-69 Checksum (modulo 10) // keplerElements.epoch = epochFromSubstring(line.substr(18, 14)); // } else { // throw ghoul::RuntimeError(fmt::format( // "File {} @ line {} does not have '1' header", filename, lineNum + 1 // )); // } // std::getline(file, line); // Get line 2 of TLE format // if (line[0] == '2') { // // Second line // // Field Columns Content // // 1 01-01 Line number // // 2 03-07 Satellite number // // 3 09-16 Inclination (degrees) // // 4 18-25 Right ascension of the ascending node (degrees) // // 5 27-33 Eccentricity (decimal point assumed) // // 6 35-42 Argument of perigee (degrees) // // 7 44-51 Mean Anomaly (degrees) // // 8 53-63 Mean Motion (revolutions per day) // // 9 64-68 Revolution number at epoch (revolutions) // // 10 69-69 Checksum (modulo 10) // std::stringstream stream; // stream.exceptions(std::ios::failbit); // // Get inclination // stream.str(line.substr(8, 8)); // stream >> keplerElements.inclination; // stream.clear(); // // Get Right ascension of the ascending node // stream.str(line.substr(17, 8)); // stream >> keplerElements.ascendingNode; // stream.clear(); // // Get Eccentricity // stream.str("0." + line.substr(26, 7)); // stream >> keplerElements.eccentricity; // stream.clear(); // // Get argument of periapsis // stream.str(line.substr(34, 8)); // stream >> keplerElements.argumentOfPeriapsis; // stream.clear(); // // Get mean anomaly // stream.str(line.substr(43, 8)); // stream >> keplerElements.meanAnomaly; // stream.clear(); // // Get mean motion // stream.str(line.substr(52, 11)); // stream >> keplerElements.meanMotion; // } else { // throw ghoul::RuntimeError(fmt::format( // "File {} @ line {} does not have '2' header", filename, lineNum + 2 // )); // } // file.close(); // // Calculate the semi major axis based on the mean motion using kepler's laws // keplerElements.semiMajorAxis = calculateSemiMajorAxis(keplerElements.meanMotion); // // Converting the mean motion (revolutions per day) to period (seconds per revolution) // using namespace std::chrono; // double period = seconds(hours(24)).count() / keplerElements.meanMotion; // //TODO: fix obv // size_t i = 0; // _orbits[i++] = KeplerTranslation::KeplerOrbit{ // keplerElements.eccentricity, // keplerElements.semiMajorAxis, // keplerElements.inclination, // keplerElements.ascendingNode, // keplerElements.argumentOfPeriapsis, // keplerElements.meanAnomaly, // period, // keplerElements.epoch // }; // /* // KeplerTranslation setKeplerElements( // keplerElements.eccentricity, // keplerElements.semiMajorAxis, // keplerElements.inclination, // keplerElements.ascendingNode, // keplerElements.argumentOfPeriapsis, // keplerElements.meanAnomaly, // period, // keplerElements.epoch // ); // */ // } // /* // * !test // */ // RenderableSatellites::~RenderableSatellites() { // } // void RenderableSatellites::initialize() { // readFromCsvFile(); // updateBuffers(); // _path.onChange([this]() { // readFromCsvFile(); // updateBuffers(); // }); // _semiMajorAxisUnit.onChange([this]() { // readFromCsvFile(); // updateBuffers(); // }); // _nSegments.onChange([this]() { // updateBuffers(); // }); // } // void RenderableSatellites::deinitialize() { // } // void RenderableSatellites::initializeGL() { // glGenVertexArrays(1, &_vertexArray); // glGenBuffers(1, &_vertexBuffer); // glGenBuffers(1, &_indexBuffer); // _programObject = SpaceModule::ProgramObjectManager.request( // ProgramName, // []() -> std::unique_ptr { // return global::renderEngine.buildRenderProgram( // ProgramName, // absPath("${MODULE_SPACE}/shaders/RenderableKeplerOrbits_vs.glsl"), // absPath("${MODULE_SPACE}/shaders/RenderableKeplerOrbits_fs.glsl") // ); // } // ); // _uniformCache.opacity = _programObject->uniformLocation("opacity"); // _uniformCache.modelView = _programObject->uniformLocation("modelViewTransform"); // _uniformCache.projection = _programObject->uniformLocation("projectionTransform"); // _uniformCache.color = _programObject->uniformLocation("color"); // _uniformCache.useLineFade = _programObject->uniformLocation("useLineFade"); // _uniformCache.lineFade = _programObject->uniformLocation("lineFade"); // setRenderBin(Renderable::RenderBin::Overlay); // } // void RenderableSatellites::deinitializeGL() { // SpaceModule::ProgramObjectManager.release(ProgramName); // glDeleteBuffers(1, &_vertexBuffer); // glDeleteBuffers(1, &_indexBuffer); // glDeleteVertexArrays(1, &_vertexArray); // } // bool RenderableSatellites::isReady() const { // return true; // } // void RenderableSatellites::update(const UpdateData&) {} // void RenderableSatellites::render(const RenderData& data, RendererTasks&) { // _programObject->activate(); // _programObject->setUniform(_uniformCache.opacity, _opacity); // glm::dmat4 modelTransform = // glm::translate(glm::dmat4(1.0), data.modelTransform.translation) * // glm::dmat4(data.modelTransform.rotation) * // glm::scale(glm::dmat4(1.0), glm::dvec3(data.modelTransform.scale)); // _programObject->setUniform( // _uniformCache.modelView, // data.camera.combinedViewMatrix() * modelTransform // ); // _programObject->setUniform(_uniformCache.projection, data.camera.projectionMatrix()); // _programObject->setUniform(_uniformCache.color, _appearance.lineColor); // //_programObject->setUniform(_uniformCache.useLineFade, _appearance.useLineFade); // /*if (_appearance.useLineFade) { // _programObject->setUniform(_uniformCache.lineFade, _appearance.lineFade); // }*/ // glDepthMask(false); // //glBlendFunc(GL_SRC_ALPHA, GL_ONE); // glBindVertexArray(_vertexArray); // glDrawElements(GL_LINES, // static_cast(_indexBufferData.size()), // GL_UNSIGNED_INT, // 0); // glBindVertexArray(0); // _programObject->deactivate(); // } // void RenderableSatellites::updateBuffers() { // const size_t nVerticesPerOrbit = _nSegments + 1; // _vertexBufferData.resize(_orbits.size() * nVerticesPerOrbit); // _indexBufferData.resize(_orbits.size() * _nSegments * 2); // size_t orbitIndex = 0; // size_t elementIndex = 0; // for (const auto& orbit : _orbits) { // KeplerTranslation keplerTranslation(orbit); // const double period = orbit.period(); // for (size_t i = 0; i <= _nSegments; ++i) { // size_t index = orbitIndex * nVerticesPerOrbit + i; // double timeOffset = period * // static_cast(i) / static_cast(_nSegments); // glm::vec3 position = // keplerTranslation.position(Time(orbit.epoch + timeOffset)); // _vertexBufferData[index].x = position.x; // _vertexBufferData[index].y = position.y; // _vertexBufferData[index].z = position.z; // _vertexBufferData[index].time = timeOffset; // if (i > 0) { // _indexBufferData[elementIndex++] = static_cast(index) - 1; // _indexBufferData[elementIndex++] = static_cast(index); // } // } // ++orbitIndex; // } // glBindVertexArray(_vertexArray); // glBindBuffer(GL_ARRAY_BUFFER, _vertexBuffer); // glBufferData(GL_ARRAY_BUFFER, // _vertexBufferData.size() * sizeof(TrailVBOLayout), // _vertexBufferData.data(), // GL_STATIC_DRAW // ); // glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, _indexBuffer); // glBufferData(GL_ELEMENT_ARRAY_BUFFER, // _indexBufferData.size() * sizeof(int), // _indexBufferData.data(), // GL_STATIC_DRAW // ); // glBindVertexArray(0); // } // void RenderableSatellites::readFromCsvFile() { // std::vector columns = { // _eccentricityColumnName, // _semiMajorAxisColumnName, // _inclinationColumnName, // _ascendingNodeColumnName, // _argumentOfPeriapsisColumnName, // _meanAnomalyAtEpochColumnName, // _epochColumnName, // }; // std::vector> data = // ghoul::loadCSVFile(_path, columns, false); // _orbits.resize(data.size()); // size_t i = 0; // for (const std::vector& line : data) { // _orbits[i++] = KeplerTranslation::KeplerOrbit{ // std::stof(line[0]), // _semiMajorAxisUnit * std::stof(line[1]) / 1000.0, // std::stof(line[2]), // std::stof(line[3]), // std::stof(line[4]), // std::stof(line[5]), // std::stof(line[6]) // }; // } // } // }