diff --git a/Auxiliary_utils/PCS_create/CMakeLists.txt b/Auxiliary_utils/PCS_create/CMakeLists.txt new file mode 100644 index 0000000..05c3d57 --- /dev/null +++ b/Auxiliary_utils/PCS_create/CMakeLists.txt @@ -0,0 +1,73 @@ +cmake_minimum_required(VERSION 3.0) +set(PROJ PCS_create) +set(MINOR_VERSION "1") +set(MID_VERSION "0") +set(MAJOR_VERSION "0") +set(VERSION "${MAJOR_VERSION}.${MID_VERSION}.${MINOR_VERSION}") + +project(${PROJ} VERSION ${PROJ_VERSION} LANGUAGES C) +#enable_language(C) + +message("VER: ${VERSION}") + +# default flags +set(CMAKE_C_FLAGS_RELEASE "") +set(CMAKE_C_FLAGS_DEBUG "") +set(CMAKE_C_FLAGS "-O2 -std=gnu99") + +set(CMAKE_COLOR_MAKEFILE ON) + +# here is one of two variants: all .c in directory or .c files in list +aux_source_directory(${CMAKE_CURRENT_SOURCE_DIR} SOURCES) + +# cmake -DDEBUG=1 -> debugging +if(DEFINED EBUG) + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wextra -Wall -Werror -W") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wextra -Wall -Werror -W") + set(CMAKE_BUILD_TYPE DEBUG) + set(CMAKE_VERBOSE_MAKEFILE "ON") + add_definitions(-DEBUG) +else() + set(CMAKE_BUILD_TYPE RELEASE) +endif() + +# find cfitsio +set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}) +find_package(CFITSIO REQUIRED) + +###### pkgconfig ###### +# pkg-config modules (for pkg-check-modules) +set(MODULES usefull_macros) + +# find packages: +find_package(PkgConfig REQUIRED) +pkg_check_modules(${PROJ} REQUIRED ${MODULES}) + +###### additional flags ###### +list(APPEND ${PROJ}_LIBRARIES "-lsofa_c") + +# change wrong behaviour with install prefix +if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT AND CMAKE_INSTALL_PREFIX MATCHES "/usr/local") +else() + message("Change default install path to /usr/local") + set(CMAKE_INSTALL_PREFIX "/usr/local") +endif() +message("Install dir prefix: ${CMAKE_INSTALL_PREFIX}") + +# exe file +add_executable(${PROJ} ${SOURCES}) +# -l +target_link_libraries(${PROJ} ${${PROJ}_LIBRARIES} ${CFITSIO_LIBRARIES}) +# -I +include_directories(${${PROJ}_INCLUDE_DIRS} ${CFITSIO_INCLUDE_DIRS}) +# -L +link_directories(${${PROJ}_LIBRARY_DIRS} ${CFITSIO_LIBRARY_DIRS}) +# -D +add_definitions(${CFLAGS} -DLOCALEDIR=\"${LOCALEDIR}\" + -DPACKAGE_VERSION=\"${VERSION}\" -DGETTEXT_PACKAGE=\"${PROJ}\" + -DMINOR_VERSION=\"${MINOR_VERSION}\" -DMID_VERSION=\"${MID_VERSION}\" + -DMAJOR_VERSION=\"${MAJOR_VERSION}\") + + +# Installation of the program +INSTALL(TARGETS ${PROJ} DESTINATION "bin") diff --git a/Auxiliary_utils/PCS_create/FindCFITSIO.cmake b/Auxiliary_utils/PCS_create/FindCFITSIO.cmake new file mode 100644 index 0000000..01dd612 --- /dev/null +++ b/Auxiliary_utils/PCS_create/FindCFITSIO.cmake @@ -0,0 +1,67 @@ +# - Try to find CFITSIO +# Once done this will define +# +# CFITSIO_FOUND - system has CFITSIO +# CFITSIO_INCLUDE_DIR - the CFITSIO include directory +# CFITSIO_LIBRARIES - Link these to use CFITSIO +# CFITSIO_VERSION_STRING - Human readable version number of cfitsio +# CFITSIO_VERSION_MAJOR - Major version number of cfitsio +# CFITSIO_VERSION_MINOR - Minor version number of cfitsio + +# Copyright (c) 2006, Jasem Mutlaq +# Based on FindLibfacile by Carsten Niehaus, +# +# Redistribution and use is allowed according to the terms of the BSD license. +# For details see the accompanying COPYING-CMAKE-SCRIPTS file. + +if (CFITSIO_INCLUDE_DIR AND CFITSIO_LIBRARIES) + + # in cache already, be quiet + set(CFITSIO_FIND_QUIETLY TRUE) + +else (CFITSIO_INCLUDE_DIR AND CFITSIO_LIBRARIES) + + # JM: Packages from different distributions have different suffixes + find_path(CFITSIO_INCLUDE_DIR fitsio.h + PATH_SUFFIXES libcfitsio3 libcfitsio0 cfitsio + PATHS + $ENV{CFITSIO} + ${_obIncDir} + ${GNUWIN32_DIR}/include + ) + + find_library(CFITSIO_LIBRARIES NAMES cfitsio + PATHS + $ENV{CFITSIO} + ${_obLinkDir} + ${GNUWIN32_DIR}/lib + ) + + if(CFITSIO_INCLUDE_DIR AND CFITSIO_LIBRARIES) + set(CFITSIO_FOUND TRUE) + else (CFITSIO_INCLUDE_DIR AND CFITSIO_LIBRARIES) + set(CFITSIO_FOUND FALSE) + endif(CFITSIO_INCLUDE_DIR AND CFITSIO_LIBRARIES) + + + if (CFITSIO_FOUND) + + # Find the version of the cfitsio header + FILE(READ "${CFITSIO_INCLUDE_DIR}/fitsio.h" FITSIO_H) + STRING(REGEX REPLACE ".*#define CFITSIO_VERSION[^0-9]*([0-9]+)\\.([0-9]+).*" "\\1.\\2" CFITSIO_VERSION_STRING "${FITSIO_H}") + STRING(REGEX REPLACE "^([0-9]+)[.]([0-9]+)" "\\1" CFITSIO_VERSION_MAJOR ${CFITSIO_VERSION_STRING}) + STRING(REGEX REPLACE "^([0-9]+)[.]([0-9]+)" "\\2" CFITSIO_VERSION_MINOR ${CFITSIO_VERSION_STRING}) + message(STATUS "found version string ${CFITSIO_VERSION_STRING}") + + if (NOT CFITSIO_FIND_QUIETLY) + message(STATUS "Found CFITSIO ${CFITSIO_VERSION_MAJOR}.${CFITSIO_VERSION_MINOR}: ${CFITSIO_LIBRARIES}") + endif (NOT CFITSIO_FIND_QUIETLY) + else (CFITSIO_FOUND) + if (CFITSIO_FIND_REQUIRED) + message(STATUS "CFITSIO not found.") + endif (CFITSIO_FIND_REQUIRED) + endif (CFITSIO_FOUND) + + mark_as_advanced(CFITSIO_INCLUDE_DIR CFITSIO_LIBRARIES) + +endif (CFITSIO_INCLUDE_DIR AND CFITSIO_LIBRARIES) diff --git a/Auxiliary_utils/PCS_create/Makefile b/Auxiliary_utils/PCS_create/Makefile deleted file mode 100644 index 363ae96..0000000 --- a/Auxiliary_utils/PCS_create/Makefile +++ /dev/null @@ -1,6 +0,0 @@ -# -all: PCS_create - -PCS_create: PCS.c - gcc PCS.c -o PCS_create -lsofa_c -lm -lcfitsio - diff --git a/Auxiliary_utils/PCS_create/PCS.c b/Auxiliary_utils/PCS_create/PCS.c index 83fa7f3..6f06c7c 100644 --- a/Auxiliary_utils/PCS_create/PCS.c +++ b/Auxiliary_utils/PCS_create/PCS.c @@ -15,22 +15,23 @@ * along with this program. If not, see . */ +#include "cmdlnopts.h" +#include "sofatools.h" + #include +#include #include -#include -#include -#include -#include -#include #include #include -#include #include +#include +static const double hpa2mm = 760. / 1013.25; // hPa->mmHg -int status = 0; +static glob_pars *G = NULL; +static int status = 0; // check fits error status -int chkstatus(){ +static int chkstatus(){ int os = status; status = 0; if(os) fits_report_error(stderr, os); @@ -44,11 +45,10 @@ int chkstatus(){ * @param key (i) - keyword * @return value (static char) or NULL */ -char *getFITSkeyval(fitsfile *fptr, const char *key){ - char card[FLEN_CARD], newcard[FLEN_CARD], comment[FLEN_COMMENT]; +static char *getFITSkeyval(fitsfile *fptr, const char *key){ + char card[FLEN_CARD], comment[FLEN_COMMENT]; static char value[FLEN_VALUE]; int status = 0; - int keytype; if(fits_read_card(fptr, key, card, &status) || !*card){ fprintf(stderr, "Keyword %s does not exist or empty\n", key); return NULL; @@ -58,17 +58,15 @@ char *getFITSkeyval(fitsfile *fptr, const char *key){ return value; } -// safely convert string to double (@return 0 if all OK) -int getdouble(double *d, const char *str){ +// safely convert string to double (@return NULL if bad or return pointer to next non-digit) +static char *getdouble(double *d, const char *str){ double res = -1.; char *endptr; - if(!str) return 1; + if(!str || *str == 0) return NULL; res = strtod(str, &endptr); - if(endptr == str || *str == '\0'){ // || *endptr != '\0'){ - return 1; - } + if(endptr == str) return NULL; if(d) *d = res; - return 0; + return endptr; } /* // convert string like DD:MM:SS or HH:MM:SS into double, ishours==1 if HH @@ -86,72 +84,31 @@ int hd2d(double *dbl, const char *str, int ishours){ return 0; } */ -char *d2s(double dbl, int ishours){ +/** + * @brief d2s - convert angle to string + * @param dbl - angle value + * @param ishours - ==0 if angle is in degrees + * @return formed string + */ +static char *d2s(double dbl, int ishours){ int d, m; char *s = ""; if(ishours) dbl /= 15.; if(dbl < 0.){ s = "-"; dbl = -dbl; - } + }else if(!ishours) s = "+"; d = (int)dbl; dbl = (dbl-d)*60.; m = (int)dbl; dbl = (dbl-m)*60.; static char res[32]; - snprintf(res, 32, "%s%02d:%02d:%04.1f", s, d, m, dbl); + // no digits after decimal point in degrees + if(ishours) snprintf(res, 32, "%s%02d:%02d:%04.1f", s, d, m, dbl); + else snprintf(res, 32, "%s%02d:%02d:%02.0f", s, d, m, dbl); return res; } -typedef struct{ - double ra; - double dec; -} polar; - -/** - * @brief J2000toJnow - convert ra/dec between epochs - * @param in - J2000 (degrees) - * @param out - Jnow (degrees) - * @return - */ -int J2000toJnow(const polar *in, polar *out){ - if(!out) return 1; - double utc1, utc2; - time_t tsec; - struct tm *ts; - tsec = time(0); // number of seconds since the Epoch, 1970-01-01 00:00:00 +0000 (UTC) - ts = gmtime(&tsec); - int result = 0; - result = iauDtf2d ( "UTC", ts->tm_year+1900, ts->tm_mon+1, ts->tm_mday, ts->tm_hour, ts->tm_min, ts->tm_sec, &utc1, &utc2 ); - if (result != 0) { - fprintf(stderr, "iauDtf2d call failed\n"); - return 1; - } - // Make TT julian date for Atci13 call - double tai1, tai2; - double tt1, tt2; - result = iauUtctai(utc1, utc2, &tai1, &tai2); - if(result){ - fprintf(stderr, "iauUtctai call failed\n"); - return 1; - } - result = iauTaitt(tai1, tai2, &tt1, &tt2); - if(result){ - fprintf(stderr, "iauTaitt call failed\n"); - return 1; - } - double pr = 0.0; // RA proper motion (radians/year; Note 2) - double pd = 0.0; // Dec proper motion (radians/year) - double px = 0.0; // parallax (arcsec) - double rv = 0.0; // radial velocity (km/s, positive if receding) - double rc = DD2R * in->ra, dc = DD2R * in->dec; // convert into radians - double ri, di, eo; - iauAtci13(rc, dc, pr, pd, px, rv, tt1, tt2, &ri, &di, &eo); - out->ra = iauAnp(ri - eo) * DR2D; - out->dec = di * DR2D; - return 0; -} - /** * @brief getDval - get double value of keyword * @param ret - double value @@ -159,11 +116,11 @@ int J2000toJnow(const polar *in, polar *out){ * @param key - keyword to search * @return 0 if found and converted */ -int getDval(double *ret, fitsfile *fptr, const char *key){ +static int getDval(double *ret, fitsfile *fptr, const char *key){ if(!ret || !key) return 2; char *val = getFITSkeyval(fptr, key); if(!val) return 1; - if(getdouble(ret, val)){ + if(!getdouble(ret, val)){ fprintf(stderr, "Wrong %s value\n", key); return 1; } @@ -171,7 +128,7 @@ int getDval(double *ret, fitsfile *fptr, const char *key){ } // run command xy2sky @fname and return stdout (up to 1024 bytes) -char *exe(char *fname){ +static char *exe(char *fname){ #define die(text) do{fprintf(stderr, "%s\n", text); return NULL; }while(0) int link[2]; pid_t pid; @@ -190,6 +147,7 @@ char *exe(char *fname){ char *ptr = ret; while(0 != (r = read(link[0], ptr, nleave))){ ptr += r; + nleave -= r; *ptr = 0; } wait(NULL); @@ -198,58 +156,158 @@ char *exe(char *fname){ #undef die } -int main(int argc, char **argv) { - if(argc != 2){ - fprintf(stderr, "USAGE: %s fitsfile\n\tCalculate data :newalpt by given plate-recognized file\n", argv[0]); - return 1; - } - double ra_center = 12., dec_center = 21., ra_scope, dec_scope; +static int parse_fits_file(char *name){ + double ra_center = 400., dec_center = 400., ra_scope, dec_scope; // get CENTER: - char *val = exe(argv[1]); + char *val = exe(name); if(!val) return 1; - char *p = strchr(val, ' '); - if(!p) return 1; - getdouble(&ra_center, val); - getdouble(&dec_center, p); + DBG("EXE gives: %s", val); +// char *p = strchr(val, ' '); +// if(!p) return 1; + val = getdouble(&ra_center, val); + if(!val) return 1; + if(!getdouble(&dec_center, val)) return 1; + DBG("J2000=%g/%g", ra_center, dec_center); // get FITS keywords fitsfile *fptr; int iomode = READONLY; - fits_open_file(&fptr, argv[1], iomode, &status); + fits_open_file(&fptr, name, iomode, &status); iomode = chkstatus(); if(iomode) return iomode; if(getDval(&ra_scope, fptr, "RA")) return 4; if(getDval(&dec_scope, fptr, "DEC")) return 5; + double uxt; + if(getDval(&uxt, fptr, "UNIXTIME")) return 55; + struct timeval tv; + tv.tv_sec = (time_t) uxt; + tv.tv_usec = 0; val = getFITSkeyval(fptr, "PIERSIDE"); if(!val) return 6; char pierside = 'W'; if(strstr(val, "East")) pierside = 'E'; - val = getFITSkeyval(fptr, "LSTEND"); - if(!val) return 7; - char *s = strchr(val, '\''); - if(!s) return 8; - char *sidtm = strdup(++s); - s = strchr(sidtm, '\''); - if(!s) return 9; - *s = 0; fits_close_file(fptr, &status); chkstatus(); - polar J2000 = {.ra = ra_center, .dec = dec_center}, Jnow; - if(J2000toJnow(&J2000, &Jnow)) return 1; + polarCrds J2000 = {.ra = DD2R * ra_center, .dec = DD2R * dec_center}, Jnow; + DBG("J2000=%g/%g", ra_center, dec_center); + DBG("J2000=%g/%g", J2000.ra/DD2R, J2000.dec/DD2R); + if(get_ObsPlace(&tv, &J2000, &Jnow, NULL)) return 1; + DBG("JNOW: RA=%g, DEC=%g, EO=%g", Jnow.ra/DD2R, Jnow.dec/DD2R, Jnow.eo/DD2R); + sMJD mjd; + if(get_MJDt(&tv, &mjd)) return 1; + double ST; + almDut adut; + if(getDUT(&adut)) return 1; + placeData *place = getPlace(); + if(!place) return 1; + if(get_LST(&mjd, adut.DUT1, place->slong, &ST)) return 1; + ST /= DD2R; // convert radians to degrees - char *sra_scope = strdup(d2s(ra_scope, 1)), - *sdec_scope = strdup(d2s(dec_scope, 0)), - *sra_center = strdup(d2s(Jnow.ra, 1)), - *sdec_center= strdup(d2s(Jnow.dec, 0)); + double ra_now = (Jnow.ra - Jnow.eo)/DD2R, dec_now = Jnow.dec/DD2R; + DBG("RA_now=%g, DEC_now=%g", ra_now, dec_now); + if(G->ha){ // print HA instead of RA + ra_scope = ST - ra_scope; + if(ra_scope < 0.) ra_scope += D2PI; + ra_now = ST - ra_now; + if(ra_now < 0.) ra_now += D2PI; + } + if(G->delta){ + ra_now -= ra_scope; + dec_now -= dec_scope; + } + int rainhrs = !G->raindeg; + char *sidtm = NULL, *sra_scope = NULL, *sdec_scope = NULL, *sra_center = NULL, *sdec_center= NULL; + if(G->crdstrings){ // string form + sidtm = strdup(d2s(ST, G->stindegr ? 0 : 1)); + sra_scope = strdup(d2s(ra_scope, rainhrs)); + sdec_scope = strdup(d2s(dec_scope, 0)); + sra_center = strdup(d2s(ra_now, rainhrs)); + sdec_center= strdup(d2s(dec_now, 0)); + } - // create line like + // for 10-micron output create line like // :newalpt10:10:34.8,-12:21:14,E,10:10:3.6,-12:12:56,11:15:12.04# // MRA MDEC MSIDE PRA PDEC SIDTIME - printf(":newalpt%s,%s,%c,%s,%s,%s#\n", sra_scope, sdec_scope, pierside, sra_center, sdec_center, sidtm); - free(sra_scope); free(sdec_scope); - free(sra_center); free(sdec_center); - //free(sidtm); - - return(0); + // MRA: HH.MM.SS.S - mount-reported RA + // MDEC: sDD:MM:SS - mount-reported DEC + // MSIDE: 'E'/'W' - pier-side + // PRA: HH:MM:SS.S - plate-solved RA + // PDEC: sDD:MM:SS - plate-solved DEC + // SIDTIME: HH:MM:SS.S - local sid.time + if(G->for10m){ + printf(":newalpt%s,%s,%c,%s,%s,%s#\n", sra_scope, sdec_scope, pierside, sra_center, sdec_center, sidtm); + }else{ + if(G->crdstrings){ + printf("%-16s%-16s %c %-18s%-19s", sra_scope, sdec_scope, pierside, sra_center, sdec_center); + if(!G->ha) printf("%-15s", sidtm); + }else{ + if(!G->raindeg){ ra_scope /= 15.; ra_now /= 15.; } + printf("%-16.8f%-16.8f %c %-18.8f%-19.8f", ra_scope, dec_scope, pierside, ra_now, dec_now); + if(!G->ha) printf("%-15.8f", G->stindegr ? ST : ST/15.); + } + printf("%-15s\n", basename(name)); + } + FREE(sra_scope); FREE(sdec_scope); + FREE(sra_center); FREE(sdec_center); + FREE(sidtm); + return 0; +} + +static void printheader(){ + printf("# Pointing data @ p=%.f %s, T=%.1f degrC\n", G->pressure*(G->pmm ? hpa2mm : 1.), G->pmm ? "mmHg" : "hPa", G->temperature); + const char *raha = G->ha ? "HA" : "RA"; + const char *raunits, *decunits; + if(G->crdstrings){ + raunits = G->raindeg ? "dms" : "hms"; + decunits = "dms"; + }else{ + raunits = G->raindeg ? "deg" : "hrs"; + decunits = "deg"; + } + const char *apparent = G->delta ? "(app-enc)" : "Apparent"; + char a[4][32]; + snprintf(a[0], 32, "Encoder %s,%s", raha, raunits); + snprintf(a[1], 32, "Encoder DEC,%s", decunits); + snprintf(a[2], 32, "%s %s,%s", apparent, raha, raunits); + snprintf(a[3], 32, "%s DEC,%s", apparent, decunits); + printf("%-16s%-16s Pier %-18s%-19s", a[0], a[1], a[2], a[3]); + if(!G->ha){ + printf("Sid. time,"); + if(G->crdstrings){ + if(G->stindegr) printf("dms"); + else printf("hms"); + }else{ + if(G->stindegr) printf("deg"); + else printf("hrs"); + } + printf(" "); + } + printf("Filename\n"); +} + +int main(int argc, char **argv) { + initial_setup(); + G = parse_args(argc, argv); + if(G->pressure < 0.) ERRX("Pressure should be greater than zero"); + if(G->temperature < -100. || G->temperature > 100.) ERRX("Temperature over the range -100..+100"); + if(G->pmm) G->pressure /= hpa2mm; + setWeath(G->pressure, G->temperature, 0.5); + if(G->for10m){ + G->crdstrings = 1; + G->raindeg = 0; + G->ha = 0; + G->stindegr = 0; + } + if(G->printhdr){ + printheader(); + if(G->nfiles < 1) return 0; + } + if(G->nfiles < 1){ + WARNX("Need at least one FITS filename"); + return 1; + } + for(int i = 0; i < G->nfiles; ++i) + if(parse_fits_file(G->infiles[i])) WARNX("Can't parse file %s", G->infiles[i]); + return 0; } diff --git a/Auxiliary_utils/PCS_create/Readme b/Auxiliary_utils/PCS_create/Readme deleted file mode 100644 index 8591cc4..0000000 --- a/Auxiliary_utils/PCS_create/Readme +++ /dev/null @@ -1 +0,0 @@ -Calculate all data required for align model calculation (Pointing Correction System). \ No newline at end of file diff --git a/Auxiliary_utils/PCS_create/cmdlnopts.c b/Auxiliary_utils/PCS_create/cmdlnopts.c new file mode 100644 index 0000000..d26b6fd --- /dev/null +++ b/Auxiliary_utils/PCS_create/cmdlnopts.c @@ -0,0 +1,80 @@ +/* + * This file is part of the ttyterm project. + * Copyright 2020 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 // assert +#include // printf +#include // memcpy +#include +#include "cmdlnopts.h" + +/* + * here are global parameters initialisation + */ +static int help; +static glob_pars G; + +// DEFAULTS +// default global parameters +glob_pars const Gdefault = { + .pressure = 780. // 580mmHg + ,.temperature = 5. // @5degrC +}; + +/* + * Define command line options by filling structure: + * name has_arg flag val type argptr help +*/ +static myoption cmdlnopts[] = { + // set 1 to param despite of its repeating number: + {"help", NO_ARGS, NULL, 'h', arg_int, APTR(&help), _("show this help")}, + {"10m", NO_ARGS, NULL, 't', arg_int, APTR(&G.for10m), _("make output suitable for 10-micron mount")}, + {"header", NO_ARGS, NULL, 'H', arg_int, APTR(&G.printhdr), _("print header for output file")}, + {"raindeg", NO_ARGS, NULL, 'd', arg_int, APTR(&G.raindeg), _("RA in degrees instead of hours")}, + {"strings", NO_ARGS, NULL, 's', arg_int, APTR(&G.crdstrings),_("coordinates in string form")}, + {"ha", NO_ARGS, NULL, 0, arg_int, APTR(&G.ha), _("print HA instead of RA")}, + {"stdeg", NO_ARGS, NULL, 0, arg_int, APTR(&G.stindegr), _("sidereal time in degrees instead of hours")}, + {"delta", NO_ARGS, NULL, 'D', arg_int, APTR(&G.delta), _("show delta: apparent-encoder instead of apparent coordinates")}, + {"pressure",NEED_ARG, NULL, 'P', arg_double, APTR(&G.pressure), _("atmospheric pressure (hPa)")}, + {"pinmm", NO_ARGS, NULL, 'm', arg_int, APTR(&G.pmm), _("pressure in mmHg instead of hPa")}, + {"temperature",NEED_ARG,NULL, 'T', arg_double, APTR(&G.temperature),_("temperature, degrC")}, + end_option +}; + +/** + * Parse command line options and return dynamically allocated structure + * to global parameters + * @param argc - copy of argc from main + * @param argv - copy of argv from main + * @return allocated structure with global parameters + */ +glob_pars *parse_args(int argc, char **argv){ + void *ptr = memcpy(&G, &Gdefault, sizeof(G)); assert(ptr); + // format of help: "Usage: progname [args]\n" + change_helpstring(_("Usage: %s [args] FITS_files\nMake PCS list for equatorial mount\n\tWhere args are:\n")); + // parse arguments + parseargs(&argc, &argv, cmdlnopts); + if(help) showhelp(-1, cmdlnopts); + G.nfiles = argc; + G.infiles = MALLOC(char*, argc); + for(int i = 0; i < argc; i++){ + G.infiles[i] = strdup(argv[i]); + } + return &G; +} + diff --git a/Auxiliary_utils/PCS_create/cmdlnopts.h b/Auxiliary_utils/PCS_create/cmdlnopts.h new file mode 100644 index 0000000..97f52ed --- /dev/null +++ b/Auxiliary_utils/PCS_create/cmdlnopts.h @@ -0,0 +1,42 @@ +/* + * This file is part of the ttyterm project. + * Copyright 2020 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 +#ifndef CMDLNOPTS_H__ +#define CMDLNOPTS_H__ + +/* + * here are some typedef's for global data + */ +typedef struct{ + int for10m; // output suitable for 10-micron mount + int printhdr; // print header + int raindeg; // RA in degrees + int crdstrings; // coordinates in string form + int ha; // print HA instead of RA + int stindegr; // sidereal time in degrees instead of hours + int delta; // show delta: apparent-encoder instead of apparent coordinates + double pressure; // atmospheric pressure (HPa or mmHg if pmm==1) + int pmm; // pressure in mmHg + double temperature; // temperature, degrC + int nfiles; // amount of input files (amount of free arguments) + char **infiles; // input file[s] name[s] +} glob_pars; + +glob_pars *parse_args(int argc, char **argv); +#endif // CMDLNOPTS_H__ diff --git a/Auxiliary_utils/PCS_create/sofatools.c b/Auxiliary_utils/PCS_create/sofatools.c new file mode 100644 index 0000000..8d78c1c --- /dev/null +++ b/Auxiliary_utils/PCS_create/sofatools.c @@ -0,0 +1,280 @@ +/* + * This file is part of the sofa project. + * Copyright 2020 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 "sofatools.h" + +static placeData *pldata = NULL; +// temporal stubs for weather/place/DUT1 data; return 0 if all OK +placeData *getPlace(){ + if(pldata) return pldata; + pldata = malloc(sizeof(placeData)); + /* Site longitude, latitude (radians) and height above the geoid (m). */ + pldata->slong = 0.7232763200; + pldata->slat = 0.7618977414; + pldata->salt = 2070.0; // altitude + return pldata; +} +static placeWeather W = {0}; +// set weather parameters: pressure, temperature and humidity +void setWeath(double P, double T, double H){ + W.php = P; + W.tc = T; + W.relhum = H; +} +int getWeath(placeWeather *w){ + if(!w) return 0; + memcpy(w, &W, sizeof(placeWeather)); + return 0; +} +int getDUT(almDut *a){ + if(!a) return 0; + a->px = a->py = 0; + a->DUT1 = -0.25080; + return 0; +} + +char *radtodeg(double r){ + static char buf[128]; + int i[4]; char pm; + r = iauAnpm(r); + iauA2af(2, r, &pm, i); + snprintf(buf, 128, "%c%02d %02d %02d.%02d", pm, i[0],i[1],i[2],i[3]); + return buf; +} + +char *radtohrs(double r){ + static char buf[128]; + int i[4]; char pm; + r = iauAnp(r); + iauA2tf(2, r, &pm, i); + snprintf(buf, 128, "%02d:%02d:%02d.%02d", i[0],i[1],i[2],i[3]); + return buf; +} + +/** + * @brief get_MJDt - calculate MJD of date from argument + * @param tval (i) - given date (or NULL for current) + * @param MJD (o) - time (or NULL just to check) + * @return 0 if all OK + */ +int get_MJDt(struct timeval *tval, sMJD *MJD){ + struct tm tms; + double tSeconds; + if(!tval){ + //DBG("MJD for current time"); + 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. */ + if(iauDtf2d("UTC", y, m, d, tms.tm_hour, tms.tm_min, tSeconds, &utc1, &utc2) < 0) return -1; + if(!MJD) return 0; + MJD->MJD = utc1 - 2400000.5 + utc2; + MJD->utc1 = utc1; + MJD->utc2 = utc2; + //DBG("UTC(m): %g, %.8f\n", utc1 - 2400000.5, utc2); + if(iauUtctai(utc1, utc2, &MJD->tai1, &MJD->tai2)) return -1; + //DBG("TAI"); + if(iauTaitt(MJD->tai1, MJD->tai2, &MJD->tt1, &MJD->tt2)) return -1; + //DBG("TT"); + return 0; +} + +/** + * @brief get_LST - calculate local siderial time + * @param mjd (i) - date/time for LST (utc1 & tt used) + * @param dUT1 - (UT1-UTC) + * @param slong - site longitude (radians) + * @param LST (o) - local sidereal time (radians) + * @return 0 if all OK + */ +int get_LST(sMJD *mjd, double dUT1, double slong, double *LST){ + double ut11, ut12; + if(iauUtcut1(mjd->utc1, mjd->utc2, dUT1, &ut11, &ut12)) return 1; + /*double era = iauEra00(ut11, ut12) + slong; + double eo = iauEe06a(mjd->tt1, mjd->tt2); + printf("ERA = %s; ", radtohrs(era)); + printf("ERA-eo = %s\n", radtohrs(era-eo));*/ + if(!LST) return 0; + double ST = iauGst06a(ut11, ut12, mjd->tt1, mjd->tt2); + ST += slong; + if(ST > D2PI) ST -= D2PI; + else if(ST < -D2PI) ST += D2PI; + *LST = ST; + return 0; +} + +/** + * @brief hor2eq - convert horizontal coordinates to polar + * @param h (i) - horizontal coordinates + * @param pc (o) - polar coordinates + * @param sidTime - sidereal time + */ +void hor2eq(horizCrds *h, polarCrds *pc, double sidTime){ + if(!h || !pc) return; + placeData *p = getPlace(); + iauAe2hd(h->az, DPI/2. - h->zd, p->slat, &pc->ha, &pc->dec); // A,H -> HA,DEC; phi - site latitude + pc->ra = sidTime - pc->ha; + pc->eo = 0.; +} + +/** + * @brief eq2horH - convert polar coordinates to horizontal + * @param pc (i) - polar coordinates (only HA used) + * @param h (o) - horizontal coordinates + * @param sidTime - sidereal time + */ +void eq2horH(polarCrds *pc, horizCrds *h){ + if(!h || !pc) return; + placeData *p = getPlace(); + double alt; + iauHd2ae(pc->ha, pc->dec, p->slat, &h->az, &alt); + h->zd = DPI/2. - alt; +} + +/** + * @brief eq2hor - convert polar coordinates to horizontal + * @param pc (i) - polar coordinates (only RA used) + * @param h (o) - horizontal coordinates + * @param sidTime - sidereal time + */ +void eq2hor(polarCrds *pc, horizCrds *h, double sidTime){ + if(!h || !pc) return; + double ha = sidTime - pc->ra + pc->eo; + placeData *p = getPlace(); + double alt; + iauHd2ae(ha, pc->dec, p->slat, &h->az, &alt); + h->zd = DPI/2. - alt; +} + +/** + * @brief get_ObsPlace - calculate observed place (without PM etc) for given date @550nm + * @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 0 if all OK + */ +int get_ObsPlace(struct timeval *tval, polarCrds *p2000, polarCrds *pnow, horizCrds *hnow){ + double pr = 0.0; // RA proper motion (radians/year; Note 2) + double pd = 0.0; // Dec proper motion (radians/year) + double px = 0.0; // parallax (arcsec) + double rv = 0.0; // radial velocity (km/s, positive if receding) + sMJD MJD; + if(get_MJDt(tval, &MJD)) return -1; + if(!p2000) return -1; + placeData *p = getPlace(); + placeWeather w; + almDut d; + if(!p) return -1; + if(getWeath(&w)) return -1; + if(getDUT(&d)) return -1; + /* Effective wavelength (microns) */ + double wl = 0.55; + /* ICRS to observed. */ + double aob, zob, hob, dob, rob, eo; + if(iauAtco13(p2000->ra, p2000->dec, + pr, pd, px, rv, + MJD.utc1, MJD.utc2, + d.DUT1, + p->slong, p->slat, p->salt, + d.px, d.py, + w.php, w.tc, w.relhum, + wl, + &aob, &zob, + &hob, &dob, &rob, &eo)) return -1; + DBG("(RA/HA/DEC) J2000: %g/%g/%g; Jnow: %g/%g/%g", p2000->ra, p2000->ha, p2000->dec, rob, hob, dob); + if(pnow){ + pnow->eo = eo; + pnow->ha = hob; + pnow->ra = rob; + pnow->dec = dob; + } + if(hnow){ + hnow->az = aob; + hnow->zd = zob; + } + return 0; +} + + + + + + +#if 0 +typedef struct{ + double ra; + double dec; +} polar; + + +/** + * @brief J2000toJnow - convert ra/dec between epochs + * @param in - J2000 (degrees) + * @param out - Jnow (degrees) + * @return + */ +int J2000toJnow(const polar *in, polar *out){ + if(!out) return 1; + double utc1, utc2; + time_t tsec; + struct tm *ts; + tsec = time(0); // number of seconds since the Epoch, 1970-01-01 00:00:00 +0000 (UTC) + ts = gmtime(&tsec); + int result = 0; + result = iauDtf2d ( "UTC", ts->tm_year+1900, ts->tm_mon+1, ts->tm_mday, ts->tm_hour, ts->tm_min, ts->tm_sec, &utc1, &utc2 ); + if (result != 0) { + fprintf(stderr, "iauDtf2d call failed\n"); + return 1; + } + // Make TT julian date for Atci13 call + double tai1, tai2; + double tt1, tt2; + result = iauUtctai(utc1, utc2, &tai1, &tai2); + if(result){ + fprintf(stderr, "iauUtctai call failed\n"); + return 1; + } + result = iauTaitt(tai1, tai2, &tt1, &tt2); + if(result){ + fprintf(stderr, "iauTaitt call failed\n"); + return 1; + } + double pr = 0.0; // RA proper motion (radians/year; Note 2) + double pd = 0.0; // Dec proper motion (radians/year) + double px = 0.0; // parallax (arcsec) + double rv = 0.0; // radial velocity (km/s, positive if receding) + double rc = DD2R * in->ra, dc = DD2R * in->dec; // convert into radians + double ri, di, eo; + iauAtci13(rc, dc, pr, pd, px, rv, tt1, tt2, &ri, &di, &eo); + out->ra = iauAnp(ri - eo) * DR2D; + out->dec = di * DR2D; + return 0; +} +#endif diff --git a/Auxiliary_utils/PCS_create/sofatools.h b/Auxiliary_utils/PCS_create/sofatools.h new file mode 100644 index 0000000..210f2d1 --- /dev/null +++ b/Auxiliary_utils/PCS_create/sofatools.h @@ -0,0 +1,86 @@ +/* + * This file is part of the sofa project. + * Copyright 2020 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 +#ifndef SOFA_H__ +#define SOFA_H__ + +#include +#include +#include +#include +#include +#include +#include + +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 +} sMJD; + +// polar coordinates & equation of origins (all in radians) +typedef struct{ + double ha; // hour angle + double dec; // declination + double ra; // right ascension + double eo; // equation of origins +} polarCrds; + +// horizontal coordinates (all in radians) +typedef struct{ + double az; // azimuth, 0 @ south, positive clockwise + double zd; // zenith distance +} horizCrds; + +// observational place coordinates and altitude; all coordinates are in radians! +typedef struct{ + double slong; // longitude + double slat; // lattitude + double salt; // altitude, m +} placeData; + +// place weather data +typedef struct{ + double relhum; // rel. humidity, 0..1 + double php; // atm. pressure (hectopascales) + double tc; // temperature, degrC +} placeWeather; + +// DUT/polar almanach data +typedef struct{ + double DUT1; // UT1-UTC, sec + double px; // polar coordinates, arcsec + double py; +} almDut; + +placeData *getPlace(); +void setWeath(double P, double T, double H); +int getWeath(placeWeather *w); +int getDUT(almDut *a); +char *radtodeg(double r); +char *radtohrs(double r); +int get_MJDt(struct timeval *tval, sMJD *MJD); +int get_LST(sMJD *mjd, double dUT1, double slong, double *LST); +void hor2eq(horizCrds *h, polarCrds *pc, double sidTime); +void eq2horH(polarCrds *pc, horizCrds *h); +void eq2hor(polarCrds *pc, horizCrds *h, double sidTime);; +int get_ObsPlace(struct timeval *tval, polarCrds *p2000, polarCrds *pnow, horizCrds *hnow); + +#endif // SOFA_H__