xypor/index.php

/*

    SPDX-FileCopyrightText: 2001 Jason Harris <jharris@30doradus.org>


    SPDX-License-Identifier: GPL-2.0-or-later

*/


#include "ksasteroid.h"


#include "dms.h"

#include "ksnumbers.h"

#include "Options.h"

#ifdef KSTARS_LITE

#include "skymaplite.h"

#else

#include "skymap.h"

#endif


#include <qdebug.h>


#include <typeinfo>


KSAsteroid::KSAsteroid(int _catN, const QString &s, const QString &imfile, long double _JD, double _a, double _e,

                       dms _i, dms _w, dms _Node, dms _M, double _H, double _G)

    : KSPlanetBase(s, imfile), catN(_catN), JD(_JD), a(_a), e(_e), i(_i), w(_w), M(_M), N(_Node), H(_H), G(_G)

{

    setType(SkyObject::ASTEROID);

    setMag(H);

    //Compute the orbital Period from Kepler's 3rd law:

    P = 365.2568984 * pow(a, 1.5); //period in days

}


KSAsteroid *KSAsteroid::clone() const

{

    Q_ASSERT(typeid(this) ==

             typeid(static_cast<const KSAsteroid *>(this))); // Ensure we are not slicing a derived class

    return new KSAsteroid(*this);

}


bool KSAsteroid::findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth)

{

    // TODO: (Valentin) TOP LEVEL CONTROL OF CALCULATION FOR ALL OBJECTS

    if(!toCalculate()){

        return false;

    }


    //determine the mean anomaly for the desired date.  This is the mean anomaly for the

    //ephemeis epoch, plus the number of days between the desired date and ephemeris epoch,

    //times the asteroid's mean daily motion (360/P):


    // All elements are in the heliocentric ecliptic J2000 reference frame.

    // Mean anomaly is supplied at the Epoch (which is JD here)


    dms m = dms(double(M.Degrees() + (lastPrecessJD - JD) * 360.0 / P)).reduce();

    double sinm, cosm;

    m.SinCos(sinm, cosm);


    //compute eccentric anomaly:

    double E = m.Degrees() + e * 180.0 / dms::PI * sinm * (1.0 + e * cosm);


    if (e > 0.005) //need more accurate approximation, iterate...

    {

        double E0;

        int iter(0);

        do

        {

            E0 = E;

            iter++;

            E = E0 -

                (E0 - e * 180.0 / dms::PI * sin(E0 * dms::DegToRad) - m.Degrees()) / (1 - e * cos(E0 * dms::DegToRad));

        } while (fabs(E - E0) > 0.001 && iter < 1000);

    }


    // Assert that the solution of the Kepler equation E = M + e sin E is accurate to about 0.1 arcsecond

    //Q_ASSERT( fabs( E - ( m.Degrees() + ( e * 180.0 / dms::PI ) * sin( E * dms::DegToRad ) ) ) < 0.10/3600.0 );


    double sinE, cosE;

    dms E1(E);

    E1.SinCos(sinE, cosE);


    double xv = a * (cosE - e);

    double yv = a * sqrt(1.0 - e * e) * sinE;


    //v is the true anomaly; r is the distance from the Sun

    double v = atan2(yv, xv) / dms::DegToRad;

    double r = sqrt(xv * xv + yv * yv);


    //vw is the sum of the true anomaly and the argument of perihelion

    dms vw(v + w.Degrees());

    double sinN, cosN, sinvw, cosvw, sini, cosi;


    N.SinCos(sinN, cosN);

    vw.SinCos(sinvw, cosvw);

    i.SinCos(sini, cosi);


    //xh, yh, zh are the heliocentric cartesian coords with the ecliptic plane congruent with zh=0.

    double xh = r * (cosN * cosvw - sinN * sinvw * cosi);

    double yh = r * (sinN * cosvw + cosN * sinvw * cosi);

    double zh = r * (sinvw * sini);


    //the spherical ecliptic coordinates:

    double ELongRad = atan2(yh, xh);

    double ELatRad  = atan2(zh, r);


    helEcPos.longitude.setRadians(ELongRad);

    helEcPos.longitude.reduceToRange(dms::ZERO_TO_2PI);

    helEcPos.latitude.setRadians(ELatRad);

    setRsun(r);


    if (Earth)

    {

        //xe, ye, ze are the Earth's heliocentric cartesian coords

        double cosBe, sinBe, cosLe, sinLe;

        Earth->ecLong().SinCos(sinLe, cosLe);

        Earth->ecLat().SinCos(sinBe, cosBe);


        double xe = Earth->rsun() * cosBe * cosLe;

        double ye = Earth->rsun() * cosBe * sinLe;

        double ze = Earth->rsun() * sinBe;


        //convert to geocentric ecliptic coordinates by subtracting Earth's coords:

        xh -= xe;

        yh -= ye;

        zh -= ze;

    }


    //the spherical geocentricecliptic coordinates:

    ELongRad  = atan2(yh, xh);

    double rr = sqrt(xh * xh + yh * yh);

    ELatRad   = atan2(zh, rr);


    ep.longitude.setRadians(ELongRad);

    ep.longitude.reduceToRange(dms::ZERO_TO_2PI);

    ep.latitude.setRadians(ELatRad);

    if (Earth)

        setRearth(Earth);


    EclipticToEquatorial(num->obliquity());


    // JM 2017-09-10: The calculations above produce J2000 RESULTS

    // So we have to precess as well

    setRA0(ra());

    setDec0(dec());

    apparentCoord(J2000, lastPrecessJD);

    //nutate(num);

    //aberrate(num);


    return true;

}


void KSAsteroid::findMagnitude(const KSNumbers *)

{

    double param     = 5.0 * log10(rsun() * rearth());

    double phase_rad = phase().radians();

    double phi1      = exp(-3.33 * pow(tan(phase_rad / 2.0), 0.63));

    double phi2      = exp(-0.187 * pow(tan(phase_rad / 2.0), 1.22));


    setMag(H + param - 2.5 * log10((1.0 - G) * phi1 + G * phi2));

}


void KSAsteroid::setPerihelion(double perihelion)

{

    q = perihelion;

}


void KSAsteroid::setEarthMOID(double earth_moid)

{

    EarthMOID = earth_moid;

}


void KSAsteroid::setAlbedo(float albedo)

{

    Albedo = albedo;

}


void KSAsteroid::setDiameter(float diam)

{

    Diameter = diam;

}


void KSAsteroid::setDimensions(QString dim)

{

    Dimensions = dim;

}


void KSAsteroid::setNEO(bool neo)

{

    NEO = neo;

}


void KSAsteroid::setOrbitClass(QString orbit_class)

{

    OrbitClass = orbit_class;

}


void KSAsteroid::setOrbitID(QString orbit_id)

{

    OrbitID = orbit_id;

}


void KSAsteroid::setPeriod(float per)

{

    Period = per;

}


bool KSAsteroid::toCalculate()

{

    // Filter by magnitude, but draw focused asteroids anyway :)

    return ((mag() < Options::magLimitAsteroid())|| (std::isnan(mag()) != 0) ||

#ifdef KSTARS_LITE

            SkyMapLite::Instance()->focusObject() == this

#else

            SkyMap::Instance()->focusObject() == this

#endif

            );

}


QDataStream &operator<<(QDataStream &out, const KSAsteroid &asteroid)

{

    out << asteroid.Name << asteroid.OrbitClass << asteroid.Dimensions << asteroid.OrbitID

        << asteroid.catN << static_cast<double>(asteroid.JD) << asteroid.a

        << asteroid.e << asteroid.i << asteroid.w << asteroid.N

        << asteroid.M << asteroid.H << asteroid.G << asteroid.q

        << asteroid.NEO << asteroid.Diameter

        << asteroid.Albedo << asteroid.RotationPeriod

        << asteroid.Period << asteroid.EarthMOID;

    return out;

}


QDataStream &operator>>(QDataStream &in, KSAsteroid *&asteroid)

{

    QString name, orbit_id, orbit_class, dimensions;

    double q, a, e, H, G, earth_moid;

    dms i, w, N, M;

    double JD;

    float diameter, albedo, rot_period, period;

    bool neo;

    int catN;


    in >> name;

    in >> orbit_class;

    in >> dimensions;

    in >> orbit_id;


    in >> catN >> JD >> a >> e >> i >> w >> N >> M >> H >> G >> q >> neo >> diameter

            >> albedo >> rot_period >> period >> earth_moid;


    asteroid = new KSAsteroid(catN, name, QString(), JD, a, e, i, w, N, M, H, G);

    asteroid->setPerihelion(q);

    asteroid->setOrbitID(orbit_id);

    asteroid->setNEO(neo);

    asteroid->setDiameter(diameter);

    asteroid->setDimensions(dimensions);

    asteroid->setAlbedo(albedo);

    asteroid->setRotationPeriod(rot_period);

    asteroid->setPeriod(period);

    asteroid->setEarthMOID(earth_moid);

    asteroid->setOrbitClass(orbit_class);

    asteroid->setPhysicalSize(diameter);


    return in;

}


void KSAsteroid::setRotationPeriod(float rot_per)

{

    RotationPeriod = rot_per;

}


//Unused virtual function from KSPlanetBase


bool KSAsteroid::loadData()

{

    return false;

}


SkyObject::UID KSAsteroid::getUID() const

{

    return solarsysUID(UID_SOL_ASTEROID) | catN;

}


KSAsteroid
A subclass of KSPlanetBase that implements asteroids.
Definition ksasteroid.h:42

KSAsteroid::getUID
SkyObject::UID getUID() const override
Return UID for object.
Definition ksasteroid.cpp:274

KSAsteroid::setPerihelion
void setPerihelion(double perihelion)
Sets the asteroid's perihelion distance.
Definition ksasteroid.cpp:160

KSAsteroid::clone
KSAsteroid * clone() const override
Create copy of object.
Definition ksasteroid.cpp:32

KSAsteroid::setPeriod
void setPeriod(float per)
Sets the asteroid's period.
Definition ksasteroid.cpp:200

KSAsteroid::setDiameter
void setDiameter(float diam)
Sets the asteroid's diameter.
Definition ksasteroid.cpp:175

KSAsteroid::setDimensions
void setDimensions(QString dim)
Sets the asteroid's dimensions.
Definition ksasteroid.cpp:180

KSAsteroid::setEarthMOID
void setEarthMOID(double earth_moid)
Sets the asteroid's earth minimum orbit intersection distance.
Definition ksasteroid.cpp:165

KSAsteroid::setNEO
void setNEO(bool neo)
Sets if the comet is a near earth object.
Definition ksasteroid.cpp:185

KSAsteroid::setOrbitID
void setOrbitID(QString orbit_id)
Sets the asteroid's orbit solution ID.
Definition ksasteroid.cpp:195

KSAsteroid::toCalculate
bool toCalculate()
toCalculate
Definition ksasteroid.cpp:205

KSAsteroid::setAlbedo
void setAlbedo(float albedo)
Sets the asteroid's albedo.
Definition ksasteroid.cpp:170

KSAsteroid::findGeocentricPosition
bool findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth=nullptr) override
Calculate the geocentric RA, Dec coordinates of the Asteroid.
Definition ksasteroid.cpp:39

KSAsteroid::loadData
bool loadData() override
This is inherited from KSPlanetBase.
Definition ksasteroid.cpp:269

KSAsteroid::KSAsteroid
KSAsteroid(int catN, const QString &s, const QString &image_file, long double JD, double a, double e, dms i, dms w, dms N, dms M, double H, double G)
Constructor.
Definition ksasteroid.cpp:22

KSAsteroid::setRotationPeriod
void setRotationPeriod(float rot_per)
Sets the asteroid's rotation period.
Definition ksasteroid.cpp:263

KSAsteroid::setOrbitClass
void setOrbitClass(QString orbit_class)
Sets the asteroid's orbit class.
Definition ksasteroid.cpp:190

KSNumbers
There are several time-dependent values used in position calculations, that are not specific to an ob...
Definition ksnumbers.h:43

KSNumbers::obliquity
const CachingDms * obliquity() const
Definition ksnumbers.h:56

KSPlanetBase
A subclass of TrailObject that provides additional information needed for most solar system objects.
Definition ksplanetbase.h:50

KSPlanetBase::solarsysUID
UID solarsysUID(UID type) const
Compute high 32-bits of UID.
Definition ksplanetbase.h:223

KSPlanetBase::setPhysicalSize
void setPhysicalSize(double size)
set the planet's physical size, in km.
Definition ksplanetbase.h:197

KSPlanetBase::UID_SOL_ASTEROID
static const UID UID_SOL_ASTEROID
Asteroids.
Definition ksplanetbase.h:218

KSPlanetBase::setRearth
void setRearth(double r)
Set the distance from Earth, in AU.
Definition ksplanetbase.h:145

KSPlanetBase::ecLat
const dms & ecLat() const
Definition ksplanetbase.h:94

KSPlanetBase::ecLong
const dms & ecLong() const
Definition ksplanetbase.h:91

KSPlanetBase::rearth
double rearth() const
Definition ksplanetbase.h:139

KSPlanetBase::EclipticToEquatorial
void EclipticToEquatorial(const CachingDms *Obliquity)
Convert Ecliptic longitude/latitude to Right Ascension/Declination.
Definition ksplanetbase.cpp:81

KSPlanetBase::phase
dms phase()
Definition ksplanetbase.h:200

KSPlanetBase::rsun
double rsun() const
Definition ksplanetbase.h:130

KSPlanetBase::setRsun
void setRsun(double r)
Set the solar distance in AU.
Definition ksplanetbase.h:136

SkyObject::setMag
void setMag(float m)
Set the object's sorting magnitude.
Definition skyobject.h:404

SkyObject::UID
qint64 UID
Type for Unique object IDenticator.
Definition skyobject.h:49

SkyObject::mag
float mag() const
Definition skyobject.h:207

SkyObject::setType
void setType(int t)
Set the object's type identifier to the argument.
Definition skyobject.h:196

SkyPoint::apparentCoord
void apparentCoord(long double jd0, long double jdf)
Computes the apparent coordinates for this SkyPoint for any epoch, accounting for the effects of prec...
Definition skypoint.cpp:700

SkyPoint::dec
const CachingDms & dec() const
Definition skypoint.h:269

SkyPoint::ra
const CachingDms & ra() const
Definition skypoint.h:263

SkyPoint::setRA0
void setRA0(dms r)
Sets RA0, the catalog Right Ascension.
Definition skypoint.h:94

SkyPoint::setDec0
void setDec0(dms d)
Sets Dec0, the catalog Declination.
Definition skypoint.h:119

dms
An angle, stored as degrees, but expressible in many ways.
Definition dms.h:38

dms::PI
static constexpr double PI
PI is a const static member; it's public so that it can be used anywhere, as long as dms....
Definition dms.h:385

dms::DegToRad
static constexpr double DegToRad
DegToRad is a const static member equal to the number of radians in one degree (dms::PI/180....
Definition dms.h:390

QDataStream

QString