libastro.cpp

Go to the documentation of this file.

/*
    libastro
 
    functions used for coordinate conversions, based on libnova
 
    Copyright (C) 2020 Chris Rowland
    Copyright (C) 2021 Jasem Mutlaq
 
    This library is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published
    by the Free Software Foundation; either version 2.1 of the License, or
    (at your option) any later version.
 
    This library is distributed in the hope that it will be useful, but
    WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
    License for more details.
 
    You should have received a copy of the GNU Lesser General Public License
    along with this library; if not, write to the Free Software Foundation,
    Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA
 
*/
 
// file libastro.c
//
// holds extensions to the libnova library to provide useful functions for INDI
//
// Chris Rowland April 2020
//
 
 
#include "libastro.h"
#include "indicom.h"
 
#include <math.h>
 
#include <libnova/precession.h>
#include <libnova/aberration.h>
#include <libnova/transform.h>
#include <libnova/nutation.h>
#include <libnova/sidereal_time.h>
 
namespace INDI
{
 
// converts the Observed (JNow) position to a J2000 catalogue position by removing
// aberration, nutation and precession using the libnova library
void ObservedToJ2000(IEquatorialCoordinates * observed, double jd, IEquatorialCoordinates * J2000pos)
{
    ln_equ_posn tempPos;
    // RA Hours --> Degrees
    struct ln_equ_posn libnova_observed = {observed->rightascension * 15.0, observed->declination};
    // remove the aberration
    ln_get_equ_aber(&libnova_observed, jd, &tempPos);
    // this conversion has added the aberration, we want to subtract it
    tempPos.ra = libnova_observed.ra - (tempPos.ra - libnova_observed.ra);
    tempPos.dec = libnova_observed.dec * 2 - tempPos.dec;
 
    // remove the nutation
    ln_get_equ_nut(&tempPos, jd, true);
 
    struct ln_equ_posn libnova_J2000Pos;
    // precess from now to J2000
    ln_get_equ_prec2(&tempPos, jd, JD2000, &libnova_J2000Pos);
 
    J2000pos->rightascension = libnova_J2000Pos.ra / 15.0;
    J2000pos->declination = libnova_J2000Pos.dec;
 
 
}
 
void J2000toObserved(IEquatorialCoordinates *J2000pos, double jd, IEquatorialCoordinates *observed)
{
    ln_equ_posn tempPosn;
    struct ln_equ_posn libnova_J2000Pos = {J2000pos->rightascension * 15.0, J2000pos->declination };
 
    // apply precession from J2000 to jd
    ln_get_equ_prec2(&libnova_J2000Pos, JD2000, jd, &tempPosn);
 
    // apply nutation
    ln_get_equ_nut(&tempPosn, jd, false);
 
    struct ln_equ_posn libnova_observed;
    // apply aberration
    ln_get_equ_aber(&tempPosn, jd, &libnova_observed);
 
    observed->rightascension = libnova_observed.ra / 15.0;
    observed->declination = libnova_observed.dec;
}
 
void ln_get_equ_nut(ln_equ_posn *posn, double jd, bool reverse)
{
    // code lifted from libnova ln_get_equ_nut
    // with the option to add or remove nutation
    struct ln_nutation nut;
    ln_get_nutation (jd, &nut);
 
    double mean_ra, mean_dec, delta_ra, delta_dec;
 
    mean_ra = DEG_TO_RAD(posn->ra);
    mean_dec = DEG_TO_RAD(posn->dec);
 
    // Equ 22.1
 
    double nut_ecliptic = DEG_TO_RAD(nut.ecliptic + nut.obliquity);
    double sin_ecliptic = sin(nut_ecliptic);
 
    double sin_ra = sin(mean_ra);
    double cos_ra = cos(mean_ra);
 
    double tan_dec = tan(mean_dec);
 
    delta_ra = (cos (nut_ecliptic) + sin_ecliptic * sin_ra * tan_dec) * nut.longitude - cos_ra * tan_dec * nut.obliquity;
    delta_dec = (sin_ecliptic * cos_ra) * nut.longitude + sin_ra * nut.obliquity;
 
    // the sign changed to remove nutation
    if (reverse)
    {
        delta_ra = -delta_ra;
        delta_dec = -delta_dec;
    }
    posn->ra += delta_ra;
    posn->dec += delta_dec;
}
 
void EquatorialToHorizontal(IEquatorialCoordinates *object, IGeographicCoordinates *observer, double JD,
                            IHorizontalCoordinates *position)
{
    // Convert from INDI standard location to libnova standard location
    struct ln_lnlat_posn libnova_location = {observer->longitude > 180 ? observer->longitude - 360 : observer->longitude, observer->latitude};
    // RA Hours --> Degrees
    struct ln_equ_posn libnova_object = {object->rightascension * 15.0, object->declination};
    struct ln_hrz_posn horizontalPos;
    // Use apparent sidereal time to match ln_get_equ_from_hrz (which already uses apparent).
    // ln_get_hrz_from_equ uses mean sidereal time — a libnova bug that introduces a ~17" RA
    // offset (equation of the equinoxes) in the round-trip.
    double apparentSidereal = ln_get_apparent_sidereal_time(JD);
    ln_get_hrz_from_equ_sidereal_time(&libnova_object, &libnova_location, apparentSidereal, &horizontalPos);
    position->azimuth = range360(180 + horizontalPos.az);
    position->altitude = horizontalPos.alt;
}
 
void HorizontalToEquatorial(IHorizontalCoordinates *object, IGeographicCoordinates *observer, double JD,
                            IEquatorialCoordinates *position)
{
    // Convert from INDI standard location to libnova standard location
    struct ln_lnlat_posn libnova_location = {observer->longitude > 180 ? observer->longitude - 360 : observer->longitude, observer->latitude};
    // Convert from INDI standard location to libnova standard location
    struct ln_hrz_posn libnova_object = {range360(object->azimuth + 180), object->altitude};
    struct ln_equ_posn equatorialPos;
    // ln_get_equ_from_hrz uses apparent sidereal time, matching EquatorialToHorizontal above.
    ln_get_equ_from_hrz(&libnova_object, &libnova_location, JD, &equatorialPos);
    // Degrees --> Hours
    position->rightascension = equatorialPos.ra / 15.0;
    position->declination = equatorialPos.dec;
 
}
 
 
}