mirror of
https://github.com/OpenSpace/OpenSpace.git
synced 2026-01-01 09:10:18 -06:00
269 lines
12 KiB
C++
269 lines
12 KiB
C++
/*****************************************************************************************
|
|
* *
|
|
* OpenSpace *
|
|
* *
|
|
* Copyright (c) 2014-2024 *
|
|
* *
|
|
* 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 <modules/gaia/tasks/readfilejob.h>
|
|
|
|
#include <openspace/util/distanceconversion.h>
|
|
#include <ghoul/format.h>
|
|
#include <ghoul/logging/logmanager.h>
|
|
#include <ghoul/misc/dictionary.h>
|
|
|
|
namespace {
|
|
constexpr std::string_view _loggerCat = "ReadFileJob";
|
|
} // namespace
|
|
|
|
namespace openspace::gaia {
|
|
|
|
ReadFileJob::ReadFileJob(std::filesystem::path filePath,
|
|
std::vector<std::string> allColumns, int firstRow, int lastRow,
|
|
size_t nDefaultCols, int nValuesPerStar,
|
|
std::shared_ptr<FitsFileReader> fitsReader)
|
|
: _inFilePath(std::move(filePath))
|
|
, _firstRow(firstRow)
|
|
, _lastRow(lastRow)
|
|
, _nDefaultCols(nDefaultCols)
|
|
, _nValuesPerStar(nValuesPerStar)
|
|
, _allColumns(std::move(allColumns))
|
|
, _fitsFileReader(std::move(fitsReader))
|
|
, _octants(8)
|
|
{}
|
|
|
|
void ReadFileJob::execute() {
|
|
// Read columns from FITS file. If rows aren't specified then full table will be read.
|
|
const std::shared_ptr<TableData<float>> table = _fitsFileReader->readTable<float>(
|
|
_inFilePath,
|
|
_allColumns,
|
|
_firstRow,
|
|
_lastRow
|
|
);
|
|
|
|
if (!table) {
|
|
throw ghoul::RuntimeError(
|
|
std::format("Failed to open Fits file '{}'", _inFilePath
|
|
));
|
|
}
|
|
|
|
const int nStars = table->readRows - _firstRow + 1;
|
|
|
|
const size_t nColumnsRead = _allColumns.size();
|
|
if (nColumnsRead != _nDefaultCols) {
|
|
LINFO(
|
|
"Additional columns will be read! Consider add column in code for "
|
|
"significant speedup"
|
|
);
|
|
}
|
|
|
|
// Copy columns to local variables.
|
|
std::unordered_map<std::string, std::vector<float>>& tableContent = table->contents;
|
|
|
|
// Default columns parameters.
|
|
//std::vector<float> l_longitude = std::move(tableContent[_allColumns[0]]);
|
|
//std::vector<float> b_latitude = std::move(tableContent[_allColumns[1]]);
|
|
std::vector<float> ra = std::move(tableContent[_allColumns[0]]);
|
|
std::vector<float> ra_err = std::move(tableContent[_allColumns[1]]);
|
|
std::vector<float> dec = std::move(tableContent[_allColumns[2]]);
|
|
std::vector<float> dec_err = std::move(tableContent[_allColumns[3]]);
|
|
std::vector<float> parallax = std::move(tableContent[_allColumns[4]]);
|
|
std::vector<float> parallax_err = std::move(tableContent[_allColumns[5]]);
|
|
std::vector<float> pmra = std::move(tableContent[_allColumns[6]]);
|
|
std::vector<float> pmra_err = std::move(tableContent[_allColumns[7]]);
|
|
std::vector<float> pmdec = std::move(tableContent[_allColumns[8]]);
|
|
std::vector<float> pmdec_err = std::move(tableContent[_allColumns[9]]);
|
|
std::vector<float> meanMagG = std::move(tableContent[_allColumns[10]]);
|
|
std::vector<float> meanMagBp = std::move(tableContent[_allColumns[11]]);
|
|
std::vector<float> meanMagRp = std::move(tableContent[_allColumns[12]]);
|
|
std::vector<float> bp_rp = std::move(tableContent[_allColumns[13]]);
|
|
std::vector<float> bp_g = std::move(tableContent[_allColumns[14]]);
|
|
std::vector<float> g_rp = std::move(tableContent[_allColumns[15]]);
|
|
std::vector<float> radial_vel = std::move(tableContent[_allColumns[16]]);
|
|
std::vector<float> radial_vel_err = std::move(tableContent[_allColumns[17]]);
|
|
|
|
|
|
// Construct data array. OBS: ORDERING IS IMPORTANT! This is where slicing happens.
|
|
for (int i = 0; i < nStars; i++) {
|
|
std::vector<float> values(_nValuesPerStar);
|
|
size_t idx = 0;
|
|
|
|
// Default order for rendering:
|
|
// Position [X, Y, Z]
|
|
// Mean G-band Magnitude
|
|
// -- Mean Bp-band Magnitude
|
|
// -- Mean Rp-band Magnitude
|
|
// Bp-Rp Color
|
|
// -- Bp-G Color
|
|
// -- G-Rp Color
|
|
// Velocity [X, Y, Z]
|
|
|
|
// Return early if star doesn't have a measured position.
|
|
if (std::isnan(ra[i]) || std::isnan(dec[i])) {
|
|
continue;
|
|
}
|
|
|
|
// Store positions. Set to a default distance if parallax doesn't exist.
|
|
float radiusInKiloParsec = 9.0;
|
|
if (!std::isnan(parallax[i])) {
|
|
// Parallax is in milliArcseconds -> distance in kiloParsecs
|
|
// https://gea.esac.esa.int/archive/documentation/GDR2/Gaia_archive/
|
|
// chap_datamodel/sec_dm_main_tables/ssec_dm_gaia_source.html
|
|
//LINFO("Parallax: " + std::to_string(parallax[i]));
|
|
radiusInKiloParsec = 1.f / parallax[i];
|
|
}
|
|
/*// Convert to Galactic Coordinates from Galactic Lon & Lat.
|
|
// https://gea.esac.esa.int/archive/documentation/GDR2/Data_processing/
|
|
// chap_cu3ast/sec_cu3ast_intro/ssec_cu3ast_intro_tansforms.html#SSS1
|
|
values[idx++] = radiusInKiloParsec * cos(glm::radians(b_latitude[i])) *
|
|
cos(glm::radians(l_longitude[i])); // Pos X
|
|
values[idx++] = radiusInKiloParsec * cos(glm::radians(b_latitude[i])) *
|
|
sin(glm::radians(l_longitude[i])); // Pos Y
|
|
values[idx++] = radiusInKiloParsec * sin(glm::radians(b_latitude[i])); // Pos Z
|
|
*/
|
|
|
|
// Convert ICRS Equatorial Ra and Dec to Galactic latitude and longitude.
|
|
const glm::mat3 aPrimG = glm::mat3(
|
|
// Col 0
|
|
glm::vec3(-0.0548755604162154, 0.4941094278755837, -0.8676661490190047),
|
|
// Col 1
|
|
glm::vec3(-0.8734370902348850, -0.4448296299600112, -0.1980763734312015),
|
|
// Col 2
|
|
glm::vec3(-0.4838350155487132, 0.7469822444972189, 0.4559837761750669)
|
|
);
|
|
const glm::vec3 rICRS = glm::vec3(
|
|
std::cos(glm::radians(ra[i])) * std::cos(glm::radians(dec[i])),
|
|
std::sin(glm::radians(ra[i])) * std::cos(glm::radians(dec[i])),
|
|
std::sin(glm::radians(dec[i]))
|
|
);
|
|
const glm::vec3 rGal = aPrimG * rICRS;
|
|
values[idx++] = radiusInKiloParsec * rGal.x; // Pos X
|
|
values[idx++] = radiusInKiloParsec * rGal.y; // Pos Y
|
|
values[idx++] = radiusInKiloParsec * rGal.z; // Pos Z
|
|
|
|
/*if (abs(rGal.x - values[0]) > 1e-5 || abs(rGal.y - values[1]) > 1e-5 ||
|
|
abs(rGal.z - values[2]) > 1e-5) {
|
|
LINFO("rGal: " + std::to_string(rGal) +
|
|
" - LB: [" + std::to_string(values[0]) + ", " + std::to_string(values[1]) +
|
|
", " + std::to_string(values[2]) + "]");
|
|
}*/
|
|
|
|
// Store magnitude render value. (Set default to high mag = low brightness)
|
|
values[idx++] = std::isnan(meanMagG[i]) ? 20.f : meanMagG[i]; // Mean G-band Mag
|
|
|
|
// Store color render value. (Default value is bluish stars)
|
|
values[idx++] = std::isnan(bp_rp[i]) ? 0.f : bp_rp[i]; // Bp-Rp Color
|
|
|
|
|
|
// Store velocity.
|
|
if (std::isnan(pmra[i])) {
|
|
pmra[i] = 0.f;
|
|
}
|
|
if (std::isnan(pmdec[i])) {
|
|
pmdec[i] = 0.f;
|
|
}
|
|
|
|
// Convert Proper Motion from ICRS [Ra,Dec] to Galactic Tanget Vector [l,b].
|
|
const glm::vec3 uICRS = glm::vec3(
|
|
-std::sin(glm::radians(ra[i])) * pmra[i] -
|
|
std::cos(glm::radians(ra[i])) * std::sin(glm::radians(dec[i])) * pmdec[i],
|
|
std::cos(glm::radians(ra[i])) * pmra[i] -
|
|
std::sin(glm::radians(ra[i])) * std::sin(glm::radians(dec[i])) * pmdec[i],
|
|
std::cos(glm::radians(dec[i])) * pmdec[i]
|
|
);
|
|
const glm::vec3 pmVecGal = aPrimG * uICRS;
|
|
|
|
// Convert to Tangential vector [m/s] from Proper Motion vector [mas/yr]
|
|
const float tanVelX = 1000.f * 4.74f * radiusInKiloParsec * pmVecGal.x;
|
|
const float tanVelY = 1000.f * 4.74f * radiusInKiloParsec * pmVecGal.y;
|
|
const float tanVelZ = 1000.f * 4.74f * radiusInKiloParsec * pmVecGal.z;
|
|
|
|
// Calculate True Space Velocity [m/s] if we have the radial velocity
|
|
if (!std::isnan(radial_vel[i])) {
|
|
// Calculate Radial Velocity in the direction of the star.
|
|
// radial_vel is given in [km/s] -> convert to [m/s].
|
|
const float radVelX = 1000.f * radial_vel[i] * rGal.x;
|
|
const float radVelY = 1000.f * radial_vel[i] * rGal.y;
|
|
const float radVelZ = 1000.f * radial_vel[i] * rGal.z;
|
|
|
|
// Use Pythagoras theorem for the final Space Velocity [m/s].
|
|
values[idx++] = static_cast<float>(
|
|
std::sqrt(std::pow(radVelX, 2) + std::pow(tanVelX, 2)) // Vel X [U]
|
|
);
|
|
values[idx++] = static_cast<float>(
|
|
std::sqrt(std::pow(radVelY, 2) + std::pow(tanVelY, 2)) // Vel Y [V]
|
|
);
|
|
values[idx++] = static_cast<float>(
|
|
std::sqrt(std::pow(radVelZ, 2) + std::pow(tanVelZ, 2)) // Vel Z [W]
|
|
);
|
|
}
|
|
// Otherwise use the vector [m/s] we got from proper motion.
|
|
else {
|
|
radial_vel[i] = 0.f;
|
|
values[idx++] = tanVelX; // Vel X [U]
|
|
values[idx++] = tanVelY; // Vel Y [V]
|
|
values[idx++] = tanVelZ; // Vel Z [W]
|
|
}
|
|
|
|
// Store additional parameters to filter by.
|
|
values[idx++] = std::isnan(meanMagBp[i]) ? 20.f : meanMagBp[i];
|
|
values[idx++] = std::isnan(meanMagRp[i]) ? 20.f : meanMagRp[i];
|
|
values[idx++] = std::isnan(bp_g[i]) ? 0.f : bp_g[i];
|
|
values[idx++] = std::isnan(g_rp[i]) ? 0.f : g_rp[i];
|
|
values[idx++] = ra[i];
|
|
values[idx++] = std::isnan(ra_err[i]) ? 0.f : ra_err[i];
|
|
values[idx++] = dec[i];
|
|
values[idx++] = std::isnan(dec_err[i]) ? 0.f : dec_err[i];
|
|
values[idx++] = std::isnan(parallax[i]) ? 0.f : parallax[i];
|
|
values[idx++] = std::isnan(parallax_err[i]) ? 0.f : parallax_err[i];
|
|
values[idx++] = pmra[i];
|
|
values[idx++] = std::isnan(pmra_err[i]) ? 0.f : pmra_err[i];
|
|
values[idx++] = pmdec[i];
|
|
values[idx++] = std::isnan(pmdec_err[i]) ? 0.f : pmdec_err[i];
|
|
values[idx++] = radial_vel[i];
|
|
values[idx++] = std::isnan(radial_vel_err[i]) ? 0.f : radial_vel_err[i];
|
|
|
|
// Read extra columns, if any. This will slow down the sorting tremendously!
|
|
for (size_t col = _nDefaultCols; col < nColumnsRead; ++col) {
|
|
std::vector<float> vecData = std::move(tableContent[_allColumns[col]]);
|
|
values[idx++] = std::isnan(vecData[col]) ? 0.f : vecData[col];
|
|
}
|
|
|
|
size_t index = 0;
|
|
if (values[0] < 0.0) {
|
|
index += 1;
|
|
}
|
|
if (values[1] < 0.0) {
|
|
index += 2;
|
|
}
|
|
if (values[2] < 0.0) {
|
|
index += 4;
|
|
}
|
|
|
|
_octants[index].insert(_octants[index].end(), values.begin(), values.end());
|
|
}
|
|
}
|
|
|
|
std::vector<std::vector<float>> ReadFileJob::product() {
|
|
return _octants;
|
|
}
|
|
|
|
} // namespace openspace::gaia
|