Fixed bug in PCS_create (xy2sky data parsing), add some parameters

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")

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 <mutlaqja@ikarustech.com>
+# Based on FindLibfacile by Carsten Niehaus, <cniehaus@gmx.de>
+#
+# 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)

Auxiliary_utils/PCS_create/Makefile

@@ -1,6 +0,0 @@
-#
-all: PCS_create
-
-PCS_create: PCS.c
-	gcc PCS.c -o PCS_create -lsofa_c -lm -lcfitsio
-

Auxiliary_utils/PCS_create/PCS.c

@@ -15,22 +15,23 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include "cmdlnopts.h"
+#include "sofatools.h"
+
 #include <fitsio.h>
+#include <libgen.h>
 #include <math.h>
-#include <sofa.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <sys/time.h>
 #include <sys/types.h>
 #include <sys/wait.h>
-#include <time.h>
 #include <unistd.h>
+#include <usefull_macros.h>
 
+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){
+static char *getFITSkeyval(fitsfile *fptr, const char *key){
-    char card[FLEN_CARD], newcard[FLEN_CARD], comment[FLEN_COMMENT];
+    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)
+// safely convert string to double (@return NULL if bad or return pointer to next non-digit)
-int getdouble(double *d, const char *str){
+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'){
+    if(endptr == str) return NULL;
-        return 1;
-    }
     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) {
+static int parse_fits_file(char *name){
-    if(argc != 2){
+    double ra_center = 400., dec_center = 400., ra_scope, dec_scope;
-        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;
     // get CENTER:
-    char *val = exe(argv[1]);
+    char *val = exe(name);
     if(!val) return 1;
-    char *p = strchr(val, ' ');
+    DBG("EXE gives: %s", val);
-    if(!p) return 1;
+//    char *p = strchr(val, ' ');
-    getdouble(&ra_center, val);
+//    if(!p) return 1;
-    getdouble(&dec_center, p);
+    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;
+    polarCrds J2000 = {.ra = DD2R * ra_center, .dec = DD2R * dec_center}, Jnow;
-    if(J2000toJnow(&J2000, &Jnow)) return 1;
+    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)),
+    double ra_now = (Jnow.ra - Jnow.eo)/DD2R, dec_now = Jnow.dec/DD2R;
-            *sdec_scope = strdup(d2s(dec_scope, 0)),
+    DBG("RA_now=%g, DEC_now=%g", ra_now, dec_now);
-            *sra_center = strdup(d2s(Jnow.ra, 1)),
+    if(G->ha){ // print HA instead of RA
-            *sdec_center= strdup(d2s(Jnow.dec, 0));
+        ra_scope = ST - ra_scope;
-
+        if(ra_scope < 0.) ra_scope += D2PI;
-    // create line like
+        ra_now = ST - ra_now;
-    // :newalpt10:10:34.8,-12:21:14,E,10:10:3.6,-12:12:56,11:15:12.04#
+        if(ra_now < 0.) ra_now += D2PI;
-    //         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);
+    if(G->delta){
-    free(sra_scope); free(sdec_scope);
+        ra_now -= ra_scope;
-    free(sra_center); free(sdec_center);
+        dec_now -= dec_scope;
-    //free(sidtm);
+    }
-
+    int rainhrs = !G->raindeg;
-    return(0);
+    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));
+    }
+
+    // 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
+    // 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;
 }
 

Auxiliary_utils/PCS_create/Readme

@@ -1 +0,0 @@
-Calculate all data required for align model calculation (Pointing Correction System).

Auxiliary_utils/PCS_create/cmdlnopts.c

@@ -0,0 +1,80 @@
+/*
+ * This file is part of the ttyterm project.
+ * Copyright 2020 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+
+#include <assert.h> // assert
+#include <stdio.h>  // printf
+#include <string.h> // memcpy
+#include <usefull_macros.h>
+#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;
+}
+

Auxiliary_utils/PCS_create/cmdlnopts.h

@@ -0,0 +1,42 @@
+/*
+ * This file is part of the ttyterm project.
+ * Copyright 2020 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+#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__

Auxiliary_utils/PCS_create/sofatools.c

@@ -0,0 +1,280 @@
+/*
+ * This file is part of the sofa project.
+ * Copyright 2020 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+#include <usefull_macros.h>
+#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(&currentTime, NULL);
+        gmtime_r(&currentTime.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

Auxiliary_utils/PCS_create/sofatools.h

@@ -0,0 +1,86 @@
+/*
+ * This file is part of the sofa project.
+ * Copyright 2020 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+#ifndef SOFA_H__
+#define SOFA_H__
+
+#include <sofa.h>
+#include <sofam.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/time.h>
+#include <time.h>
+
+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__