Celestia

1
// stardbbuilder.cpp
2
//
3
// Copyright (C) 2001-2024, the Celestia Development Team
4
// Original version by Chris Laurel <claurel@gmail.com>
5
//
6
// This program is free software; you can redistribute it and/or
7
// modify it under the terms of the GNU General Public License
8
// as published by the Free Software Foundation; either version 2
9
// of the License, or (at your option) any later version.
10

11
#include "stardbbuilder.h"
12

13
#include <algorithm>
14
#include <array>
15
#include <cassert>
16
#include <cmath>
17
#include <cstdint>
18
#include <istream>
19
#include <iterator>
20
#include <string_view>
21
#include <type_traits>
22
#include <utility>
23

24
#include <boost/smart_ptr/intrusive_ptr.hpp>
25

26
#include <Eigen/Geometry>
27

28
#include <fmt/format.h>
29

30
#include <celastro/astro.h>
31
#include <celmath/geomutil.h>
32
#include <celmath/mathlib.h>
33
#include <celutil/binaryread.h>
34
#include <celutil/fsutils.h>
35
#include <celutil/gettext.h>
36
#include <celutil/logger.h>
37
#include <celutil/timer.h>
38
#include <celutil/tokenizer.h>
39
#include "hash.h"
40
#include "meshmanager.h"
41
#include "octreebuilder.h"
42
#include "parser.h"
43
#include "stardb.h"
44
#include "stellarclass.h"
45
#include "value.h"
46

47
using namespace std::string_view_literals;
48

49
namespace astro = celestia::astro;
50
namespace engine = celestia::engine;
51
namespace ephem = celestia::ephem;
52
namespace math = celestia::math;
53
namespace util = celestia::util;
54

55
using util::GetLogger;
56

57
struct StarDatabaseBuilder::StcHeader
58
{
59
    explicit StcHeader(const fs::path&);
60
    explicit StcHeader(fs::path&&) = delete;
61

62
    const fs::path* path;
63
    int lineNumber{ 0 };
64
    DataDisposition disposition{ DataDisposition::Add };
65
    bool isStar{ true };
66
    AstroCatalog::IndexNumber catalogNumber{ AstroCatalog::InvalidIndex };
67
    std::vector<std::string> names;
68
};
69

70
StarDatabaseBuilder::StcHeader::StcHeader(const fs::path& _path) :
71
    path(&_path)
72
{
73
}
74

75
template<>
76
struct fmt::formatter<StarDatabaseBuilder::StcHeader> : formatter<std::string_view>
77
{
78
    format_context::iterator format(const StarDatabaseBuilder::StcHeader& header, format_context& ctx)
79
    {
80
        fmt::basic_memory_buffer<char> data;
81
        fmt::format_to(std::back_inserter(data), "line {}", header.lineNumber);
82
        if (header.catalogNumber <= Star::MaxTychoCatalogNumber)
83
            fmt::format_to(std::back_inserter(data), " - HIP {}", header.catalogNumber);
84
        if (!header.names.empty())
85
            fmt::format_to(std::back_inserter(data), " - {}", header.names.front());
86
        return formatter<std::string_view>::format(std::string_view(data.data(), data.size()), ctx);
87
    }
88
};
89

90
namespace
91
{
92

93
// In testing, changing SPLIT_THRESHOLD from 100 to 50 nearly
94
// doubled the number of nodes in the tree, but provided only between a
95
// 0 to 5 percent frame rate improvement.
96
constexpr engine::OctreeObjectIndex StarOctreeSplitThreshold = 75;
97

98
// The octree node into which a star is placed is dependent on two properties:
99
// its obsPosition and its luminosity--the fainter the star, the deeper the node
100
// in which it will reside.  Each node stores an absolute magnitude; no child
101
// of the node is allowed contain a star brighter than this value, making it
102
// possible to determine quickly whether or not to cull subtrees.
103

104
struct StarOctreeTraits
105
{
106
    using ObjectType = Star;
107
    using PrecisionType = float;
108

109
    static Eigen::Vector3f getPosition(const ObjectType&);
110
    static float getRadius(const ObjectType&);
111
    static float getMagnitude(const ObjectType&);
112
    static float applyDecay(float);
113
};
114

115
inline Eigen::Vector3f
116
StarOctreeTraits::getPosition(const ObjectType& obj)
117
{
118
    return obj.getPosition();
119
}
120

121
inline float
122
StarOctreeTraits::getRadius(const ObjectType& obj)
123
{
124
    return obj.getOrbitalRadius();
125
}
126

127
inline float
128
StarOctreeTraits::getMagnitude(const ObjectType& obj)
129
{
130
    return obj.getAbsoluteMagnitude();
131
}
132

133
inline float
134
StarOctreeTraits::applyDecay(float factor)
135
{
136
    // Decrease in luminosity by factor of 4
137
    // -2.5 * log10(1.0 / 4.0) = 1.50515 (nearest float)
138
    return factor + 1.50515f;
139
}
140

141
constexpr float STAR_OCTREE_MAGNITUDE = 6.0f;
142

143
// We can't compute the intrinsic brightness of the star from
144
// the apparent magnitude if the star is within a few AU of the
145
// origin.
146
constexpr float VALID_APPMAG_DISTANCE_THRESHOLD = 1e-5f;
147

148
constexpr std::string_view STARSDAT_MAGIC = "CELSTARS"sv;
149
constexpr std::uint16_t StarDBVersion     = 0x0100;
150

151
#pragma pack(push, 1)
152

153
// stars.dat header structure
154
struct StarsDatHeader
155
{
156
    StarsDatHeader() = delete;
157
    char magic[8]; //NOSONAR
158
    std::uint16_t version;
159
    std::uint32_t counter;
160
};
161

162
// stars.dat record structure
163
struct StarsDatRecord
164
{
165
    StarsDatRecord() = delete;
166
    AstroCatalog::IndexNumber catNo;
167
    float x;
168
    float y;
169
    float z;
170
    std::int16_t absMag;
171
    std::uint16_t spectralType;
172
};
173

174
#pragma pack(pop)
175

176
static_assert(std::is_standard_layout_v<StarsDatHeader>);
177
static_assert(std::is_standard_layout_v<StarsDatRecord>);
178

179
bool
180
parseStarsDatHeader(std::istream& in, std::uint32_t& nStarsInFile)
181
{
182
    std::array<char, sizeof(StarsDatHeader)> header;
183
    if (!in.read(header.data(), header.size()).good()) /* Flawfinder: ignore */
184
        return false;
185

186
    // Verify the magic string
187
    if (auto magic = std::string_view(header.data() + offsetof(StarsDatHeader, magic), STARSDAT_MAGIC.size());
188
        magic != STARSDAT_MAGIC)
189
    {
190
        return false;
191
    }
192

193
    // Verify the version
194
    if (auto version = util::fromMemoryLE<std::uint16_t>(header.data() + offsetof(StarsDatHeader, version));
195
        version != StarDBVersion)
196
    {
197
        return false;
198
    }
199

200
    // Read the star count
201
    nStarsInFile = util::fromMemoryLE<std::uint32_t>(header.data() + offsetof(StarsDatHeader, counter));
202
    return true;
203
}
204

205
inline void
206
stcError(const StarDatabaseBuilder::StcHeader& header, std::string_view msg)
207
{
208
    GetLogger()->error(_("Error in .stc file ({}): {}\n"), header, msg);
209
}
210

211
inline void
212
stcWarn(const StarDatabaseBuilder::StcHeader& header, std::string_view msg)
213
{
214
    GetLogger()->warn(_("Warning in .stc file ({}): {}\n"), header, msg);
215
}
216

217
bool
218
parseStcHeader(Tokenizer& tokenizer, StarDatabaseBuilder::StcHeader& header)
219
{
220
    header.lineNumber = tokenizer.getLineNumber();
221

222
    header.isStar = true;
223

224
    // Parse the disposition--either Add, Replace, or Modify. The disposition
225
    // may be omitted. The default value is Add.
226
    header.disposition = DataDisposition::Add;
227
    if (auto tokenValue = tokenizer.getNameValue(); tokenValue.has_value())
228
    {
229
        if (*tokenValue == "Modify")
230
        {
231
            header.disposition = DataDisposition::Modify;
232
            tokenizer.nextToken();
233
        }
234
        else if (*tokenValue == "Replace")
235
        {
236
            header.disposition = DataDisposition::Replace;
237
            tokenizer.nextToken();
238
        }
239
        else if (*tokenValue == "Add")
240
        {
241
            header.disposition = DataDisposition::Add;
242
            tokenizer.nextToken();
243
        }
244
    }
245

246
    // Parse the object type--either Star or Barycenter. The object type
247
    // may be omitted. The default is Star.
248
    if (auto tokenValue = tokenizer.getNameValue(); tokenValue.has_value())
249
    {
250
        if (*tokenValue == "Star")
251
        {
252
            header.isStar = true;
253
        }
254
        else if (*tokenValue == "Barycenter")
255
        {
256
            header.isStar = false;
257
        }
258
        else
259
        {
260
            stcError(header, _("unrecognized object type"));
261
            return false;
262
        }
263
        tokenizer.nextToken();
264
    }
265

266
    // Parse the catalog number; it may be omitted if a name is supplied.
267
    header.catalogNumber = AstroCatalog::InvalidIndex;
268
    if (auto tokenValue = tokenizer.getNumberValue(); tokenValue.has_value())
269
    {
270
        header.catalogNumber = static_cast<AstroCatalog::IndexNumber>(*tokenValue);
271
        tokenizer.nextToken();
272
    }
273

274
    header.names.clear();
275
    if (auto tokenValue = tokenizer.getStringValue(); tokenValue.has_value())
276
    {
277
        for (std::string_view remaining = *tokenValue; !remaining.empty();)
278
        {
279
            auto pos = remaining.find(':');
280
            if (std::string_view name = remaining.substr(0, pos);
281
                !name.empty() && std::find(header.names.cbegin(), header.names.cend(), name) == header.names.cend())
282
            {
283
                header.names.emplace_back(name);
284
            }
285

286
            if (pos == std::string_view::npos || header.names.size() == StarDatabase::MAX_STAR_NAMES)
287
                break;
288

289
            remaining = remaining.substr(pos + 1);
290
        }
291

292
        tokenizer.nextToken();
293
    }
294
    else if (header.catalogNumber == AstroCatalog::InvalidIndex)
295
    {
296
        stcError(header, _("entry missing name and catalog number"));
297
        return false;
298
    }
299

300
    return true;
301
}
302

303
bool
304
checkSpectralType(const StarDatabaseBuilder::StcHeader& header,
305
                  const AssociativeArray* starData,
306
                  const Star* star,
307
                  boost::intrusive_ptr<StarDetails>& newDetails)
308
{
309
    const std::string* spectralType = starData->getString("SpectralType");
310
    if (!header.isStar)
311
    {
312
        if (spectralType != nullptr)
313
            stcWarn(header, _("ignoring SpectralType on Barycenter"));
314
        newDetails = StarDetails::GetBarycenterDetails();
315
    }
316
    else if (spectralType != nullptr)
317
    {
318
        newDetails = StarDetails::GetStarDetails(StellarClass::parse(*spectralType));
319
        if (newDetails == nullptr)
320
        {
321
            stcError(header, _("invalid SpectralType"));
322
            return false;
323
        }
324
    }
325
    else if (header.disposition != DataDisposition::Modify || star->isBarycenter())
326
    {
327
        stcError(header, _("missing SpectralType on Star"));
328
        return false;
329
    }
330

331
    return true;
332
}
333

334
bool
335
checkPolarCoordinates(const StarDatabaseBuilder::StcHeader& header,
336
                      const AssociativeArray* starData,
337
                      const Star* star,
338
                      std::optional<Eigen::Vector3f>& position)
339
{
340
    constexpr unsigned int has_ra = 1;
341
    constexpr unsigned int has_dec = 2;
342
    constexpr unsigned int has_distance = 4;
343
    constexpr unsigned int has_all = has_ra | has_dec | has_distance;
344

345
    unsigned int status = 0;
346

347
    auto raValue = starData->getAngle<double>("RA", astro::DEG_PER_HRA, 1.0);
348
    auto decValue = starData->getAngle<double>("Dec");
349
    auto distanceValue = starData->getLength<double>("Distance", astro::KM_PER_LY<double>);
350
    status = (static_cast<unsigned int>(raValue.has_value()) * has_ra)
351
           | (static_cast<unsigned int>(decValue.has_value()) * has_dec)
352
           | (static_cast<unsigned int>(distanceValue.has_value()) * has_distance);
353

354
    if (status == 0)
355
        return true;
356

357
    if (status == has_all)
358
    {
359
        position = astro::equatorialToCelestialCart(*raValue, *decValue, *distanceValue).cast<float>();
360
        return true;
361
    }
362

363
    if (header.disposition != DataDisposition::Modify)
364
    {
365
        stcError(header, _("incomplete set of coordinates RA/Dec/Distance specified"));
366
        return false;
367
    }
368

369
    // Partial modification of polar coordinates
370
    assert(star != nullptr);
371

372
    // Convert from Celestia's coordinate system
373
    const Eigen::Vector3f& p = star->getPosition();
374
    Eigen::Vector3d v = math::XRotation(math::degToRad(astro::J2000Obliquity)) * Eigen::Vector3f(p.x(), -p.z(), p.y()).cast<double>();
375
    // Disable Sonar on the below: suggests using value-or which would eagerly-evaluate the replacement value
376
    double distance = distanceValue.has_value() ? *distanceValue : v.norm(); //NOSONAR
377
    double ra = raValue.has_value() ? *raValue : (math::radToDeg(std::atan2(v.y(), v.x())) / astro::DEG_PER_HRA); //NOSONAR
378
    double dec = decValue.has_value() ? *decValue : math::radToDeg(std::asin(std::clamp(v.z(), -1.0, 1.0))); //NOSONAR
379

380
    position = astro::equatorialToCelestialCart(ra, dec, distance).cast<float>();
381
    return true;
382
}
383

384
bool
385
checkMagnitudes(const StarDatabaseBuilder::StcHeader& header,
386
                const AssociativeArray* starData,
387
                const Star* star,
388
                float distance,
389
                std::optional<float>& absMagnitude,
390
                std::optional<float>& extinction)
391
{
392
    assert(header.disposition != DataDisposition::Modify || star != nullptr);
393
    absMagnitude = starData->getNumber<float>("AbsMag");
394
    auto appMagnitude = starData->getNumber<float>("AppMag");
395

396
    if (!header.isStar)
397
    {
398
        if (absMagnitude.has_value())
399
            stcWarn(header, _("AbsMag ignored on Barycenter"));
400
        if (appMagnitude.has_value())
401
            stcWarn(header, _("AppMag ignored on Barycenter"));
402
        absMagnitude = 30.0f;
403
        return true;
404
    }
405

406
    extinction = starData->getNumber<float>("Extinction");
407
    if (extinction.has_value() && distance < VALID_APPMAG_DISTANCE_THRESHOLD)
408
    {
409
        stcWarn(header, _("Extinction ignored for stars close to the origin"));
410
        extinction = std::nullopt;
411
    }
412

413
    if (absMagnitude.has_value())
414
    {
415
        if (appMagnitude.has_value())
416
            stcWarn(header, _("AppMag ignored when AbsMag is supplied"));
417
    }
418
    else if (appMagnitude.has_value())
419
    {
420
        if (distance < VALID_APPMAG_DISTANCE_THRESHOLD)
421
        {
422
            stcError(header, _("AppMag cannot be used close to the origin"));
423
            return false;
424
        }
425

426
        float extinctionValue = 0.0;
427
        if (extinction.has_value())
428
            extinctionValue = *extinction;
429
        else if (header.disposition == DataDisposition::Modify)
430
            extinctionValue = star->getExtinction() * distance;
431

432
        absMagnitude = astro::appToAbsMag(*appMagnitude, distance) - extinctionValue;
433
    }
434
    else if (header.disposition != DataDisposition::Modify || star->isBarycenter())
435
    {
436
        stcError(header, _("no magnitude defined for star"));
437
        return false;
438
    }
439

440
    return true;
441
}
442

443
void
444
mergeStarDetails(boost::intrusive_ptr<StarDetails>& existingDetails,
445
                 const boost::intrusive_ptr<StarDetails>& referenceDetails)
446
{
447
    if (referenceDetails == nullptr)
448
        return;
449

450
    if (existingDetails->shared())
451
    {
452
        // If there are no extended information values set, just
453
        // use the new reference details object
454
        existingDetails = referenceDetails;
455
    }
456
    else
457
    {
458
        // There are custom details: copy the new data into the
459
        // existing record
460
        existingDetails->mergeFromStandard(referenceDetails.get());
461
    }
462
}
463

464
void
465
applyTemperatureBoloCorrection(const StarDatabaseBuilder::StcHeader& header,
466
                               const AssociativeArray* starData,
467
                               boost::intrusive_ptr<StarDetails>& details)
468
{
469
    auto bolometricCorrection = starData->getNumber<float>("BoloCorrection");
470
    if (bolometricCorrection.has_value())
471
    {
472
        if (!header.isStar)
473
            stcWarn(header, _("BoloCorrection is ignored on Barycenters"));
474
        else
475
            StarDetails::setBolometricCorrection(details, *bolometricCorrection);
476
    }
477

478
    if (auto temperature = starData->getNumber<float>("Temperature"); temperature.has_value())
479
    {
480
        if (!header.isStar)
481
        {
482
            stcWarn(header, _("Temperature is ignored on Barycenters"));
483
        }
484
        else if (*temperature > 0.0)
485
        {
486
            StarDetails::setTemperature(details, *temperature);
487
            if (!bolometricCorrection.has_value())
488
            {
489
                // if we change the temperature, recalculate the bolometric
490
                // correction using formula from formula for main sequence
491
                // stars given in B. Cameron Reed (1998), "The Composite
492
                // Observational-Theoretical HR Diagram", Journal of the Royal
493
                // Astronomical Society of Canada, Vol 92. p36.
494

495
                double logT = std::log10(static_cast<double>(*temperature)) - 4.0;
496
                double bc = -8.499 * std::pow(logT, 4) + 13.421 * std::pow(logT, 3)
497
                            - 8.131 * logT * logT - 3.901 * logT - 0.438;
498

499
                StarDetails::setBolometricCorrection(details, static_cast<float>(bc));
500
            }
501
        }
502
        else
503
        {
504
            stcWarn(header, _("Temperature value must be greater than zero"));
505
        }
506
    }
507
}
508

509
void
510
applyCustomDetails(const StarDatabaseBuilder::StcHeader& header,
511
                   const AssociativeArray* starData,
512
                   boost::intrusive_ptr<StarDetails>& details)
513
{
514
    if (const auto* mesh = starData->getString("Mesh"); mesh != nullptr)
515
    {
516
        if (!header.isStar)
517
        {
518
            stcWarn(header, _("Mesh is ignored on Barycenters"));
519
        }
520
        else if (auto meshPath = util::U8FileName(*mesh); meshPath.has_value())
521
        {
522
            using engine::GeometryInfo;
523
            using engine::GetGeometryManager;
524
            ResourceHandle geometryHandle = GetGeometryManager()->getHandle(GeometryInfo(*meshPath,
525
                                                                                         *header.path,
526
                                                                                         Eigen::Vector3f::Zero(),
527
                                                                                         1.0f,
528
                                                                                         true));
529
            StarDetails::setGeometry(details, geometryHandle);
530
        }
531
        else
532
        {
533
            stcError(header, _("invalid filename in Mesh"));
534
        }
535
    }
536

537
    if (const auto* texture = starData->getString("Texture"); texture != nullptr)
538
    {
539
        if (!header.isStar)
540
            stcWarn(header, _("Texture is ignored on Barycenters"));
541
        else if (auto texturePath = util::U8FileName(*texture); texturePath.has_value())
542
            StarDetails::setTexture(details, MultiResTexture(*texturePath, *header.path));
543
        else
544
            stcError(header, _("invalid filename in Texture"));
545
    }
546

547
    if (auto rotationModel = CreateRotationModel(starData, *header.path, 1.0); rotationModel != nullptr)
548
    {
549
        if (!header.isStar)
550
            stcWarn(header, _("Rotation is ignored on Barycenters"));
551
        else
552
            StarDetails::setRotationModel(details, rotationModel);
553
    }
554

555
    if (auto semiAxes = starData->getLengthVector<float>("SemiAxes"); semiAxes.has_value())
556
    {
557
        if (!header.isStar)
558
            stcWarn(header, _("SemiAxes is ignored on Barycenters"));
559
        else if (semiAxes->minCoeff() >= 0.0)
560
            StarDetails::setEllipsoidSemiAxes(details, *semiAxes);
561
        else
562
            stcWarn(header, _("SemiAxes must be greater than zero"));
563
    }
564

565
    if (auto radius = starData->getLength<float>("Radius"); radius.has_value())
566
    {
567
        if (!header.isStar)
568
            stcWarn(header, _("Radius is ignored on Barycenters"));
569
        else if (*radius >= 0.0)
570
            StarDetails::setRadius(details, *radius);
571
        else
572
            stcWarn(header, _("Radius must be greater than zero"));
573
    }
574

575
    applyTemperatureBoloCorrection(header, starData, details);
576

577
    if (const auto* infoUrl = starData->getString("InfoURL"); infoUrl != nullptr)
578
        StarDetails::setInfoURL(details, *infoUrl);
579
}
580

581
} // end unnamed namespace
582

583
StarDatabaseBuilder::~StarDatabaseBuilder() = default;
584

585
bool
586
StarDatabaseBuilder::loadBinary(std::istream& in)
587
{
588
    Timer timer;
589
    std::uint32_t nStarsInFile;
590
    if (!parseStarsDatHeader(in, nStarsInFile))
591
        return false;
592

593
    constexpr std::uint32_t BUFFER_RECORDS = UINT32_C(4096) / sizeof(StarsDatRecord);
594
    std::vector<char> buffer(sizeof(StarsDatRecord) * BUFFER_RECORDS);
595
    std::uint32_t nStarsRemaining = nStarsInFile;
596
    while (nStarsRemaining > 0)
597
    {
598
        std::uint32_t recordsToRead = std::min(BUFFER_RECORDS, nStarsRemaining);
599
        if (!in.read(buffer.data(), sizeof(StarsDatRecord) * recordsToRead).good()) /* Flawfinder: ignore */
600
            return false;
601

602
        const char* ptr = buffer.data();
603
        for (std::uint32_t i = 0; i < recordsToRead; ++i)
604
        {
605
            auto catNo = util::fromMemoryLE<AstroCatalog::IndexNumber>(ptr + offsetof(StarsDatRecord, catNo));
606
            Eigen::Vector3f position(util::fromMemoryLE<float>(ptr + offsetof(StarsDatRecord, x)),
607
                                     util::fromMemoryLE<float>(ptr + offsetof(StarsDatRecord, y)),
608
                                     util::fromMemoryLE<float>(ptr + offsetof(StarsDatRecord, z)));
609
            auto absMag = util::fromMemoryLE<std::int16_t>(ptr + offsetof(StarsDatRecord, absMag));
610
            auto spectralType = util::fromMemoryLE<std::uint16_t>(ptr + offsetof(StarsDatRecord, spectralType));
611

612
            boost::intrusive_ptr<StarDetails> details = nullptr;
613
            if (StellarClass sc; sc.unpackV1(spectralType))
614
                details = StarDetails::GetStarDetails(sc);
615

616
            if (details == nullptr)
617
            {
618
                GetLogger()->error(_("Bad spectral type in star database, star #{}\n"), catNo);
619
                continue;
620
            }
621

622
            Star& star = unsortedStars.emplace_back(catNo, details);
623
            star.setPosition(position);
624
            star.setAbsoluteMagnitude(static_cast<float>(absMag) / 256.0f);
625

626
            ptr += sizeof(StarsDatRecord);
627
        }
628

629
        nStarsRemaining -= recordsToRead;
630
    }
631

632
    if (in.bad())
633
        return false;
634

635
    auto loadTime = timer.getTime();
636

637
    GetLogger()->debug("StarDatabase::read: nStars = {}, time = {} ms\n", nStarsInFile, loadTime);
638
    GetLogger()->info(_("{} stars in binary database\n"), unsortedStars.size());
639

640
    // Create the temporary list of stars sorted by catalog number; this
641
    // will be used to lookup stars during file loading. After loading is
642
    // complete, the stars are sorted into an octree and this list gets
643
    // replaced.
644
    binFileCatalogNumberIndex.reserve(unsortedStars.size());
645
    for (Star& star : unsortedStars)
646
        binFileCatalogNumberIndex.push_back(&star);
647

648
    std::sort(binFileCatalogNumberIndex.begin(), binFileCatalogNumberIndex.end(),
649
                [](const Star* star0, const Star* star1) { return star0->getIndex() < star1->getIndex(); });
650

651
    return true;
652
}
653

654
/*! Load an STC file with star definitions. Each definition has the form:
655
 *
656
 *  [disposition] [object type] [catalog number] [name]
657
 *  {
658
 *      [properties]
659
 *  }
660
 *
661
 *  Disposition is either Add, Replace, or Modify; Add is the default.
662
 *  Object type is either Star or Barycenter, with Star the default
663
 *  It is an error to omit both the catalog number and the name.
664
 *
665
 *  The dispositions are slightly more complicated than suggested by
666
 *  their names. Every star must have an unique catalog number. But
667
 *  instead of generating an error, Adding a star with a catalog
668
 *  number that already exists will actually replace that star. Here
669
 *  are how all of the possibilities are handled:
670
 *
671
 *  <name> or <number> already exists:
672
 *  Add <name>        : new star
673
 *  Add <number>      : replace star
674
 *  Replace <name>    : replace star
675
 *  Replace <number>  : replace star
676
 *  Modify <name>     : modify star
677
 *  Modify <number>   : modify star
678
 *
679
 *  <name> or <number> doesn't exist:
680
 *  Add <name>        : new star
681
 *  Add <number>      : new star
682
 *  Replace <name>    : new star
683
 *  Replace <number>  : new star
684
 *  Modify <name>     : error
685
 *  Modify <number>   : error
686
 */
687
bool
688
StarDatabaseBuilder::load(std::istream& in, const fs::path& resourcePath)
689
{
690
    Tokenizer tokenizer(&in);
691
    Parser parser(&tokenizer);
692

693
#ifdef ENABLE_NLS
694
    std::string domain = resourcePath.string();
695
    const char *d = domain.c_str();
696
    bindtextdomain(d, d); // domain name is the same as resource path
697
#else
698
    std::string domain;
699
#endif
700

701
    StcHeader header(resourcePath);
702
    while (tokenizer.nextToken() != Tokenizer::TokenEnd)
703
    {
704
        if (!parseStcHeader(tokenizer, header))
705
            return false;
706

707
        // now goes the star definition
708
        tokenizer.pushBack();
709
        const Value starDataValue = parser.readValue();
710
        const Hash* starData = starDataValue.getHash();
711
        if (starData == nullptr)
712
        {
713
            GetLogger()->error(_("Bad star definition at line {}.\n"), tokenizer.getLineNumber());
714
            return false;
715
        }
716

717
        if (header.disposition != DataDisposition::Add && header.catalogNumber == AstroCatalog::InvalidIndex)
718
            header.catalogNumber = starDB->namesDB->findCatalogNumberByName(header.names.front(), false);
719

720
        Star* star = findWhileLoading(header.catalogNumber);
721
        if (star == nullptr)
722
        {
723
            if (header.disposition == DataDisposition::Modify)
724
            {
725
                GetLogger()->error(_("Modify requested for nonexistent star.\n"));
726
                continue;
727
            }
728

729
            if (header.catalogNumber == AstroCatalog::InvalidIndex)
730
            {
731
                header.catalogNumber = nextAutoCatalogNumber;
732
                --nextAutoCatalogNumber;
733
            }
734
        }
735

736
        if (createOrUpdateStar(header, starData, star))
737
        {
738
            loadCategories(header, starData, domain);
739

740
            if (!header.names.empty())
741
            {
742
                starDB->namesDB->erase(header.catalogNumber);
743
                for (const auto& name : header.names)
744
                    starDB->namesDB->add(header.catalogNumber, name);
745
            }
746
        }
747
    }
748

749
    return true;
750
}
751

752
void
753
StarDatabaseBuilder::setNameDatabase(std::unique_ptr<StarNameDatabase>&& nameDB)
754
{
755
    starDB->namesDB = std::move(nameDB);
756
}
757

758
std::unique_ptr<StarDatabase>
759
StarDatabaseBuilder::finish()
760
{
761
    GetLogger()->info(_("Total star count: {}\n"), unsortedStars.size());
762

763
    buildOctree();
764
    buildIndexes();
765

766
    // Resolve all barycenters; this can't be done before star sorting. There's
767
    // still a bug here: final orbital radii aren't available until after
768
    // the barycenters have been resolved, and these are required when building
769
    // the octree.  This will only rarely cause a problem, but it still needs
770
    // to be addressed.
771
    for (const auto [starIdx, barycenterIdx] : barycenters)
772
    {
773
        Star* star = starDB->find(starIdx);
774
        Star* barycenter = starDB->find(barycenterIdx);
775
        assert(star != nullptr);
776
        assert(barycenter != nullptr);
777
        if (star != nullptr && barycenter != nullptr)
778
        {
779
            StarDetails::setOrbitBarycenter(star->details, barycenter);
780
            StarDetails::addOrbitingStar(barycenter->details, star);
781
        }
782
    }
783

784
    for (const auto& [catalogNumber, category] : categories)
785
    {
786
        Star* star = starDB->find(catalogNumber);
787
        UserCategory::addObject(star, category);
788
    }
789

790
    return std::move(starDB);
791
}
792

793
/*! Load star data from a property list into a star instance.
794
 */
795
bool
796
StarDatabaseBuilder::createOrUpdateStar(const StcHeader& header,
797
                                        const AssociativeArray* starData,
798
                                        Star* star)
799
{
800
    boost::intrusive_ptr<StarDetails> newDetails = nullptr;
801
    if (!checkSpectralType(header, starData, star, newDetails))
802
        return false;
803

804
    std::optional<Eigen::Vector3f> position = std::nullopt;
805
    std::optional<AstroCatalog::IndexNumber> barycenterNumber = std::nullopt;
806
    std::shared_ptr<const ephem::Orbit> orbit = nullptr;
807
    if (!checkStcPosition(header, starData, star, position, barycenterNumber, orbit))
808
        return false;
809

810
    std::optional<float> absMagnitude = std::nullopt;
811
    std::optional<float> extinction = std::nullopt;
812
    float distance;
813
    if (position.has_value())
814
    {
815
        distance = position->norm();
816
    }
817
    else
818
    {
819
        assert(star != nullptr);
820
        distance = star->getPosition().norm();
821
    }
822

823
    if (!checkMagnitudes(header, starData, star, distance, absMagnitude, extinction))
824
        return false;
825

826
    if (star == nullptr)
827
    {
828
        assert(newDetails != nullptr);
829
        star = &unsortedStars.emplace_back(header.catalogNumber, newDetails);
830
        stcFileCatalogNumberIndex[header.catalogNumber] = star;
831
    }
832
    else if (header.disposition == DataDisposition::Modify)
833
    {
834
        mergeStarDetails(star->details, newDetails);
835
    }
836
    else
837
    {
838
        assert(newDetails != nullptr);
839
        star->details = newDetails;
840
    }
841

842
    if (position.has_value())
843
        star->setPosition(*position);
844

845
    if (absMagnitude.has_value())
846
        star->setAbsoluteMagnitude(*absMagnitude);
847

848
    if (extinction.has_value())
849
        star->setExtinction(*extinction / distance);
850

851
    if (barycenterNumber == AstroCatalog::InvalidIndex)
852
        barycenters.erase(header.catalogNumber);
853
    else if (barycenterNumber.has_value())
854
        barycenters[header.catalogNumber] = *barycenterNumber;
855

856
    if (orbit != nullptr)
857
        StarDetails::setOrbit(star->details, orbit);
858

859
    applyCustomDetails(header, starData, star->details);
860
    return true;
861
}
862

863
bool
864
StarDatabaseBuilder::checkStcPosition(const StarDatabaseBuilder::StcHeader& header,
865
                                      const AssociativeArray* starData,
866
                                      const Star* star,
867
                                      std::optional<Eigen::Vector3f>& position,
868
                                      std::optional<AstroCatalog::IndexNumber>& barycenterNumber,
869
                                      std::shared_ptr<const ephem::Orbit>& orbit) const
870
{
871
    position = std::nullopt;
872
    barycenterNumber = std::nullopt;
873

874
    if (!checkPolarCoordinates(header, starData, star, position))
875
        return false;
876

877
    if (auto positionValue = starData->getLengthVector<float>("Position", astro::KM_PER_LY<double>);
878
        positionValue.has_value())
879
    {
880
        if (position.has_value())
881
            stcWarn(header, _("ignoring RA/Dec/Distance in favor of Position"));
882
        position = *positionValue;
883
    }
884

885
    if (!checkBarycenter(header, starData, position, barycenterNumber))
886
        return false;
887

888
    // we consider a star to have a barycenter if it has an OrbitBarycenter defined
889
    // or the star is modified without overriding its position, and it has no other
890
    // position overrides.
891
    bool hasBarycenter = (barycenterNumber.has_value() && *barycenterNumber != AstroCatalog::InvalidIndex)
892
                         || (header.disposition == DataDisposition::Modify
893
                             && !position.has_value()
894
                             && barycenters.find(header.catalogNumber) != barycenters.end());
895

896
    if (auto newOrbit = CreateOrbit(Selection(), starData, *header.path, true); newOrbit != nullptr)
897
    {
898
        if (hasBarycenter)
899
            orbit = std::move(newOrbit);
900
        else
901
            stcWarn(header, _("ignoring orbit for object without OrbitBarycenter"));
902
    }
903
    else if (hasBarycenter && star != nullptr && star->getOrbit() == nullptr)
904
    {
905
        stcError(header, _("no orbit specified for star with OrbitBarycenter"));
906
        return false;
907
    }
908

909
    return true;
910
}
911

912
bool
913
StarDatabaseBuilder::checkBarycenter(const StarDatabaseBuilder::StcHeader& header,
914
                                     const AssociativeArray* starData,
915
                                     std::optional<Eigen::Vector3f>& position,
916
                                     std::optional<AstroCatalog::IndexNumber>& barycenterNumber) const
917
{
918
    // If we override RA/Dec/Position, remove the barycenter
919
    if (position.has_value())
920
        barycenterNumber = AstroCatalog::InvalidIndex;
921

922
    const Value* orbitBarycenterValue = starData->getValue("OrbitBarycenter");
923
    if (orbitBarycenterValue == nullptr)
924
        return true;
925

926
    if (auto bcNumber = orbitBarycenterValue->getNumber(); bcNumber.has_value())
927
    {
928
        barycenterNumber = static_cast<AstroCatalog::IndexNumber>(*bcNumber);
929
    }
930
    else if (auto bcName = orbitBarycenterValue->getString(); bcName != nullptr)
931
    {
932
        barycenterNumber = starDB->namesDB->findCatalogNumberByName(*bcName, false);
933
    }
934
    else
935
    {
936
        stcError(header, _("OrbitBarycenter should be either a string or an integer"));
937
        return false;
938
    }
939

940
    if (*barycenterNumber == header.catalogNumber)
941
    {
942
        stcError(header, _("OrbitBarycenter cycle detected"));
943
        return false;
944
    }
945

946
    if (const Star* barycenter = findWhileLoading(*barycenterNumber); barycenter != nullptr)
947
    {
948
        if (position.has_value())
949
            stcWarn(header, "ignoring stellar coordinates in favor of OrbitBarycenter");
950
        position = barycenter->getPosition();
951
    }
952
    else
953
    {
954
        stcError(header, _("OrbitBarycenter refers to nonexistent star"));
955
        return false;
956
    }
957

958
    for (auto it = barycenters.find(*barycenterNumber); it != barycenters.end(); it = barycenters.find(it->second))
959
    {
960
        if (it->second == header.catalogNumber)
961
        {
962
            stcError(header, _("OrbitBarycenter cycle detected"));
963
            return false;
964
        }
965
    }
966

967
    return true;
968
}
969

970
void
971
StarDatabaseBuilder::loadCategories(const StcHeader& header,
972
                                    const AssociativeArray* starData,
973
                                    const std::string& domain)
974
{
975
    if (header.disposition == DataDisposition::Replace)
976
        categories.erase(header.catalogNumber);
977

978
    const Value* categoryValue = starData->getValue("Category");
979
    if (categoryValue == nullptr)
980
        return;
981

982
    if (const std::string* categoryName = categoryValue->getString(); categoryName != nullptr)
983
    {
984
        if (categoryName->empty())
985
            return;
986

987
        addCategory(header.catalogNumber, *categoryName, domain);
988
        return;
989
    }
990

991
    const ValueArray *arr = categoryValue->getArray();
992
    if (arr == nullptr)
993
        return;
994

995
    for (const auto& it : *arr)
996
    {
997
        const std::string* categoryName = it.getString();
998
        if (categoryName == nullptr || categoryName->empty())
999
            continue;
1000

1001
        addCategory(header.catalogNumber, *categoryName, domain);
1002
    }
1003
}
1004

1005
void
1006
StarDatabaseBuilder::addCategory(AstroCatalog::IndexNumber catalogNumber,
1007
                                 const std::string& name,
1008
                                 const std::string& domain)
1009
{
1010
    auto category = UserCategory::findOrAdd(name, domain);
1011
    if (category == UserCategoryId::Invalid)
1012
        return;
1013

1014
    auto [start, end] = categories.equal_range(catalogNumber);
1015
    if (start == end)
1016
    {
1017
        categories.emplace(catalogNumber, category);
1018
        return;
1019
    }
1020

1021
    if (std::any_of(start, end, [category](const auto& it) { return it.second == category; }))
1022
        return;
1023

1024
    categories.emplace_hint(end, catalogNumber, category);
1025
}
1026

1027
/*! While loading the star catalogs, this function must be called instead of
1028
 *  find(). The final catalog number index for stars cannot be built until
1029
 *  after all stars have been loaded. During catalog loading, there are two
1030
 *  separate indexes: one for the binary catalog and another index for stars
1031
 *  loaded from stc files. They binary catalog index is a sorted array, while
1032
 *  the stc catalog index is an STL map. Since the binary file can be quite
1033
 *  large, we want to avoid creating a map with as many nodes as there are
1034
 *  stars. Stc files should collectively contain many fewer stars, and stars
1035
 *  in an stc file may reference each other (barycenters). Thus, a dynamic
1036
 *  structure like a map is both practical and essential.
1037
 */
1038
Star*
1039
StarDatabaseBuilder::findWhileLoading(AstroCatalog::IndexNumber catalogNumber) const
1040
{
1041
    if (catalogNumber == AstroCatalog::InvalidIndex)
1042
        return nullptr;
1043

1044
    // First check for stars loaded from the binary database
1045
    if (auto it = std::lower_bound(binFileCatalogNumberIndex.cbegin(), binFileCatalogNumberIndex.cend(),
1046
                                   catalogNumber,
1047
                                   [](const Star* star, AstroCatalog::IndexNumber catNum) { return star->getIndex() < catNum; });
1048
        it != binFileCatalogNumberIndex.cend() && (*it)->getIndex() == catalogNumber)
1049
    {
1050
        return *it;
1051
    }
1052

1053
    // Next check for stars loaded from an stc file
1054
    if (auto it = stcFileCatalogNumberIndex.find(catalogNumber); it != stcFileCatalogNumberIndex.end())
1055
        return it->second;
1056

1057
    // Star not found
1058
    return nullptr;
1059
}
1060

1061
void
1062
StarDatabaseBuilder::buildOctree()
1063
{
1064
    // This should only be called once for the database
1065
    GetLogger()->debug("Sorting stars into octree . . .\n");
1066
    auto starCount = static_cast<engine::OctreeObjectIndex>(unsortedStars.size());
1067

1068
    float absMag = astro::appToAbsMag(STAR_OCTREE_MAGNITUDE,
1069
                                      StarDatabase::STAR_OCTREE_ROOT_SIZE * celestia::numbers::sqrt3_v<float>);
1070

1071
    auto root = engine::makeDynamicOctree<StarOctreeTraits>(std::move(unsortedStars),
1072
                                                            Eigen::Vector3f(1000.0f, 1000.0f, 1000.0f),
1073
                                                            StarDatabase::STAR_OCTREE_ROOT_SIZE,
1074
                                                            absMag,
1075
                                                            StarOctreeSplitThreshold);
1076

1077
    GetLogger()->debug("Spatially sorting stars for improved locality of reference . . .\n");
1078
    starDB->octreeRoot = root->build();
1079

1080
    GetLogger()->debug("{} stars total\nOctree has {} nodes and {} stars.\n",
1081
                       starCount,
1082
                       starDB->octreeRoot->nodeCount(),
1083
                       starDB->octreeRoot->size());
1084

1085
    unsortedStars.clear();
1086
}
1087

1088
void
1089
StarDatabaseBuilder::buildIndexes()
1090
{
1091
    // This should only be called once for the database
1092
    // assert(catalogNumberIndexes[0] == nullptr);
1093

1094
    GetLogger()->info("Building catalog number indexes . . .\n");
1095

1096
    auto nStars = starDB->octreeRoot->size();
1097

1098
    starDB->catalogNumberIndex.clear();
1099
    starDB->catalogNumberIndex.reserve(nStars);
1100
    for (std::uint32_t i = 0; i < nStars; ++i)
1101
        starDB->catalogNumberIndex.push_back(i);
1102

1103
    const auto& octreeRoot = *starDB->octreeRoot;
1104
    std::sort(starDB->catalogNumberIndex.begin(), starDB->catalogNumberIndex.end(),
1105
              [&octreeRoot](std::uint32_t idx0, std::uint32_t idx1)
1106
              {
1107
                  return octreeRoot[idx0].getIndex() < octreeRoot[idx1].getIndex();
1108
              });
1109
}
1110