diff --git a/Auxiliary_utils/libastrotools/CMakeLists.txt b/Auxiliary_utils/libastrotools/CMakeLists.txt new file mode 100644 index 0000000..3ced0b2 --- /dev/null +++ b/Auxiliary_utils/libastrotools/CMakeLists.txt @@ -0,0 +1,96 @@ +cmake_minimum_required(VERSION 3.30) +set(PROJ astrotools) +set(MINOR_VERSION "1") +set(MID_VERSION "0") +set(MAJOR_VERSION "0") +set(VERSION "${MAJOR_VERSION}.${MID_VERSION}.${MINOR_VERSION}") + +project(${PROJ} VERSION ${VERSION} LANGUAGES C) + +# default flags +set(CMAKE_C_FLAGS "${CFLAGS} -O2") +set(CMAKE_C_FLAGS_RELEASE "${CMAKE_C_FLAGS}") +set(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS} -Wextra -Wall -Werror -W") +set(CMAKE_COLOR_MAKEFILE ON) + +option(DEBUG "Compile in debug mode" OFF) +option(EXAMPLES "Compile also some examples" ON) + +# cmake -DDEBUG=on -> debugging +if(NOT CMAKE_BUILD_TYPE STREQUAL "Debug") + if(DEBUG) + set(CMAKE_BUILD_TYPE "Debug") + else() + set(CMAKE_BUILD_TYPE "Release") + endif() +endif() + +if(CMAKE_BUILD_TYPE STREQUAL "Debug") + add_definitions(-DEBUG) + set(CMAKE_VERBOSE_MAKEFILE true) + if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) + message("install to ${CMAKE_CURRENT_SOURCE_DIR}/install ") + set(CMAKE_INSTALL_PREFIX ${CMAKE_CURRENT_SOURCE_DIR}/install) + endif() + set(CMAKE_C_FLAGS ${CMAKE_C_FLAGS_DEBUG}) +else() + set(CMAKE_C_FLAGS ${CMAKE_C_FLAGS_RELEASE}) +endif() + +message("Build type: ${CMAKE_BUILD_TYPE}, cflags: ${CMAKE_C_FLAGS}") + +aux_source_directory(${CMAKE_CURRENT_SOURCE_DIR} SOURCES) + +###### pkgconfig ###### +# pkg-config modules (for pkg-check-modules) +#set(MODULES cfitsio fftw3) +# find packages: +#find_package(PkgConfig REQUIRED) +#pkg_check_modules(${PROJ} REQUIRED ${MODULES}) + +# external modules like OpenMP: +include(FindOpenMP) +if(OPENMP_FOUND) + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}") + add_definitions(-DOMP_FOUND) +endif() +###### additional flags ###### +#list(APPEND ${PROJ}_LIBRARIES "-lfftw3_threads") + +# exe file +#add_executable(${PROJ} ${SOURCES}) +# library +add_library(${PROJ} SHARED ${SOURCES}) +# library header files +set(LIBHEADER "astrotools.h") +# -I +include_directories(${${PROJ}_INCLUDE_DIRS}) +# -L +link_directories(${${PROJ}_LIBRARY_DIRS}) +# -D +add_definitions( + -DPACKAGE_VERSION=\"${VERSION}\" -DMINOR_VERSION=\"${MINOR_VERSION}\" + -DMID_VERSION=\"${MID_VERSION}\" -DMAJOR_VERSION=\"${MAJOR_VESION}\" +) + +# -l +target_link_libraries(${PROJ} ${${PROJ}_LIBRARIES}) + +set(PCFILE "${CMAKE_BINARY_DIR}/${PROJ}.pc") +configure_file("${PROJ}.pc.in" ${PCFILE} @ONLY) + +set_target_properties(${PROJ} PROPERTIES VERSION ${VERSION}) +set_target_properties(${PROJ} PROPERTIES PUBLIC_HEADER ${LIBHEADER}) + +# Installation of the program +include(GNUInstallDirs) +#install(TARGETS ${PROJ} DESTINATION "bin") +install(TARGETS ${PROJ} LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) +install(FILES ${PCFILE} DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/pkgconfig) + +# EXAMPLES +if(EXAMPLES) + add_subdirectory(examples) +endif() diff --git a/Auxiliary_utils/libastrotools/Changelog b/Auxiliary_utils/libastrotools/Changelog new file mode 100644 index 0000000..158da70 --- /dev/null +++ b/Auxiliary_utils/libastrotools/Changelog @@ -0,0 +1,11 @@ +2025/01/27: Ver. 0.0.1 +Base operations: +- set/get place, weather, DUT1 and polar move +- dynamic string and conversion of time and angle to human-readable string +- get julian date with all need components (including MJD) by timeval, UNIX time and full julian date; get local sidereal time +- convert RADEC <> AZ +- get ha/ra by ra/ha and LST +- get observed place for given location/weather/polar data of place or star (with pm etc) +- get catalog mean position of coordinates for given epoch and observed place + + diff --git a/Auxiliary_utils/libastrotools/Readme.md b/Auxiliary_utils/libastrotools/Readme.md new file mode 100644 index 0000000..a0d4d51 --- /dev/null +++ b/Auxiliary_utils/libastrotools/Readme.md @@ -0,0 +1,4 @@ +Astrotools +========== + +Simplest wrapper over ERFA (or/and other libraries). diff --git a/Auxiliary_utils/libastrotools/astrotools.cflags b/Auxiliary_utils/libastrotools/astrotools.cflags new file mode 100644 index 0000000..85d51b3 --- /dev/null +++ b/Auxiliary_utils/libastrotools/astrotools.cflags @@ -0,0 +1 @@ +-std=c17 diff --git a/Auxiliary_utils/libastrotools/astrotools.config b/Auxiliary_utils/libastrotools/astrotools.config new file mode 100644 index 0000000..33039ab --- /dev/null +++ b/Auxiliary_utils/libastrotools/astrotools.config @@ -0,0 +1,6 @@ +// Add predefined macros for your project here. For example: +// #define THE_ANSWER 42 +#define EBUG +#define _POSIX_C_SOURCE +#define PACKAGE_VERSION "0.0.1" + diff --git a/Auxiliary_utils/libastrotools/astrotools.creator b/Auxiliary_utils/libastrotools/astrotools.creator new file mode 100644 index 0000000..e94cbbd --- /dev/null +++ b/Auxiliary_utils/libastrotools/astrotools.creator @@ -0,0 +1 @@ +[General] diff --git a/Auxiliary_utils/libastrotools/astrotools.creator.user b/Auxiliary_utils/libastrotools/astrotools.creator.user new file mode 100644 index 0000000..002b8fd --- /dev/null +++ b/Auxiliary_utils/libastrotools/astrotools.creator.user @@ -0,0 +1,184 @@ + + + + + + EnvironmentId + {cf63021e-ef53-49b0-b03b-2f2570cdf3b6} + + + ProjectExplorer.Project.ActiveTarget + 0 + + + ProjectExplorer.Project.EditorSettings + + true + false + true + + Cpp + + CppGlobal + + + + QmlJS + + QmlJSGlobal + + + 2 + KOI8-R + false + 4 + false + 0 + 80 + true + true + 1 + 0 + false + false + false + 1 + true + true + 0 + 8 + true + false + 1 + true + true + true + *.md, *.MD, Makefile + true + true + true + + + + ProjectExplorer.Project.PluginSettings + + + true + false + true + true + true + true + + false + + + 0 + true + + true + true + Builtin.DefaultTidyAndClazy + 4 + true + + + + true + + + + + ProjectExplorer.Project.Target.0 + + Desktop + Desktop + Desktop + {91347f2c-5221-46a7-80b1-0a054ca02f79} + 0 + 0 + 0 + + /home/eddy/Docs/SAO/10micron/C-sources/erfa_functions + + + + all + + true + GenericProjectManager.GenericMakeStep + + 1 + Сборка + Сборка + ProjectExplorer.BuildSteps.Build + + + + + clean + + true + GenericProjectManager.GenericMakeStep + + 1 + Очистка + Очистка + ProjectExplorer.BuildSteps.Clean + + 2 + false + + false + + По умолчанию + GenericProjectManager.GenericBuildConfiguration + + 1 + + + 0 + Развёртывание + Развёртывание + ProjectExplorer.BuildSteps.Deploy + + 1 + + false + ProjectExplorer.DefaultDeployConfiguration + + 1 + + true + true + 0 + true + + + 2 + + false + -e cpu-cycles --call-graph dwarf,4096 -F 250 + + ProjectExplorer.CustomExecutableRunConfiguration + + false + true + true + + 1 + + + + ProjectExplorer.Project.TargetCount + 1 + + + ProjectExplorer.Project.Updater.FileVersion + 22 + + + Version + 22 + + diff --git a/Auxiliary_utils/libastrotools/astrotools.creator.user.7bd84e3 b/Auxiliary_utils/libastrotools/astrotools.creator.user.7bd84e3 new file mode 100644 index 0000000..e594907 --- /dev/null +++ b/Auxiliary_utils/libastrotools/astrotools.creator.user.7bd84e3 @@ -0,0 +1,165 @@ + + + + + + EnvironmentId + {7bd84e39-ca37-46d3-be9d-99ebea85bc0d} + + + ProjectExplorer.Project.ActiveTarget + 0 + + + ProjectExplorer.Project.EditorSettings + + true + false + true + + Cpp + + CppGlobal + + + + QmlJS + + QmlJSGlobal + + + 2 + KOI8-R + false + 4 + false + 80 + true + true + 1 + true + false + 0 + true + true + 0 + 8 + true + 1 + true + false + true + false + + + + ProjectExplorer.Project.PluginSettings + + + true + + + + ProjectExplorer.Project.Target.0 + + Desktop + Desktop + {65a14f9e-e008-4c1b-89df-4eaa4774b6e3} + 0 + 0 + 0 + + /Big/Data/00__Small_tel/C-sources/erfa + + + + all + + false + + + false + true + Сборка + + GenericProjectManager.GenericMakeStep + + 1 + Сборка + + ProjectExplorer.BuildSteps.Build + + + + + clean + + false + + + false + true + Сборка + + GenericProjectManager.GenericMakeStep + + 1 + Очистка + + ProjectExplorer.BuildSteps.Clean + + 2 + false + + По умолчанию + По умолчанию + GenericProjectManager.GenericBuildConfiguration + + 1 + + + 0 + Развёртывание + + ProjectExplorer.BuildSteps.Deploy + + 1 + Конфигурация развёртывания + + ProjectExplorer.DefaultDeployConfiguration + + 1 + + + 2 + + + Особая программа + + ProjectExplorer.CustomExecutableRunConfiguration + + 3768 + false + true + false + false + true + + + + 1 + + + + ProjectExplorer.Project.TargetCount + 1 + + + ProjectExplorer.Project.Updater.FileVersion + 22 + + + Version + 22 + + diff --git a/Auxiliary_utils/libastrotools/astrotools.cxxflags b/Auxiliary_utils/libastrotools/astrotools.cxxflags new file mode 100644 index 0000000..6435dfc --- /dev/null +++ b/Auxiliary_utils/libastrotools/astrotools.cxxflags @@ -0,0 +1 @@ +-std=c++17 \ No newline at end of file diff --git a/Auxiliary_utils/libastrotools/astrotools.files b/Auxiliary_utils/libastrotools/astrotools.files new file mode 100644 index 0000000..4e072c9 --- /dev/null +++ b/Auxiliary_utils/libastrotools/astrotools.files @@ -0,0 +1,7 @@ +CMakeLists.txt +astrotools.h +erfa.c +erfa.c +examples/CMakeLists.txt +examples/coordstest.c +examples/transform.c diff --git a/Auxiliary_utils/libastrotools/astrotools.h b/Auxiliary_utils/libastrotools/astrotools.h new file mode 100644 index 0000000..87fffe7 --- /dev/null +++ b/Auxiliary_utils/libastrotools/astrotools.h @@ -0,0 +1,127 @@ +/* + * This file is part of the astrotools project. + * Copyright 2025 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include +#include +#include +#include +#include + +#ifndef FALSE +#define FALSE (0) +#endif + +#ifndef TRUE +#define TRUE (1) +#endif + +#ifndef ERFA_DR2S +#define ERFA_DR2S 1.3750987083139757010431557155385240879777313391975e4 +#endif + +typedef struct{ + double utc1; double utc2; // UTC JD, commonly used MJD = utc1+utc2-2400000.5 + double MJD; + double tai1; double tai2; // TAI JD + double tt1; double tt2; // TT JD +} at_MJD_t; + +// equatorial coordinates & equation of origins (all in radians) +typedef struct{ + double ha; // hour angle + double ra; // right ascencion + double dec; // declination + double eo; // equation of origins +} at_equat_t; + +// horizontal coordinates (all in radians) +typedef struct{ + double az; // azimuth, 0 @ south, positive clockwise + double zd; // zenith distance +} at_horiz_t; + +// observational place coordinates and altitude; all angle coordinates are in radians! +typedef struct{ + double longitude; // longitude + double latitude; // latitude + double altitude; // altitude, m +} at_place_t; + +// place weather data +typedef struct{ + double relhum; // rel. humidity, 0..1 + double phpa; // atm. pressure (hectopascales) + double tdegc; // temperature, degrC +} at_weather_t; + +// DUT/polar almanach data +typedef struct{ + double DUT1; // UT1-UTC, sec + double px; // polar coordinates, arcsec + double py; +} at_dut_t; + +typedef struct{ + char *str; + size_t maxlen; + size_t strlen; +} at_string_t; + +typedef struct{ + double pmRA; // proper motion (radians/year) + double pmDec; + double parallax; // (arcsec) + double radvel; // radial velocity (km/s, positive if receding) +} at_star_t; + +int at_getPlace(at_place_t *p); +int at_setPlace(at_place_t *p); +int at_getWeath(at_weather_t *w); +int at_setWeath(at_weather_t *w); +int at_getDUT(at_dut_t *a); +int at_setDUT(at_dut_t *a); +void at_setEffWvl(double w); + +// strings +int at_chkstring(at_string_t *s, size_t minlen); +at_string_t *at_newstring(size_t maxlen); +void at_delstring(at_string_t **s); +const char *at_libversion(); + +// human readable +int at_radtoHdeg(double r, at_string_t *s); +int at_radtoHtime(double r, at_string_t *s); + +// time-date +int at_get_MJDt(struct timeval *tval, at_MJD_t *mjd); +int at_get_MJDu(double UNIXtime, at_MJD_t *mjd); +int at_get_MJDj(double JD, at_MJD_t *mjd); +double at_get_LST(at_MJD_t *mjd); + +// coordinates +void at_hor2eq(at_horiz_t *h, at_equat_t *pc, double LST); +void at_eq2hor(at_equat_t *pc, at_horiz_t *h); +double at_getRA(at_equat_t *pc, double LST); +double at_getHA(at_equat_t *pc, double LST); +int at_get_ObsPlace(at_MJD_t *MJD, at_equat_t *p2000, at_equat_t *pnow, at_horiz_t *hnow); +int at_get_ObsPlaceStar(at_MJD_t *MJD, at_equat_t *p2000, at_star_t *star, at_equat_t *pnow, at_horiz_t *hnow);; +int at_get_mean(at_MJD_t *MJD, at_equat_t *pnow, at_equat_t *p2000); +int at_obs2catP(at_MJD_t *MJD, at_equat_t *pnow, at_equat_t *p2000); +int at_obs2catA(at_MJD_t *MJD, at_horiz_t *hnow, at_equat_t *p2000); diff --git a/Auxiliary_utils/libastrotools/astrotools.includes b/Auxiliary_utils/libastrotools/astrotools.includes new file mode 100644 index 0000000..263898c --- /dev/null +++ b/Auxiliary_utils/libastrotools/astrotools.includes @@ -0,0 +1,2 @@ +. +.. diff --git a/Auxiliary_utils/libastrotools/astrotools.pc.in b/Auxiliary_utils/libastrotools/astrotools.pc.in new file mode 100644 index 0000000..6f40abc --- /dev/null +++ b/Auxiliary_utils/libastrotools/astrotools.pc.in @@ -0,0 +1,10 @@ +prefix=@CMAKE_INSTALL_PREFIX@ +exec_prefix=${prefix} +libdir=${exec_prefix}/lib +includedir=${prefix}/include + +Name: @PROJ@ +Description: simplest astrotools @ ERFA and NOVA +Version: @VERSION@ +Libs: -L${libdir} -l@PROJ@ +Cflags: -I${includedir} diff --git a/Auxiliary_utils/libastrotools/erfa.c b/Auxiliary_utils/libastrotools/erfa.c new file mode 100644 index 0000000..cb09ea6 --- /dev/null +++ b/Auxiliary_utils/libastrotools/erfa.c @@ -0,0 +1,481 @@ +/* + * This file is part of the astrotools project. + * Copyright 2025 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +#include +#include + +#include "astrotools.h" + +#define ERFA(f, ...) do{if(f(__VA_ARGS__)){WARNX("Error in " #f); return FALSE;}}while(0) + +/*********** default parameters (could/should be changed by setters) ***********/ +// default - BTA +static at_place_t s_pldata = { + .longitude = 0.7232763200, + .latitude = 0.7618977414, + .altitude = 2070.0 +}; +static pthread_mutex_t s_pldata_m = PTHREAD_MUTEX_INITIALIZER; +// default weather +static at_weather_t s_weather = { + .relhum = 0.5, + .phpa = 800., + .tdegc = 0. +}; +static pthread_mutex_t s_weather_m = PTHREAD_MUTEX_INITIALIZER; +// default DUT1 and polar motion +static at_dut_t s_dut = { + .DUT1 = -0.01697, + .px = 0., + .py = 0. +}; +static pthread_mutex_t s_dut_m = PTHREAD_MUTEX_INITIALIZER; + +// eff. wavelength for observed place calculation +static double s_effwavelength = 0.5; + +/** + * @brief at_getPlace - get stored place data + * @param p (o) - data + * @return FALSE if !p + */ +int at_getPlace(at_place_t *p){ + if(!p) return FALSE; + pthread_mutex_lock(&s_pldata_m); + *p = s_pldata; + pthread_mutex_unlock(&s_pldata_m); + return TRUE; +} +/** + * @brief at_setPlace - change stored place data + * @param p (i) - new place + * @return FALSE if !p or p is wrong + */ +int at_setPlace(at_place_t *p){ + if(!p) return FALSE; + if(p->longitude > ERFA_DPI || p->longitude < -ERFA_DPI + || p->latitude > ERFA_DPI/2 || p->latitude < -ERFA_DPI/2) return FALSE; + pthread_mutex_lock(&s_pldata_m); + s_pldata = *p; + pthread_mutex_unlock(&s_pldata_m); + return TRUE; +} + +/** + * @brief at_getWeath - get last stored weather data + * @param w (o) - weather + * @return FALSE if !w + */ +int at_getWeath(at_weather_t *w){ + if(!w) return FALSE; + pthread_mutex_lock(&s_weather_m); + *w = s_weather; + pthread_mutex_unlock(&s_weather_m); + return TRUE; +} +/** + * @brief at_setWeath - set new weather data + * @param w (i) - new weather + * @return FALSE if !w or w is wrong + */ +int at_setWeath(at_weather_t *w){ + if(!w) return FALSE; + if(w->phpa < 0 || w->phpa > 2000 || + w->relhum < 0. || w->relhum > 1. || + w->tdegc < -273.15 || w->tdegc > 100.) return FALSE; + pthread_mutex_lock(&s_weather_m); + s_weather = *w; + pthread_mutex_unlock(&s_weather_m); + return TRUE; +} + +/** + * @brief at_getDUT - get last stored DUT1 and polar motion data + * @param a (o) - data + * @return FALSE if !a + */ +int at_getDUT(at_dut_t *a){ + if(!a) return FALSE; + pthread_mutex_lock(&s_dut_m); + *a = s_dut; + pthread_mutex_unlock(&s_dut_m); + return TRUE; +} +/** + * @brief at_setDUT - set new DUT1 and PM + * @param a (i) - data + * @return FALSE if !a or a is wrong + */ +int at_setDUT(at_dut_t *a){ + if(!a) return FALSE; + if(a->DUT1 < -1. || a->DUT1 > 1. || + a->px < -1000. || a->px > 1000. || + a->py < -1000. || a->py > 1000.) return FALSE; + pthread_mutex_lock(&s_dut_m); + s_dut = *a; + pthread_mutex_unlock(&s_dut_m); + return TRUE; +} + +/** + * @brief at_chkstring - check if string is large enough + * @param s - string + * @param minlen - minimal length (or 0 if not need to check) + * @return TRUE if all OK + */ +int at_chkstring(at_string_t *s, size_t minlen){ + if(!s || !s->str) return FALSE; + if(s->maxlen < minlen) s->str = realloc(s->str, minlen+1); + return TRUE; +} +/** + * @brief at_newstring + * @param maxlen + * @return + */ +at_string_t *at_newstring(size_t maxlen){ + at_string_t *s = MALLOC(at_string_t, 1); + s->str = MALLOC(char, maxlen+1); + s->maxlen = maxlen; + return s; +} +void at_delstring(at_string_t **s){ + if(!s || !*s) return; + FREE((*s)->str); + FREE(*s); +} +const char *at_libversion(){ + const char *v = PACKAGE_VERSION; + return v; +} + +/** + * @brief at_radtoHdeg - reformat angle in radians to human readable string + * @param r - angle + * @param s - string + * @return + */ +int at_radtoHdeg(double r, at_string_t *s){ +#define FMTSTRLEN 13 + if(!at_chkstring(s, FMTSTRLEN)) return FALSE; + int i[4]; char pm; + r = eraAnpm(r); + eraA2af(2, r, &pm, i); + snprintf(s->str, FMTSTRLEN+1, "%c%03d:%02d:%02d.%02d", pm, i[0],i[1],i[2],i[3]); + s->strlen = FMTSTRLEN; +#undef FMTSTRLEN + return TRUE; +} +/** + * @brief at_radtoHtime - reformat time from radians into human readable string + * @param r - time / hour angle / etc. + * @param s - string + * @return + */ +int at_radtoHtime(double r, at_string_t *s){ +#define FMTSTRLEN 11 + if(!at_chkstring(s, FMTSTRLEN)) return FALSE; + int i[4]; char pm; + r = eraAnp(r); + eraA2tf(2, r, &pm, i); + snprintf(s->str, FMTSTRLEN+1, "%02d:%02d:%02d.%02d", i[0],i[1],i[2],i[3]); + s->strlen = FMTSTRLEN; +#undef FMTSTRLEN + return TRUE; +} + +/** + * @brief at_get_MJDt - calculate MJD of `struct timeval` from argument + * @param tval (i) - given date (or NULL for current) + * @param MJD (o) - time (or NULL just to check) + * @return TRUE if all OK + */ +int at_get_MJDt(struct timeval *tval, at_MJD_t *mjd){ + struct tm tms; + double tSeconds; + if(!tval){ + struct timeval currentTime; + gettimeofday(¤tTime, NULL); + gmtime_r(¤tTime.tv_sec, &tms); + tSeconds = tms.tm_sec + ((double)currentTime.tv_usec)/1e6; + }else{ + gmtime_r(&tval->tv_sec, &tms); + tSeconds = tms.tm_sec + ((double)tval->tv_usec)/1e6; + } + int y, m, d; + y = 1900 + tms.tm_year; + m = tms.tm_mon + 1; + d = tms.tm_mday; + double utc1, utc2; + /* UTC date. */ + ERFA(eraDtf2d, "UTC", y, m, d, tms.tm_hour, tms.tm_min, tSeconds, &utc1, &utc2); + at_MJD_t MJD; + MJD.MJD = utc1 - ERFA_DJM0 + utc2; + MJD.utc1 = utc1; + MJD.utc2 = utc2; + ERFA(eraUtctai, utc1, utc2, &MJD.tai1, &MJD.tai2); + ERFA(eraTaitt, MJD.tai1, MJD.tai2, &MJD.tt1, &MJD.tt2); + if(mjd) *mjd = MJD; + return TRUE; +} +// calculate by UNIX-time +int at_get_MJDu(double UNIXtime, at_MJD_t *mjd){ + struct timeval t; + t.tv_sec = (time_t) UNIXtime; + t.tv_usec = (suseconds_t) ((UNIXtime-t.tv_sec) * 1e6); + return at_get_MJDt(&t, mjd); +} +// get full date by JulianDate +int at_get_MJDj(double JD, at_MJD_t *mjd){ + int y, m, d, H, M; + double fd, utc1, utc2; + ERFA(eraJd2cal, JD, 0., &y, &m, &d, &fd); + fd *= 24.; + H = (int)fd; + fd = (fd - H)*60.; + M = (int)fd; + fd = (fd - M)*60.; + ERFA(eraDtf2d, "UTC", y, m, d, H, M, fd, &utc1, &utc2); + at_MJD_t MJD; + MJD.MJD = utc1 - ERFA_DJM0 + utc2; + MJD.utc1 = utc1; + MJD.utc2 = utc2; + ERFA(eraUtctai, utc1, utc2, &MJD.tai1, &MJD.tai2); + ERFA(eraTaitt, MJD.tai1, MJD.tai2, &MJD.tt1, &MJD.tt2); + if(mjd) *mjd = MJD; + return TRUE; +} + +/** + * @brief at_get_LST - calculate local siderial time for current dut and place (change them if need with setters) + * @param mjd (i) - date/time for LST (utc1 & tt used) + * @param dUT1 - (UT1-UTC) + * @param longitude - site longitude (radians) + * @param LST (o) - local sidereal time (radians) + * @return lst (radians) or -1 in case of error + */ +double at_get_LST(at_MJD_t *mjd){ + double ut11, ut12; + if(eraUtcut1(mjd->utc1, mjd->utc2, s_dut.DUT1, &ut11, &ut12)){ + WARNX("error in eraUtcut1"); + return -1.; + } + double ST = eraGst06a(ut11, ut12, mjd->tt1, mjd->tt2); + ST += s_pldata.longitude; + if(ST > ERFA_D2PI) ST -= ERFA_D2PI; + return ST; +} + +/** + * @brief at_hor2eq - convert horizontal coordinates to polar + * @param h (i) - horizontal coordinates + * @param pc (o) - polar coordinates + * @param LST - local sidereal time + */ +void at_hor2eq(at_horiz_t *h, at_equat_t *pc, double LST){ + if(!h || !pc) return; + eraAe2hd(h->az, ERFA_DPI/2. - h->zd, s_pldata.latitude, &pc->ha, &pc->dec); // A,H -> HA,DEC; phi - site latitude + pc->eo = 0.; + if(LST >= 0. && LST < ERFA_D2PI) at_getRA(pc, LST); +} + +/** + * @brief at_eq2hor - convert polar coordinates to horizontal + * @param pc (i) - polar coordinates (only HA used) + * @param h (o) - horizontal coordinates + * @param sidTime - sidereal time + */ +void at_eq2hor(at_equat_t *pc, at_horiz_t *h){ + if(!h || !pc) return; + double alt; + eraHd2ae(pc->ha, pc->dec, s_pldata.latitude, &h->az, &alt); + h->zd = ERFA_DPI/2. - alt; +} + +/** + * @brief at_getRA - get RA from HA and LST (and change pc->ra) + * @param pc - polar coordinates + * @param LST - local sidereal time + * @return right ascension (radians) + * WARNING: this function doesn't check arguments (e.g. for 0 < LST < 2pi) + */ +double at_getRA(at_equat_t *pc, double LST){ + double ra = eraAnp(LST - pc->ha + pc->eo); + pc->ra = ra; + return ra; +} +/** + * @brief at_getHA - get HA from RA and LST (and change pc->ha) + * @param pc - polar coordinates + * @param LST - local sidereal time + * @return hour angle (radians) + * WARNING: this function doesn't check arguments (e.g. for 0 < LST < 2pi) + */ +double at_getHA(at_equat_t *pc, double LST){ + double ha = eraAnpm(LST - pc->ra + pc->eo); + pc->ha = ha; + return ha; +} + +/** + * @brief at_setEffWvl - set effective wavelength for el_get_ObsPlace + * @param w - wavelength in um + * WITHOUT CHECKING! + */ +void at_setEffWvl(double w){ + s_effwavelength = w; +} + +/** + * @brief at_get_ObsPlaceStar - calculate observed place for given star @MJD + * @param MJD (i) - date + * @param p2000 (i) - catalog coordinates + * @param star (i) - star with proper motion etc. + * @param pnow (o) - observed place in polar coordinates + * @param hnow (o) - observed AZ + * @return FALSE if failed + */ +int at_get_ObsPlaceStar(at_MJD_t *MJD, at_equat_t *p2000, at_star_t *star, at_equat_t *pnow, at_horiz_t *hnow){ + if(!MJD || !p2000 || !star) return -1; + double aob, zob, hob, dob, rob, eo; + /*DBG("p2000->ra=%g, p2000->dec=%g," + "star->pmRA=%g, star->pmDec=%g, star->parallax=%g, star->radvel=%g," + "MJD->utc1=%g, MJD->utc2=%g," + "s_dut.DUT1=%g," + "s_pldata.longitude=%g, s_pldata.latitude=%g, s_pldata.altitude=%g," + "s_dut.px=%g, s_dut.py=%g," + "s_weather.phpa=%g, s_weather.tdegc=%g, s_weather.relhum=%g," + "s_effwavelength=%g", p2000->ra, p2000->dec, + star->pmRA, star->pmDec, star->parallax, star->radvel, + MJD->utc1, MJD->utc2, + s_dut.DUT1, + s_pldata.longitude, s_pldata.latitude, s_pldata.altitude, + s_dut.px, s_dut.py, + s_weather.phpa, s_weather.tdegc, s_weather.relhum, + s_effwavelength);*/ + ERFA(eraAtco13,p2000->ra, p2000->dec, + star->pmRA, star->pmDec, star->parallax, star->radvel, + MJD->utc1, MJD->utc2, + s_dut.DUT1, + s_pldata.longitude, s_pldata.latitude, s_pldata.altitude, + s_dut.px, s_dut.py, + s_weather.phpa, s_weather.tdegc, s_weather.relhum, + s_effwavelength, + &aob, &zob, + &hob, &dob, &rob, &eo); + if(pnow){ + pnow->eo = eo; + pnow->ha = hob; + pnow->ra = rob; + pnow->dec = dob; + } + if(hnow){ + hnow->az = aob; + hnow->zd = zob; + } + return TRUE; +} + +/** + * @brief at_get_ObsPlace - calculate observed place (without PM etc) for given date + * @param tval (i) - time + * @param p2000 (i) - polar coordinates for J2000 (only ra/dec used), ICRS (catalog) + * @param pnow (o) - polar coordinates for given epoch (or NULL) + * @param hnow (o) - horizontal coordinates for given epoch (or NULL) + * @return FALSE if failed + */ +int at_get_ObsPlace(at_MJD_t *MJD, at_equat_t *p2000, at_equat_t *pnow, at_horiz_t *hnow){ + at_star_t star = {0}; + return at_get_ObsPlaceStar(MJD, p2000, &star, pnow, hnow); +} + +/** + * @brief + * @param + */ +/** + * @brief at_get_mean - convert apparent coordinates (for MJD) to mean (JD2000) + * @param MJD (i) - given date + * @param pnow (i) - polar coordinates for given date + * @param p2000 (o) - catalog coordinates + * @return TRUE if OK + */ +int at_get_mean(at_MJD_t *MJD, at_equat_t *pnow, at_equat_t *p2000){ + if(!MJD || !pnow) return FALSE; + double ra, dec, eo, ri; + eraAtic13(pnow->ra, pnow->dec, MJD->tt1, MJD->tt2, &ri, &dec, &eo); + ra = eraAnp(ri + eo); + if(p2000){ + p2000->ra = ra; + p2000->dec = dec; + p2000->eo = eo; + at_getHA(p2000, 0.); + } + return TRUE; +} + +/** + * @brief at_obs2catE - convert observed place in Ra-Dec coordinates into catalog J2000 + * @param MJD (i) - epoch for `pnow` + * @param pnow (i) - coordinates + * @param p2000 (o) - catalog coordinates + * @return TRUE if OK + */ +int at_obs2catP(at_MJD_t *MJD, at_equat_t *pnow, at_equat_t *p2000){ + double ra, dec; + ERFA(eraAtoc13, "R", pnow->ra, pnow->dec, MJD->utc1, MJD->utc2, + s_dut.DUT1, s_pldata.longitude, s_pldata.latitude, s_pldata.altitude, + s_dut.px, s_dut.py, + s_weather.phpa, s_weather.tdegc, s_weather.relhum, + s_effwavelength, + &ra, &dec); + if(p2000){ + p2000->ra = ra; + p2000->dec = dec; + p2000->eo = 0.; + at_getHA(p2000, 0.); + } + return TRUE; +} + +/** + * @brief at_obs2catAconvert - observed place in A-Z coordinates into catalog J2000 + * @param MJD (i) - epoch for `hnow` + * @param hnow (i) - coordinates + * @param p2000 (o) - catalog coordinates + * @return TRUE if OK + */ +int at_obs2catA(at_MJD_t *MJD, at_horiz_t *hnow, at_equat_t *p2000){ + double ra, dec; + ERFA(eraAtoc13, "A", hnow->az, hnow->zd, MJD->utc1, MJD->utc2, + s_dut.DUT1, s_pldata.longitude, s_pldata.latitude, s_pldata.altitude, + s_dut.px, s_dut.py, + s_weather.phpa, s_weather.tdegc, s_weather.relhum, + s_effwavelength, + &ra, &dec); + if(p2000){ + p2000->ra = ra; + p2000->dec = dec; + p2000->eo = 0.; + at_getHA(p2000, 0.); + } + return TRUE; +} diff --git a/Auxiliary_utils/libastrotools/examples/CMakeLists.txt b/Auxiliary_utils/libastrotools/examples/CMakeLists.txt new file mode 100644 index 0000000..11343c0 --- /dev/null +++ b/Auxiliary_utils/libastrotools/examples/CMakeLists.txt @@ -0,0 +1,9 @@ +project(examples) + +# common includes & library +include_directories(../) +link_libraries(astrotools usefull_macros erfa) + +# exe list +add_executable(coordstest coordstest.c) +add_executable(transform transform.c) diff --git a/Auxiliary_utils/libastrotools/examples/coordstest.c b/Auxiliary_utils/libastrotools/examples/coordstest.c new file mode 100644 index 0000000..ebc47c3 --- /dev/null +++ b/Auxiliary_utils/libastrotools/examples/coordstest.c @@ -0,0 +1,84 @@ +/* + * This file is part of the ERFA project. + * Copyright 2025 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +#include "astrotools.h" + +static at_string_t *S = NULL; + +static const char *radtohrs(double r){ + if(!S) S = at_newstring(256); + if(!at_radtoHtime(r, S)) ERRX("at_radtoHtime"); + return S->str; +} +static const char *radtodeg(double r){ + if(!S) S = at_newstring(256); + if(!at_radtoHdeg(r, S)) ERRX("at_radtoHdeg"); + return S->str; +} + + +int main(){ + initial_setup(); + at_MJD_t mjd; + if(!at_get_MJDu(time(NULL), &mjd)) ERRX("at_get_MJDu"); + printf("MJD=%g; TAI=%g/%g, TT=%g/%g, UTC=%g/%g\n", mjd.MJD, mjd.tai1, mjd.tai2, mjd.tt1, mjd.tt2, mjd.utc1, mjd.utc2); + double ST = at_get_LST(&mjd); + if(ST < 0.) ERRX("at_get_LST"); + printf("ST = %s\n", radtohrs(ST)); + at_horiz_t htest = {.az = ERFA_DD2R*91., .zd = ERFA_DPI/4.}; // Z=45degr, A=1degr from south to west + printf("hzd=%g\n", htest.zd); + at_equat_t ptest; + at_hor2eq(&htest, &ptest, ST); + printf("A=%s, ", radtodeg(htest.az)); + printf("Z=%s; ", radtodeg(htest.zd)); + printf("HOR->EQ: HA=%s, ", radtohrs(ptest.ha)); + printf("RA=%s, ", radtohrs(ptest.ra)); + printf("DEC=%s\n", radtodeg(ptest.dec)); + at_horiz_t h2; + at_eq2hor(&ptest, &h2); + printf("Back conversion EQ->HOR: A=%s, ", radtodeg(h2.az)); + printf("Z=%s\n", radtodeg(h2.zd)); + at_equat_t pnow; + if(!at_get_ObsPlace(&mjd, &ptest, &pnow, &h2)) ERRX("at_get_ObsPlace"); + printf("\nApparent place, RA=%s, ", radtohrs(pnow.ra-pnow.eo)); + printf("HA=%s, ", radtohrs(pnow.ha)); + printf("ST-RA=%s, ", radtohrs(ST-pnow.ra+pnow.eo)); + printf("DEC=%s; ", radtodeg(pnow.dec)); + printf("A=%s, ", radtodeg(h2.az)); + printf("Z=%s\n", radtodeg(h2.zd)); + at_hor2eq(&h2, &ptest, ST); + printf("\tHOR->EQ: RA=%s, ", radtohrs(ptest.ra-ptest.eo)); + printf("HA=%s, ", radtohrs(ptest.ha)); + printf("ST-RA=%s, ", radtohrs(ST-ptest.ra+ptest.eo)); + printf("DEC=%s\n", radtodeg(ptest.dec)); + at_eq2hor(&pnow, &h2); + printf("\tEQ->HOR: A=%s, ", radtodeg(h2.az)); + printf("Z=%s\n", radtodeg(h2.zd)); + if(!at_get_mean(&mjd, &pnow, &ptest)) ERRX("at_get_mean"); + printf("\nBack conversion pnow to mean place, "); + printf("RA=%s, ", radtohrs(ptest.ra)); + printf("Dec=%s\n", radtodeg(ptest.dec)); + if(!at_obs2catP(&mjd, &pnow, &ptest)) ERRX("at_obs2catP"); + printf("And back to J2000 by observed pnow: "); + printf("RA=%s, ", radtohrs(ptest.ra)); + printf("Dec=%s\n", radtodeg(ptest.dec)); + return 0; +} diff --git a/Auxiliary_utils/libastrotools/examples/transform.c b/Auxiliary_utils/libastrotools/examples/transform.c new file mode 100644 index 0000000..47f4a3b --- /dev/null +++ b/Auxiliary_utils/libastrotools/examples/transform.c @@ -0,0 +1,135 @@ +/* + * This file is part of the astrotools project. + * Copyright 2025 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include "astrotools.h" + +// convert coordinates from given epoch to J2000 + +typedef struct{ + int help; + int obsplace; + double ra; + double dec; + double JD; + double unixtime; + double longitude; + double latitude; + double altitude; + double relhum; + double phpa; + double tdegc; + double DUT1; + double px; + double py; +} parameters; + +static parameters G = { + .ra = -100., + .dec = -100., + .JD = -1., + .unixtime = -1., + .longitude = -1000., + .latitude = -1000., + .altitude = -1000., + .relhum = -1., + .phpa = -1., + .tdegc = -300., + .DUT1 = -100., + .px = -10000., + .py = -10000. +}; + +static myoption cmdlnopts[] = { + {"help", NO_ARGS, NULL, 'h', arg_int, APTR(&G.help), "show this help"}, + {"obsplace", NO_ARGS, NULL, 'O', arg_int, APTR(&G.obsplace), "input RA/Dec is observed place"}, + {"JD", NEED_ARG, NULL, 'J', arg_double, APTR(&G.JD), "Julian date"}, + {"unixtime", NEED_ARG, NULL, 'u', arg_double, APTR(&G.unixtime), "UNIX-time (seconds)"}, + {"ra", NEED_ARG, NULL, 'R', arg_double, APTR(&G.ra), "Right Ascention for given date (degrees)"}, + {"dec", NEED_ARG, NULL, 'D', arg_double, APTR(&G.dec), "Declination for given date (degrees)"}, + {"longitude", NEED_ARG, NULL, 'o', arg_double, APTR(&G.longitude), "site longitude (degr)"}, + {"latitude", NEED_ARG, NULL, 'l', arg_double, APTR(&G.latitude), "site latitude (degr)"}, + {"altitude", NEED_ARG, NULL, 'a', arg_double, APTR(&G.altitude), "site altitude (meters)"}, + {"relhum", NEED_ARG, NULL, 'H', arg_double, APTR(&G.relhum), "relative humidity (0..1)"}, + {"phpa", NEED_ARG, NULL, 'P', arg_double, APTR(&G.phpa), "atm. pressure (hPa)"}, + {"tdegc", NEED_ARG, NULL, 'T', arg_double, APTR(&G.tdegc), "ambient temperature (degC)"}, + {"DUT1", NEED_ARG, NULL, 'd', arg_double, APTR(&G.DUT1), "DUT1 (seconds)"}, + {"px", NEED_ARG, NULL, 'x', arg_double, APTR(&G.px), "polar X (m)"}, + {"py", NEED_ARG, NULL, 'y', arg_double, APTR(&G.py), "polar Y (m)"}, + end_option +}; + + +int main(int argc, char **argv){ + initial_setup(); + parseargs(&argc, &argv, cmdlnopts); + if(G.help) showhelp(-1, cmdlnopts); + at_MJD_t MJD; + G.ra *= ERFA_DD2R; + G.dec *= ERFA_DD2R; + if(G.ra < 0. || G.ra > ERFA_D2PI || G.dec < -ERFA_DPI/2. || G.dec > ERFA_DPI/2.) + ERRX("Need RA (0..360 degr) and Dec (-90..90 degr)"); + if(G.JD < 0. && G.unixtime < 0.) + ERRX("Need JD or unixtime"); + at_place_t place; + at_getPlace(&place); + at_weather_t weather; + at_getWeath(&weather); + at_dut_t dut; + at_getDUT(&dut); + G.longitude *= ERFA_DD2R; + G.latitude *= ERFA_DD2R; + if(G.longitude >= -ERFA_DPI && G.longitude <= ERFA_DPI) place.longitude = G.longitude; + if(G.latitude >= -ERFA_DPI/2 && G.latitude <= ERFA_DPI/2) place.latitude = G.latitude; + if(G.altitude > -100. && G.altitude < 12000.) place.altitude = G.altitude; + if(G.relhum >= 0. && G.relhum <= 1.) weather.relhum = G.relhum; + if(G.phpa >= 0. && G.phpa <= 1300.) weather.phpa = G.phpa; + if(G.tdegc > -273.15 && G.tdegc < 100.) weather.tdegc = G.tdegc; + if(G.DUT1 > -1. && G.DUT1 < 1.) dut.DUT1 = G.DUT1; + if(G.px > -1000. && G.px < 1000.) dut.px = G.px; + if(G.py > -1000. && G.py < 1000.) dut.py = G.py; + at_setPlace(&place); + DBG("Place: long=%g, lat=%g, alt=%g", + place.longitude*ERFA_DR2D, place.latitude*ERFA_DR2D, place.altitude); + at_setWeath(&weather); + DBG("Weather: P=%g hPa, rho=%g%%, T=%g degrC", + weather.phpa, weather.relhum*100., weather.tdegc); + at_setDUT(&dut); + DBG("DUT1=%g, px=%g, py=%g", dut.DUT1, dut.px, dut.py); + if(G.JD > 0. && !at_get_MJDj(G.JD, &MJD)) ERRX("Bad julian date"); + if(G.unixtime > 0. && !at_get_MJDu(G.unixtime, &MJD)) ERRX("Bad UNIX time"); + DBG("Julian: MJD=%g, TT=%.2f+%g", MJD.MJD, MJD.tt1, MJD.tt2); + at_equat_t p2000, pnow = {.ra = G.ra, .dec = G.dec, .eo = 0.}; + at_getHA(&pnow, at_get_LST(&MJD)); + DBG("in: ra=%.10f, dec=%.10f, ha=%.10f", pnow.ra*ERFA_DR2D, pnow.dec*ERFA_DR2D, pnow.ha*ERFA_DR2D); + if(G.obsplace){ // ra/dec is observed place + if(!at_obs2catP(&MJD, &pnow, &p2000)) ERRX("at_obs2catP"); + DBG("Observed"); + }else{ // ra/dec is catalog for given epoch + if(!at_get_mean(&MJD, &pnow, &p2000)) ERRX("at_get_mean"); + DBG("Catalog"); + } + at_string_t *s = at_newstring(128); + at_radtoHtime(p2000.ra, s); + printf("RA(h:m:s)=%s, ", s->str); + at_radtoHdeg(p2000.dec, s); + printf("Dec(d:m:s)=%s\n", s->str); + printf("RA(degr)=%g, Dec(degr)=%g\n", p2000.ra*ERFA_DR2D, p2000.dec*ERFA_DR2D); + return 0; +}